diff --git a/dune/contact/solvers/filtermultigridsolver.cc b/dune/contact/solvers/filtermultigridsolver.cc index a2043cb30d4fcbda59fc7a03cbcb9e22e9ed922a..51eef118e9e09dfeaa0f2795211a2ff1c586c62a 100644 --- a/dune/contact/solvers/filtermultigridsolver.cc +++ b/dune/contact/solvers/filtermultigridsolver.cc @@ -212,8 +212,20 @@ void FilterMultigridContactSolver<ProblemType>::solve() " and infeasibility "<<newInfeasibility << std::endl; field_type infNorm = this->correction_.infinity_norm(); + + // if we want to properly decrease the trust-region according to the inf-norm, + // then we have to take into account the possible scaling + field_type scaledInfNorm = infNorm; + if (trScaling_) { // Changed from not trScaling! + const auto& trScaling = problem_->contactAssembler()->getColNorms(); + scaledInfNorm = 0; + for (size_t j = 0 ; j<this->correction_.size();j++) + for (int k=0; k < dim; k++) + scaledInfNorm = std::max(scaledInfNorm, std::fabs(this->correction_[j][k]/trScaling[j][k])); + } + if (!std::isfinite(newEnergy)) { - trustRegion = 0.1*std::min(trustRegion, infNorm); + trustRegion = 0.1*std::min(trustRegion, scaledInfNorm); correction_ = dirichletValues_; std::cout<<"Energy is not finite, reduce trust--region "<<trustRegion<<std::endl; continue; @@ -248,22 +260,26 @@ void FilterMultigridContactSolver<ProblemType>::solve() field_type eps = 1e-16; if (std::abs(modelDecrease) <eps or std::abs(energy_-newEnergy)<eps or - std::abs(relativeModelDecrease)<eps or std::abs(relativeEnergyDecrease)<eps - or std::abs(newInfeasibility) <eps) { + std::abs(relativeModelDecrease)<eps or std::abs(relativeEnergyDecrease)<eps) { + // or std::abs(newInfeasibility) <eps) { if (this->verbosity_ == Solver::FULL) std::cout << "Computing trust--region ratio might lead to cancellation...exiting." << std::endl; // if the trust-region norm is very small then increase it a bit to allow full steps - if (std::abs(infNorm-trustRegion)<eps) - trustRegion = std::min(initTR_, 2*infNorm); + if (std::abs(scaledInfNorm-trustRegion)<eps) + trustRegion = std::min(initTR_, 2*scaledInfNorm); break; } // if iterate is not acceptable to the filter // then reject it and reduce the trust-region - if (!filter_.acceptable(std::make_pair(infeasibility_, energy_), newInfeasibility, newEnergy)) { + bool acceptable{true}; + if (it == 0) + acceptable = (it == 0) ? filter_.acceptable(newInfeasibility, newEnergy) + : filter_.acceptable({infeasibility_, energy_}, newInfeasibility, newEnergy); + if (not acceptable) { // some heuristic to reduce the trust--region - trustRegion = 0.25*infNorm; + trustRegion = 0.25*std::min(trustRegion, scaledInfNorm); correction_ = dirichletValues_; if (this->verbosity_== Solver::FULL) @@ -280,7 +296,7 @@ void FilterMultigridContactSolver<ProblemType>::solve() if (suffDecrease && (ratio < 0.1)) { // some heuristic to reduce the trust--region - trustRegion = 0.25*infNorm; + trustRegion = 0.25*std::min(trustRegion, scaledInfNorm); correction_ = dirichletValues_; if (this->verbosity_ == Solver::FULL) @@ -305,7 +321,7 @@ void FilterMultigridContactSolver<ProblemType>::solve() std::cout << "Very successful step!\n"; // some heuristic to enlage the trust-region - trustRegion = std::min(maxTR_, std::fmax(2.5*infNorm, trustRegion)); + trustRegion = std::min(maxTR_, std::fmax(2.5*scaledInfNorm, trustRegion)); break; // Approximation is acceptable