multistagefvassembler.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_EXPERIMENTAL_MULTISTAGE_FV_ASSEMBLER_HH
#define DUMUX_EXPERIMENTAL_MULTISTAGE_FV_ASSEMBLER_HH
 
#include <vector>
#include <deque>
#include <type_traits>
#include <memory>
 
#include <dune/istl/matrixindexset.hh>
 
#include <dumux/common/properties.hh>
#include <dumux/common/gridcapabilities.hh>
#include <dumux/common/typetraits/vector.hh>
 
#include <dumux/discretization/method.hh>
#include <dumux/linear/parallelhelpers.hh>
#include <dumux/linear/dunevectors.hh>
 
#include <dumux/assembly/coloring.hh>
#include <dumux/assembly/jacobianpattern.hh>
#include <dumux/assembly/diffmethod.hh>
 
#include <dumux/parallel/multithreading.hh>
#include <dumux/parallel/parallel_for.hh>
 
#include <dumux/experimental/assembly/cvfelocalassembler.hh>
#include <dumux/experimental/assembly/cclocalassembler.hh>
#include <dumux/experimental/assembly/multistagefvlocaloperator.hh>
 
#include <dumux/experimental/timestepping/multistagemethods.hh>
#include <dumux/experimental/timestepping/multistagetimestepper.hh>
 
namespace Dumux::Experimental::Detail {
 
template<class DiscretizationMethod>
struct LocalAssemblerChooser;
 
template<class DM>
struct LocalAssemblerChooser<DiscretizationMethods::CVFE<DM>>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod>
    using type = Experimental::CVFELocalAssembler<TypeTag, Impl, diffMethod>;
};
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::CCMpfa>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod>
    using type = Experimental::CCLocalAssembler<TypeTag, Impl, diffMethod>;
};
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::CCTpfa>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod>
    using type = Experimental::CCLocalAssembler<TypeTag, Impl, diffMethod>;
};
 
template<class TypeTag, class Impl, DiffMethod diffMethod>
using LocalAssemblerChooser_t = typename LocalAssemblerChooser<
    typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
>::template type<TypeTag, Impl, diffMethod>;
 
} // end namespace Dumux::Detail
 
namespace Dumux::Experimental {

template<class TypeTag, DiffMethod diffMethod>
class MultiStageFVAssembler
{
    using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GridGeo::GridView;
    using LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>;
    using Element = typename GridView::template Codim<0>::Entity;
    using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed;
 
    static constexpr bool isBox = GridGeo::discMethod == DiscretizationMethods::box;
 
    using ThisType = MultiStageFVAssembler<TypeTag, diffMethod>;
    using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod>;
 
public:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using StageParams = Experimental::MultiStageParams<Scalar>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
    using ResidualType = typename Dumux::Detail::NativeDuneVectorType<SolutionVector>::type;
 
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
 
    using GridGeometry = GridGeo;
    using Problem = GetPropType<TypeTag, Properties::Problem>;

    MultiStageFVAssembler(std::shared_ptr<const Problem> problem,
                          std::shared_ptr<const GridGeometry> gridGeometry,
                          std::shared_ptr<GridVariables> gridVariables,
                          std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> msMethod,
                          const SolutionVector& prevSol)
    : timeSteppingMethod_(msMethod)
    , problem_(problem)
    , gridGeometry_(gridGeometry)
    , gridVariables_(gridVariables)
    , prevSol_(&prevSol)
    {
        enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
            && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }

    template<class PartialReassembler = DefaultPartialReassembler>
    void assembleJacobianAndResidual(const SolutionVector& curSol, const PartialReassembler* partialReassembler = nullptr)
    {
        resetJacobian_(partialReassembler);
 
        resetResidual_();
        spatialOperatorEvaluations_.back() = 0.0;
        temporalOperatorEvaluations_.back() = 0.0;
 
        if (stageParams_->size() != spatialOperatorEvaluations_.size())
            DUNE_THROW(Dune::InvalidStateException, "Wrong number of residuals");
 
        assemble_([&](const Element& element)
        {
            LocalAssembler localAssembler(*this, element, curSol);
            localAssembler.assembleJacobianAndResidual(
                *jacobian_, *residual_, *gridVariables_,
                *stageParams_,
                temporalOperatorEvaluations_.back(),
                spatialOperatorEvaluations_.back(),
                constrainedDofs_,
                partialReassembler
            );
        });
 
        // Autonomous spatial operator: finalize the previous stage's stored evaluation for
        // free. On the first assembly of this stage, curSol is the converged previous-stage
        // solution and the value just assembled into the current slot equals R/M at that
        // solution; time-independence makes the stage time irrelevant, so a copy suffices
        // (no extra evaluation, no dependence on the solver re-assembling at convergence).
        const auto curStage = stageParams_->size() - 1;
        if (curStage > 1 && autonomousSpatialOperator_ && firstStageAssembly_)
        {
            spatialOperatorEvaluations_[curStage-1] = spatialOperatorEvaluations_.back();
            temporalOperatorEvaluations_[curStage-1] = temporalOperatorEvaluations_.back();
        }
        firstStageAssembly_ = false;
 
        // assemble the full residual for the time integration stage
        auto constantResidualComponent = (*residual_);
        constantResidualComponent = 0.0;
        for (std::size_t k = 0; k < stageParams_->size()-1; ++k)
        {
            if (!stageParams_->skipTemporal(k))
                constantResidualComponent.axpy(stageParams_->temporalWeight(k), temporalOperatorEvaluations_[k]);
            if (!stageParams_->skipSpatial(k))
                constantResidualComponent.axpy(stageParams_->spatialWeight(k), spatialOperatorEvaluations_[k]);
        }
 
        // masked summation of constant residual component onto this stage's resiudal component
        for (std::size_t i = 0; i < constantResidualComponent.size(); ++i)
            for (std::size_t ii = 0; ii < constantResidualComponent[i].size(); ++ii)
                (*residual_)[i][ii] += constrainedDofs_[i][ii] > 0.5 ? 0.0 : constantResidualComponent[i][ii];
    }

    void assembleResidual(const SolutionVector& curSol)
    { DUNE_THROW(Dune::NotImplemented, "residual"); }

    void setLinearSystem()
    {
        jacobian_ = std::make_shared<JacobianMatrix>();
        jacobian_->setBuildMode(JacobianMatrix::random);
        residual_ = std::make_shared<ResidualType>();
 
        setResidualSize_(*residual_);
        setJacobianPattern_();
    }

    void updateAfterGridAdaption()
    {
        setResidualSize_(*residual_);
        setJacobianPattern_();
        maybeComputeColors_();
    }

    std::size_t numDofs() const
    { return gridGeometry_->numDofs(); }

    const Problem& problem() const
    { return *problem_; }

    const GridGeometry& gridGeometry() const
    { return *gridGeometry_; }

    const GridView& gridView() const
    { return gridGeometry().gridView(); }

    GridVariables& gridVariables()
    { return *gridVariables_; }

    const GridVariables& gridVariables() const
    { return *gridVariables_; }

    JacobianMatrix& jacobian()
    { return *jacobian_; }

    ResidualType& residual()
    { return *residual_; }

    const SolutionVector& prevSol() const
    { return *prevSol_; }

    MultiStageFVLocalOperator<LocalResidual> localResidual() const
    { return { LocalResidual{problem_.get(), nullptr} }; }

    void updateGridVariables(const SolutionVector &cursol)
    { gridVariables().update(cursol); }

    void resetTimeStep(const SolutionVector &cursol)
    {
        gridVariables().resetTimeStep(cursol);
        this->clearStages();
    }
 
    void clearStages()
    {
        spatialOperatorEvaluations_.clear();
        temporalOperatorEvaluations_.clear();
        stageParams_.reset();
    }
 
    template<class StageParams>
    void prepareStage(SolutionVector& x, StageParams params)
    {
        stageParams_ = std::move(params);
        const auto curStage = stageParams_->size() - 1;
        firstStageAssembly_ = true;
 
        // in the first stage, also assemble the residual
        // at the previous time level (stage 0 residual)
        if (curStage == 1)
        {
            // update time in variables?
            setProblemTime_(*problem_, stageParams_->timeAtStage(0));
 
            resetResidual_(); // residual resized and zero
 
            assert(spatialOperatorEvaluations_.size() >= 0);
            if (spatialOperatorEvaluations_.size() == 0)
            {
                spatialOperatorEvaluations_.push_back(*residual_);
                temporalOperatorEvaluations_.push_back(*residual_);
 
                // assemble stage 0 residuals
                assemble_([&](const auto& element)
                {
                    LocalAssembler localAssembler(*this, element, *prevSol_);
                    localAssembler.localResidual().spatialWeight(1.0);
                    localAssembler.localResidual().temporalWeight(1.0);
                    localAssembler.assembleCurrentResidual(spatialOperatorEvaluations_.back(), temporalOperatorEvaluations_.back());
                });
            }
 
            // we don't delete the first stage so it can be reused in a restarted
            // time integration step. The evaluations are only deleted
            // when explicitly requested by calling clearStages().
            // So if here the vector is non-empty, we don't need to evaluate again
            // (this should only occur if we are restarting time integration, e.g.
            // with a different time step size)
            else if (spatialOperatorEvaluations_.size() > 0)
            {
                updateGridVariables(x);
                spatialOperatorEvaluations_.resize(1);
                temporalOperatorEvaluations_.resize(1);
            }
        }
 
        // Non-autonomous spatial operator: the constant residual component reuses each previous
        // stage's operators at its *converged* solution. During the stage solve the residual
        // assembly stored them at the solver's last iterate — correct for solvers that re-assemble
        // at the converged solution (e.g. NewtonSolver), but stale for single-assembly solvers
        // (e.g. LinearPDESolver). So we re-evaluate here at the converged solution and the previous
        // stage's time level (which matters for time-dependent operators), making multi-stage
        // schemes correct for any solver. This is one extra residual (no Jacobian) evaluation per
        // stage; it is skipped entirely for an autonomous operator, where the value is instead
        // reused for free in assembleJacobianAndResidual (see setAutonomousSpatialOperator).
        if (curStage > 1 && !autonomousSpatialOperator_)
        {
            setProblemTime_(*problem_, stageParams_->timeAtStage(curStage-1));
            spatialOperatorEvaluations_.back() = 0.0;
            temporalOperatorEvaluations_.back() = 0.0;
            assemble_([&](const auto& element)
            {
                LocalAssembler localAssembler(*this, element, x);
                localAssembler.localResidual().spatialWeight(1.0);
                localAssembler.localResidual().temporalWeight(1.0);
                localAssembler.assembleCurrentResidual(spatialOperatorEvaluations_.back(), temporalOperatorEvaluations_.back());
            });
        }
 
        // update time in variables?
        setProblemTime_(*problem_, stageParams_->timeAtStage(curStage));
 
        resetResidual_(); // residual resized and zero
 
        if (spatialOperatorEvaluations_.size() != curStage)
            DUNE_THROW(Dune::InvalidStateException,
                "Invalid state. Maybe you forgot to call clearStages()");
 
        // allocate memory for this stage
        spatialOperatorEvaluations_.push_back(*residual_);
        temporalOperatorEvaluations_.push_back(*residual_);
    }

    bool isStationaryProblem() const
    { return false; }
 
    bool isImplicit() const
    { return timeSteppingMethod_->implicit(); }

    void setAutonomousSpatialOperator(bool a)
    { autonomousSpatialOperator_ = a; }
 
private:
    void setJacobianPattern_()
    {
        // resize the jacobian and the residual
        const auto numDofs = this->numDofs();
        jacobian_->setSize(numDofs, numDofs);
 
        // create occupation pattern of the jacobian
        if (timeSteppingMethod_->implicit())
            getJacobianPattern<true>(gridGeometry()).exportIdx(*jacobian_);
        else
            getJacobianPattern<false>(gridGeometry()).exportIdx(*jacobian_);
    }

    void setResidualSize_(ResidualType& res)
    { res.resize(numDofs()); }

    void maybeComputeColors_()
    {
        if (enableMultithreading_)
            elementSets_ = computeColoring(gridGeometry()).sets;
    }
 
    // reset the residual vector to 0.0
    void resetResidual_()
    {
        if(!residual_)
        {
            residual_ = std::make_shared<ResidualType>();
            setResidualSize_(*residual_);
        }
 
        setResidualSize_(constrainedDofs_);
 
        (*residual_) = 0.0;
        constrainedDofs_ = 0.0;
    }
 
    // reset the Jacobian matrix to 0.0
    template <class PartialReassembler = DefaultPartialReassembler>
    void resetJacobian_(const PartialReassembler *partialReassembler = nullptr)
    {
        if(!jacobian_)
        {
            jacobian_ = std::make_shared<JacobianMatrix>();
            jacobian_->setBuildMode(JacobianMatrix::random);
            setJacobianPattern_();
        }
 
        if (partialReassembler)
            partialReassembler->resetJacobian(*this);
        else
            *jacobian_ = 0.0;
    }

    template<typename AssembleElementFunc>
    void assemble_(AssembleElementFunc&& assembleElement) const
    {
        // a state that will be checked on all processes
        bool succeeded = false;
 
        // try assembling using the local assembly function
        try
        {
            if (enableMultithreading_)
            {
                assert(elementSets_.size() > 0);
 
                // 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
                for (const auto& elements : elementSets_)
                {
                    Dumux::parallelFor(elements.size(), [&](const std::size_t i)
                    {
                        const auto element = gridView().grid().entity(elements[i]);
                        assembleElement(element);
                    });
                }
            }
            else
                for (const auto& element : elements(gridView()))
                    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().comm().rank()
                      << " caught an exception while assembling:" << e.what()
                      << "\n";
            succeeded = false;
        }
 
        // make sure everything worked well on all processes
        if (gridView().comm().size() > 1)
            succeeded = gridView().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");
    }
 
    template<class P>
    void setProblemTime_(const P& p, const Scalar t)
    {
        if constexpr (requires { p.setTime(t); })
            p.setTime(t);
        else
            static_assert(!requires (P& q) { q.setTime(Scalar{}); },
                "The multi-stage assembler sets the stage time through a const problem: setTime has to be const.");
    }
 
    std::shared_ptr<const Experimental::MultiStageMethod<Scalar>> timeSteppingMethod_;
    std::vector<ResidualType> spatialOperatorEvaluations_;
    std::vector<ResidualType> temporalOperatorEvaluations_;
    ResidualType constrainedDofs_;
    std::shared_ptr<const StageParams> stageParams_;

    std::shared_ptr<const Problem> problem_;

    std::shared_ptr<const GridGeometry> gridGeometry_;

    std::shared_ptr<GridVariables> gridVariables_;

    const SolutionVector* prevSol_ = nullptr;

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

    bool autonomousSpatialOperator_ = false;

    bool firstStageAssembly_ = false;

    bool enableMultithreading_ = false;
    std::deque<std::vector<ElementSeed>> elementSets_;
};
 
} // namespace Dumux
 
#endif