towards clean solver construction

dune/tectonic/spatial-solving/functionalfactory.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_GLOBALPROBLEM_HH
+
+#include <vector>
+#include <memory>
+
+#include "tnnmg/functional.hh"
+
+template <class NBodyAssembler, class Nonlinearity, class MatrixType, class VectorType>
+class FunctionalFactory {
+    using Vector = VectorType;
+    using Matrix = MatrixType;
+
+    const static int dims = Vector::block_type::dimension;
+
+public:
+    using Functional = Functional<Matrix&, Vector&, Nonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
+
+    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<Matrix>& matrices, const std::vector<Vector>& rhs) :
+        nBodyAssembler_(nBodyAssembler),
+        nonlinearity_(nonlinearity),
+        rhs_(rhs) {
+
+        const auto nBodies = nBodyAssembler_.nGrids();
+
+        matrices_.resize(nBodies);
+        for (int i=0; i<nBodies; i++) {
+            matrices_[i] = &matrices[i];
+        }
+    }
+
+    FunctionalFactory(const NBodyAssembler& nBodyAssembler, const Nonlinearity& nonlinearity, const std::vector<std::shared_ptr<Matrix>>& matrices, const std::vector<Vector>& rhs) :
+        nBodyAssembler_(nBodyAssembler),
+        nonlinearity_(nonlinearity),
+        rhs_(rhs) {
+
+        const auto nBodies = nBodyAssembler_.nGrids();
+
+        matrices_.resize(nBodies);
+        for (size_t i=0; i<nBodies; i++) {
+            matrices_[i] = matrices[i].get();
+      }
+    }
+
+    void build() {
+        // assemble full global contact problem
+        nBodyAssembler_.assembleJacobian(matrices_, bilinearForm_);
+        nBodyAssembler_.assembleRightHandSide(rhs_, totalRhs_);
+
+        // get lower and upper obstacles
+        const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+        lower_.resize(totalObstacles.size());
+        upper_.resize(totalObstacles.size());
+
+        for (size_t j=0; j<totalObstacles.size(); ++j) {
+            const auto& totalObstaclesj = totalObstacles[j];
+            auto& lowerj = lower_[j];
+            auto& upperj = upper_[j];
+          for (size_t d=0; d<dims; ++d) {
+              lowerj[d] = totalObstaclesj[d][0];
+              upperj[d] = totalObstaclesj[d][1];
+          }
+        }
+
+        // set up functional
+        functional_ = std::make_shared<Functional>(bilinearForm_, totalRhs_, nonlinearity_, lower_, upper_);
+    }
+
+    std::shared_ptr<Functional> functional() const {
+        return functional_;
+    }
+
+    const Matrix& bilinearForm() const {
+        return bilinearForm_;
+    }
+
+    const Vector& rhs() const {
+        return totalRhs_;
+    }
+
+    const Vector& lower() const {
+        return lower_;
+    }
+
+    const Vector& upper() const {
+        return upper_;
+    }
+
+private:
+    const NBodyAssembler nBodyAssembler_;
+    const Nonlinearity& nonlinearity_;
+    std::vector<const Matrix*> matrices_;
+    const std::vector<Vector>& rhs_;
+
+    Matrix bilinearForm_;
+    Vector totalRhs_;
+    Vector lower_;
+    Vector upper_;
+
+    std::shared_ptr<Functional> functional_;
+};
+#endif

dune/tectonic/spatial-solving/makelinearsolver.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_MAKELINEARSOLVER_HH
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "preconditioners/multilevelpatchpreconditioner.hh"
+
+template <class ContactNetwork, class VectorType>
+auto makeLinearSolver(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+    // make linear solver for linear correction in TNNMGStep
+    using Vector = VectorType;
+    using Matrix = typename ContactNetwork::Matrix;
+
+        /* old, pre multi threading, was unused
+        using Norm =  EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        const auto& preconditionerParset = parset_.sub("solver.tnnmg.linear.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        return std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.linear.tolerance"), norm, Solver::QUIET);
+        */
+
+        // patch preconditioner only needs to be computed once per advance()
+        // make linear solver for linear correction in TNNMGStep
+        using Norm = EnergyNorm<Matrix, Vector>;
+        using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        /*const auto& preconditionerParset = parset_.sub("solver.tnnmg.preconditioner");
+
+        Dune::BitSetVector<1> activeLevels(contactNetwork_.nLevels(), true);
+        Preconditioner preconditioner(preconditionerParset, contactNetwork_, activeLevels);
+        preconditioner.setPatchDepth(preconditionerParset.template get<size_t>("patchDepth"));
+        preconditioner.build();
+
+        auto cgStep = std::make_shared<Dune::Solvers::CGStep<Matrix, Vector>>();
+        cgStep->setPreconditioner(preconditioner);
+
+        Norm norm(*cgStep);
+
+        auto linearSolver = std::make_shared<LinearSolver>(cgStep, parset_.template get<int>("solver.tnnmg.main.multi"), parset_.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+        */
+
+        // set multigrid solver
+        auto smoother = TruncatedBlockGSStep<Matrix, Vector>();
+
+        // transfer operators need to be recomputed on change due to setDeformation()
+        using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+        TransferOperators transfer(contactNetwork.nLevels()-1);
+        for (size_t i=0; i<transfer.size(); i++) {
+            transfer[i] = std::make_shared<TransferOperator>();
+            transfer[i]->setup(contactNetwork, i, i+1);
+        }
+
+        // Remove any recompute filed so that initially the full transferoperator is assembled
+        for (size_t i=0; i<transfer.size(); i++)
+            std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+        linearMultigridStep->setMGType(1, 3, 3);
+        linearMultigridStep->setSmoother(smoother);
+        linearMultigridStep->setTransferOperators(transfer);
+
+        Norm norm(*linearMultigridStep);
+
+        return std::make_shared<LinearSolver>(linearMultigridStep, parset.template get<int>("solver.tnnmg.main.multi"), parset.template get<double>("solver.tnnmg.preconditioner.basesolver.tolerance"), norm, Solver::QUIET);
+    }
+#endif

dune/tectonic/spatial-solving/nonlinearsolver.hh

+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_SOLVENONLINEARPROBLEM_HH
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
+#include "../data-structures/enums.hh"
+#include "../data-structures/enumparser.hh"
+
+
+
+#include "../utils/tobool.hh"
+#include "../utils/debugutils.hh"
+
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/contact/nbodycontacttransfer.hh"
+
+#include "solverfactory.hh"
+#include "solverfactory.cc"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "fixedpointiterator.hh"
+
+
+template <class Functional, class BitVector>
+class NonlinearSolver {
+protected:
+    using Matrix = typename Functional::Matrix;
+    using Vector = typename Functional::Vector;
+    using Factory = SolverFactory<Functional, BitVector>;
+
+    using Norm = EnergyNorm<Matrix, Vector>;
+    using SolverType = Dune::Solvers::LoopSolver<Vector>;
+
+    const static int dims = Vector::block_type::dimension;
+
+public:
+    template <class LinearSolver>
+    NonlinearSolver(const Dune::ParameterTree& tnnmgParset, std::shared_ptr<LinearSolver> linearSolver, std::shared_ptr<Functional> functional, const BitVector& dirichletNodes) :
+        functional_(functional),
+        norm_(functional_->quadraticPart()) {
+
+        // set up TNMMG step
+        solverFactory_ = std::make_shared<Factory>(tnnmgParset, functional_, dirichletNodes);
+        solverFactory_->build(linearSolver);
+    }
+
+    auto solve(const Dune::ParameterTree& solverParset, Vector& x) {
+        auto tnnmgStep = solverFactory_->step();
+
+        SolverType solver(*tnnmgStep.get(), solverParset.get<size_t>("maximumIterations"),
+            solverParset.get<double>("tolerance"), norm_,
+            solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+        const auto& lower = functional_->lowerObstacle();
+        const auto& upper = functional_->upperObstacle();
+
+        // project in onto admissible set
+        for (size_t i=0; i<x.size(); i++) {
+            for (size_t j=0; j<dims; j++) {
+                if (x[i][j] < lower[i][j]) {
+                    x[i][j] = lower[i][j];
+                }
+
+                if (x[i][j] > upper[i][j]) {
+                    x[i][j] = upper[i][j];
+                }
+            }
+        }
+
+        solverFactory_->setProblem(x);
+
+        solver.preprocess();
+        solver.solve();
+
+        return solver.getResult();
+    }
+
+private:
+    std::shared_ptr<Functional> functional_;
+    std::shared_ptr<Factory> solverFactory_;
+
+    Norm norm_;
+};
+
+#endif