Add some options to RiemannianTrustRegionSolver

Add new setup method, option to choose different norms (infinity-norm and H1-seminorm) for stopping criterion and read instrumented-path from parameterFile.

dune/gfe/riemanniantrsolver.cc

+#include <filesystem>
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/timer.hh>
 
@@ -28,45 +29,96 @@
 #include <dune/gfe/parallel/vectorcommunicator.hh>
 #endif
 
-template <class Basis, class TargetSpace, class Assembler>
-void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::
-setup(const GridType& grid,
-      const Assembler* assembler,
-         const SolutionType& x,
-         const Dune::BitSetVector<blocksize>& dirichletNodes,
-         double tolerance,
-         int maxTrustRegionSteps,
-         double initialTrustRegionRadius,
-         int multigridIterations,
-         double mgTolerance,
-         int mu,
-         int nu1,
-         int nu2,
-         int baseIterations,
-         double baseTolerance,
-         bool instrumented)
-{
-    int rank = grid.comm().rank();
+    template <class Basis, class TargetSpace, class Assembler>
+    void RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
+        setup(const GridType& grid,
+            const Assembler* assembler,
+            const SolutionType& x,
+            const Dune::BitSetVector<blocksize>& dirichletNodes,
+            const Dune::ParameterTree& parameterSet)
+    {
+
+        const double tolerance                = parameterSet.get<double>("tolerance");
+        const int maxTrustRegionSteps         = parameterSet.get<int>("maxTrustRegionSteps");
+        const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
+        const int multigridIterations         = parameterSet.get<int>("numIt");
+        const double mgTolerance              = parameterSet.get<double>("mgTolerance");
+        const int nu1                         = parameterSet.get<int>("nu1");
+        const int nu2                         = parameterSet.get<int>("nu2");
+        const int mu                          = parameterSet.get<int>("mu");
+        const int baseIterations              = parameterSet.get<int>("baseIt");
+        const double baseTolerance            = parameterSet.get<double>("baseTolerance");
+        const bool instrumented               = parameterSet.get<bool>("instrumented", 0);
+
+        norm_                                 = parameterSet.get("norm", "infinity");
+        instrumentedPath_                     = parameterSet.get("instrumentedPath", "/tmp");
+
+        //create intrumented-folder and mgHistory subfolder if it does not exist.
+        if (!(std::filesystem::exists(instrumentedPath_)))
+        {
+            std::filesystem::create_directory(instrumentedPath_);
+            std::filesystem::create_directory(instrumentedPath_ + "/mgHistory");
+        }
+
+        setup(grid,
+              assembler,
+              x,
+              dirichletNodes,
+              tolerance,
+              maxTrustRegionSteps,
+              initialTrustRegionRadius,
+              multigridIterations,
+              mgTolerance,
+              mu, nu1, nu2,
+              baseIterations,
+              baseTolerance,
+              instrumented);
+
+
+    }
+
+
+
+    template <class Basis, class TargetSpace, class Assembler>
+    void RiemannianTrustRegionSolver<Basis, TargetSpace, Assembler>::
+        setup(const GridType& grid,
+              const Assembler* assembler,
+              const SolutionType& x,
+              const Dune::BitSetVector<blocksize>& dirichletNodes,
+              double tolerance,
+              int maxTrustRegionSteps,
+              double initialTrustRegionRadius,
+              int multigridIterations,
+              double mgTolerance,
+              int mu,
+              int nu1,
+              int nu2,
+              int baseIterations,
+              double baseTolerance,
+              bool instrumented)
+    {
+        int rank = grid.comm().rank();
+
+        grid_ = &grid;
+        assembler_ = assembler;
+        x_ = x;
+        this->tolerance_ = tolerance;
+        maxTrustRegionSteps_ = maxTrustRegionSteps;
+        initialTrustRegionRadius_ = initialTrustRegionRadius;
+        innerIterations_ = multigridIterations;
+        innerTolerance_ = mgTolerance;
+        instrumented_ = instrumented;
+        ignoreNodes_ = &dirichletNodes;
+
+        int numLevels = grid_->maxLevel() + 1;
 
-    grid_                     = &grid;
-    assembler_                = assembler;
-    x_                        = x;
-    this->tolerance_          = tolerance;
-    maxTrustRegionSteps_      = maxTrustRegionSteps;
-    initialTrustRegionRadius_ = initialTrustRegionRadius;
-    innerIterations_          = multigridIterations;
-    innerTolerance_           = mgTolerance;
-    instrumented_             = instrumented;
-    ignoreNodes_              = &dirichletNodes;
 
-    int numLevels = grid_->maxLevel()+1;
 
 // The VectorCommunicator and MatrixCommunicator work only for GRID_DIM == WORLD_DIM == 2 or GRID_DIM == WORLD_DIM == 3
 #if HAVE_MPI && (!defined(GRID_DIM) or (defined(GRID_DIM) && GRID_DIM < 3)) && (!defined(WORLD_DIM) or (defined(WORLD_DIM) && WORLD_DIM < 3))
     //////////////////////////////////////////////////////////////////
     //  Create global numbering for matrix and vector transfer
     //////////////////////////////////////////////////////////////////
-
     globalMapper_ = std::make_unique<GlobalMapper>(grid_->leafGridView());
 #endif
 
@@ -171,7 +223,7 @@ setup(const GridType& grid,
     // Write all intermediate solutions, if requested
     if (instrumented_
         && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
-        dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get())->historyBuffer_ = "tmp/mgHistory";
+        dynamic_cast<IterativeSolver<CorrectionType> *>(innerSolver_.get())->historyBuffer_ = instrumentedPath_ + "/mgHistory";
 
     // ////////////////////////////////////////////////////////////
     //    Create Hessian matrix and its occupation structure
@@ -482,9 +534,16 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
         corr = corr_global;
 #endif
 
-        // Make infinity norm of corr_global known on all processors
-        double corrNorm = corr.infinity_norm();
-        double corrGlobalInfinityNorm = grid_->comm().max(corrNorm);
+        double corrNorm;
+        if (norm_ == "infinity")
+            corrNorm = corr.infinity_norm();
+        else if (norm_ == "H1-Semi")
+            corrNorm = h1SemiNorm_->operator()(corr);
+        else
+            DUNE_THROW(Dune::Exception, "Unknown norm type for stopping criterion!");
+
+        // Make norm of corr_global known on all processors
+        double corrGlobalNorm = grid_->comm().max(corrNorm);
 
         if (instrumented_) {
 
@@ -512,10 +571,10 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
                 // read iteration from file
                 CorrectionType intermediateSol(grid_->size(gridDim));
                 intermediateSol = 0;
-                char iSolFilename[100];
-                sprintf(iSolFilename, "tmp/mgHistory/intermediatesolution_%04d", j);
+                char iSolFilename[200];
+                sprintf(iSolFilename, (instrumentedPath_ + "/mgHistory/intermediatesolution_%04d").c_str(), j);
 
-                FILE* fpInt = fopen(iSolFilename, "rb");
+                FILE *fpInt = fopen(iSolFilename, "rb");
                 if (!fpInt)
                     DUNE_THROW(Dune::IOError, "Couldn't open intermediate solution");
                 for (size_t k=0; k<intermediateSol.size(); k++)
@@ -560,10 +619,12 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
             std::cout << i+1 << " trust-region steps were taken, the maximum was reached." << std::endl << "Total solver time: " << totalSolverTime << " sec., total assembly time: " << totalAssemblyTime << " sec." << std::endl;
 
         if (solved) {
+
             if (this->verbosity_ == NumProc::FULL and rank==0)
-                std::cout << "Infinity norm of the correction: " << corrGlobalInfinityNorm << std::endl;
+                std::cout << norm_ + "Norm of the correction: " << corrNorm << std::endl;
+
+            if (corrGlobalNorm < this->tolerance_ && corrNorm < trustRegion.radius()*smallestScalingParameter) {
 
-            if (corrGlobalInfinityNorm < this->tolerance_ && corrGlobalInfinityNorm < trustRegion.radius()*smallestScalingParameter) {
                 if (this->verbosity_ == NumProc::FULL and rank==0)
                     std::cout << "CORRECTION IS SMALL ENOUGH" << std::endl;
 
@@ -681,10 +742,10 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
 
         if (instrumented_) {
 
-            char iFilename[100];
-            sprintf(iFilename, "tmp/intermediateSolution_%04ld", i);
+            char iFilename[200];
+            sprintf(iFilename, (instrumentedPath_ + "/intermediateSolution_%04ld").c_str(), i);
 
-            FILE* fpIterate = fopen(iFilename, "wb");
+            FILE *fpIterate = fopen(iFilename, "wb");
             if (!fpIterate)
                 DUNE_THROW(SolverError, "Couldn't open file " << iFilename << " for writing");
 
@@ -692,7 +753,6 @@ void RiemannianTrustRegionSolver<Basis,TargetSpace,Assembler>::solve()
                 fwrite(&x_[j], sizeof(TargetSpace), 1, fpIterate);
 
             fclose(fpIterate);
-
         }
 
         if (rank==0)

[Sander, Oliver]

No need to copy all these parameters to variables just to use them once later. Use parameterSet.get directly as argument to setup below.

[Sander, Oliver]

Ideally, this would be a separate commit, because it adds new functionality.

[Sander, Oliver]

Omit these empty lines.

[Sander, Oliver]

Please keep the indentation here.

[Sander, Oliver]

Omit this spurious empty line.

[Sander, Oliver]

This is yet another new feature, and could also be a separate commit.

[Sander, Oliver]

Keep the formatting as it is.

[Sander, Oliver]

Spurious empty line.