From 83b6a702ea43a8b31295fe091e0f02427db84be4 Mon Sep 17 00:00:00 2001
From: Lisa Julia Nebel <lisa_julia.nebel@tu-dresden.de>
Date: Fri, 4 Sep 2020 18:49:53 +0200
Subject: [PATCH] WHITESPACES: Wrap energy calculation and solve-call in the
Mixed-Riemannian-TR solver in a try-catch statement
---
dune/gfe/mixedriemanniantrsolver.cc | 180 ++++++++++++++--------------
1 file changed, 90 insertions(+), 90 deletions(-)
diff --git a/dune/gfe/mixedriemanniantrsolver.cc b/dune/gfe/mixedriemanniantrsolver.cc
index 2fce2ad0..2cd74a7e 100644
--- a/dune/gfe/mixedriemanniantrsolver.cc
+++ b/dune/gfe/mixedriemanniantrsolver.cc
@@ -415,30 +415,30 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
Dune::Timer solutionTimer;
int ii = 0;
try {
- for (; ii<innerIterations_; ii++) {
- residual[_0] = rhs_global[_0];
- stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
- mmgStep0->setRhs(residual[_0]);
- mmgStep0->iterate();
-
- residual[_1] = rhs_global[_1];
- stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
- mmgStep1->setRhs(residual[_1]);
- mmgStep1->iterate();
-
- // Compute energy
- CorrectionType tmp(corr_global);
- stiffnessMatrix.mv(corr_global,tmp);
- double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
-
- if (energy > oldEnergy)
- //DUNE_THROW(Dune::Exception, "energy increase!");
- std::cout << "Warning: energy increase!" << std::endl;
-
- oldEnergy = energy;
- if (corr_global.infinity_norm() < innerTolerance_)
- break;
- }
+ for (; ii<innerIterations_; ii++) {
+ residual[_0] = rhs_global[_0];
+ stiffnessMatrix[_0][_1].mmv(corr_global[_1], residual[_0]);
+ mmgStep0->setRhs(residual[_0]);
+ mmgStep0->iterate();
+
+ residual[_1] = rhs_global[_1];
+ stiffnessMatrix[_1][_0].mmv(corr_global[_0], residual[_1]);
+ mmgStep1->setRhs(residual[_1]);
+ mmgStep1->iterate();
+
+ // Compute energy
+ CorrectionType tmp(corr_global);
+ stiffnessMatrix.mv(corr_global,tmp);
+ double energy = 0.5 * (tmp*corr_global) - (rhs_global*corr_global);
+
+ if (energy > oldEnergy)
+ //DUNE_THROW(Dune::Exception, "energy increase!");
+ std::cout << "Warning: energy increase!" << std::endl;
+
+ oldEnergy = energy;
+ if (corr_global.infinity_norm() < innerTolerance_)
+ break;
+ }
} catch (Dune::Exception &e) {
std::cerr << "Error while solving: " << e << std::endl;
solvedByInnerSolver = false;
@@ -462,78 +462,78 @@ void MixedRiemannianTrustRegionSolver<GridType,Basis,Basis0,TargetSpace0,Basis1,
SolutionType newIterate = x_;
if (solvedByInnerSolver) {
- if (this->verbosity_ == NumProc::FULL)
- std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
-
- if (corr_global.infinity_norm() < tolerance_) {
- if (verbosity_ == NumProc::FULL and rank==0)
- std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
-
- if (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[_0].size(); j++)
- newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
-
- for (size_t j=0; j<newIterate[_1].size(); j++)
- newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
-
- try {
- energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
- } 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) {
- energy = mpiHelper.getCollectiveCommunication().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);
- hessianMatrix_->mv(corr,tmp);
- modelDecrease = rhs*corr - 0.5 * (corr*tmp);
- modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
-
- double relativeModelDecrease = modelDecrease / std::fabs(energy);
-
- if (verbosity_ == NumProc::FULL and rank==0) {
- std::cout << "Absolute model decrease: " << modelDecrease
- << ", functional decrease: " << oldEnergy - energy << std::endl;
- std::cout << "Relative model decrease: " << relativeModelDecrease
- << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
- }
+ if (this->verbosity_ == NumProc::FULL)
+ std::cout << "Infinity norm of the correction: " << corr.infinity_norm() << std::endl;
+
+ if (corr_global.infinity_norm() < tolerance_) {
+ if (verbosity_ == NumProc::FULL and rank==0)
+ std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+ if (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[_0].size(); j++)
+ newIterate[_0][j] = TargetSpace0::exp(newIterate[_0][j], corr[_0][j]);
+
+ for (size_t j=0; j<newIterate[_1].size(); j++)
+ newIterate[_1][j] = TargetSpace1::exp(newIterate[_1][j], corr[_1][j]);
+
+ try {
+ energy = assembler_->computeEnergy(newIterate[_0],newIterate[_1]);
+ } 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) {
+ energy = mpiHelper.getCollectiveCommunication().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);
+ hessianMatrix_->mv(corr,tmp);
+ modelDecrease = rhs*corr - 0.5 * (corr*tmp);
+ modelDecrease = mpiHelper.getCollectiveCommunication().sum(modelDecrease);
+
+ double relativeModelDecrease = modelDecrease / std::fabs(energy);
+
+ if (verbosity_ == NumProc::FULL and rank==0) {
+ std::cout << "Absolute model decrease: " << modelDecrease
+ << ", functional decrease: " << oldEnergy - energy << std::endl;
+ std::cout << "Relative model decrease: " << relativeModelDecrease
+ << ", functional decrease: " << (oldEnergy - energy)/energy << std::endl;
+ }
- assert(modelDecrease >= 0);
+ assert(modelDecrease >= 0);
- if (energy >= oldEnergy and rank==0) {
- if (this->verbosity_ == NumProc::FULL)
- printf("Richtung ist keine Abstiegsrichtung!\n");
- }
+ if (energy >= oldEnergy and rank==0) {
+ if (this->verbosity_ == NumProc::FULL)
+ printf("Richtung ist keine Abstiegsrichtung!\n");
+ }
- if (energy >= oldEnergy &&
- (std::abs((oldEnergy-energy)/energy) < 1e-9 || relativeModelDecrease < 1e-9)) {
- if (this->verbosity_ == NumProc::FULL and rank==0)
- std::cout << "Suspecting rounding problems" << std::endl;
+ if (energy >= oldEnergy &&
+ (std::abs((oldEnergy-energy)/energy) < 1e-9 || 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;
+ if (this->verbosity_ != NumProc::QUIET and rank==0)
+ std::cout << i+1 << " trust-region steps were taken." << std::endl;
- x_ = newIterate;
- break;
- }
+ x_ = newIterate;
+ break;
+ }
- }
+ }
}
// //////////////////////////////////////////////
// Check for acceptance of the step
--
GitLab