Merge branch 'feature/new_tnnmg_allow_linearSolver' into 'feature/new_tnnmg'

Allow TNNMG to use Solvers::LinearSolver object

See merge request !5

dune/tnnmg/iterationsteps/tnnmgstep.hh

@@ -12,6 +12,7 @@
 #include "dune/solvers/iterationsteps/iterationstep.hh"
 #include "dune/solvers/iterationsteps/lineariterationstep.hh"
 #include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
 
 namespace Dune {
 namespace TNNMG {
@@ -33,11 +34,36 @@ class TNNMGStep :
 public:
 
   using Vector = typename F::Vector;
+  using ConstrainedVector = typename Linearization::ConstrainedVector;
+  using ConstrainedMatrix = typename Linearization::ConstrainedMatrix;
   using BitVector = typename Base::BitVector;
+  using ConstrainedBitVector = typename Linearization::ConstrainedBitvector;
   using Functional = F;
-  using LinearSolver = Solvers::IterativeSolver<typename Linearization::ConstrainedVector, Solvers::DefaultBitVector_t<Vector> >;
+  using IterativeSolver = Solvers::IterativeSolver< ConstrainedVector, Solvers::DefaultBitVector_t<ConstrainedVector> >;
+  using LinearSolver = Solvers::LinearSolver< ConstrainedMatrix,  ConstrainedVector >;
 
-  /** \brief Constructor with a solver object for the linear correction
+  /** \brief Constructor with an iterative solver object for the linear correction
+   * \param iterativeSolver This is a callback used to solve the constrained linearized system
+   * \param projection This is a callback used to compute a projection into a defect-admissible set
+   * \param lineSolver This is a callback used to minimize a directional restriction of the functional
+   *        for computing a damping parameter
+   */
+  TNNMGStep(const Functional& f,
+            Vector& x,
+            std::shared_ptr<IterationStep<Vector,BitVector> > nonlinearSmoother,
+            std::shared_ptr<IterativeSolver> iterativeSolver,
+            const DefectProjection& projection,
+            const LineSearchSolver& lineSolver)
+  : Base(x),
+    f_(&f),
+    nonlinearSmoother_(nonlinearSmoother),
+    preSmoothingSteps_(1),
+    iterativeSolver_(iterativeSolver),
+    projection_(projection),
+    lineSolver_(lineSolver)
+  {}
+
+  /** \brief Constructor with a linear solver object for the linear correction
    * \param linearSolver This is a callback used to solve the constrained linearized system
    * \param projection This is a callback used to compute a projection into a defect-admissible set
    * \param lineSolver This is a callback used to minimize a directional restriction of the functional
@@ -67,7 +93,7 @@ public:
   TNNMGStep(const Functional& f,
             Vector& x,
             std::shared_ptr<Solvers::IterationStep<Vector,BitVector> > nonlinearSmoother,
-            std::shared_ptr<Solvers::LinearIterationStep<typename Linearization::Matrix,typename Linearization::Vector> > linearIterationStep,
+            std::shared_ptr<Solvers::LinearIterationStep<ConstrainedMatrix,ConstrainedVector> > linearIterationStep,
             unsigned int noOfLinearIterationSteps,
             const DefectProjection& projection,
             const LineSearchSolver& lineSolver)
@@ -125,7 +151,10 @@ public:
     Solvers::resizeInitializeZero(correction_, x);
     Solvers::resizeInitializeZero(constrainedCorrection_, r);
 
-    auto emptyIgnore = ignore;
+    // TNNMGStep assumes that the linearization and the solver for the
+    // linearized problem will not use the ignoreNodes field
+    auto emptyIgnore = ConstrainedBitVector();
+
     Solvers::resizeInitialize(emptyIgnore, constrainedCorrection_, false);
 
     // Do the linear correction with a LinearIterationStep object
@@ -138,33 +167,39 @@ public:
       for (unsigned int i=0; i<noOfLinearIterationSteps_; i++)
         linearIterationStep_->iterate();
     }
-    else  // Do the linear correction with a Solver object
+    else if (iterativeSolver_)  // Do the linear correction with an iterative Solver object
     {
-    // Hand the linear problem to the linear solver.
-    // Currently we only support IterativeSolvers.  The IterationStep member
-    // needs to be a LinearIterationStep, so we can give it the matrix.
-    using LinearIterationStepType = Solvers::LinearIterationStep<std::decay_t<decltype(A)>,
-                                                                 typename Linearization::ConstrainedVector,
-                                                                 decltype(emptyIgnore) >;
-
-    LinearIterationStepType* linearIterationStep;
-
-    auto iterativeSolver = std::dynamic_pointer_cast<Solvers::IterativeSolver<typename Linearization::ConstrainedVector> >(linearSolver_);
-    if (iterativeSolver)
+      // Hand the linear problem to the iterative solver.
+      // The IterationStep member needs to be a LinearIterationStep,
+      // so we can give it the matrix.
+      using LinearIterationStepType = Solvers::LinearIterationStep<std::decay_t<decltype(A)>,
+                                                                  typename Linearization::ConstrainedVector,
+                                                                  decltype(emptyIgnore) >;
+
+      LinearIterationStepType* linearIterationStep;
+
+      auto iterativeSolver = std::dynamic_pointer_cast<Solvers::IterativeSolver<typename Linearization::ConstrainedVector> >(iterativeSolver_);
+      if (iterativeSolver)
+      {
+        iterativeSolver->getIterationStep().setIgnore(emptyIgnore);
+        linearIterationStep = dynamic_cast<LinearIterationStepType*>(&(iterativeSolver->getIterationStep()));
+      } else
+        DUNE_THROW(Exception, "Linear solver has to be an IterativeSolver!");
+
+      if (linearIterationStep)
+        linearIterationStep->setProblem(A, constrainedCorrection_, r);
+      else
+        DUNE_THROW(Exception, "Linear solver does not accept matrices!");
+
+
+      iterativeSolver_->preprocess();
+      iterativeSolver_->solve();
+    }
+    else // Do the linear correction with a linear Solver object
     {
-      iterativeSolver->getIterationStep().setIgnore(emptyIgnore);
-      linearIterationStep = dynamic_cast<LinearIterationStepType*>(&(iterativeSolver->getIterationStep()));
-    } else
-      DUNE_THROW(Exception, "Linear solver has to be an IterativeSolver!");
-
-    if (linearIterationStep)
-      linearIterationStep->setProblem(A, constrainedCorrection_, r);
-    else
-      DUNE_THROW(Exception, "Linear solver does not accept matrices!");
-
-
-    linearSolver_->preprocess();
-    linearSolver_->solve();
+      linearSolver_->setProblem(A,constrainedCorrection_, r);
+      linearSolver_->preprocess();
+      linearSolver_->solve();
     }
 
     linearization_->extendCorrection(constrainedCorrection_, correction_);
@@ -209,12 +244,15 @@ private:
   std::shared_ptr<Linearization> linearization_;
 
   // The following members cannot all be used at once:
-  // either we have a Dune::Solvers::Solver that implements the linear correction...
+  // either we have a Dune::Solvers::IterativeSolver that implements the linear correction...
+  std::shared_ptr<IterativeSolver> iterativeSolver_;
+
+  // or we have a Dune::Solvers::LinearSolver that implements the linear correction...
   std::shared_ptr<LinearSolver> linearSolver_;
 
   // ... or we have a mere LinearIterationStep, together with a number of times
-  // it is supposed to be called.  You cannot have both at once.
-  std::shared_ptr<LinearIterationStep<typename Linearization::Matrix,typename Linearization::Vector> > linearIterationStep_;
+  // it is supposed to be called.  You cannot have more than one at once.
+  std::shared_ptr<LinearIterationStep<typename Linearization::ConstrainedMatrix,typename Linearization::ConstrainedVector> > linearIterationStep_;
   unsigned int noOfLinearIterationSteps_;
 
   typename Linearization::ConstrainedVector constrainedCorrection_;