From bb0c60491ffe69e7b84e8f7ee02a72babd59e7bb Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Wed, 25 Mar 2020 13:53:56 +0100
Subject: [PATCH] =?UTF-8?q?Add=20try{}catch{}=20around=20the=20energy=20ca?=
=?UTF-8?q?lculation=20and=20the=20solve()=20call=20inside=20the=20trustre?=
=?UTF-8?q?gion=20solver=20Throw=20an=20exception=20in=20the=20Mooney-Rivl?=
=?UTF-8?q?in-Energy-Calculation=20if=20det(=E2=88=87=CF=86)=20<=200=20The?=
=?UTF-8?q?n:=20=20=20=20=20In=20case=20an=20excpetion=20is=20thrown=20in?=
=?UTF-8?q?=20either=20of=20these=20two=20calls,=20the=20trustregion=20sol?=
=?UTF-8?q?ver=20will=20not=20accept=20the=20step=20and=20continue=20with?=
=?UTF-8?q?=20a=20reduced=20radius?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
dune/elasticity/common/trustregionsolver.cc | 180 ++++++++++--------
.../materials/mooneyrivlinenergy.hh | 2 +
2 files changed, 102 insertions(+), 80 deletions(-)
diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc
index f7ed1ab..1440bde 100644
--- a/dune/elasticity/common/trustregionsolver.cc
+++ b/dune/elasticity/common/trustregionsolver.cc
@@ -316,6 +316,7 @@ void TrustRegionSolver<BasisType,VectorType>::solve()
CorrectionType corr(rhs.size());
corr = 0;
+ bool solved = 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();
@@ -394,97 +395,116 @@ 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;
+ try {
+ innerSolver_->solve();
+ } catch (Dune::Exception &e) {
+ std::cerr << "Error while solving: " << e << std::endl;
+ solved = false;
+ corr = 0;
+ }
+ double energy = 0;
+ double totalModelDecrease = 0;
+ SolutionType newIterate = x_;
+ if (solved) {
#if ! HAVE_DUNE_PARMG
- if (mgStep)
- corr = mgStep->getSol();
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl;
+ if (mgStep)
+ corr = mgStep->getSol();
- auto correctionInfinityNorm = corr.infinity_norm();
+ auto correctionInfinityNorm = corr.infinity_norm();
#else
- corr = iterationStep->getSol();
+ std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds and " << step << " steps." << std::endl;
+ corr = iterationStep->getSol();
- auto correctionInfinityNorm = parallelInfinityNorm(corr);
+ auto correctionInfinityNorm = parallelInfinityNorm(corr);
#endif
- //std::cout << "Correction: " << std::endl << corr << std::endl;
-
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
-
- if (correctionInfinityNorm < this->tolerance_) {
- if (this->verbosity_ == NumProc::FULL and rank==0)
- std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
- if (this->verbosity_ != NumProc::QUIET and rank==0)
- std::cout << i+1 << " trust-region steps were taken." << std::endl;
- break;
- }
-
- 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;
- }
-
- // ////////////////////////////////////////////////////
- // Check whether trust-region step can be accepted
- // ////////////////////////////////////////////////////
-
- 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 (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;
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
+
+ if (correctionInfinityNorm < this->tolerance_) {
+ if (this->verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (this->verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
+ break;
+ }
+
+ 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;
+ }
+
+ // ////////////////////////////////////////////////////
+ // Check whether trust-region step can be accepted
+ // ////////////////////////////////////////////////////
+
+ for (size_t j=0; j<newIterate.size(); j++)
+ newIterate[j] += corr[j];
+
+ 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_;
+ solved = false;
+ energy = oldEnergy;
+ }
+ if (solved) {
+ 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 (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 (solved && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
// very successful iteration
x_ = newIterate;
@@ -495,7 +515,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 (solved && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
// successful iteration
x_ = newIterate;
diff --git a/dune/elasticity/materials/mooneyrivlinenergy.hh b/dune/elasticity/materials/mooneyrivlinenergy.hh
index 508112e..6f58961 100644
--- a/dune/elasticity/materials/mooneyrivlinenergy.hh
+++ b/dune/elasticity/materials/mooneyrivlinenergy.hh
@@ -152,6 +152,8 @@ energy(const Entity& element,
for (int j = i+1; j < dim; j++)
c2Tilde += sigmaSquared[j]*sigmaSquared[i];
}
+ if (detF < 0)
+ DUNE_THROW( Dune::Exception, "det(∇φ) < 0, so it is not possible to calculate the MooneyRivlinEnergy.");
field_type trCTildeMinus3 = normFSquared/pow(detF, 2.0/dim) - 3;
// \tilde{D} = \frac{1}{\det{F}^{4/3}}D -> divide by det(F)^(4/3)= det(C)^(2/3)
--
GitLab