Merge branch 'fix/minor-improvements' into 'master'

Fix/minor improvements

See merge request agnumpde/dune-elasticity!35

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -66,7 +66,12 @@ energy(const LocalView& localView,
     for (size_t i=0; i<localConfiguration.size(); i++)
         localAConfiguration[i] <<= localConfiguration[i];
 
-    energy = localEnergy_->energy(localView,localAConfiguration);
+    try {
+        energy = localEnergy_->energy(localView,localAConfiguration);
+    } catch (Dune::Exception &e) {
+        trace_off(rank);
+        throw e;
+    }
 
     energy >>= pureEnergy;
 

dune/elasticity/common/trustregionsolver.cc

@@ -383,6 +383,8 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
     CorrectionType rhs;
     MatrixType stiffnessMatrix;
 
+    Dune::Timer problemTimer;
+
     for (int i=0; i<maxTrustRegionSteps_; i++) {
 
         Dune::Timer totalTimer;
@@ -440,6 +442,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
         CorrectionType corr(rhs.size());
         corr = 0;
+        bool solvedByInnerSolver = true;
 
         //Take the obstacles on the finest grid and give them to the multigrid solver, it will create a hierarchy for all coarser grids
         trustRegionObstacles = trustRegion.obstacles();
@@ -519,97 +522,121 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
         std::cout << "Solve quadratic problem..." << std::endl;
 
         Dune::Timer solutionTimer;
-        innerSolver_->solve();
-        std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
 
-#if ! HAVE_DUNE_PARMG
-        if (mgStep)
-            corr = mgStep->getSol();
+        try {
+            innerSolver_->solve();
+        } catch (Dune::Exception &e) {
+            std::cerr << "Error while solving: " << e << std::endl;
+            solvedByInnerSolver = false;
+            corr = 0;
+        }
+        double energy = 0;
+        double totalModelDecrease = 0;
+        SolutionType newIterate = x_;
 
-        auto correctionInfinityNorm = corr.infinity_norm();
-#else
-        corr = iterationStep->getSol();
+        if (solvedByInnerSolver) {
 
-        auto correctionInfinityNorm = parallelInfinityNorm(corr);
-#endif
+            std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
 
-        //std::cout << "Correction: " << std::endl << corr << std::endl;
+    #if ! HAVE_DUNE_PARMG
+            if (mgStep)
+                corr = mgStep->getSol();
 
-        if (this->verbosity_ == NumProc::FULL)
-            std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
+            auto correctionInfinityNorm = corr.infinity_norm();
+    #else
+            corr = iterationStep->getSol();
 
-        if (correctionInfinityNorm < this->tolerance_) {
-            if (this->verbosity_ == NumProc::FULL and rank==0)
-                std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+            auto correctionInfinityNorm = parallelInfinityNorm(corr);
+    #endif
 
-            if (this->verbosity_ != NumProc::QUIET and rank==0)
-                std::cout << i+1 << " trust-region steps were taken." << std::endl;
-            break;
-        }
+            //std::cout << "Correction: " << std::endl << corr << std::endl;
 
-        if (trustRegion.radius() < this->tolerance_) {
-          if (this->verbosity_ == NumProc::FULL and rank==0)
-                std::cout << "TRUST REGION RADIUS IS TOO SMALL, SMALLER THAN TOLERANCE, STOPPING NOW" << std::endl;
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
 
-            if (this->verbosity_ != NumProc::QUIET and rank==0)
-                std::cout << i+1 << " trust-region steps were taken." << std::endl;
-            break;
-        }
+            if (correctionInfinityNorm < this->tolerance_) {
+                if (this->verbosity_ == NumProc::FULL and rank==0)
+                    std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
 
-        // ////////////////////////////////////////////////////
-        //   Check whether trust-region step can be accepted
-        // ////////////////////////////////////////////////////
+                if (this->verbosity_ != NumProc::QUIET and rank==0)
+                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
+                break;
+            }
 
-        SolutionType newIterate = x_;
-        for (size_t j=0; j<newIterate.size(); j++)
-            newIterate[j] += corr[j];
-
-        double energy    = assembler_->computeEnergy(newIterate);
-        energy = grid_->comm().sum(energy);
-
-        // compute the model decrease
-        // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
-        // Note that rhs = -g
-        CorrectionType tmp(corr.size());
-        tmp = 0;
-        hessianMatrix_->umv(corr, tmp);
-        double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
-        double totalModelDecrease = grid_->comm().sum(modelDecrease);
-
-        assert(totalModelDecrease >= 0);
-
-        double relativeModelDecrease = totalModelDecrease / std::fabs(energy);
-        double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
-
-        if (this->verbosity_ == NumProc::FULL and rank==0) {
-            std::cout << "Absolute model decrease: " << totalModelDecrease << std::endl;
-            std::cout << "Functional decrease: " << oldEnergy - energy << std::endl;
-            std::cout << "oldEnergy = " << oldEnergy << " and new energy = " << energy << std::endl;
-            std::cout << "Relative model decrease: " << relativeModelDecrease << std::endl;
-            std::cout << "Relative functional decrease: " << relativeFunctionalDecrease << std::endl;
-        }
+            if (trustRegion.radius() < this->tolerance_) {
+              if (this->verbosity_ == NumProc::FULL and rank==0)
+                    std::cout << "TRUST REGION RADIUS IS TOO SMALL, SMALLER THAN TOLERANCE, STOPPING NOW" << std::endl;
 
+                if (this->verbosity_ != NumProc::QUIET and rank==0)
+                    std::cout << i+1 << " trust-region steps were taken." << std::endl;
+                break;
+            }
 
-        if (energy >= oldEnergy) {
-            if (this->verbosity_ == NumProc::FULL and rank==0)
-                printf("Richtung ist keine Abstiegsrichtung!\n");
-        }
+            // ////////////////////////////////////////////////////
+            //   Check whether trust-region step can be accepted
+            // ////////////////////////////////////////////////////
 
-        if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
-            if (this->verbosity_ == NumProc::FULL and rank==0)
-                std::cout << "Suspecting rounding problems" << std::endl;
+            for (size_t j=0; j<newIterate.size(); j++)
+                newIterate[j] += corr[j];
 
-            if (this->verbosity_ != NumProc::QUIET and rank==0)
-                std::cout << i+1 << " trust-region steps were taken." << std::endl;
+            try {
+                energy  = assembler_->computeEnergy(newIterate);
+            } catch (Dune::Exception &e) {
+                std::cerr << "Error while computing the energy of the new Iterate: " << e << std::endl;
+                std::cerr << "Redoing trust region step with smaller radius..." << std::endl;
+                newIterate = x_;
+                solvedByInnerSolver = false;
+                energy = oldEnergy;
+            }
+            if (solvedByInnerSolver) {
 
-            x_ = newIterate;
-            break;
+                energy = grid_->comm().sum(energy);
+
+                // compute the model decrease
+                // It is $ m(x) - m(x+s) = -<g,s> - 0.5 <s, Hs>
+                // Note that rhs = -g
+                CorrectionType tmp(corr.size());
+                tmp = 0;
+                hessianMatrix_->umv(corr, tmp);
+                double modelDecrease = (rhs*corr) - 0.5 * (corr*tmp);
+                totalModelDecrease = grid_->comm().sum(modelDecrease);
+
+                assert(totalModelDecrease >= 0);
+
+                double relativeModelDecrease = totalModelDecrease / std::fabs(energy);
+                double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
+
+                if (this->verbosity_ == NumProc::FULL and rank==0) {
+                    std::cout << "Absolute model decrease: " << totalModelDecrease << std::endl;
+                    std::cout << "Functional decrease: " << oldEnergy - energy << std::endl;
+                    std::cout << "oldEnergy = " << oldEnergy << " and new energy = " << energy << std::endl;
+                    std::cout << "Relative model decrease: " << relativeModelDecrease << std::endl;
+                    std::cout << "Relative functional decrease: " << relativeFunctionalDecrease << std::endl;
+                }
+
+
+                if (energy >= oldEnergy) {
+                    if (this->verbosity_ == NumProc::FULL and rank==0)
+                        printf("Richtung ist keine Abstiegsrichtung!\n");
+                }
+
+                if (energy >= oldEnergy && (std::fabs(relativeFunctionalDecrease) < 1e-9 || std::fabs(relativeModelDecrease) < 1e-9)) {
+                    if (this->verbosity_ == NumProc::FULL and rank==0)
+                        std::cout << "Suspecting rounding problems" << std::endl;
+
+                    if (this->verbosity_ != NumProc::QUIET and rank==0)
+                        std::cout << i+1 << " trust-region steps were taken." << std::endl;
+
+                    x_ = newIterate;
+                    break;
+                }
+            }
         }
 
         // //////////////////////////////////////////////
         //   Check for acceptance of the step
         // //////////////////////////////////////////////
-        if (energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
+        if (solvedByInnerSolver && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
             // very successful iteration
 
             x_ = newIterate;
@@ -620,7 +647,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
             recomputeGradientHessian = true;
 
-        } else if ( ( (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
+        } else if (solvedByInnerSolver && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
             // successful iteration
             x_ = newIterate;
 
@@ -642,5 +669,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
 
         std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
     }
+    if (rank==0)
+        std::cout << "The whole trust-region step took " << problemTimer.elapsed() << " sec." << std::endl;
 
 }

dune/elasticity/materials/mooneyrivlindensity.hh

@@ -19,23 +19,26 @@ public:
   */
   MooneyRivlinDensity(const Dune::ParameterTree& parameters)
   {
-    // Mooneyrivlin constants
-    mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
-    mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
-    mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
-
-    mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
-    mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
-    mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
-    mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02");
-    mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11");
-    mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30");
-    mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
-    mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
-    mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
-    mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
-
     mooneyrivlin_energy = parameters.get<std::string>("mooneyrivlin_energy");
+    // Mooneyrivlin constants
+    if (mooneyrivlin_energy == "ciarlet") {
+      mooneyrivlin_a = parameters.get<double>("mooneyrivlin_a");
+      mooneyrivlin_b = parameters.get<double>("mooneyrivlin_b");
+      mooneyrivlin_c = parameters.get<double>("mooneyrivlin_c");
+    } else if (mooneyrivlin_energy == "log" or mooneyrivlin_energy == "square") {
+      mooneyrivlin_10 = parameters.get<double>("mooneyrivlin_10");
+      mooneyrivlin_01 = parameters.get<double>("mooneyrivlin_01");
+      mooneyrivlin_20 = parameters.get<double>("mooneyrivlin_20");
+      mooneyrivlin_02 = parameters.get<double>("mooneyrivlin_02");
+      mooneyrivlin_11 = parameters.get<double>("mooneyrivlin_11");
+      mooneyrivlin_30 = parameters.get<double>("mooneyrivlin_30");
+      mooneyrivlin_03 = parameters.get<double>("mooneyrivlin_03");
+      mooneyrivlin_21 = parameters.get<double>("mooneyrivlin_21");
+      mooneyrivlin_12 = parameters.get<double>("mooneyrivlin_12");
+      mooneyrivlin_k = parameters.get<double>("mooneyrivlin_k");
+    } else {
+      DUNE_THROW(Exception, "Error: Selected mooneyrivlin implementation (" << mooneyrivlin_energy << ") not available!");
+    }
   }
 
   /** \brief Evaluation with the deformation gradient
@@ -67,6 +70,9 @@ public:
     field_type normFSquared = gradient.frobenius_norm2();
     field_type detF = gradient.determinant();
 
+    if (detF < 0)
+      DUNE_THROW( Dune::Exception, "det(F) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
+
     field_type normFinvSquared = 0;
 
     field_type c2Tilde = 0;
@@ -81,11 +87,14 @@ public:
     // \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
     c2Tilde /= pow(detF, 4.0/dim);
     field_type c2TildeMinus3 = c2Tilde - 3;
-    field_type detFMinus1 = detF - 1;
 
     // Add the local energy density
-    auto strainEnergyCiarlet = (mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF));
-    auto strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
+    field_type strainEnergy = 0;
+
+    if (mooneyrivlin_energy == "ciarlet")
+      return mooneyrivlin_a*normFSquared + mooneyrivlin_b*normFinvSquared*detF + mooneyrivlin_c*detF*detF - ((dim-1)*mooneyrivlin_a + mooneyrivlin_b + 2*mooneyrivlin_c)*std::log(detF);
+    else {
+      strainEnergy = mooneyrivlin_10 * trCTildeMinus3 +
                         mooneyrivlin_01 * c2TildeMinus3  +
                         mooneyrivlin_20 * trCTildeMinus3 * trCTildeMinus3 +
                         mooneyrivlin_02 * c2TildeMinus3  * c2TildeMinus3  +
@@ -94,20 +103,14 @@ public:
                         mooneyrivlin_21 * trCTildeMinus3 * trCTildeMinus3 * c2TildeMinus3  +
                         mooneyrivlin_12 * trCTildeMinus3 * c2TildeMinus3  * c2TildeMinus3  +
                         mooneyrivlin_03 * c2TildeMinus3  * c2TildeMinus3  * c2TildeMinus3;
-    field_type logDetF = std::log(detF);
-    auto strainEnergyWithLog =  ( strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF );
-    auto strainEnergyWithSquare =  ( strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1);
-
-
-    std::cout << std::scientific;
-    if (mooneyrivlin_energy == "log") {
-      return strainEnergyWithLog;
-    } else if (mooneyrivlin_energy == "square") {
-      return strainEnergyWithSquare;
+      if (mooneyrivlin_energy == "log") {
+        field_type logDetF = std::log(detF);
+        return strainEnergy + 0.5 * mooneyrivlin_k* logDetF * logDetF;
+      } else if (mooneyrivlin_energy == "square") {
+        field_type detFMinus1 = detF - 1;
+        return strainEnergy + mooneyrivlin_k* detFMinus1 * detFMinus1;
+      }
     }
-    std::cout << std::fixed;
-
-    return strainEnergyCiarlet;
   }
 
   /** \brief Lame constants */

src/finite-strain-elasticity.cc

@@ -223,6 +223,7 @@ int main (int argc, char *argv[]) try
   vtkWriter.addVertexData(localDisplacementFunction, VTK::FieldInfo("displacement", VTK::FieldInfo::Type::vector, dim));
   vtkWriter.write(resultPath + "finite-strain_homotopy_0");
 
+  Dune::Timer homotopyTimer;
   for (int i=0; i<numHomotopySteps; i++)
   {
     double homotopyParameter = (i+1)*(1.0/numHomotopySteps);
@@ -354,6 +355,9 @@ int main (int argc, char *argv[]) try
     vtkWriter.write(resultPath + "finite-strain_homotopy_" + std::to_string(i+1));
   }
 
+  if (mpiHelper.rank()==0)
+    std::cout << "Complete duration: " << homotopyTimer.elapsed() << " sec." << std::endl;
+
 } catch (Exception& e) {
     std::cout << e.what() << std::endl;
 }