multidomain/couplingmanager.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_COUPLING_MANAGER_HH
#define DUMUX_MULTIDOMAIN_COUPLING_MANAGER_HH
 
#include <memory>
#include <tuple>
#include <vector>
#include <dune/common/exceptions.hh>
#include <dune/common/indices.hh>
#include <dune/common/shared_ptr.hh>
#include <dune/common/hybridutilities.hh>
#include <dune/istl/multitypeblockvector.hh>
 
#include <dumux/assembly/numericepsilon.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/typetraits/typetraits.hh>
 
namespace Dumux {
 
namespace Detail {
 
// helper to create a multitype vector of references to solution vectors
template<class... Args, std::size_t ...Is>
auto toRef(const std::tuple<Args...>& v, std::index_sequence<Is...> indices)
{
    return Dune::MultiTypeBlockVector<std::add_lvalue_reference_t<typename Args::element_type>...>(*std::get<Is>(v)...);
}
 
} // end namespace Detail
 
template<class Traits>
class CouplingManager
{
    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 GridView = typename GetPropType<SubDomainTypeTag<id>, Properties::GridGeometry>::GridView;
    template<std::size_t id> using Element = typename GridView<id>::template Codim<0>::Entity;
    template<std::size_t id> using Problem = GetPropType<SubDomainTypeTag<id>, Properties::Problem>;
    template<std::size_t id> using ProblemPtr = const Problem<id> *;
    using ProblemPtrs = typename Traits::template Tuple<ProblemPtr>;
 
    template<std::size_t id>
    using SubSolutionVector
        = std::decay_t<decltype(std::declval<typename Traits::SolutionVector>()[Dune::index_constant<id>()])>;
 
public:
    template<std::size_t i, std::size_t j>
    using CouplingStencilType = std::vector<std::size_t>;
 
    using SolutionVector = typename Traits::SolutionVector;
 
protected:
    using SolutionVectorStorage = typename Traits::template TupleOfSharedPtr<SubSolutionVector>;
 
public:
    CouplingManager()
    {
        using namespace Dune::Hybrid;
        forEach(problems_, [](auto& problem){
            problem = nullptr;
        });
 
        forEach(curSols_, [](auto& solutionVector){
            solutionVector = std::make_shared<typename std::decay_t<decltype(solutionVector)>::element_type>();
        });
    }
 
    // \{
 
    template<std::size_t i, std::size_t j>
    const CouplingStencilType<i, j>& couplingStencil(Dune::index_constant<i> domainI,
                                                     const Element<i>& elementI,
                                                     Dune::index_constant<j> domainJ) const
    {
        static_assert(i != j, "Domain i cannot be coupled to itself!");
        static_assert(AlwaysFalse<Dune::index_constant<i>>::value,
                      "The coupling manager does not implement the couplingStencil() function");
    }
 
    template<std::size_t id, class JacobianPattern>
    void extendJacobianPattern(Dune::index_constant<id> domainI, JacobianPattern& pattern) const
    {}
 
    // \}
 
    // \{
 
    template<std::size_t i, class Assembler>
    void bindCouplingContext(Dune::index_constant<i> domainI,
                             const Element<i>& elementI,
                             const Assembler& assembler)
    {}
 
 
    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)
    {
        curSol(domainJ)[dofIdxGlobalJ][pvIdxJ] = priVarsJ[pvIdxJ];
    }
 
    template<std::size_t i, class LocalAssemblerI, class UpdatableElementVolVars, class UpdatableFluxVarCache>
    void updateCoupledVariables(Dune::index_constant<i> domainI,
                                const LocalAssemblerI& localAssemblerI,
                                UpdatableElementVolVars& elemVolVars,
                                UpdatableFluxVarCache& elemFluxVarsCache)
    {}
 
    void updateSolution(const SolutionVector& curSol)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(curSols_)), [&](const auto id)
        {
            // copy external solution into object stored in this class
            *std::get<id>(curSols_) = curSol[id];
        });
    }
 
    // \}
 
    template<std::size_t i, std::size_t j, class LocalAssemblerI>
    decltype(auto) evalCouplingResidual(Dune::index_constant<i> domainI,
                                        const LocalAssemblerI& localAssemblerI,
                                        Dune::index_constant<j> domainJ,
                                        std::size_t dofIdxGlobalJ) const
    {
        return localAssemblerI.evalLocalResidual();
    }
 
    template<std::size_t i, class LocalAssemblerI, class JacobianMatrixDiagBlock, class GridVariables>
    void evalAdditionalDomainDerivatives(Dune::index_constant<i> domainI,
                                         const LocalAssemblerI& localAssemblerI,
                                         const typename LocalAssemblerI::LocalResidual::ElementResidualVector& origResiduals,
                                         JacobianMatrixDiagBlock& A,
                                         GridVariables& gridVariables)
    {}
 
    template<std::size_t i>
    decltype(auto) numericEpsilon(Dune::index_constant<i>,
                                  const std::string& paramGroup) const
    {
        constexpr auto numEq = PrimaryVariables<i>::dimension;
        return NumericEpsilon<typename Traits::Scalar, numEq>(paramGroup);
    }
 
    template<typename... SubProblems>
    void setSubProblems(const std::tuple<std::shared_ptr<SubProblems>...>& problems)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(size(problems_)), [&](const auto i)
        { setSubProblem(std::get<i>(problems), i); });
    }
 
    template<class SubProblem, std::size_t i>
    void setSubProblem(std::shared_ptr<SubProblem> problem, Dune::index_constant<i> domainIdx)
    { std::get<i>(problems_) = problem.get(); }
 
    template<std::size_t i>
    const Problem<i>& problem(Dune::index_constant<i> domainIdx) const
    {
        const Problem<i>* p = std::get<i>(problems_);
        assert(p && "The problem pointer is invalid. Use setSubProblems() before calling this function");
        return *p;
    }
 
protected:
    void attachSolution(SolutionVectorStorage& curSol)
    {
        using namespace Dune::Hybrid;
        forEach(integralRange(Dune::Hybrid::size(curSols_)), [&](const auto id)
        {
            // do not take ownership of the external pointer's object
            std::get<id>(curSols_) = Dune::stackobject_to_shared_ptr(*std::get<id>(curSol));
        });
    }
 
    template<std::size_t i>
    SubSolutionVector<i>& curSol(Dune::index_constant<i> domainIdx)
    { return *std::get<i>(curSols_); }
 
    template<std::size_t i>
    const SubSolutionVector<i>& curSol(Dune::index_constant<i> domainIdx) const
    { return *std::get<i>(curSols_); }
 
    [[deprecated("This function returns a Dune::MultiTypeBlockVector<SubDomainSolutionVector&, ....> (i.e. storing references). "
                 "Use curSol(domainIdx) to get a reference to the corresponding subdomain solution vector. Will be removed after 3.5")]]
    decltype(auto) curSol()
    {
        return Detail::toRef(curSols_, std::make_index_sequence<Traits::numSubDomains>());
    }
 
    [[deprecated("This function returns a Dune::MultiTypeBlockVector<SubDomainSolutionVector&, ....> (i.e. storing references). "
                 "Use curSol(domainIdx) to get a reference to the corresponding subdomain solution vector. Will be removed after 3.5")]]
    decltype(auto) curSol() const
    {
        return Detail::toRef(curSols_, std::make_index_sequence<Traits::numSubDomains>());
    }
 
private:
    SolutionVectorStorage curSols_;
 
    ProblemPtrs problems_;
};
 
} // end namespace Dumux
 
#endif