assembly/fvassembler.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_FV_ASSEMBLER_HH
#define DUMUX_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/timeloop.hh>
#include <dumux/discretization/method.hh>
#include <dumux/linear/parallelhelpers.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 "boxlocalassembler.hh"
#include "cclocalassembler.hh"
#include "fclocalassembler.hh"
 
namespace Dumux::Detail {
 
template<class DiscretizationMethod>
struct LocalAssemblerChooser;
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::Box>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
    using type = BoxLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::CCMpfa>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
    using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::CCTpfa>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
    using type = CCLocalAssembler<TypeTag, Impl, diffMethod, isImplicit>;
};
 
template<>
struct LocalAssemblerChooser<DiscretizationMethods::FCStaggered>
{
    template<class TypeTag, class Impl, DiffMethod diffMethod, bool isImplicit>
    using type = FaceCenteredLocalAssembler<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>;
 
} // end namespace Dumux::Detail
 
namespace Dumux {
 
template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
class FVAssembler
{
    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;
    using TimeLoop = TimeLoopBase<GetPropType<TypeTag, Properties::Scalar>>;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
 
    static constexpr bool isBox = GridGeo::discMethod == DiscretizationMethods::box;
 
    using ThisType = FVAssembler<TypeTag, diffMethod, isImplicit>;
    using LocalAssembler = typename Detail::LocalAssemblerChooser_t<TypeTag, ThisType, diffMethod, isImplicit>;
 
public:
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
    using GridGeometry = GridGeo;
    using Problem = GetPropType<TypeTag, Properties::Problem>;
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
 
    using ResidualType = SolutionVector;
 
    FVAssembler(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
            && !Multithreading::isSerial()
            && getParam<bool>("Assembly.Multithreading", true);
 
        maybeComputeColors_();
    }
 
    FVAssembler(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
            && !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_);
    }
 
    void assembleJacobian(const SolutionVector& curSol)
    {
        checkAssemblerState_();
        resetJacobian_();
 
        assemble_([&](const Element& element)
        {
            LocalAssembler localAssembler(*this, element, curSol);
            localAssembler.assembleJacobianAndResidual(*jacobian_, *gridVariables_);
        });
    }
 
    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);
        });
    }
 
    Scalar residualNorm(const SolutionVector& curSol) const
    {
        ResidualType residual(numDofs());
        assembleResidual(residual, curSol);
 
        // issue a warning if the caluclation is used in parallel with overlap
        static bool warningIssued = false;
 
        if (gridView().comm().size() > 1 && gridView().overlapSize(0) == 0)
        {
            if constexpr (isBox)
            {
                using DM = typename GridGeometry::VertexMapper;
                using PVHelper = ParallelVectorHelper<GridView, DM, GridView::dimension>;
 
                PVHelper vectorHelper(gridView(), gridGeometry_->vertexMapper());
 
                vectorHelper.makeNonOverlappingConsistent(residual);
            }
        }
        else if (!warningIssued)
        {
            if (gridView().comm().size() > 1 && gridView().comm().rank() == 0)
                std::cout << "\nWarning: norm calculation adds entries corresponding to\n"
                << "overlapping entities multiple times. Please use the norm\n"
                << "function provided by a linear solver instead." << std::endl;
 
            warningIssued = true;
        }
 
        // calculate the square norm of the residual
        Scalar result2 = residual.two_norm2();
        if (gridView().comm().size() > 1)
            result2 = gridView().comm().sum(result2);
 
        using std::sqrt;
        return sqrt(result2);
    }
 
    void setLinearSystem(std::shared_ptr<JacobianMatrix> A,
                         std::shared_ptr<SolutionVector> 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<SolutionVector>();
 
        setResidualSize_();
        setJacobianPattern_();
    }
 
    void updateAfterGridAdaption()
    {
        setResidualSize_();
        setJacobianPattern_();
        maybeComputeColors_();
    }
 
    [[deprecated("Use updateAfterGridAdaption. Will be removed after release 3.5.")]]
    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_);
 
        // maybe recompute colors
        maybeComputeColors_();
    }
 
    [[deprecated("Use updateAfterGridAdaption. Will be removed after release 3.5.")]]
    void setResidualSize()
    { residual_->resize(numDofs()); }
 
    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_; }
 
    SolutionVector& 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<SolutionVector>();
            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 ocurred
        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 GG> std::enable_if_t<GG::discMethod == DiscretizationMethods::box, void>
    enforcePeriodicConstraints_(JacobianMatrix& jac, SolutionVector& res, const SolutionVector& curSol, const GG& gridGeometry)
    {
        for (const auto& m : gridGeometry.periodicVertexMap())
        {
            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];
                for (auto it = jac[m.second].begin(); it != end; ++it)
                    (*it) = it.index() == m.second ? 1.0 : it.index() == m.first ? -1.0 : 0.0;
            }
        }
    }
 
    template<class GG> std::enable_if_t<GG::discMethod != DiscretizationMethods::box, void>
    enforcePeriodicConstraints_(JacobianMatrix& jac, SolutionVector& res, const SolutionVector& curSol, const GG& gridGeometry) {}
 
    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<SolutionVector> residual_;
 
    bool enableMultithreading_ = false;
    std::deque<std::vector<ElementSeed>> elementSets_;
};
 
} // namespace Dumux
 
#endif