Improve trust-region conditions a bit including the scaled inf-norm

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