From 9ad6f95ca90991c30faf06925a9214fc5ca310ca Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 3 Sep 2014 15:29:23 +0000
Subject: [PATCH] Monitor the L^2-norm of the gradient

[[Imported from SVN: r9878]]
---
 dune/gfe/riemanniantrsolver.cc | 28 +++++++++++++++++++++++++---
 dune/gfe/riemanniantrsolver.hh |  3 +++
 2 files changed, 28 insertions(+), 3 deletions(-)

diff --git a/dune/gfe/riemanniantrsolver.cc b/dune/gfe/riemanniantrsolver.cc
index 7f57ff7a..1767d62c 100644
--- a/dune/gfe/riemanniantrsolver.cc
+++ b/dune/gfe/riemanniantrsolver.cc
@@ -136,6 +136,19 @@ setup(const GridType& grid,
                                                                                                    h1SemiNorm_,
                                                                                                  Solver::QUIET));
 
+    // //////////////////////////////////////////////////////////////////////////////////////
+    //   Assemble a mass matrix to create a norm that's equivalent to the L2-norm
+    //   This will be used to monitor the gradient
+    // //////////////////////////////////////////////////////////////////////////////////////
+
+    MassAssembler<GridType, typename BasisType::LocalFiniteElement, typename BasisType::LocalFiniteElement> massStiffness;
+    ScalarMatrixType localMassMatrix;
+
+    operatorAssembler.assemble(massStiffness, localMassMatrix);
+
+    ScalarMatrixType* massMatrix = new ScalarMatrixType(matrixComm.reduceAdd(localMassMatrix));
+    l2Norm_ = std::make_shared<H1SemiNorm<CorrectionType> >(*massMatrix);
+
     // Write all intermediate solutions, if requested
     if (instrumented_
         && dynamic_cast<IterativeSolver<CorrectionType>*>(innerSolver_.get()))
@@ -309,15 +322,24 @@ void RiemannianTrustRegionSolver<GridType,TargetSpace>::solve()
 
             rhs *= -1;        // The right hand side is the _negative_ gradient
 
+            // Transfer vector data
+            rhs_global = vectorComm.reduceAdd(rhs);
+
+            CorrectionType gradient = rhs_global;
+            for (size_t j=0; j<gradient.size(); j++)
+              for (int k=0; k<gradient[j].size(); k++)
+                if ((*ignoreNodes_)[j][k])
+                  gradient[j][k] = 0;
+
+            if (this->verbosity_ == Solver::FULL and rank==0)
+              std::cout << "Gradient norm: " << l2Norm_->operator()(gradient) << std::endl;
+
             if (this->verbosity_ == Solver::FULL)
               std::cout << "Assembly took " << gradientTimer.elapsed() << " sec." << std::endl;
 
             // Transfer matrix data
             stiffnessMatrix = matrixComm.reduceAdd(*hessianMatrix_);
 
-            // Transfer vector data
-            rhs_global = vectorComm.reduceAdd(rhs);
-
             recomputeGradientHessian = false;
 
         }
diff --git a/dune/gfe/riemanniantrsolver.hh b/dune/gfe/riemanniantrsolver.hh
index 44b2a7b2..1bb0f998 100644
--- a/dune/gfe/riemanniantrsolver.hh
+++ b/dune/gfe/riemanniantrsolver.hh
@@ -136,6 +136,9 @@ protected:
     /** \brief The norm used to measure multigrid convergence */
     H1SemiNorm<CorrectionType>* h1SemiNorm_;
 
+    /** \brief An L2-norm, really.  The H1SemiNorm class is badly named */
+    std::shared_ptr<H1SemiNorm<CorrectionType> > l2Norm_;
+
     /** \brief If set to true we log convergence speed and other stuff */
     bool instrumented_;
 
-- 
GitLab