stokes_solver.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_LINEAR_STOKES_SOLVER_HH
#define DUMUX_LINEAR_STOKES_SOLVER_HH
 
#include <type_traits>
#include <memory>
#include <tuple>
 
#include <dune/common/parametertree.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/common/exceptions.hh>
 
#include <dune/istl/matrixindexset.hh>
#include <dune/istl/preconditioner.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh>
 
#include <dumux/common/math.hh>
#include <dumux/common/exceptions.hh>
#include <dumux/linear/solver.hh>
#include <dumux/linear/preconditioners.hh>
#include <dumux/linear/linearsolverparameters.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/linearsolvertraits.hh>
#include <dumux/linear/parallelmatrixadapter.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/linear/symmetrize_constraints.hh>
#include <dumux/assembly/jacobianpattern.hh>
 
namespace Dumux::Detail {
 
template<class M, class X, class Y, int l = 2>
class StokesPreconditioner : public Dune::Preconditioner<X,Y>
{
    static_assert(Dumux::isMultiTypeBlockMatrix<M>::value && M::M() == 2 && M::N() == 2, "Expects a 2x2 MultiTypeBlockMatrix.");
    static_assert(l == 2, "StokesPreconditioner expects a block level of 2 for Matrix type M.");
 
    using A = std::decay_t<decltype(std::declval<M>()[Dune::Indices::_0][Dune::Indices::_0])>;
    using U = std::decay_t<decltype(std::declval<X>()[Dune::Indices::_0])>;
 
    using P = std::decay_t<decltype(std::declval<M>()[Dune::Indices::_1][Dune::Indices::_1])>;
    using V = std::decay_t<decltype(std::declval<X>()[Dune::Indices::_1])>;
 
    enum class Mode { symmetric, triangular,  diagonal };
 
public:
    using matrix_type = M;
    using domain_type = X;
    using range_type = Y;
    using field_type = typename X::field_type;
    using scalar_field_type = Dune::Simd::Scalar<field_type>;
    using PressureLinearOperator = Dune::MatrixAdapter<P,V,V>;
 
    StokesPreconditioner(
        const std::shared_ptr<const Dune::AssembledLinearOperator<M,X,Y>>& fullLinearOperator,
        const std::shared_ptr<const PressureLinearOperator>& pressureLinearOperator,
        const Dune::ParameterTree& params
    )
    : matrix_(fullLinearOperator->getmat())
    , pmatrix_(pressureLinearOperator->getmat())
    , verbosity_(params.get<int>("verbosity"))
    , paramGroup_(params.get<std::string>("ParameterGroup"))
    {
        const auto mode = getParamFromGroup<std::string>(
            paramGroup_, "LinearSolver.Preconditioner.Mode", "Diagonal"
        );
 
        if (mode == "Symmetric")
            mode_ = Mode::symmetric;
        else if (mode == "Triangular")
            mode_ = Mode::triangular;
        else
            mode_ = Mode::diagonal;
 
        initPreconditioner_(params);
    }
 
    void pre(X& update, Y& currentDefect) override {}
 
    void apply(X& update, const Y& currentDefect) override
    {
        using namespace Dune::Indices;
 
        if (mode_ == Mode::symmetric)
        {
            update = applyRightBlock_(currentDefect);
            update = applyDiagBlock_(update);
            update = applyLeftBlock_(update);
        }
        else if (mode_ == Mode::triangular)
        {
            update = applyRightBlock_(currentDefect);
            update = applyDiagBlock_(update);
        }
        else
            update = applyDiagBlock_(currentDefect);
    }
 
    void post(X& update) override {}
 
    Dune::SolverCategory::Category category() const override
    {
        return Dune::SolverCategory::sequential;
    }
 
private:
    X applyRightBlock_(const Y& d)
    {
        using namespace Dune::Indices;
 
        // right bit of LDU decomposition
        // applied to d
        //
        // | I         0 |
        // | -C*inv(A) I |
        //
 
        auto dTmp0 = d[_0];
        auto vTmp = d; vTmp = 0.0;
 
        // invert velocity block (apply to first block of d)
        applyPreconditionerForA_(vTmp[_0], dTmp0);
 
        // then multiply with C
        matrix_[_1][_0].mv(vTmp[_0], vTmp[_1]);
 
        // and subtract from d
        auto v = d;
        v[_0] = vTmp[_0]; // already do A^-1 d of the diagonal block because we already computed it here
        v[_1] -= vTmp[_1];
        return v;
    }
 
    X applyDiagBlock_(const Y& d)
    {
        using namespace Dune::Indices;
 
        // diagonal middle bit of LDU decomposition
        auto dTmp = d;
        auto v = d; v = 0.0;
 
        // invert velocity block
        if (mode_ == Mode::diagonal)
            applyPreconditionerForA_(v[_0], dTmp[_0]);
 
        // reuse the already computed A^-1 d (see applyRightBlock_)
        else
            v[_0] = dTmp[_0];
 
        // invert pressure block
        applyPreconditionerForP_(v[_1], dTmp[_1]);
 
        return v;
    }
 
    X applyLeftBlock_(const Y& d)
    {
        using namespace Dune::Indices;
 
        // left bit of LDU decomposition
        // applied to d
        //
        // | I  -inv(A)*B |
        // | 0       I    |
        //
 
        auto dTmp = d;
        auto vTmp = d; vTmp = 0.0;
 
        // multiply with B
        matrix_[_0][_1].umv(dTmp[_1], vTmp[_0]);
 
        // invert velocity block (apply to first block of d)
        auto vTmp0 = vTmp[_0]; vTmp0 = 0.0;
        applyPreconditionerForA_(vTmp0, vTmp[_0]);
 
        // and subtract from d
        auto v = d;
        v[_0] -= vTmp0;
 
        return v;
    }
 
    void initPreconditioner_(const Dune::ParameterTree& params)
    {
        using namespace Dune::Indices;
 
        if (getParamFromGroup<bool>(paramGroup_, "LinearSolver.DirectSolverForVelocity", false))
        {
#if HAVE_UMFPACK
            directSolverForA_ = std::make_shared<Dune::UMFPack<A>>(matrix_[_0][_0], verbosity_);
            using Wrap = Dune::InverseOperator2Preconditioner<Dune::InverseOperator<U, U>>;
            preconditionerForA_ = std::make_shared<Wrap>(*directSolverForA_);
#else
            DUNE_THROW(Dune::InvalidStateException, "Selected direct solver but UMFPack is not available.");
#endif
        }
        else
        {
            using VelLinearOperator = Dune::MatrixAdapter<A, U, U>;
            auto lopV = std::make_shared<VelLinearOperator>(matrix_[_0][_0]);
            preconditionerForA_ = std::make_shared<
                Dune::Amg::AMG<VelLinearOperator, U, Dumux::ParMTSSOR<A,U,U>>
            >(lopV, params);
        }
 
        if (getParamFromGroup<bool>(paramGroup_, "LinearSolver.DirectSolverForPressure", false))
        {
#if HAVE_UMFPACK
            directSolverForP_ = std::make_shared<Dune::UMFPack<P>>(pmatrix_, verbosity_);
            using Wrap = Dune::InverseOperator2Preconditioner<Dune::InverseOperator<V, V>>;
            preconditionerForP_ = std::make_shared<Wrap>(*directSolverForP_);
#else
            DUNE_THROW(Dune::InvalidStateException, "Selected direct solver but UMFPack is not available.");
#endif
        }
        else
        {
            using PressJacobi = Dumux::ParMTJac<P, V, V>;
            std::size_t numIterations = pmatrix_.nonzeroes() == pmatrix_.N() ? 1 : 10;
            preconditionerForP_ = std::make_shared<PressJacobi>(pmatrix_, numIterations, 1.0);
        }
    }
 
    template<class Sol, class Rhs>
    void applyPreconditionerForA_(Sol& sol, Rhs& rhs) const
    {
        preconditionerForA_->pre(sol, rhs);
        preconditionerForA_->apply(sol, rhs);
        preconditionerForA_->post(sol);
    }
 
    template<class Sol, class Rhs>
    void applyPreconditionerForP_(Sol& sol, Rhs& rhs) const
    {
        preconditionerForP_->pre(sol, rhs);
        preconditionerForP_->apply(sol, rhs);
        preconditionerForP_->post(sol);
    }
 
    const M& matrix_;
    const P& pmatrix_;
    const int verbosity_;
 
    std::shared_ptr<Dune::Preconditioner<U, U>> preconditionerForA_;
    std::shared_ptr<Dune::Preconditioner<V, V>> preconditionerForP_;
    std::shared_ptr<Dune::InverseOperator<U, U>> directSolverForA_;
    std::shared_ptr<Dune::InverseOperator<V, V>> directSolverForP_;
 
    const std::string paramGroup_;
    Mode mode_;
};
 
} // end namespace Detail
 
namespace Dumux {
 
template<class Matrix, class Vector, class VelocityGG, class PressureGG>
class StokesSolver
: public LinearSolver
{
    using Preconditioner = Detail::StokesPreconditioner<Matrix, Vector, Vector>;
public:
    StokesSolver(std::shared_ptr<const VelocityGG> vGridGeometry,
                 std::shared_ptr<const PressureGG> pGridGeometry,
                 const Vector& dirichletDofs,
                 const std::string& paramGroup = "")
    : LinearSolver(paramGroup)
    , vGridGeometry_(std::move(vGridGeometry))
    , pGridGeometry_(std::move(pGridGeometry))
    , dirichletDofs_(dirichletDofs)
    {
        params_ = LinearSolverParameters<LinearSolverTraits<VelocityGG>>::createParameterTree(this->paramGroup());
        density_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidDensity");
 
        if (hasParamInGroup(this->paramGroup(), "Component.LiquidDynamicViscosity"))
            viscosity_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidDynamicViscosity");
        else if (hasParamInGroup(this->paramGroup(), "Component.LiquidKinematicViscosity"))
            viscosity_ = getParamFromGroup<double>(this->paramGroup(), "Component.LiquidKinematicViscosity") * density_;
        else
        {
            const std::string group = this->paramGroup() == "" ? "Component" : this->paramGroup() + ".Component";
            DUNE_THROW(ParameterException, "Stokes solver requires the parameters"
                           << " LiquidDynamicViscosity or LiquidKinematicViscosity"
                           << " in parameter group [" << group << "]");
        }
 
        weight_ = getParamFromGroup<double>(this->paramGroup(), "LinearSolver.Preconditioner.MassMatrixWeight", 1.0);
        solverType_ = getParamFromGroup<std::string>(this->paramGroup(), "LinearSolver.Type", "gmres");
        scalarProduct_ = std::make_shared<Dune::ScalarProduct<Vector>>();
    }
 
    bool solve(const Matrix& A, Vector& x, const Vector& b)
    {
        auto bTmp = b;
        auto ATmp = A;
 
        return applyIterativeSolver_(ATmp, x, bTmp);
    }
 
    Scalar norm(const Vector& b) const
    {
        return scalarProduct_->norm(b);
    }
 
    std::string name() const
    {
        return "Block-preconditioned Stokes solver";
    }
 
    const Dune::InverseOperatorResult& result() const
    {
        return result_;
    }
 
private:
    bool applyIterativeSolver_(Matrix& A, Vector& x, Vector& b)
    {
        // make Dirichlet boundary conditions symmetric
        if (getParamFromGroup<bool>(this->paramGroup(), "LinearSolver.SymmetrizeDirichlet", true))
            symmetrizeConstraints(A, b, dirichletDofs_);
 
        // make Matrix symmetric on the block-scale
        using namespace Dune::Indices;
        A[_1] *= -1.0/density_;
        b[_1] *= -1.0/density_;
 
        auto op = std::make_shared<Dumux::ParallelMultiTypeMatrixAdapter<Matrix, Vector, Vector>>(A);
        auto pop = makePressureLinearOperator_<typename Preconditioner::PressureLinearOperator>();
        auto preconditioner = std::make_shared<Preconditioner>(op, pop, params_.sub("preconditioner"));
        params_["verbose"] = pGridGeometry_->gridView().comm().rank() == 0 ? params_["verbose"] : "0";
 
        // defaults to restarted GMRes
        std::unique_ptr<Dune::InverseOperator<Vector, Vector>> solver;
        if (solverType_ == "minres")
            solver = std::make_unique<Dune::MINRESSolver<Vector>>(op, scalarProduct_, preconditioner, params_);
        else if (solverType_ == "bicgstab")
            solver = std::make_unique<Dune::BiCGSTABSolver<Vector>>(op, scalarProduct_, preconditioner, params_);
        else if (solverType_ == "gmres")
            solver = std::make_unique<Dune::RestartedGMResSolver<Vector>>(op, scalarProduct_, preconditioner, params_);
        else
            DUNE_THROW(Dune::NotImplemented, "Solver choice " << solverType_ << " is not implemented");
 
        solver->apply(x, b, result_);
 
        return result_.converged;
    }
 
    template<class LinearOperator>
    std::shared_ptr<LinearOperator> makePressureLinearOperator_()
    {
        using M = typename LinearOperator::matrix_type;
        auto massMatrix = createMassMatrix_<M>();
        return std::make_shared<LinearOperator>(massMatrix);
    }
 
    template<class M, class DM = typename PressureGG::DiscretizationMethod>
    std::shared_ptr<M> createMassMatrix_()
    {
        auto massMatrix = std::make_shared<M>();
        massMatrix->setBuildMode(M::random);
        const auto numDofs = pGridGeometry_->numDofs();
 
        if constexpr (DM{} == DiscretizationMethods::cctpfa || DM{} == DiscretizationMethods::ccmpfa)
        {
            Dune::MatrixIndexSet pattern;
            pattern.resize(numDofs, numDofs);
            for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
                pattern.add(globalI, globalI);
            pattern.exportIdx(*massMatrix);
 
            const auto& gv = pGridGeometry_->gridView();
            auto fvGeometry = localView(*pGridGeometry_);
            for (const auto& element : elements(gv))
            {
                fvGeometry.bindElement(element);
                for (const auto& scv : scvs(fvGeometry))
                {
                    using Extrusion = Extrusion_t<PressureGG>;
                    const auto dofIndex = scv.dofIndex();
                    if (element.partitionType() == Dune::GhostEntity) // do not modify ghosts
                        (*massMatrix)[dofIndex][dofIndex] = 1.0;
                    else
                        (*massMatrix)[dofIndex][dofIndex] += weight_*Extrusion::volume(fvGeometry, scv)/(2.0*viscosity_);
                }
            }
 
            return massMatrix;
        }
        else if constexpr (DM{} == DiscretizationMethods::box || DM{} == DiscretizationMethods::fcdiamond)
        {
            getJacobianPattern<true>(*pGridGeometry_).exportIdx(*massMatrix);
            (*massMatrix) = 0.0;
 
            const auto& gv = pGridGeometry_->gridView();
            auto fvGeometry = localView(*pGridGeometry_);
            std::vector<Dune::FieldVector<double, 1>> values;
 
            for (const auto& element : elements(gv))
            {
                const auto geometry = element.geometry();
                const auto& localBasis = pGridGeometry_->feCache().get(element.type()).localBasis();
                const auto numLocalDofs = localBasis.size();
                values.resize(numLocalDofs);
 
                fvGeometry.bindElement(element);
                for (const auto& scvJ : scvs(fvGeometry))
                {
                    // Use mid-point rule (only works for linear ansatz functions)
                    const auto globalJ = scvJ.dofIndex();
                    const auto qWeightJ = PressureGG::Extrusion::volume(fvGeometry, scvJ);
                    const auto quadPos = geometry.local(scvJ.center());
                    localBasis.evaluateFunction(quadPos, values);
 
                    for (const auto& scvI : scvs(fvGeometry))
                    {
                        const auto valueIJ = values[scvI.localDofIndex()]*qWeightJ/(2.0*viscosity_);
                        (*massMatrix)[scvI.dofIndex()][globalJ][0][0] += valueIJ;
                    }
                }
            }
 
            return massMatrix;
        }
 
        DUNE_THROW(Dune::NotImplemented, "Mass matrix for discretization method not implemented");
    }
 
    double density_, viscosity_, weight_;
    Dune::InverseOperatorResult result_;
    Dune::ParameterTree params_;
    std::shared_ptr<const VelocityGG> vGridGeometry_;
    std::shared_ptr<const PressureGG> pGridGeometry_;
    const Vector& dirichletDofs_;
    std::string solverType_;
 
    std::shared_ptr<Dune::ScalarProduct<Vector>> scalarProduct_;
};
 
} // end namespace Dumux
 
#endif