multidomain/assembler.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-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
//
#ifndef DUMUX_MULTIDOMAIN_ASSEMBLER_HH
#define DUMUX_MULTIDOMAIN_ASSEMBLER_HH
 
#include <algorithm>
#include <cassert>
#include <iostream>
#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/typetraits/periodic.hh>
#include <dumux/common/gridcapabilities.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/discretization/method.hh>
#include <dumux/common/concepts/variables_.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 "subdomaincvfelocalassembler_.hh"
#include "subdomaincvfelocalassembler.hh"
#include "subdomainfclocalassembler.hh"
#include "assemblerview.hh"
#include "assemblytraits.hh"
 
#include <dumux/discretization/method.hh>
#if HAVE_DUMUX_OLD_STAGGERED
#include <dumux/multidomain/subdomainstaggeredlocalassembler.hh>
#endif
 
namespace Dumux::Experimental {

template<class MDTraits, class CMType, DiffMethod diffMethod, bool useImplicitAssembly = true>
class MultiDomainAssembler
{
    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 = typename MDTraits::template SubDomain<id>::LocalResidual;
 
    template<std::size_t id>
    using GridVariables = typename MDTraits::template SubDomain<id>::GridVariables;
 
    template<std::size_t id>
    using GridDiscretization = GetPropType<SubDomainTypeTag<id>, Properties::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 = typename MDTraits::ResidualVector;
 
    using CouplingManager = CMType;

    static constexpr bool isImplicit()
    { return useImplicitAssembly; }
 
private:
 
    using ProblemTuple = typename MDTraits::template TupleOfSharedPtrConst<Problem>;
    using GridDiscretizationTuple = typename MDTraits::template TupleOfSharedPtrConst<GridDiscretization>;
    using GridVariablesTuple = typename MDTraits::template TupleOfSharedPtr<GridVariables>;
 
    using TimeLoop = TimeLoopBase<Scalar>;
    using ThisType = MultiDomainAssembler<MDTraits, CouplingManager, diffMethod, isImplicit()>;
 
    template<std::size_t id>
    using SubDomainAssemblerView = MultiDomainAssemblerSubDomainView<ThisType, id>;
 
    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>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::CCMpfa, id>
    {
        using type = SubDomainCCLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id, class DM>
    struct SubDomainAssemblerType<DiscretizationMethods::CVFE<DM>, id>
    {
    private:
        using GV = GridVariables<id>;
        static constexpr bool usesGeneralGridVariables =
            Dumux::Concept::GridVariables<GV> && !Dumux::Concept::FVGridVariables<GV>;
    public:
        using type = std::conditional_t<usesGeneralGridVariables,
                                        Dumux::Experimental::SubDomainCVFELocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>,
                                        Dumux::SubDomainCVFELocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>>;
    };
 
#if HAVE_DUMUX_OLD_STAGGERED
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::Staggered, id>
    {
        using type = SubDomainStaggeredLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
    };
#endif
 
    template<std::size_t id>
    struct SubDomainAssemblerType<DiscretizationMethods::FCStaggered, id>
    {
        using type = SubDomainFaceCenteredLocalAssembler<id, SubDomainTypeTag<id>, SubDomainAssemblerView<id>, diffMethod, isImplicit()>;
    };
 
    template<std::size_t id>
    using SubDomainAssembler = typename SubDomainAssemblerType<typename GridDiscretization<id>::DiscretizationMethod, id>::type;
 
public:

    MultiDomainAssembler(ProblemTuple problem,
                         GridDiscretizationTuple gridDiscretization,
                         GridVariablesTuple gridVariables,
                         std::shared_ptr<CouplingManager> couplingManager)
    : couplingManager_(couplingManager)
    , problemTuple_(std::move(problem))
    , gridDiscretizationTuple_(std::move(gridDiscretization))
    , 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(gridDiscretizationTuple_)
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }

    MultiDomainAssembler(ProblemTuple problem,
                         GridDiscretizationTuple gridDiscretization,
                         GridVariablesTuple gridVariables,
                         std::shared_ptr<CouplingManager> couplingManager,
                         std::shared_ptr<const TimeLoop> timeLoop,
                         const SolutionVector& prevSol)
    : couplingManager_(couplingManager)
    , problemTuple_(std::move(problem))
    , gridDiscretizationTuple_(std::move(gridDiscretization))
    , 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(gridDiscretizationTuple_)
            && !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];
            const auto& subSol = curSol[domainId];
            this->assembleJacobianAndResidual_(domainId, jacRow, subRes, curSol);
 
            const auto gridDiscretization = std::get<domainId>(gridDiscretizationTuple_);
            enforcePeriodicConstraints_(domainId, jacRow, subRes, *gridDiscretization, subSol);
            enforceProblemConstraints_(domainId, jacRow, subRes, *gridDiscretization, subSol);
        });
    }

    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];
            const auto& subSol = curSol[domainId];
            this->assembleResidual_(domainId, subRes, curSol);
 
            const auto gridDiscretization = std::get<domainId>(gridDiscretizationTuple_);
            enforcePeriodicConstraints_(domainId, subRes, *gridDiscretization, subSol);
            enforceProblemConstraints_(domainId, subRes, *gridDiscretization, subSol);
        });
    }

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

    void setLinearSystem()
    {
        jacobian_ = std::make_shared<JacobianMatrix>();
        residual_ = std::make_shared<ResidualType>();
 
        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()
    {
        setJacobianPattern_(*jacobian_);
        setResidualSize_(*residual_);
        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>(gridDiscretizationTuple_)->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& gridDiscretization(Dune::index_constant<i> domainId) const
    { return *std::get<domainId>(gridDiscretizationTuple_); }

    template<std::size_t i>
    const auto& gridView(Dune::index_constant<i> domainId) const
    { return gridDiscretization(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_; }

    ResidualType& 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_(ResidualType& 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<ResidualType>();
            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)
        {
            MultiDomainAssemblerSubDomainView view{*this, domainId};
            SubDomainAssembler<i> subDomainAssembler(view, 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)
        {
            MultiDomainAssemblerSubDomainView view{*this, domainId};
            SubDomainAssembler<i> subDomainAssembler(view, 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 = gridDiscretization(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, gridDiscretization(domainI),
                                                        domainJ, gridDiscretization(domainJ));
    }
 
    // build periodic constraints into the system matrix
    template<std::size_t i, class JacRow, class Res, class GG, class Sol>
    void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Res& res, const GG& gridDiscretization, const Sol& curSol) const
    {
        if constexpr (Dumux::Detail::hasPeriodicDofMap<GG>())
        {
            for (const auto& m : gridDiscretization.periodicDofMap())
            {
                const auto& periodicallyMappedDofs = Dumux::Deprecated::ensureRangeOfPeriodicDofs(m.second);
                if (std::ranges::all_of(periodicallyMappedDofs, [&](auto second){ return m.first < second; }))
                {
                    for (const auto& second : periodicallyMappedDofs)
                    {
                        auto& jac = jacRow[domainI];
 
                        // add the second row to the first
                        res[m.first] += res[second];
 
                        const auto end = jac[second].end();
                        for (auto it = jac[second].begin(); it != end; ++it)
                            jac[m.first][it.index()] += (*it);
 
 
                        // enforce the solution of the first periodic DOF to the second one
                        res[second] = curSol[second] - curSol[m.first];
 
                        // set derivatives accordingly in jacobian, i.e. id for second and -id for m.first
                        auto setMatrixBlock = [] (auto& matrixBlock, double diagValue)
                        {
                            for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
                                matrixBlock[eIdx][eIdx] = diagValue;
                        };
 
                        for (auto it = jac[second].begin(); it != end; ++it)
                        {
                            auto& matrixBlock = *it;
                            matrixBlock = 0.0;
 
                            assert(matrixBlock.N() == matrixBlock.M());
                            if(it.index() == second)
                                setMatrixBlock(matrixBlock, 1.0);
 
                            if(it.index() == m.first)
                                setMatrixBlock(matrixBlock, -1.0);
 
                        }
 
                        using namespace Dune::Hybrid;
                        forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
                        {
                            auto& jacCoupling = jacRow[couplingDomainId];
 
                            for (auto it = jacCoupling[second].begin(); it != jacCoupling[second].end(); ++it)
                                jacCoupling[m.first][it.index()] += (*it);
 
                            for (auto it = jacCoupling[second].begin(); it != jacCoupling[second].end(); ++it)
                                (*it) = 0.0;
                        });
                    }
                }
            }
        }
    }
 
    // enforce global constraints into the system matrix
    template<std::size_t i, class JacRow, class Res, class GG, class Sol>
    void enforceProblemConstraints_(Dune::index_constant<i> domainI, JacRow& jacRow, Res& res, const GG& gridDiscretization, const Sol& curSol) const
    {
        if constexpr (Dumux::Detail::hasSubProblemGlobalConstraints<Problem<domainI>>())
        {
            auto& jac = jacRow[domainI];
            auto applyDirichletConstraint = [&] (const auto& dofIdx,
                                                 const auto& values,
                                                 const auto eqIdx,
                                                 const auto pvIdx)
            {
                (res)[dofIdx][eqIdx] = curSol[dofIdx][pvIdx] - values[pvIdx];
 
                auto& row = jac[dofIdx];
                for (auto col = row.begin(); col != row.end(); ++col)
                    row[col.index()][eqIdx] = 0.0;
 
                jac[dofIdx][dofIdx][eqIdx][pvIdx] = 1.0;
 
                using namespace Dune::Hybrid;
                forEach(makeIncompleteIntegerSequence<JacRow::size(), domainI>(), [&](const auto couplingDomainId)
                {
                    auto& jacCoupling = jacRow[couplingDomainId];
 
                    auto& rowCoupling = jacCoupling[dofIdx];
                    for (auto colCoupling = rowCoupling.begin(); colCoupling != rowCoupling.end(); ++colCoupling)
                        rowCoupling[colCoupling.index()][eqIdx] = 0.0;
                });
 
            };
 
            for (const auto& constraintData : this->problem(domainI).constraints())
            {
                const auto& constraintInfo = constraintData.constraintInfo();
                const auto& values = constraintData.values();
                const auto dofIdx = constraintData.dofIndex();
                // set the constraint in residual and jacobian
                for (int eqIdx = 0; eqIdx < constraintInfo.size(); ++eqIdx)
                {
                    if (constraintInfo.isConstraintEquation(eqIdx))
                    {
                        const auto pvIdx = constraintInfo.eqToPriVarIndex(eqIdx);
                        assert(0 <= pvIdx && pvIdx < constraintInfo.size());
                        applyDirichletConstraint(dofIdx, values, eqIdx, pvIdx);
                    }
                }
            }
        }
    }
 
    // build periodic constraints in residual only
    template<std::size_t i, class Res, class GG, class Sol>
    void enforcePeriodicConstraints_(Dune::index_constant<i> domainI, Res& res, const GG& gridDiscretization, const Sol& curSol) const
    {
        if constexpr (Dumux::Detail::hasPeriodicDofMap<GG>())
        {
            for (const auto& m : gridDiscretization.periodicDofMap())
            {
                const auto& periodicallyMappedDofs = Dumux::Deprecated::ensureRangeOfPeriodicDofs(m.second);
                if (std::ranges::all_of(periodicallyMappedDofs, [&](auto second){ return m.first < second; }))
                {
                    for (const auto& second : periodicallyMappedDofs)
                    {
                        // add the second row to the first
                        res[m.first] += res[second];
 
                        // enforce the solution of the first periodic DOF to the second one
                        res[second] = curSol[second] - curSol[m.first];
                    }
                }
            }
        }
    }
 
    // enforce global constraints in residual only
    template<std::size_t i, class Res, class GG, class Sol>
    void enforceProblemConstraints_(Dune::index_constant<i> domainI, Res& res, const GG& gridDiscretization, const Sol& curSol) const
    {
        if constexpr (Dumux::Detail::hasSubProblemGlobalConstraints<Problem<domainI>>())
        {
            for (const auto& constraintData : this->problem(domainI).constraints())
            {
                const auto& constraintInfo = constraintData.constraintInfo();
                const auto& values = constraintData.values();
                const auto dofIdx = constraintData.dofIndex();
                for (int eqIdx = 0; eqIdx < constraintInfo.size(); ++eqIdx)
                {
                    if (constraintInfo.isConstraintEquation(eqIdx))
                    {
                        const auto pvIdx = constraintInfo.eqToPriVarIndex(eqIdx);
                        assert(0 <= pvIdx && pvIdx < constraintInfo.size());
                        res[dofIdx][eqIdx] = curSol[dofIdx][pvIdx] - values[pvIdx];
                    }
                }
            }
        }
    }

    ProblemTuple problemTuple_;

    GridDiscretizationTuple gridDiscretizationTuple_;

    GridVariablesTuple gridVariablesTuple_;

    std::shared_ptr<const TimeLoop> timeLoop_;

    const SolutionVector* prevSol_ = nullptr;

    bool isStationaryProblem_;

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

    mutable bool warningIssued_;

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