couplingmanager_cvfe.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_FREEFLOW_COUPLING_MANAGER_CVFE_HH
#define DUMUX_MULTIDOMAIN_FREEFLOW_COUPLING_MANAGER_CVFE_HH
 
#include <memory>
#include <tuple>
#include <vector>
#include <deque>
#include <iostream>
 
#include <dune/common/exceptions.hh>
#include <dune/common/indices.hh>
#include <dune/common/float_cmp.hh>
#include <dune/geometry/referenceelements.hh>
 
#include <dumux/common/properties.hh>
#include <dumux/common/deprecated.hh>
#include <dumux/common/concepts/variables_.hh>
#include <dumux/common/typetraits/typetraits.hh>
#include <dumux/discretization/cvfe/localdof.hh>
 
#include <dumux/discretization/method.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/discretization/cvfe/interpolationpointdata.hh>
 
#include <dumux/multidomain/couplingmanager.hh>
#include <dumux/multidomain/fvassembler.hh>
 
#include <dumux/parallel/parallel_for.hh>
#include <dumux/assembly/coloring.hh>
 
#include "typetraits.hh"
 
namespace Dumux {

template<class Traits>
class CVFEFreeFlowCouplingManager
: public CouplingManager<Traits>
{
    using ParentType = CouplingManager<Traits>;
public:
    static constexpr auto freeFlowMomentumIndex = typename Traits::template SubDomain<0>::Index();
    static constexpr auto freeFlowMassIndex = typename Traits::template SubDomain<1>::Index();
 
    // this can be used if the coupling manager is used inside a meta-coupling manager (e.g. multi-binary)
    // to manager the solution vector storage outside this class
    using SolutionVectorStorage = typename ParentType::SolutionVectorStorage;
private:
    template<std::size_t id> using SubDomainTypeTag = typename Traits::template SubDomain<id>::TypeTag;
    template<std::size_t id> using PrimaryVariables = GetPropType<SubDomainTypeTag<id>, Properties::PrimaryVariables>;
    template<std::size_t id> using GridGeometry = GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>;
    template<std::size_t id> using GridView = typename GridGeometry<id>::GridView;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using ElementSeed = typename GridView<id>::Grid::template Codim<0>::EntitySeed;
    template<std::size_t id> using FVElementGeometry = typename GridGeometry<id>::LocalView;
    template<std::size_t id> using SubControlVolume = typename FVElementGeometry<id>::SubControlVolume;
    template<std::size_t id> using SubControlVolumeFace = typename FVElementGeometry<id>::SubControlVolumeFace;
    template<std::size_t id> using GridVariables = typename Traits::template SubDomain<id>::GridVariables;
    template<std::size_t id> using GridVariablesCache = Concept::GridVariablesCache_t<GridVariables<id>>;
    template<std::size_t id> using ElementVariables = typename GridVariablesCache<id>::LocalView;
    template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
    template<std::size_t id> using Variables = Concept::Variables_t<GridVariables<id>>;
 
    using Scalar = typename Traits::Scalar;
    using SolutionVector = typename Traits::SolutionVector;
 
    using CouplingStencilType = std::vector<std::size_t>;
 
    using GridVariablesTuple = typename Traits::template TupleOfSharedPtr<GridVariables>;
 
    using FluidSystem = typename Variables<freeFlowMassIndex>::FluidSystem;
 
    using GlobalPosition = typename SubControlVolumeFace<freeFlowMassIndex>::GlobalPosition;
    using VelocityVector = GlobalPosition;
    using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
 
    static_assert(std::is_same_v<VelocityVector, typename SubControlVolumeFace<freeFlowMomentumIndex>::GlobalPosition>);
 
    template<class ElementSolution>
    struct MomentumCouplingContextNoCaching
    {
        MomentumCouplingContextNoCaching(ElementSolution&& elemSol)
        : elemSol_(std::move(elemSol)) {}
 
        template<class GridVarsCache, class FvElementGeometry, class SubControlVolume>
        auto vars(const GridVarsCache& gridVarsCache, const FvElementGeometry& fvGeometry, const SubControlVolume& scv) const
        {
            const auto& problem = gridVarsCache.problem();
            Variables<freeFlowMassIndex> variables;
            if constexpr (Concept::FVGridVariables<GridVariables<freeFlowMassIndex>>)
                variables.update(elemSol_, problem, fvGeometry.element(), scv);
            else
                variables.update(elemSol_, problem, fvGeometry, ipData(fvGeometry, scv));
            return variables;
        }
 
        ElementSolution elemSol_;
    };
 
    struct MomentumCouplingContextGlobalCaching
    {
        template<class GridVarsCache, class FvElementGeometry, class SubControlVolume>
        const auto& vars(const GridVarsCache& gridVarsCache, const FvElementGeometry& fvGeometry, const SubControlVolume& scv) const
        {
            if constexpr (requires { gridVarsCache.volVars(scv); })
                return gridVarsCache.volVars(scv);
            else
                return gridVarsCache.variables(scv);
        }
    };
 
    using MomentumDiscretizationMethod = typename GridGeometry<freeFlowMomentumIndex>::DiscretizationMethod;
    using MassDiscretizationMethod = typename GridGeometry<freeFlowMassIndex>::DiscretizationMethod;
 
    template<std::size_t id> using IpData
        = Dumux::CVFE::InterpolationPointData<typename GridView<id>::template Codim<0>::Entity::Geometry::LocalCoordinate,
                                              typename GridView<id>::template Codim<0>::Entity::Geometry::GlobalCoordinate>;
 
public:
 
    static constexpr auto pressureIdx = Variables<freeFlowMassIndex>::Indices::pressureIdx;

    // \{

    //! use as regular coupling manager
    void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem,
              std::shared_ptr<Problem<freeFlowMassIndex>> massProblem,
              GridVariablesTuple&& gridVariables,
              const SolutionVector& curSol)
    {
        this->setSubProblems(std::make_tuple(momentumProblem, massProblem));
        gridVariables_ = gridVariables;
        this->updateSolution(curSol);
        prevSol_ = nullptr;
 
        computeCouplingStencils_();
    }

    //! use as regular coupling manager in a transient setting
    void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem,
              std::shared_ptr<Problem<freeFlowMassIndex>> massProblem,
              GridVariablesTuple&& gridVariables,
              const SolutionVector& curSol,
              const SolutionVector& prevSol)
    {
        init(momentumProblem, massProblem, std::forward<GridVariablesTuple>(gridVariables), curSol);
        prevSol_ = &prevSol;
    }

    //! use as binary coupling manager in multi model context
    void init(std::shared_ptr<Problem<freeFlowMomentumIndex>> momentumProblem,
              std::shared_ptr<Problem<freeFlowMassIndex>> massProblem,
              GridVariablesTuple&& gridVariables,
              const typename ParentType::SolutionVectorStorage& curSol)
    {
        this->setSubProblems(std::make_tuple(momentumProblem, massProblem));
        gridVariables_ = gridVariables;
        this->attachSolution(curSol);
        prevSol_ = nullptr;
 
        computeCouplingStencils_();
    }
 
    // \}

    // \{

    [[deprecated("This method will be removed after release (3.11). Use pressure(..., ipData) instead!")]]
    Scalar pressure(const Element<freeFlowMomentumIndex>& element,
                    const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                    const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
                    const bool considerPreviousTimeStep = false) const
    {
        const auto& globalPos = scvf.ipGlobal();
        const auto& localPos = element.geometry().local(globalPos);
        return this->pressure(element, fvGeometry, IpData<freeFlowMassIndex>(localPos, globalPos), considerPreviousTimeStep);
    }

    [[deprecated("This method will be removed after release (3.11). Use pressure(..., ipData) instead!")]]
    Scalar pressure(const Element<freeFlowMomentumIndex>& element,
                    const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                    const SubControlVolume<freeFlowMomentumIndex>& scv,
                    const bool considerPreviousTimeStep = false) const
    {
        return this->pressure(element, fvGeometry, ipData(fvGeometry, scv), considerPreviousTimeStep);
    }

    template <class IpData>
    Scalar pressure(const Element<freeFlowMomentumIndex>& element,
                    const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                    const IpData& ipData,
                    const bool considerPreviousTimeStep = false) const
    {
        assert(!(considerPreviousTimeStep && !isTransient_()));
        const auto& gg = Deprecated::gridGeometry(this->problem(freeFlowMassIndex));
        const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex]
                                                   :  this->curSol(freeFlowMassIndex);
        const auto elemSol = elementSolution(element, sol, gg);
        return evalSolutionAtLocalPos(element, element.geometry(), gg, elemSol, ipData.local())[pressureIdx];
    }

    [[deprecated("This method will be removed after release (3.11). Use density(..., ipData) instead!")]]
    Scalar density(const Element<freeFlowMomentumIndex>& element,
                   const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                   const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
                   const bool considerPreviousTimeStep = false) const
    {
        const auto& globalPos = scvf.ipGlobal();
        const auto& localPos = element.geometry().local(globalPos);
        return this->density(element, fvGeometry,  IpData<freeFlowMassIndex>(localPos, globalPos), considerPreviousTimeStep);
    }

    [[deprecated("This method will be removed after release (3.11). Use density(..., ipData) instead!")]]
    Scalar density(const Element<freeFlowMomentumIndex>& element,
                   const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                   const SubControlVolume<freeFlowMomentumIndex>& scv,
                   const bool considerPreviousTimeStep = false) const
    {
        return this->density(element, fvGeometry, ipData(fvGeometry, scv), considerPreviousTimeStep);
    }

    template <class IpData>
    Scalar density(const Element<freeFlowMomentumIndex>& element,
                   const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                   const IpData& ipData,
                   const bool considerPreviousTimeStep = false) const
    {
        assert(!(considerPreviousTimeStep && !isTransient_()));
 
        const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex]
                                                   :  this->curSol(freeFlowMassIndex);
 
        auto massFvGeometry = localView(Deprecated::gridGeometry(this->problem(freeFlowMassIndex)));
        massFvGeometry.bind(element);
        const auto context = makeMomentumCouplingContext_(massFvGeometry, sol);
        const auto& gridVarsCache = subDomainGridVars_(Dune::index_constant<freeFlowMassIndex>{}, considerPreviousTimeStep);
 
        if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa)
        {
            const auto eIdx = fvGeometry.elementIndex();
            const auto& scv = massFvGeometry.scv(eIdx);
 
            const auto& variables = context.vars(gridVarsCache, massFvGeometry, scv);
 
            return variables.density();
        }
        else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box
                           || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond)
        {
            // TODO: cache the shape values
            using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
            const auto& localBasis = massFvGeometry.feLocalBasis();
            std::vector<ShapeValue> shapeValues;
            localBasis.evaluateFunction(ipData.local(), shapeValues);
 
            Scalar rho = 0.0;
            for (const auto& scv : scvs(massFvGeometry))
            {
                const auto& variables = context.vars(gridVarsCache, massFvGeometry, scv);
                rho += variables.density()*shapeValues[scv.localDofIndex()][0];
            }
 
            return rho;
        }
        else
            DUNE_THROW(Dune::NotImplemented,
                "Density interpolation for discretization scheme " << MassDiscretizationMethod{}
            );
    }

    [[deprecated("This method will be removed after release (3.11). Use effectiveViscosity(..., ipData) instead!")]]
    Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element,
                              const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                              const SubControlVolumeFace<freeFlowMomentumIndex>& scvf,
                              const bool considerPreviousTimeStep = false) const
    {
        const auto& globalPos = scvf.ipGlobal();
        const auto& localPos = element.geometry().local(globalPos);
        return this->effectiveViscosity(element, fvGeometry, IpData<freeFlowMassIndex>(localPos, globalPos), considerPreviousTimeStep);
    }

    [[deprecated("This method will be removed after release (3.11). Use effectiveViscosity(..., ipData) instead!")]]
    Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element,
                              const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                              const SubControlVolume<freeFlowMomentumIndex>& scv,
                              const bool considerPreviousTimeStep = false) const
    {
        return this->effectiveViscosity(element, fvGeometry, ipData(fvGeometry, scv), considerPreviousTimeStep);
    }

    template <class IpData>
    Scalar effectiveViscosity(const Element<freeFlowMomentumIndex>& element,
                              const FVElementGeometry<freeFlowMomentumIndex>& fvGeometry,
                              const IpData& ipData,
                              const bool considerPreviousTimeStep = false) const
    {
        assert(!(considerPreviousTimeStep && !isTransient_()));
 
        const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMassIndex]
                                                   :  this->curSol(freeFlowMassIndex);
 
        auto massFvGeometry = localView(Deprecated::gridGeometry(this->problem(freeFlowMassIndex)));
        massFvGeometry.bind(element);
        const auto context = makeMomentumCouplingContext_(massFvGeometry, sol);
        const auto& gridVarsCache = subDomainGridVars_(Dune::index_constant<freeFlowMassIndex>{}, considerPreviousTimeStep);
 
        if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa)
        {
            const auto eIdx = fvGeometry.elementIndex();
            const auto& scv = massFvGeometry.scv(eIdx);
 
            const auto& variables = context.vars(gridVarsCache, massFvGeometry, scv);
 
            return variables.viscosity();
        }
        else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box
                           || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond)
        {
            // TODO: cache the shape values
            using ShapeValue = typename Dune::FieldVector<Scalar, 1>;
            const auto& localBasis = massFvGeometry.feLocalBasis();
            std::vector<ShapeValue> shapeValues;
            localBasis.evaluateFunction(ipData.local(), shapeValues);
 
            Scalar mu = 0.0;
            for (const auto& scv : scvs(massFvGeometry))
            {
                const auto& variables = context.vars(gridVarsCache, massFvGeometry, scv);
                mu += variables.viscosity()*shapeValues[scv.localDofIndex()][0];
            }
 
            return mu;
        }
        else
            DUNE_THROW(Dune::NotImplemented,
                "Viscosity interpolation for discretization scheme " << MassDiscretizationMethod{}
            );
    }

     */
    VelocityVector faceVelocity(const Element<freeFlowMassIndex>& element,
                                const SubControlVolumeFace<freeFlowMassIndex>& scvf) const
    {
        // TODO: optimize this function for tpfa where the scvf ip coincides with the dof location
        auto fvGeometry = localView(Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex)));
        fvGeometry.bindElement(element);
 
        const auto& localBasis = fvGeometry.feLocalBasis();
 
        std::vector<ShapeValue> shapeValues;
        const auto ipLocal = element.geometry().local(scvf.ipGlobal());
        localBasis.evaluateFunction(ipLocal, shapeValues);
 
        // interpolate velocity at scvf
        VelocityVector velocity(0.0);
        for (const auto& localDof : localDofs(fvGeometry))
            velocity.axpy(shapeValues[localDof.index()][0], this->curSol(freeFlowMomentumIndex)[localDof.dofIndex()]);
 
        return velocity;
    }

     */
    VelocityVector elementVelocity(const FVElementGeometry<freeFlowMassIndex>& fvGeometry) const
    {
        auto momentumFvGeometry = localView(Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex)));
        momentumFvGeometry.bindElement(fvGeometry.element());
 
        const auto& localBasis = momentumFvGeometry.feLocalBasis();
 
        // interpolate velocity at scvf
        VelocityVector velocity(0.0);
        std::vector<ShapeValue> shapeValues;
        localBasis.evaluateFunction(referenceElement(fvGeometry.element()).position(0,0), shapeValues);
 
        for (const auto& localDof : localDofs(momentumFvGeometry))
            velocity.axpy(shapeValues[localDof.index()][0], this->curSol(freeFlowMomentumIndex)[localDof.dofIndex()]);
 
        return velocity;
    }

    template <class IpData>
    VelocityVector velocity(const FVElementGeometry<freeFlowMassIndex>& fvGeometry,
                            const IpData& ipData,
                            const bool considerPreviousTimeStep = false) const
    {
        assert(!(considerPreviousTimeStep && !isTransient_()));
 
        const auto& element = fvGeometry.element();
        const auto& gg = Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex));
        auto momentumFvGeometry = localView(gg);
        momentumFvGeometry.bindElement(fvGeometry.element());
 
        const auto& sol = considerPreviousTimeStep ? (*prevSol_)[freeFlowMomentumIndex]
                                                   :  this->curSol(freeFlowMomentumIndex);
 
        const auto elemSol = elementSolution(element, sol, gg);
        return evalSolutionAtLocalPos(element, element.geometry(), gg, elemSol, ipData.local());
    }

    template<std::size_t j>
    const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
                                               const Element<freeFlowMomentumIndex>& elementI,
                                               const SubControlVolume<freeFlowMomentumIndex>& scvI,
                                               Dune::index_constant<j> domainJ) const
    { return emptyStencil_; }

     */
    const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMassIndex> domainI,
                                              const Element<freeFlowMassIndex>& elementI,
                                              Dune::index_constant<freeFlowMomentumIndex> domainJ) const
    {
        const auto eIdx = Deprecated::gridGeometry(this->problem(freeFlowMassIndex)).elementMapper().index(elementI);
        return massAndEnergyToMomentumStencils_[eIdx];
    }

     */
    const CouplingStencilType& couplingStencil(Dune::index_constant<freeFlowMomentumIndex> domainI,
                                               const Element<freeFlowMomentumIndex>& elementI,
                                               Dune::index_constant<freeFlowMassIndex> domainJ) const
    {
        const auto eIdx = Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex)).elementMapper().index(elementI);
        return momentumToMassAndEnergyStencils_[eIdx];
    }
 
    // \}

    // \{

    template<std::size_t i, std::size_t j, class LocalAssemblerI>
    void updateCouplingContext(Dune::index_constant<i> domainI,
                               const LocalAssemblerI& localAssemblerI,
                               Dune::index_constant<j> domainJ,
                               std::size_t dofIdxGlobalJ,
                               const PrimaryVariables<j>& priVarsJ,
                               int pvIdxJ)
    {
        this->curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
 
        if constexpr (GridVariablesCache<freeFlowMassIndex>::cachingEnabled)
        {
            if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::cctpfa)
            {
                if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex)
                {
                    const auto& problem = this->problem(domainJ);
                    const auto& deflectedElement = Deprecated::gridGeometry(problem).element(dofIdxGlobalJ);
                    const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), Deprecated::gridGeometry(problem));
                    auto fvGeometry = localView(Deprecated::gridGeometry(problem));
                    fvGeometry.bind(deflectedElement);
                    const auto& scv = fvGeometry.scv(dofIdxGlobalJ);
 
                    if constexpr (Concept::FVGridVariables<GridVariables<freeFlowMassIndex>>)
                        subDomainVariables_(Dune::index_constant<freeFlowMassIndex>{}, /*current*/true, scv).update(std::move(elemSol), problem, deflectedElement, scv);
                    else
                        subDomainVariables_(Dune::index_constant<freeFlowMassIndex>{}, /*current*/true, ipData(fvGeometry, scv)).update(std::move(elemSol), problem, deflectedElement, scv);
                }
            }
            else if constexpr (MassDiscretizationMethod{} == DiscretizationMethods::box
                            || MassDiscretizationMethod{} == DiscretizationMethods::fcdiamond)
            {
                if constexpr (domainI == freeFlowMomentumIndex && domainJ == freeFlowMassIndex)
                {
                    const auto& problem = this->problem(domainJ);
                    const auto deflectedElementIdx = Deprecated::gridGeometry(problem).elementMapper().index(localAssemblerI.element());
                    const auto& deflectedElement = Deprecated::gridGeometry(problem).element(deflectedElementIdx);
                    const auto elemSol = elementSolution(deflectedElement, this->curSol(domainJ), Deprecated::gridGeometry(problem));
                    auto fvGeometry = localView(Deprecated::gridGeometry(problem));
                    fvGeometry.bind(deflectedElement);
 
                    // ToDo: Replace once all mass models are also working with local dofs
                    for (const auto& scv : scvs(fvGeometry))
                    {
                        if (scv.dofIndex() == dofIdxGlobalJ)
                        {
                            if constexpr (Concept::FVGridVariables<GridVariables<freeFlowMassIndex>>)
                                this->subDomainVariables_(Dune::index_constant<freeFlowMassIndex>{}, /*current*/true, scv).update(std::move(elemSol), problem, deflectedElement, scv);
                            else
                                this->subDomainVariables_(Dune::index_constant<freeFlowMassIndex>{}, /*current*/true, scv).update(std::move(elemSol), problem, fvGeometry, ipData(fvGeometry, scv));
                        }
                    }
                }
            }
            else
                DUNE_THROW(Dune::NotImplemented,
                    "Context update for discretization scheme " << MassDiscretizationMethod{}
                );
        }
    }
 
    // \}

     */
    void computeColorsForAssembly()
    {
        if constexpr (MomentumDiscretizationMethod{} == DiscretizationMethods::fcdiamond)
        {
            // use coloring of the mass discretization for both domains
            // the diamond coloring is a subset (minimum amount of colors) of cctpfa/box coloring
            elementSets_ = computeColoring(Deprecated::gridGeometry(this->problem(freeFlowMassIndex))).sets;
        }
        else
        {
            // use coloring of the momentum discretization for both domains
            elementSets_ = computeColoring(Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex))).sets;
        }
    }

    template<std::size_t i, class AssembleElementFunc>
    void assembleMultithreaded(Dune::index_constant<i> domainId, AssembleElementFunc&& assembleElement) const
    {
        if (elementSets_.empty())
            DUNE_THROW(Dune::InvalidStateException, "Call computeColorsForAssembly before assembling in parallel!");
 
        // make this element loop run in parallel
        // for this we have to color the elements so that we don't get
        // race conditions when writing into the global matrix
        // each color can be assembled using multiple threads
        const auto& grid = Deprecated::gridGeometry(this->problem(freeFlowMassIndex)).gridView().grid();
        for (const auto& elements : elementSets_)
        {
            Dumux::parallelFor(elements.size(), [&](const std::size_t eIdx)
            {
                const auto element = grid.entity(elements[eIdx]);
                assembleElement(element);
            });
        }
    }
 
private:
    template<std::size_t i>
    const auto& subDomainGridVars_(Dune::index_constant<i> domainIdx, bool current) const
    {
        if constexpr (Concept::FVGridVariables<GridVariables<i>>)
            return current ? gridVars_(domainIdx).curGridVolVars()
                           : gridVars_(domainIdx).prevGridVolVars();
        else
            return current ? gridVars_(domainIdx).curGridVars()
                           : gridVars_(domainIdx).prevGridVars();
    }
 
    template<std::size_t i>
    auto& subDomainGridVars_(Dune::index_constant<i> domainIdx, bool current)
    {
        if constexpr (Concept::FVGridVariables<GridVariables<i>>)
            return current ? gridVars_(domainIdx).curGridVolVars()
                           : gridVars_(domainIdx).prevGridVolVars();
        else
            return current ? gridVars_(domainIdx).curGridVars()
                           : gridVars_(domainIdx).prevGridVars();
    }
 
    template<std::size_t i, class Scv>
    auto& subDomainVariables_(Dune::index_constant<i> domainIdx, bool current, const Scv& scv)
    {
        if constexpr (Concept::FVGridVariables<GridVariables<i>>)
            return subDomainGridVars_(domainIdx, current).volVars(scv);
        else
            return subDomainGridVars_(domainIdx, current).variables(scv);
    }
 
    template<std::size_t i, class Scv>
    const auto& subDomainVariables_(Dune::index_constant<i> domainIdx, bool current, const Scv& scv) const
    {
        if constexpr (Concept::FVGridVariables<GridVariables<i>>)
            return subDomainGridVars_(domainIdx, current).volVars(scv);
        else
            return subDomainGridVars_(domainIdx, current).variables(scv);
    }
 
    template<class SolutionVector>
    auto makeMomentumCouplingContext_(const FVElementGeometry<freeFlowMassIndex>& fvGeometry,
                                      const SolutionVector& sol) const
    {
        if constexpr (GridVariablesCache<freeFlowMassIndex>::cachingEnabled)
            return MomentumCouplingContextGlobalCaching{};
        else
            return MomentumCouplingContextNoCaching{elementSolution(fvGeometry.element(), sol, Deprecated::gridGeometry(fvGeometry))};
    }

    template<std::size_t i>
    const GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx) const
    {
        if (std::get<i>(gridVariables_))
            return *std::get<i>(gridVariables_);
        else
            DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function");
    }

    template<std::size_t i>
    GridVariables<i>& gridVars_(Dune::index_constant<i> domainIdx)
    {
        if (std::get<i>(gridVariables_))
            return *std::get<i>(gridVariables_);
        else
            DUNE_THROW(Dune::InvalidStateException, "The gridVariables pointer was not set. Use setGridVariables() before calling this function");
    }
 
    void computeCouplingStencils_()
    {
        const auto& momentumGridGeometry = Deprecated::gridGeometry(this->problem(freeFlowMomentumIndex));
        const auto& massGridGeometry = Deprecated::gridGeometry(this->problem(freeFlowMassIndex));
        auto momentumFvGeometry = localView(momentumGridGeometry);
        auto massFvGeometry = localView(massGridGeometry);
 
        massAndEnergyToMomentumStencils_.clear();
        massAndEnergyToMomentumStencils_.resize(massGridGeometry.gridView().size(0));
 
        momentumToMassAndEnergyStencils_.clear();
        momentumToMassAndEnergyStencils_.resize(momentumGridGeometry.gridView().size(0));
 
        assert(massAndEnergyToMomentumStencils_.size() == momentumToMassAndEnergyStencils_.size());
        bool hasCubeElements = false;
 
        for (const auto& element : elements(momentumGridGeometry.gridView()))
        {
            if(element.type().isCube())
                hasCubeElements = true;
 
            momentumFvGeometry.bindElement(element);
            massFvGeometry.bindElement(element);
            const auto eIdx = momentumFvGeometry.elementIndex();
 
            for (const auto& localDof : localDofs(momentumFvGeometry))
                massAndEnergyToMomentumStencils_[eIdx].push_back(localDof.dofIndex());
 
            // ToDo: Replace once all mass models are also working with local dofs
            for (const auto& scv : scvs(massFvGeometry))
                momentumToMassAndEnergyStencils_[eIdx].push_back(scv.dofIndex());
        }
 
        // Print warning for pq1bubble scheme on cube elements if not using the hybrid variant
        if constexpr ((MomentumDiscretizationMethod{} == DiscretizationMethods::pq1bubble))
        {
            if(hasCubeElements && !GridGeometry<freeFlowMomentumIndex>::enableHybridCVFE)
            {
                std::cerr << "Warning: Coupled Navier-Stokes problem on cube elements uses non-hybrid pq1bubble. "
                          << "The hybrid variant is recommended because it implements two bubble functions for stability reasons."
                          << std::endl;
            }
        }
    }
 
    CouplingStencilType emptyStencil_;
    std::vector<CouplingStencilType> momentumToMassAndEnergyStencils_;
    std::vector<CouplingStencilType> massAndEnergyToMomentumStencils_;

    GridVariablesTuple gridVariables_;

    bool isTransient_() const
    { return prevSol_ != nullptr; }
 
    const SolutionVector* prevSol_ = nullptr;
 
    std::deque<std::vector<ElementSeed<freeFlowMomentumIndex>>> elementSets_;
};
 
namespace Detail {
 
// declaration (specialize for different discretization types)
template<class Traits, class DiscretizationMethod = typename Detail::MomentumDiscretizationMethod<Traits>::type>
struct CouplingManagerSupportsMultithreadedAssemblySelector;
 
// multi-threading is not supported because we have only one coupling context instance and a mutable cache
template<class Traits, class D>
struct CouplingManagerSupportsMultithreadedAssemblySelector<Traits, DiscretizationMethods::CVFE<D>>
{ using type = std::true_type; };
 
} // end namespace Detail

template<class T>
struct CouplingManagerSupportsMultithreadedAssembly<CVFEFreeFlowCouplingManager<T>>
: public Detail::CouplingManagerSupportsMultithreadedAssemblySelector<T>::type
{};
 
} // end namespace Dumux
 
#endif