multidomain/fvassembler.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:
/*****************************************************************************
 *   See the file COPYING for full copying permissions.                      *
 *                                                                           *
 *   This program is free software: you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation, either version 3 of the License, or       *
 *   (at your option) any later version.                                     *
 *                                                                           *
 *   This program is distributed in the hope that it will be useful,         *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
 *   GNU General Public License for more details.                            *
 *                                                                           *
 *   You should have received a copy of the GNU General Public License       *
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
 *****************************************************************************/
#ifndef DUMUX_MULTIDOMAIN_FV_ASSEMBLER_HH
#define DUMUX_MULTIDOMAIN_FV_ASSEMBLER_HH
 
#include <type_traits>
#include <tuple>
 
#include <dune/common/hybridutilities.hh>
#include <dune/istl/matrixindexset.hh>
 
#include <dumux/common/exceptions.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/typetraits/utility.hh>
#include <dumux/common/gridcapabilities.hh>
#include <dumux/discretization/method.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/assembly/jacobianpattern.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/parallel/multithreading.hh>
 
#include "couplingjacobianpattern.hh"
#include "subdomaincclocalassembler.hh"
#include "subdomainboxlocalassembler.hh"
#include "subdomainstaggeredlocalassembler.hh"
#include "subdomainfclocalassembler.hh"
#include "subdomainfcdiamondlocalassembler.hh"
#include "subdomainpq1bubblelocalassembler.hh"
 
#include <dumux/discretization/method.hh>
 
namespace Dumux {
 
namespace Grid::Capabilities {
 
namespace Detail {
// helper for multi-domain models
template<class T, std::size_t... I>
bool allGridsSupportsMultithreadingImpl(const T& gridGeometries, std::index_sequence<I...>)
{
    return (... && supportsMultithreading(std::get<I>(gridGeometries)->gridView()));
}
} // end namespace Detail
 
// helper for multi-domain models (all grids have to support multithreading)
template<class... GG>
bool allGridsSupportsMultithreading(const std::tuple<GG...>& gridGeometries)
{
    return Detail::allGridsSupportsMultithreadingImpl<std::tuple<GG...>>(gridGeometries, std::make_index_sequence<sizeof...(GG)>());
}
 
} // end namespace Grid::Capabilities

template<class CM>
struct CouplingManagerSupportsMultithreadedAssembly : public std::false_type
{};

template<class MDTraits, class CMType, DiffMethod diffMethod, bool useImplicitAssembly = true>
class MultiDomainFVAssembler
{
    template<std::size_t id>
    using SubDomainTypeTag = typename MDTraits::template SubDomain<id>::TypeTag;
 
public:
    using Traits = MDTraits;
 
    using Scalar = typename MDTraits::Scalar;

    template<std::size_t id>
    using LocalResidual = GetPropType<SubDomainTypeTag<id>, Properties::LocalResidual>;
 
    template<std::size_t id>
    using GridVariables = typename MDTraits::template SubDomain<id>::GridVariables;
 
    template<std::size_t id>
    using GridGeometry = typename MDTraits::template SubDomain<id>::GridGeometry;
 
    template<std::size_t id>
    using Problem = typename MDTraits::template SubDomain<id>::Problem;
 
    using JacobianMatrix = typename MDTraits::JacobianMatrix;
    using SolutionVector = typename MDTraits::SolutionVector;
    using ResidualType = SolutionVector;
 
    using CouplingManager = CMType;

    static constexpr bool isImplicit()
    { return useImplicitAssembly; }
 
private:
 
    using ProblemTuple = typename MDTraits::template TupleOfSharedPtrConst<Problem>;
    using GridGeometryTuple = typename MDTraits::template TupleOfSharedPtrConst<GridGeometry>;
    using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
 
    using TimeLoop = TimeLoopBase<Scalar>;
    using ThisType = MultiDomainFVAssembler<MDTraits, CouplingManager, diffMethod, isImplicit()>;
 
    template<class DiscretizationMethod, std::size_t id>
    struct SubDomainAssemblerType;
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::CCTpfa, id>
    {
        using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::CCMpfa, id>
    {
        using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::Box, id>
    {
        using type = SubDomainBoxLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::Staggered, id>
    {
        using type = SubDomainStaggeredLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::FCStaggered, id>
    {
        using type = SubDomainFaceCenteredLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::FCDiamond, id>
    {
        using type = SubDomainFaceCenteredDiamondLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::PQ1Bubble, id>
    {
        using type = SubDomainPQ1BubbleLocalAssembler<id, SubDomainTypeTag<id>, ThisType, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    using SubDomainAssembler = typename SubDomainAssemblerType<typename GridGeometry<id>::DiscretizationMethod, id>::type;
 
public:
 

    MultiDomainFVAssembler(ProblemTuple problem,
                           GridGeometryTuple gridGeometry,
                           GridVariablesTuple gridVariables,
                           std::shared_ptr<CouplingManager> couplingManager)
    : couplingManager_(couplingManager)
    , problemTuple_(std::move(problem))
    , gridGeometryTuple_(std::move(gridGeometry))
    , gridVariablesTuple_(std::move(gridVariables))
    , timeLoop_()
    , isStationaryProblem_(true)
    , warningIssued_(false)
    {
        static_assert(isImplicit(), "Explicit assembler for stationary problem doesn't make sense!");
        std::cout << "Instantiated assembler for a stationary problem." << std::endl;
 
        enableMultithreading_ = CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value
            && Grid::Capabilities::allGridsSupportsMultithreading(gridGeometryTuple_)
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }

    MultiDomainFVAssembler(ProblemTuple problem,
                           GridGeometryTuple gridGeometry,
                           GridVariablesTuple gridVariables,
                           std::shared_ptr<CouplingManager> couplingManager,
                           std::shared_ptr<const TimeLoop> timeLoop,
                           const SolutionVector& prevSol)
    : couplingManager_(couplingManager)
    , problemTuple_(std::move(problem))
    , gridGeometryTuple_(std::move(gridGeometry))
    , gridVariablesTuple_(std::move(gridVariables))
    , timeLoop_(timeLoop)
    , prevSol_(&prevSol)
    , isStationaryProblem_(false)
    , warningIssued_(false)
    {
        std::cout << "Instantiated assembler for an instationary problem." << std::endl;
 
        enableMultithreading_ = CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value
            && Grid::Capabilities::allGridsSupportsMultithreading(gridGeometryTuple_)
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }

    void assembleJacobianAndResidual(const SolutionVector& curSol)
    {
        checkAssemblerState_();
        resetJacobian_();
        resetResidual_();
 
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainId)
        {
            auto& jacRow = (*jacobian_)[domainId];
            auto& subRes = (*residual_)[domainId];
            this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
 
            const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
            enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridGeometry, curSol[domainId]);
        });
    }

    void assembleResidual(const SolutionVector& curSol)
    {
        resetResidual_();
        assembleResidual(*residual_, curSol);
    }

    void assembleResidual(ResidualType& r, const SolutionVector& curSol)
    {
        r = 0.0;
 
        checkAssemblerState_();
 
        // update the grid variables for the case of active caching
        updateGridVariables(curSol);
 
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(r)), [&](const auto domainId)
        {
            auto& subRes = r[domainId];
            this->assembleResidual_(domainId, subRes, curSol);
        });
    }

    Scalar residualNorm(const SolutionVector& curSol)
    {
        ResidualType residual;
        setResidualSize_(residual);
        assembleResidual(residual, curSol);
 
        // calculate the squared norm of the residual
        Scalar resultSquared = 0.0;
 
        // for box communicate the residual with the neighboring processes
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(residual)), [&](const auto domainId)
        {
            const auto gridGeometry = std::get<domainId>(gridGeometryTuple_);
            const auto& gridView = gridGeometry->gridView();
 
            if (gridView.comm().size() > 1 && gridView.overlapSize(0) == 0)
            {
                if constexpr (GridGeometry<domainId>::discMethod == DiscretizationMethods::box)
                {
                    using GV = typename GridGeometry<domainId>::GridView;
                    using DM = typename GridGeometry<domainId>::VertexMapper;
                    using PVHelper = ParallelVectorHelper<GV, DM, GV::dimension>;
 
                    PVHelper vectorHelper(gridView, gridGeometry->vertexMapper());
 
                    vectorHelper.makeNonOverlappingConsistent(residual[domainId]);
                }
            }
            else if (!warningIssued_)
            {
                if (gridView.comm().size() > 1 && gridView.comm().rank() == 0)
                    std::cout << "\nWarning: norm calculation adds entries corresponding to\n"
                              << "overlapping entities multiple times. Please use the norm\n"
                              << "function provided by a linear solver instead." << std::endl;
 
                warningIssued_ = true;
            }
 
            Scalar localNormSquared = residual[domainId].two_norm2();
 
            if (gridView.comm().size() > 1)
            {
                localNormSquared = gridView.comm().sum(localNormSquared);
            }
 
            resultSquared += localNormSquared;
        });
 
        using std::sqrt;
        return sqrt(resultSquared);
    }

    void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
                         std::shared_ptr<SolutionVector> r)
    {
        jacobian_ = A;
        residual_ = r;
 
        setJacobianBuildMode(*jacobian_);
        setJacobianPattern_(*jacobian_);
        setResidualSize_(*residual_);
    }

    void setLinearSystem()
    {
        jacobian_ = std::make_shared<JacobianMatrix>();
        residual_ = std::make_shared<SolutionVector>();
 
        setJacobianBuildMode(*jacobian_);
        setJacobianPattern_(*jacobian_);
        setResidualSize_(*residual_);
    }

    void setJacobianBuildMode(JacobianMatrix& jac) const
    {
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto i)
        {
            forEach(jac[i], [&](auto& jacBlock)
            {
                using BlockType = std::decay_t<decltype(jacBlock)>;
                if (jacBlock.buildMode() == BlockType::BuildMode::unknown)
                    jacBlock.setBuildMode(BlockType::BuildMode::random);
                else if (jacBlock.buildMode() != BlockType::BuildMode::random)
                    DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
            });
        });
    }

    void updateAfterGridAdaption()
    {
        setResidualSize_();
        setJacobianPattern_();
        maybeComputeColors_();
    }

    void updateGridVariables(const SolutionVector& curSol)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
        { this->gridVariables(domainId).update(curSol[domainId]); });
    }

    void resetTimeStep(const SolutionVector& curSol)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(gridVariablesTuple_)), [&](const auto domainId)
        { this->gridVariables(domainId).resetTimeStep(curSol[domainId]); });
    }

    template<std::size_t i>
    std::size_t numDofs(Dune::index_constant<i> domainId) const
    { return std::get<domainId>(gridGeometryTuple_)->numDofs(); }

    template<std::size_t i>
    const auto& problem(Dune::index_constant<i> domainId) const
    { return *std::get<domainId>(problemTuple_); }

    template<std::size_t i>
    const auto& gridGeometry(Dune::index_constant<i> domainId) const
    { return *std::get<domainId>(gridGeometryTuple_); }

    template<std::size_t i>
    const auto& gridView(Dune::index_constant<i> domainId) const
    { return gridGeometry(domainId).gridView(); }

    template<std::size_t i>
    GridVariables<i>& gridVariables(Dune::index_constant<i> domainId)
    { return *std::get<domainId>(gridVariablesTuple_); }

    template<std::size_t i>
    const GridVariables<i>& gridVariables(Dune::index_constant<i> domainId) const
    { return *std::get<domainId>(gridVariablesTuple_); }

    const CouplingManager& couplingManager() const
    { return *couplingManager_; }

    JacobianMatrix& jacobian()
    { return *jacobian_; }

    SolutionVector& residual()
    { return *residual_; }

    const SolutionVector& prevSol() const
    { return *prevSol_; }

    void setTimeManager(std::shared_ptr<const TimeLoop> timeLoop)
    { timeLoop_ = timeLoop; isStationaryProblem_ = !(static_cast<bool>(timeLoop)); }

    void setPreviousSolution(const SolutionVector& u)
    { prevSol_ = &u; }

    bool isStationaryProblem() const
    { return isStationaryProblem_; }

    template<std::size_t i>
    LocalResidual<i> localResidual(Dune::index_constant<i> domainId) const
    { return LocalResidual<i>(std::get<domainId>(problemTuple_).get(), timeLoop_.get()); }
 
protected:
    std::shared_ptr<CouplingManager> couplingManager_;
 
private:
    void setJacobianPattern_(JacobianMatrix& jac) const
    {
        using namespace Dune::Hybrid;
        forEach(std::make_index_sequence<JacobianMatrix::N()>(), [&](const auto domainI)
        {
            forEach(integralRange(Dune::Hybrid::size(jac[domainI])), [&](const auto domainJ)
            {
                const auto pattern = this->getJacobianPattern_(domainI, domainJ);
                pattern.exportIdx(jac[domainI][domainJ]);
            });
        });
    }

    void setResidualSize_(SolutionVector& res) const
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(res)), [&](const auto domainId)
        { res[domainId].resize(this->numDofs(domainId)); });
    }
 
    // reset the residual vector to 0.0
    void resetResidual_()
    {
        if(!residual_)
        {
            residual_ = std::make_shared<SolutionVector>();
            setResidualSize_(*residual_);
        }
 
        (*residual_) = 0.0;
    }
 
    // reset the jacobian vector to 0.0
    void resetJacobian_()
    {
        if(!jacobian_)
        {
            jacobian_ = std::make_shared<JacobianMatrix>();
            setJacobianBuildMode(*jacobian_);
            setJacobianPattern_(*jacobian_);
        }
 
       (*jacobian_)  = 0.0;
    }

    void maybeComputeColors_()
    {
        if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
            if (enableMultithreading_)
                couplingManager_->computeColorsForAssembly();
    }
 
    // check if the assembler is in a correct state for assembly
    void checkAssemblerState_() const
    {
        if (!isStationaryProblem_ && !prevSol_)
            DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
 
        if (isStationaryProblem_ && prevSol_)
            DUNE_THROW(Dune::InvalidStateException, "Assembling stationary problem but a previous solution was set."
                                                    << " Did you forget to set the timeLoop to make this problem instationary?");
    }
 
    template<std::size_t i, class JacRow, class SubRes>
    void assembleJacobianAndResidual_(Dune::index_constant<i> domainId, JacRow& jacRow, SubRes& subRes,
                                      const SolutionVector& curSol)
    {
        assemble_(domainId, [&](const auto& element)
        {
            SubDomainAssembler<i> subDomainAssembler(*this, element, curSol, *couplingManager_);
            subDomainAssembler.assembleJacobianAndResidual(jacRow, subRes, gridVariablesTuple_);
        });
    }
 
    template<std::size_t i, class SubRes>
    void assembleResidual_(Dune::index_constant<i> domainId, SubRes& subRes,
                           const SolutionVector& curSol)
    {
        assemble_(domainId, [&](const auto& element)
        {
            SubDomainAssembler<i> subDomainAssembler(*this, element, curSol, *couplingManager_);
            subDomainAssembler.assembleResidual(subRes);
        });
    }

    template<std::size_t i, class AssembleElementFunc>
    void assemble_(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
    {
        // a state that will be checked on all processes
        bool succeeded = false;
 
        // try assembling using the local assembly function
        try
        {
            if constexpr (CouplingManagerSupportsMultithreadedAssembly<CouplingManager>::value)
            {
                if (enableMultithreading_)
                {
                    couplingManager_->assembleMultithreaded(
                        domainId, std::forward<AssembleElementFunc>(assembleElement)
                    );
                    return;
                }
            }
 
            // fallback for coupling managers that don't support multithreaded assembly (yet)
            // or if multithreaded assembly is disabled
            // let the local assembler add the element contributions
            for (const auto& element : elements(gridView(domainId)))
                assembleElement(element);
 
            // if we get here, everything worked well on this process
            succeeded = true;
        }
        // throw exception if a problem occurred
        catch (NumericalProblem &e)
        {
            std::cout << "rank " << gridView(domainId).comm().rank()
                      << " caught an exception while assembling:" << e.what()
                      << "\n";
            succeeded = false;
        }
 
        // make sure everything worked well on all processes
        if (gridView(domainId).comm().size() > 1)
            succeeded = gridView(domainId).comm().min(succeeded);
 
        // if not succeeded rethrow the error on all processes
        if (!succeeded)
            DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
    }
 
    // get diagonal block pattern
    template<std::size_t i, std::size_t j, typename std::enable_if_t<(i==j), int> = 0>
    Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
                                             Dune::index_constant<j> domainJ) const
    {
        const auto& gg = gridGeometry(domainI);
        auto pattern = getJacobianPattern<isImplicit()>(gg);
        couplingManager_->extendJacobianPattern(domainI, pattern);
        return pattern;
    }
 
    // get coupling block pattern
    template<std::size_t i, std::size_t j, typename std::enable_if_t<(i!=j), int> = 0>
    Dune::MatrixIndexSet getJacobianPattern_(Dune::index_constant<i> domainI,
                                             Dune::index_constant<j> domainJ) const
    {
        return getCouplingJacobianPattern<isImplicit()>(*couplingManager_,
                                                        domainI, gridGeometry(domainI),
                                                        domainJ, gridGeometry(domainJ));
    }
 
    // build periodic constraints into the system matrix
    template<std::size_t i, class JacRow, class Sol, class GG>
    void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Sol& res, const GG& gridGeometry, const Sol& curSol)
    {
        if constexpr (GG::discMethod == DiscretizationMethods::box || GG::discMethod == DiscretizationMethods::fcstaggered)
        {
            for (const auto& m : gridGeometry.periodicVertexMap())
            {
                if (m.first < m.second)
                {
                    auto& jac = jacRow[domainI];
 
                    // add the second row to the first
                    res[m.first] += res[m.second];
 
                    // enforce the solution of the first periodic DOF to the second one
                    res[m.second] = curSol[m.second] - curSol[m.first];
 
                    const auto end = jac[m.second].end();
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                        jac[m.first][it.index()] += (*it);
 
                    // enforce constraint in second row
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                        (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
 
                    using namespace Dune::Hybrid;
                    forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
                    {
                        auto& jacCoupling = jacRow[couplingDomainId];
 
                        for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
                            jacCoupling[m.first][it.index()] += (*it);
 
                        for (auto it = jacCoupling[m.second].begin(); it != jacCoupling[m.second].end(); ++it)
                            (*it) = 0.0;
                    });
                }
            }
        }
    }

    ProblemTuple problemTuple_;

    GridGeometryTuple gridGeometryTuple_;

    GridVariablesTuple gridVariablesTuple_;

    std::shared_ptr<const TimeLoop> timeLoop_;

    const SolutionVector* prevSol_ = nullptr;

    bool isStationaryProblem_;

    std::shared_ptr<JacobianMatrix> jacobian_;
    std::shared_ptr<SolutionVector> residual_;

    bool warningIssued_;

    bool enableMultithreading_ = false;
};
 
} // end namespace Dumux
 
#endif