diff --git a/dune/elasticity/common/trustregionsolver.cc b/dune/elasticity/common/trustregionsolver.cc index 321e9e56bdb3b54bd447e65a0af17ebf109521ac..0f88e6d546e36f3fb4eba7091c9740c868f6a1fa 100644 --- a/dune/elasticity/common/trustregionsolver.cc +++ b/dune/elasticity/common/trustregionsolver.cc @@ -536,101 +536,101 @@ void TrustRegionSolver<BasisType,VectorType>::solve() if (solvedByInnerSolver) { - std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl; - -#if ! HAVE_DUNE_PARMG - if (mgStep) - corr = mgStep->getSol(); - - auto correctionInfinityNorm = corr.infinity_norm(); -#else - corr = iterationStep->getSol(); - - 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 - // //////////////////////////////////////////////////// - - 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_; - solvedByInnerSolver = false; - energy = oldEnergy; - } - if (solvedByInnerSolver) { - - 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; - } - } + std::cout << "Solving the quadratic problem took " << solutionTimer.elapsed() << " seconds." << std::endl; + + #if ! HAVE_DUNE_PARMG + if (mgStep) + corr = mgStep->getSol(); + + auto correctionInfinityNorm = corr.infinity_norm(); + #else + corr = iterationStep->getSol(); + + 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 + // //////////////////////////////////////////////////// + + 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_; + solvedByInnerSolver = false; + energy = oldEnergy; + } + if (solvedByInnerSolver) { + + 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; + } + } } // //////////////////////////////////////////////