Divide the regularization parameter through the scaling parameter in the PN...

Divide the regularization parameter through the scaling parameter in the PN solver instead of multiplying

In the PN solver, we add identity*regularization/scaling to the stiffness matrix. This has two effects:

• make the system matrix positive definite
• penalize large corrections to ensure an energy decrease in each step

The regularization parameter is adapted if either the system matrix is not positive definite if there is an energy increase. This works similar to the trust-region radius in the TR solver. The correction must be inside the trust-region, this ensures there is an energy decrease in each step. For each component of the correction, the trust-region radius as well as the regularization parameter is multiplied with a scaling factor, because the correction for the displacement works depends on the grid size where the correction for the rotation does not.

To align the effect of the scaling parameter, we need to divide through the scaling parameter in the PN solver and multiply with it it in the TR solver.

dune/gfe/riemannianpnsolver.cc

@@ -281,7 +281,7 @@ void RiemannianProximalNewtonSolver<Basis,TargetSpace,Assembler>::solve()
             // Add the regularization - Identity Matrix for now
             for (std::size_t i=0; i<stiffnessMatrix.N(); i++)
               for(int j=0; j<blocksize; j++)
-                stiffnessMatrix[i][i][j][j] += regularization*scaling_[j];
+                stiffnessMatrix[i][i][j][j] += regularization/scaling_[j];
 
             innerSolver_->setProblem(stiffnessMatrix,corr_global,rhs_global);
             innerSolver_->preprocess();