Skip to content
Snippets Groups Projects

PARMG cleanup, write out deformationfunction to a file, Continue the trustregionsolver with a smaller trustregion radius in case of an exception

Closed PARMG cleanup, write out deformationfunction to a file, Continue the trustregionsolver with a smaller trustregion radius in case of an exception
(removed):master into master
1 unresolved thread
Closed lisa_julia.nebel_at_tu-dresden.de requested to merge (removed):master into master
1 unresolved thread
Viewing commit bb0c6049
Show latest version
2 files
+ 102
− 80
Compare changes
  • Side-by-side
  • Inline
Files
2
dune/elasticity/common/trustregionsolver.cc
+ 100
− 80
@@ -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;