Add a simple proximal newton solver

dune/elasticity/common/pnsolver.cc

+#include <dune/common/bitsetvector.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/io.hh>
+#include <dune/istl/scaledidmatrix.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+
+#include <dune/fufem/functionspacebases/dunefunctionsbasis.hh>
+#include <dune/fufem/assemblers/istlbackend.hh>
+
+#if DUNE_VERSION_GTE(DUNE_SOLVERS, 2, 8)
+// Using a cholmod solver as the inner solver, available only since 2.8
+#include <dune/solvers/solvers/cholmodsolver.hh>
+#else
+// Using a umfpack solver as the inner solver
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#endif
+
+template <class BasisType, class VectorType>
+void ProximalNewtonSolver<BasisType,VectorType>::
+setup(const typename BasisType::GridView::Grid& grid,
+         std::conditional_t< // do we have a dune-functions basis? -> choose right assembler type
+           Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+             const Dune::Elasticity::FEAssembler<BasisType,VectorType>*,
+             const Dune::FEAssembler<DuneFunctionsBasis<BasisType>,VectorType>* > assembler,
+         const SolutionType& x,
+         const Dune::BitSetVector<blocksize>& dirichletNodes,
+         double tolerance,
+         int maxProximalNewtonSteps,
+         double initialRegularization)
+
+{
+    grid_                     = &grid;
+    assembler_                = assembler;
+    x_                        = x;
+    this->tolerance_          = tolerance;
+    maxProximalNewtonSteps_   = maxProximalNewtonSteps;
+    initialRegularization_    = initialRegularization;
+    ignoreNodes_              = std::make_shared<Dune::BitSetVector<blocksize>>(dirichletNodes);
+
+    const auto dim = VectorType::value_type::dimension;
+
+    //////////////////////////////////////////////////////////////////
+    //   Create the inner solver using a cholmod (or umfpack) solver
+    //////////////////////////////////////////////////////////////////
+#if DUNE_VERSION_GTE(DUNE_SOLVERS, 2, 8)
+    innerSolver_ = std::make_shared<Dune::Solvers::CholmodSolver<MatrixType,CorrectionType> >();
+#else
+    std::cout << "using umfpacksolver" << std::endl;
+    innerSolver_ = std::make_shared<Dune::Solvers::UMFPackSolver<MatrixType,CorrectionType> >();
+#endif
+    auto globalDirichletNodes = new Dune::BitSetVector<blocksize>(dirichletNodes);
+    innerSolver_->setIgnore(*globalDirichletNodes);
+
+    // ////////////////////////////////////////////////////////////
+    //    Create pointer to Hessian matrix and set its occupation structure
+    // ////////////////////////////////////////////////////////////
+    std::conditional_t< // do we have a dune-functions basis?
+      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+        BasisType,
+        DuneFunctionsBasis<BasisType> > basis(assembler_->basis_);
+    std::conditional_t< // do we have a dune-functions basis?
+      Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>(),
+        Dune::Fufem::DuneFunctionsOperatorAssembler<BasisType,BasisType>,
+        OperatorAssembler<DuneFunctionsBasis<BasisType>,DuneFunctionsBasis<BasisType>,Dune::Partitions::Interior> > operatorAssembler(basis,basis);
+
+    hessianMatrix_ = std::make_shared<MatrixType>();
+    if constexpr ( Dune::models<Dune::Functions::Concept::GlobalBasis<GridView>, BasisType>() )
+    {
+      auto hessianBackend = Dune::Fufem::istlMatrixBackend(*hessianMatrix_);
+      auto hessianPatternBuilder = hessianBackend.patternBuilder();
+      operatorAssembler.assembleBulkPattern(hessianPatternBuilder);
+      hessianPatternBuilder.setupMatrix();
+    }
+    else
+    {
+      Dune::MatrixIndexSet indices(grid_->size(1), grid_->size(1));
+      assembler_->getNeighborsPerVertex(indices);
+      indices.exportIdx(*hessianMatrix_);
+    }
+}
+
+
+template <class BasisType, class VectorType>
+void ProximalNewtonSolver<BasisType,VectorType>::solve()
+{
+    double oldEnergy = assembler_->computeEnergy(x_);
+
+    bool recomputeGradientHessian = true;
+    CorrectionType rhs;
+    MatrixType stiffnessMatrix;
+
+    Dune::Timer problemTimer;
+
+    double regularization = initialRegularization_;
+
+    for (int i=0; i<maxProximalNewtonSteps_; i++) {
+
+        Dune::Timer totalTimer;
+        if (this->verbosity_ == Solver::FULL) {
+            std::cout << "----------------------------------------------------" << std::endl;
+            std::cout << "      ProximalNewton Step Number: " << i
+                      << ",     regularization parameter: " << regularization
+                      << ",     energy: " << oldEnergy << std::endl;
+            std::cout << "----------------------------------------------------" << std::endl;
+        }
+
+        Dune::Timer gradientTimer;
+
+        if (recomputeGradientHessian) {
+
+            assembler_->assembleGradientAndHessian(x_,
+                                                   rhs,
+                                                   *hessianMatrix_,
+                                                   i==0    // assemble occupation pattern only for the first call
+                                                   );
+
+            rhs *= -1;        // The right hand side is the _negative_ gradient
+
+            // Compute gradient norm to monitor convergence
+            CorrectionType gradient = rhs;
+            for (size_t j=0; j<gradient.size(); j++)
+              for (size_t k=0; k<gradient[j].size(); k++)
+                if ((*ignoreNodes_)[j][k])
+                  gradient[j][k] = 0;
+
+            if (this->verbosity_ == Solver::FULL)
+              std::cout << "Gradient infinity norm: " << gradient.infinity_norm() << std::endl;
+            if (this->verbosity_ == Solver::FULL)
+              std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
+
+            // Transfer matrix data
+            stiffnessMatrix = *hessianMatrix_;
+
+            recomputeGradientHessian = false;
+        }
+
+        CorrectionType corr(rhs.size());
+        corr = 0;
+        bool solvedByInnerSolver = true;
+
+        // 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;
+
+        innerSolver_->setProblem(stiffnessMatrix,corr,rhs);
+        innerSolver_->preprocess();
+
+        ///////////////////////////////
+        //    Solve !
+        ///////////////////////////////
+
+        std::cout << "Solve quadratic problem..." << std::endl;
+
+        Dune::Timer solutionTimer;
+
+        try {
+            innerSolver_->solve();
+        } catch (Dune::Exception &e) {
+            std::cerr << "Error while solving: " << e << std::endl;
+            solvedByInnerSolver = false;
+            corr = 0;
+        }
+        double energy = 0;
+        double totalModelDecrease = 0;
+
+        double correctionInfinityNorm = corr.infinity_norm();
+        if (std::isnan(correctionInfinityNorm))
+            solvedByInnerSolver = false;
+
+        SolutionType newIterate = x_;
+
+        if (this->verbosity_ == NumProc::FULL)
+            std::cout << "Infinity norm of the correction: " << correctionInfinityNorm << std::endl;
+
+        if (correctionInfinityNorm < this->tolerance_ && this->tolerance_*100000 < 1.0/regularization ) {
+            if (this->verbosity_ == NumProc::FULL)
+              std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
+
+            if (this->verbosity_ != NumProc::QUIET)
+                std::cout << i+1 << " proximal-newton steps were taken." << std::endl;
+            break;
+        } else if (correctionInfinityNorm < this->tolerance_ )
+          std::cout << "regularization overflow!" << std::endl;
+
+        // ////////////////////////////////////////////////////
+        //   Check whether proximal-newton 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 proximal-newton step with higher regularization parameter..." << std::endl;
+            newIterate = x_;
+            solvedByInnerSolver = false;
+            energy = oldEnergy;
+        }
+
+        if (!solvedByInnerSolver) {
+            energy = oldEnergy;
+            newIterate = x_;
+        } else {
+            // 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);
+
+            assert(modelDecrease >= 0);
+
+            double relativeModelDecrease = modelDecrease / std::fabs(energy);
+            double relativeFunctionalDecrease = (oldEnergy - energy)/std::fabs(energy);
+
+            if (this->verbosity_ == NumProc::FULL) {
+                std::cout << "Absolute model decrease: " << modelDecrease << 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)
+                    printf("Correction caused an energy increase!\n");
+            }
+        }
+
+        // //////////////////////////////////////////////
+        //   Check for acceptance of the step
+        // //////////////////////////////////////////////
+        if (solvedByInnerSolver && energy < oldEnergy && (oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.9) {
+            // very successful iteration
+
+            x_ = newIterate;
+            regularization *= 0.5;
+
+            // current energy becomes 'oldEnergy' for the next iteration
+            oldEnergy = energy;
+
+            recomputeGradientHessian = true;
+
+        } else if (solvedByInnerSolver && ((oldEnergy-energy) / std::fabs(totalModelDecrease) > 0.01 || std::fabs(oldEnergy-energy) < 1e-12)) {
+            // successful iteration
+            x_ = newIterate;
+
+            // current energy becomes 'oldEnergy' for the next iteration
+            oldEnergy = energy;
+
+            recomputeGradientHessian = true;
+
+        } else {
+
+            // unsuccessful iteration
+
+            // Increase the regularization parameter
+            regularization *= 2;
+
+            if (this->verbosity_ == NumProc::FULL)
+                std::cout << "Unsuccessful iteration!" << std::endl;
+        }
+
+        std::cout << "iteration took " << totalTimer.elapsed() << " sec." << std::endl;
+    }
+    std::cout << "The whole proximal-newton step took " << problemTimer.elapsed() << " sec." << std::endl;
+
+}