compositional/primaryvariableswitch.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_PRIMARY_VARIABLE_SWITCH_HH
#define DUMUX_PRIMARY_VARIABLE_SWITCH_HH
 
#include <iostream>
 
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dumux/discretization/method.hh>
#include <dumux/discretization/elementsolution.hh>
 
namespace Dumux {
 
class NoPrimaryVariableSwitch
{
public:
    template<typename... Args>
    NoPrimaryVariableSwitch(Args&&...) {}
 
    template<typename... Args> void reset(Args&&...) {}
    template<typename... Args> bool wasSwitched(Args&&...) const { return false; }
    template<typename... Args> bool update(Args&&...) { return false; }
    template<typename... Args> void updateSwitchedVolVars(Args&&...) {}
    template<typename... Args> void updateSwitchedFluxVarsCache(Args&&...) {}
    template<typename... Args> void updateDirichletConstraints(Args&&...) {}
};
 
template<class Implementation>
class PrimaryVariableSwitch
{
public:
    PrimaryVariableSwitch(int verbosity = 1)
    : verbosity_(verbosity)
    {}
 
    bool wasSwitched(std::size_t dofIdxGlobal) const
    {
        return wasSwitched_[dofIdxGlobal];
    }
 
    void reset(const std::size_t numDofs)
    {
        wasSwitched_.assign(numDofs, false);
    }
 
    template<class SolutionVector, class GridVariables, class Problem>
    bool update(SolutionVector& curSol,
                GridVariables& gridVariables,
                const Problem& problem,
                const typename GridVariables::GridGeometry& gridGeometry)
    {
        bool switched = false;
        bool successfulUpdate = false;
        const auto& comm = gridGeometry.gridView().comm();
 
        try {
            visited_.assign(wasSwitched_.size(), false);
            std::size_t countSwitched = 0;
 
            auto fvGeometry = localView(gridGeometry);
            auto elemVolVars = localView(gridVariables.curGridVolVars());
 
            for (const auto& element : elements(gridGeometry.gridView()))
            {
                // make sure FVElementGeometry is bound to the element
                fvGeometry.bindElement(element);
                elemVolVars.bindElement(element, fvGeometry, curSol);
 
                const auto curElemSol = elementSolution(element, curSol, gridGeometry);
                for (auto&& scv : scvs(fvGeometry))
                {
                    if (!asImp_().skipDof_(element, fvGeometry, scv, problem))
                    {
                        const auto dofIdxGlobal = scv.dofIndex();
                        // Note this implies that volume variables don't differ
                        // in any sub control volume associated with the dof!
                        visited_[dofIdxGlobal] = true;
                        // Compute volVars on which grounds we decide
                        // if we need to switch the primary variables
                        auto& volVars = getVolVarAccess_(gridVariables.curGridVolVars(), elemVolVars, scv);
                        volVars.update(curElemSol, problem, element, scv);
 
                        if (asImp_().update_(curSol[dofIdxGlobal], volVars, dofIdxGlobal, scv.dofPosition()))
                        {
                            switched = true;
                            ++countSwitched;
                        }
                    }
                }
            }
 
            if (verbosity_ > 0 && countSwitched > 0)
                std::cout << "Switched primary variables at " << countSwitched << " dof locations on processor "
                        << comm.rank() << "." << std::endl;
 
            successfulUpdate = true;
        }
        catch (const Dumux::NumericalProblem& e)
        {
            if (verbosity_ > 0)
                std::cout << "Primary variable switch caught: \"" << e.what() << "\"\n";
 
            successfulUpdate = false;
        }
 
        // make sure all processes were successful
        int successfulUpdateRemote = static_cast<int>(successfulUpdate);
        if (comm.size() > 1)
            successfulUpdateRemote = comm.min(static_cast<int>(successfulUpdate));
 
        if (!successfulUpdate)
            DUNE_THROW(NumericalProblem, "Primary variable switch caught exception on process " << comm.rank());
 
        else if (!successfulUpdateRemote)
            DUNE_THROW(NumericalProblem, "Primary variable switch caught exception on a remote process");
 
        // make sure that if there was a variable switch in an
        // other partition we will also set the switch flag for our partition.
        if (comm.size() > 1)
            switched = comm.max(switched);
 
        return switched;
    }
 
    template<class Problem, class GridVariables, class SolutionVector>
    void updateSwitchedVolVars(const Problem& problem,
                               const typename GridVariables::GridGeometry::GridView::template Codim<0>::Entity& element,
                               const typename GridVariables::GridGeometry& gridGeometry,
                               GridVariables& gridVariables,
                               const SolutionVector& sol)
    {
        if constexpr (GridVariables::GridVolumeVariables::cachingEnabled)
        {
            // make sure FVElementGeometry is bound to the element
            const auto fvGeometry = localView(gridGeometry).bindElement(element);
 
            // update the secondary variables if global caching is enabled
            for (auto&& scv : scvs(fvGeometry))
            {
                const auto dofIdxGlobal = scv.dofIndex();
                if (asImp_().wasSwitched(dofIdxGlobal))
                {
                    const auto elemSol = elementSolution(element, sol, gridGeometry);
                    auto& volVars = gridVariables.curGridVolVars().volVars(scv);
                    volVars.update(elemSol, problem, element, scv);
                }
            }
        }
    }
 
    template<class Problem, class GridVariables, class SolutionVector>
    void updateSwitchedFluxVarsCache(const Problem& problem,
                                     const typename GridVariables::GridGeometry::GridView::template Codim<0>::Entity& element,
                                     const typename GridVariables::GridGeometry& gridGeometry,
                                     GridVariables& gridVariables,
                                     const SolutionVector& sol)
    {
        if constexpr (GridVariables::GridFluxVariablesCache::cachingEnabled
                      && GridVariables::GridGeometry::discMethod != DiscretizationMethods::box)
        {
            // update the flux variables if global caching is enabled
            const auto dofIdxGlobal = gridGeometry.dofMapper().index(element);
 
            if (asImp_().wasSwitched(dofIdxGlobal))
            {
                // make sure FVElementGeometry and the volume variables are bound
                const auto fvGeometry = localView(gridGeometry).bind(element);
                const auto curElemVolVars = localView(gridVariables.curGridVolVars()).bind(element, fvGeometry, sol);
                gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
            }
        }
    }
 
    template<class Problem, class GridVariables, class SolutionVector>
    void updateDirichletConstraints(const Problem& problem,
                                    const typename GridVariables::GridGeometry& gridGeometry,
                                    GridVariables& gridVariables,
                                    SolutionVector& sol)
    {
        if constexpr (GridVariables::GridGeometry::discMethod == DiscretizationMethods::box || Problem::enableInternalDirichletConstraints())
        {
            std::vector<bool> stateChanged(sol.size(), false);
            std::size_t countChanged = 0;
 
            for (const auto& element : elements(gridGeometry.gridView()))
            {
                const auto fvGeometry = localView(gridGeometry).bindElement(element);
 
                // skip if the element is not at a boundary or if no internal Dirichlet constraints are set
                if (!Problem::enableInternalDirichletConstraints() && !fvGeometry.hasBoundaryScvf())
                    continue;
 
                const auto elemVolVars = localView(gridVariables.curGridVolVars()).bindElement(element, fvGeometry, sol);
 
                for (const auto& scv : scvs(fvGeometry))
                {
                    const auto dofIdx = scv.dofIndex();
 
                    if (stateChanged[dofIdx])
                        continue;
 
                    if constexpr (GridVariables::GridGeometry::discMethod == DiscretizationMethods::box)
                    {
                        if (gridGeometry.dofOnBoundary(dofIdx))
                        {
                            if (handleDirichletBoundaryCondition_(problem, element, scv, sol))
                            {
                                stateChanged[dofIdx] = true;
                                ++countChanged;
                                continue;
                            }
                        }
                    }
 
                    if constexpr (Problem::enableInternalDirichletConstraints())
                    {
                        if (handleInternalDirichletConstraint_(problem, element, scv, sol))
                        {
                                stateChanged[dofIdx] = true;
                                ++countChanged;
                        }
                    }
                }
 
                // update the volVars if caching is enabled
                if constexpr (GridVariables::GridVolumeVariables::cachingEnabled)
                {
                    if (countChanged > 0)
                    {
                        const auto curElemSol = elementSolution(element, sol, gridGeometry);
                        for (const auto& scv : scvs(fvGeometry))
                        {
                            if (stateChanged[scv.dofIndex()])
                            {
                                auto& volVars = getVolVarAccess_(gridVariables.curGridVolVars(), elemVolVars, scv);
                                volVars.update(curElemSol, problem, element, scv);
                            }
                        }
                    }
                }
            }
 
            if (verbosity_ > 0 && countChanged > 0)
                std::cout << "Changed primary variable states and solution values to Dirichlet states and values at " << countChanged << " dof locations on processor "
                          << gridGeometry.gridView().comm().rank() << "." << std::endl;
        }
    }
 
    int verbosity() const
    { return verbosity_; }
 
protected:
 
    Implementation &asImp_()
    { return *static_cast<Implementation*>(this); }
 
    const Implementation &asImp_() const
    { return *static_cast<const Implementation*>(this); }
 
    // Perform variable switch at a degree of freedom location
    template<class VolumeVariables, class GlobalPosition>
    bool update_(typename VolumeVariables::PrimaryVariables& priVars,
                 const VolumeVariables& volVars,
                 std::size_t dofIdxGlobal,
                 const GlobalPosition& globalPos)
    {
        // evaluate if the primary variable switch would switch
        // to be implemented by the deriving class
        DUNE_THROW(Dune::NotImplemented, "This model seems to use a primary variable switch but none is implemented!");
    }
 
    // Maybe skip the degree of freedom (do not switch variables at this dof)
    template<class Geometry, class Problem>
    bool skipDof_(const typename Geometry::GridGeometry::GridView::template Codim<0>::Entity& element,
                  const Geometry& fvGeometry,
                  const typename Geometry::SubControlVolume& scv,
                  const Problem& problem)
    {
        if (visited_[scv.dofIndex()])
            return true;
 
        if (isConstrainedDof_(element, fvGeometry, scv, problem))
            return true;
 
        return false;
    }
 
    template<class Geometry, class Problem>
    bool isConstrainedDof_(const typename Geometry::GridGeometry::GridView::template Codim<0>::Entity& element,
                           const Geometry& fvGeometry,
                           const typename Geometry::SubControlVolume& scv,
                           const Problem& problem)
    {
        // Dofs can be only constrained when using the Box method or when imposing internal Dirichlet constraints
        if constexpr (Geometry::GridGeometry::discMethod != DiscretizationMethods::box && !Problem::enableInternalDirichletConstraints())
            return false;
 
        // check for internally constrained Dofs
        const auto isInternallyConstrainedDof = [&]()
        {
            if constexpr (!Problem::enableInternalDirichletConstraints())
                return false;
            else
            {
                const auto internalDirichletConstraints = problem.hasInternalDirichletConstraint(element, scv);
                return internalDirichletConstraints.any();
            }
        };
 
        if (isInternallyConstrainedDof())
            return true;
 
        // check for a Dirichlet BC when using the Box method
        if constexpr (Geometry::GridGeometry::discMethod == DiscretizationMethods::box)
        {
            if (!fvGeometry.hasBoundaryScvf())
                return false;
 
            const auto dofIdx = scv.dofIndex();
            if (!fvGeometry.gridGeometry().dofOnBoundary(dofIdx))
                return false;
 
            const auto bcTypes = problem.boundaryTypes(element, scv);
            if (bcTypes.hasDirichlet())
                return true;
        }
 
        return false;
    }
 
    template<class Problem, class Element, class SubControlVolume, class SolutionVector>
    bool handleDirichletBoundaryCondition_(const Problem& problem,
                                           const Element& element,
                                           const SubControlVolume& scv,
                                           SolutionVector& sol)
    {
        bool changed = false;
        const auto bcTypes = problem.boundaryTypes(element, scv);
        if (bcTypes.hasDirichlet())
        {
            const auto dirichletValues = problem.dirichlet(element, scv);
            const auto dofIdx = scv.dofIndex();
 
            if (sol[dofIdx].state() != dirichletValues.state())
            {
                if (verbosity() > 1)
                    std::cout << "Changing primary variable state at boundary (" << sol[dofIdx].state()
                                << ") to the one given by the Dirichlet condition (" << dirichletValues.state() << ") at dof " << dofIdx
                                << ", coordinates: " << scv.dofPosition()
                                << std::endl;
 
                // make sure the solution vector has the right state (given by the Dirichlet BC)
                sol[dofIdx].setState(dirichletValues.state());
                changed = true;
 
                // overwrite initial with Dirichlet values
                for (int eqIdx = 0; eqIdx < SolutionVector::block_type::dimension; ++eqIdx)
                {
                    if (bcTypes.isDirichlet(eqIdx))
                    {
                        const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
                        sol[dofIdx][pvIdx] = dirichletValues[pvIdx];
                    }
                }
            }
        }
 
        return changed;
    }
 
    template<class Problem, class Element, class SubControlVolume, class SolutionVector>
    bool handleInternalDirichletConstraint_(const Problem& problem,
                                            const Element& element,
                                            const SubControlVolume& scv,
                                            SolutionVector& sol)
    {
        bool changed = false;
 
        const auto internalDirichletConstraints = problem.hasInternalDirichletConstraint(element, scv);
        if (internalDirichletConstraints.none())
            return changed;
 
        const auto dirichletValues = problem.internalDirichlet(element, scv);
        const auto dofIdx = scv.dofIndex();
 
        if (sol[dofIdx].state() != dirichletValues.state())
        {
            if (verbosity() > 1)
                std::cout << "Changing primary variable state at internal DOF (" << sol[dofIdx].state()
                            << ") to the one given by the internal Dirichlet constraint (" << dirichletValues.state() << ") at dof " << dofIdx
                            << ", coordinates: " << scv.dofPosition()
                            << std::endl;
 
            // make sure the solution vector has the right state (given by the Dirichlet constraint)
            sol[dofIdx].setState(dirichletValues.state());
            changed = true;
 
            // overwrite initial with Dirichlet values
            for (int pvIdx = 0; pvIdx < SolutionVector::block_type::dimension; ++pvIdx)
            {
                if (internalDirichletConstraints[pvIdx])
                    sol[dofIdx][pvIdx] = dirichletValues[pvIdx];
            }
        }
 
        return changed;
    }
 
    std::vector<bool> wasSwitched_;
    std::vector<bool> visited_;
 
private:
    template<class GridVolumeVariables, class ElementVolumeVariables, class SubControlVolume>
    typename GridVolumeVariables::VolumeVariables&
    getVolVarAccess_(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv) const
    {
        if constexpr (GridVolumeVariables::cachingEnabled)
            return gridVolVars.volVars(scv);
        else
            return elemVolVars[scv];
    }
 
    int verbosity_; 
};
 
} // end namespace dumux
 
#endif