assembly/assembler.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_ASSEMBLER_HH
#define DUMUX_ASSEMBLER_HH
 
#include <cassert>
#include <iostream>
#include <vector>
#include <deque>
#include <memory>
#include <utility>
 
#include <dune/common/std/type_traits.hh>
#include <dune/grid/common/rangegenerators.hh>
 
#include <dumux/common/exceptions.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/common/gridcapabilities.hh>
#include <dumux/common/typetraits/periodic.hh>
 
#include <dumux/discretization/method.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 "cvfelocalassembler_.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, bool isImplicit>
    using type = Dumux::Experimental::CVFELocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
 
template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
using LocalAssemblerChooser_t = typename LocalAssemblerChooser<
    typename GetPropType<TypeTag, Properties::GridGeometry>::DiscretizationMethod
>::template type<TypeTag, Impl, diffMethod, isImplicit>;

template<class P>
using ProblemConstraintsDetector = decltype(std::declval<P>().constraints());
 
template<class P>
constexpr inline bool hasGlobalConstraints()
{ return Dune::Std::is_detected<ProblemConstraintsDetector, P>::value; }
 
} // end namespace Dumux::Experimental::Detail
 
namespace Dumux::Experimental {

template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true, class LocalResidual = GetPropType<TypeTag, Properties::LocalResidual>>
class Assembler
{
    using GridGeo = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GridGeo::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using ElementSeed = typename GridView::Grid::template Codim<0>::EntitySeed;
    using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
 
    using ThisType = Assembler<TypeTag, diffMethod, isImplicit, LocalResidual>;
    using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod, isImplicit>;
 
public:
    using Scalar = GetPropType<TypeTag, Properties::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>;

    Assembler(std::shared_ptr<const Problem> problem,
              std::shared_ptr<const GridGeometry> gridGeometry,
              std::shared_ptr<GridVariables> gridVariables)
    : problem_(problem)
    , gridGeometry_(gridGeometry)
    , gridVariables_(gridVariables)
    , timeLoop_()
    , isStationaryProblem_(true)
    {
        static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
        enableMultithreading_ = SupportsColoring<typename GridGeometry::DiscretizationMethod>::value
            && Grid::Capabilities::supportsMultithreading(gridGeometry_->gridView())
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }

    Assembler(std::shared_ptr<const Problem> problem,
              std::shared_ptr<const GridGeometry> gridGeometry,
              std::shared_ptr<GridVariables> gridVariables,
              std::shared_ptr<const TimeLoop> timeLoop,
              const SolutionVector& prevSol)
    : problem_(problem)
    , gridGeometry_(gridGeometry)
    , gridVariables_(gridVariables)
    , timeLoop_(timeLoop)
    , prevSol_(&prevSol)
    , isStationaryProblem_(!timeLoop)
    {
        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)
    {
        checkAssemblerState_();
        resetJacobian_(partialReassembler);
        resetResidual_();
 
        assemble_([&](const Element& element)
        {
            LocalAssembler localAssembler(*this, element, curSol);
            localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_, partialReassembler);
        });
 
        enforcePeriodicConstraints_(*jacobian_, *residual_, curSol, *gridGeometry_);
 
        auto applyDirichletConstraint = [&] (const auto& dofIdx,
                                             const auto& values,
                                             const auto eqIdx,
                                             const auto pvIdx)
        {
            (*residual_)[dofIdx][eqIdx] = curSol[dofIdx][pvIdx] - values[pvIdx];
 
            auto& row = (*jacobian_)[dofIdx];
            for (auto col = row.begin(); col != row.end(); ++col)
                row[col.index()][eqIdx] = 0.0;
 
            (*jacobian_)[dofIdx][dofIdx][eqIdx][pvIdx] = 1.0;
        };
        enforceProblemConstraints_(*problem_, *gridGeometry_, applyDirichletConstraint);
    }

    void assembleJacobian(const SolutionVector& curSol)
    {
        checkAssemblerState_();
        resetJacobian_();
 
        assemble_([&](const Element& element)
        {
            LocalAssembler localAssembler(*this, element, curSol);
            localAssembler.assembleJacobian(*jacobian_, *gridVariables_);
        });
 
        enforcePeriodicConstraints_(*jacobian_, curSol, *gridGeometry_);
 
        auto applyDirichletConstraint = [&] (const auto& dofIdx,
                                             const auto& values,
                                             const auto eqIdx,
                                             const auto pvIdx)
        {
            auto& row = (*jacobian_)[dofIdx];
            for (auto col = row.begin(); col != row.end(); ++col)
                row[col.index()][eqIdx] = 0.0;
 
            (*jacobian_)[dofIdx][dofIdx][eqIdx][pvIdx] = 1.0;
        };
        enforceProblemConstraints_(*problem_, *gridGeometry_, applyDirichletConstraint);
    }

    void assembleResidual(const SolutionVector& curSol)
    {
        resetResidual_();
        assembleResidual(*residual_, curSol);
    }

    void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
    {
        checkAssemblerState_();
 
        assemble_([&](const Element& element)
        {
            LocalAssembler localAssembler(*this, element, curSol);
            localAssembler.assembleResidual(r);
        });
 
        enforcePeriodicConstraints_(r, curSol, *gridGeometry_);
 
        auto applyDirichletConstraint = [&] (const auto& dofIdx,
                                             const auto& values,
                                             const auto eqIdx,
                                             const auto pvIdx)
        {
            r[dofIdx][eqIdx] = curSol[dofIdx][pvIdx] - values[pvIdx];
        };
        enforceProblemConstraints_(*problem_, *gridGeometry_, applyDirichletConstraint);
    }

    void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
                         std::shared_ptr<ResidualType> r)
    {
        jacobian_ = A;
        residual_ = r;
 
        // check and/or set the BCRS matrix's build mode
        if (jacobian_->buildMode() == JacobianMatrix::BuildMode::unknown)
            jacobian_->setBuildMode(JacobianMatrix::random);
        else if (jacobian_->buildMode() != JacobianMatrix::BuildMode::random)
            DUNE_THROW(Dune::NotImplemented, "Only BCRS matrices with random build mode are supported at the moment");
 
        setResidualSize_();
        setJacobianPattern_();
    }

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

    void updateAfterGridAdaption()
    {
        setResidualSize_();
        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_; }

    void setTimeLoop(std::shared_ptr<const TimeLoop> timeLoop)
    { timeLoop_ = timeLoop; isStationaryProblem_ = !static_cast<bool>(timeLoop); }

    void setPreviousSolution(const SolutionVector& u)
    { prevSol_ = &u;  }

    bool isStationaryProblem() const
    { return isStationaryProblem_; }

    LocalResidual localResidual() const
    { return LocalResidual(problem_.get(), timeLoop_.get()); }

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

    void resetTimeStep(const SolutionVector& curSol)
    {
        gridVariables().resetTimeStep(curSol);
    }
 
private:
    void setJacobianPattern_()
    {
        // resize the jacobian and the residual
        const auto numDofs = this->numDofs();
        jacobian_->setSize(numDofs, numDofs);
 
        // create occupation pattern of the jacobian
        const auto occupationPattern = getJacobianPattern<isImplicit>(gridGeometry());
 
        // export pattern to jacobian
        occupationPattern.exportIdx(*jacobian_);
    }

    void setResidualSize_()
    { residual_->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_) = 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;
    }
 
    // check if the assembler is in a correct state for assembly
    void checkAssemblerState_() const
    {
        if (!isStationaryProblem_ && !prevSol_)
            DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
    }

    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 (const 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 GG>
    void enforcePeriodicConstraints_(JacobianMatrix& jac, ResidualType& res, const SolutionVector& curSol, const GG& gridGeometry) const
    {
        if constexpr (Dumux::Detail::hasPeriodicDofMap<GG>())
        {
            for (const auto& m : gridGeometry.periodicDofMap())
            {
                if (m.first < m.second)
                {
                    // add the second row to the first
                    res[m.first] += res[m.second];
                    const auto end = jac[m.second].end();
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                        jac[m.first][it.index()] += (*it);
 
                    // enforce constraint in second row
                    res[m.second] = curSol[m.second] - curSol[m.first];
 
                    // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first
                    auto setMatrixBlock = [] (auto& matrixBlock, double diagValue)
                    {
                        for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
                            matrixBlock[eIdx][eIdx] = diagValue;
                    };
 
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                    {
                        auto& matrixBlock = *it;
                        matrixBlock = 0.0;
 
                        assert(matrixBlock.N() == matrixBlock.M());
                        if(it.index() == m.second)
                            setMatrixBlock(matrixBlock, 1.0);
 
                        if(it.index() == m.first)
                            setMatrixBlock(matrixBlock, -1.0);
 
                    }
                }
            }
        }
    }
 
    template<class GG>
    void enforcePeriodicConstraints_(JacobianMatrix& jac, const SolutionVector&, const GG& gridGeometry) const
    {
        if constexpr (Dumux::Detail::hasPeriodicDofMap<GG>())
        {
            for (const auto& m : gridGeometry.periodicDofMap())
            {
                if (m.first < m.second)
                {
                    // add the second row to the first
                    const auto end = jac[m.second].end();
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                        jac[m.first][it.index()] += (*it);
 
                    // set derivatives accordingly in jacobian, i.e. id for m.second and -id for m.first
                    auto setMatrixBlock = [] (auto& matrixBlock, double diagValue)
                    {
                        for (int eIdx = 0; eIdx < matrixBlock.N(); ++eIdx)
                            matrixBlock[eIdx][eIdx] = diagValue;
                    };
 
                    for (auto it = jac[m.second].begin(); it != end; ++it)
                    {
                        auto& matrixBlock = *it;
                        matrixBlock = 0.0;
 
                        assert(matrixBlock.N() == matrixBlock.M());
                        if(it.index() == m.second)
                            setMatrixBlock(matrixBlock, 1.0);
 
                        if(it.index() == m.first)
                            setMatrixBlock(matrixBlock, -1.0);
 
                    }
                }
            }
        }
    }
 
    template<class GG>
    void enforcePeriodicConstraints_(ResidualType& res, const SolutionVector& curSol, const GG& gridGeometry) const
    {
        if constexpr (Dumux::Detail::hasPeriodicDofMap<GG>())
        {
            for (const auto& m : gridGeometry.periodicDofMap())
            {
                if (m.first < m.second)
                {
                    // add the second row to the first
                    res[m.first] += res[m.second];
 
                    // enforce constraint in second row
                    res[m.second] = curSol[m.second] - curSol[m.first];
                }
            }
        }
    }
 
    template<class Problem, class GG, typename ApplyFunction>
    void enforceProblemConstraints_(const Problem& problem, const GG&, const ApplyFunction& applyDirichletConstraint) const
    {
        if constexpr (Detail::hasGlobalConstraints<Problem>())
        {
            for (const auto& constraintData : problem.constraints())
            {
                const auto& constraintInfo = constraintData.constraintInfo();
                const auto& values = constraintData.values();
                const auto dofIdx = constraintData.dofIndex();
                // set the constraint in residual and jacobian
                for (int eqIdx = 0; eqIdx < constraintInfo.size(); ++eqIdx)
                {
                    if (constraintInfo.isConstraintEquation(eqIdx))
                    {
                        const auto pvIdx = constraintInfo.eqToPriVarIndex(eqIdx);
                        assert(0 <= pvIdx && pvIdx < constraintInfo.size());
                        applyDirichletConstraint(dofIdx, values, eqIdx, pvIdx);
                    }
                }
            }
        }
    }

    std::shared_ptr<const Problem> problem_;

    std::shared_ptr<const GridGeometry> gridGeometry_;

    std::shared_ptr<GridVariables> gridVariables_;

    std::shared_ptr<const TimeLoop> timeLoop_;

    const SolutionVector* prevSol_ = nullptr;

    bool isStationaryProblem_;

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

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