assembly/localassemblerbase.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_LOCAL_ASSEMBLER_BASE_HH
#define DUMUX_LOCAL_ASSEMBLER_BASE_HH
 
#include <type_traits>
#include <utility>
 
#include <dune/grid/common/gridenums.hh> // for GhostEntity
 
#include <dumux/common/properties.hh>
 
namespace Dumux::Experimental {

template<class TypeTag, class Assembler, class Implementation, bool useImplicitAssembly>
class LocalAssemblerBase
{
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using SolutionVector = typename Assembler::SolutionVector;
    using ElementDiscretization = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
    using GridVariablesCache = typename GridVariables::GridVariablesCache;
    using ElementVariables = typename GridVariablesCache::LocalView;
    using Element = typename GridView::template Codim<0>::Entity;
 
public:
    using LocalResidual = std::decay_t<decltype(std::declval<Assembler>().localResidual())>;
    using ElementResidualVector = typename LocalResidual::ElementResidualVector;

    explicit LocalAssemblerBase(const Assembler& assembler,
                                const Element& element,
                                const SolutionVector& curSol)
    : LocalAssemblerBase(assembler,
                         element,
                         curSol,
                         localView(assembler.gridDiscretization()),
                         localView(assembler.gridVariables().curGridVars()),
                         localView(assembler.gridVariables().prevGridVars()),
                         assembler.localResidual(),
                         element.partitionType() == Dune::GhostEntity)
    {}

    explicit LocalAssemblerBase(const Assembler& assembler,
                                const Element& element,
                                const SolutionVector& curSol,
                                const ElementDiscretization& elemDisc,
                                const ElementVariables& curElemVars,
                                const ElementVariables& prevElemVars,
                                const LocalResidual& localResidual,
                                const bool elementIsGhost)
    : assembler_(assembler)
    , element_(element)
    , curSol_(curSol)
    , elemDisc_(elemDisc)
    , curElemVars_(curElemVars)
    , prevElemVars_(prevElemVars)
    , localResidual_(localResidual)
    , elementIsGhost_(elementIsGhost)
    {}

    static constexpr bool isImplicit()
    { return useImplicitAssembly; }

    ElementResidualVector evalLocalResidual() const
    {
        if (!isImplicit())
            if (this->assembler().isStationaryProblem())
                DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
 
        if (elementIsGhost())
            return ElementResidualVector(0.0);
 
        return isImplicit() ? evalLocalResidual(curElemVars())
                            : evalLocalResidual(prevElemVars());
    }

    ElementResidualVector evalLocalResidual(const ElementVariables& elemVars) const
    {
        if constexpr (requires (const Assembler& a) { a.currentStageWeights(); })
        {
            const auto [temporalWeight, spatialWeight] = assembler().currentStageWeights();
            return evalLocalResidualForStage(elemVars, temporalWeight, spatialWeight);
        }
        else if (!assembler().isStationaryProblem())
        {
            ElementResidualVector residual = evalLocalFluxAndSourceResidual(elemVars);
            residual += evalLocalStorageResidual();
            return residual;
        }
        else
            return evalLocalFluxAndSourceResidual(elemVars);
    }

    template<class Scalar>
    ElementResidualVector evalLocalResidualForStage(const ElementVariables& elemVars,
                                                     Scalar temporalWeight,
                                                     Scalar spatialWeight) const
    {
        auto spatial = evalLocalFluxAndSourceResidual(elemVars);
        for (auto& r : spatial) r *= spatialWeight;
 
        auto temporal = localResidual_.evalStorageCurrentLevel(element_, elemDisc_, elemVars);
        for (auto& r : temporal) r *= temporalWeight;
 
        spatial += temporal;
        return spatial;
    }

    ElementResidualVector evalLocalFluxAndSourceResidual() const
    {
        return isImplicit() ? evalLocalFluxAndSourceResidual(curElemVars())
                            : evalLocalFluxAndSourceResidual(prevElemVars());
    }

    ElementResidualVector evalLocalFluxAndSourceResidual(const ElementVariables& elemVars) const
    {
        return localResidual_.evalFluxAndSource(element_, elemDisc_, elemVars);
    }

    ElementResidualVector evalLocalStorageResidual() const
    {
        return localResidual_.evalStorage(element_, elemDisc_, prevElemVars_, curElemVars_);
    }

    void bindLocalViews()
    {
        // get some references for convenience
        const auto& element = this->element();
        const auto& curSol = this->curSol();
        const auto& prevSol = this->assembler().prevSol();
        auto&& elemDisc = this->elemDisc();
        auto&& curElemVars = this->curElemVars();
        auto&& prevElemVars = this->prevElemVars();
 
        // bind the caches
        elemDisc.bind(element);
 
        if (isImplicit())
        {
            curElemVars.bind(element, elemDisc, curSol);
            if (!this->assembler().isStationaryProblem())
                prevElemVars.bindElement(element, elemDisc, this->assembler().prevSol());
        }
        else
        {
            curElemVars.bindElement(element, elemDisc, curSol);
            prevElemVars.bind(element, elemDisc, prevSol);
        }
    }

    const Problem& problem() const
    { return assembler_.problem(); }

    const Assembler& assembler() const
    { return assembler_; }

    const Element& element() const
    { return element_; }

    bool elementIsGhost() const
    { return elementIsGhost_; }

    const SolutionVector& curSol() const
    { return curSol_; }

    ElementDiscretization& elemDisc()
    { return elemDisc_; }

    ElementVariables& curElemVars()
    { return curElemVars_; }

    ElementVariables& prevElemVars()
    { return prevElemVars_; }

    LocalResidual& localResidual()
    { return localResidual_; }

    const ElementDiscretization& elemDisc() const
    { return elemDisc_; }

    const ElementVariables& curElemVars() const
    { return curElemVars_; }

    const ElementVariables& prevElemVars() const
    { return prevElemVars_; }

    const LocalResidual& localResidual() const
    { return localResidual_; }
 
protected:
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }
 
    const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }
 
private:
 
    const Assembler& assembler_; 
    const Element& element_; 
    const SolutionVector& curSol_; 
 
    ElementDiscretization elemDisc_;
    ElementVariables curElemVars_;
    ElementVariables prevElemVars_;
 
    LocalResidual localResidual_; 
    bool elementIsGhost_; 
};
 
} // end namespace Dumux::Experimental
 
#endif