gridadaptindicator.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_TWOP_ADAPTION_INDICATOR_HH
#define DUMUX_TWOP_ADAPTION_INDICATOR_HH
 
#include <memory>
#include <dune/common/exceptions.hh>
#include <dune/grid/common/partitionset.hh>
 
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/gridcapabilities.hh>
#include <dumux/discretization/elementsolution.hh>
#include <dumux/discretization/evalsolution.hh>
#include <dumux/parallel/vectorcommdatahandle.hh>
 
namespace Dumux {

template<class TypeTag>
class TwoPGridAdaptIndicator
{
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
    using Element = typename GridView::template Codim<0>::Entity;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
 
    enum { saturationIdx = Indices::saturationIdx };
 
public:
    TwoPGridAdaptIndicator(std::shared_ptr<const GridGeometry> gridGeometry, const std::string& paramGroup = "")
    : gridGeometry_(gridGeometry)
    , refineBound_(std::numeric_limits<Scalar>::max())
    , coarsenBound_(std::numeric_limits<Scalar>::lowest())
    , maxSaturationDelta_(gridGeometry_->gridView().size(0), 0.0)
    , minLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MinLevel", 0))
    , maxLevel_(getParamFromGroup<std::size_t>(paramGroup, "Adaptive.MaxLevel", 0))
    {}

    void setMinLevel(std::size_t minLevel)
    {
        minLevel_ = minLevel;
    }

    void setMaxLevel(std::size_t maxLevel)
    {
        maxLevel_ = maxLevel;
    }

    void setLevels(std::size_t minLevel, std::size_t maxLevel)
    {
        minLevel_ = minLevel;
        maxLevel_ = maxLevel;
    }

    void calculate(const SolutionVector& sol,
                   Scalar refineTol = 0.05,
                   Scalar coarsenTol = 0.001)
    {
        const auto& gridView = gridGeometry_->gridView();

        refineBound_ = std::numeric_limits<Scalar>::max();
        coarsenBound_ = std::numeric_limits<Scalar>::lowest();
        maxSaturationDelta_.assign(gridView.size(0), 0.0);

        if (minLevel_ >= maxLevel_)
            return;

        if (coarsenTol > refineTol)
            DUNE_THROW(Dune::InvalidStateException, "Refine tolerance must be higher than coarsen tolerance");

        Scalar globalMax = std::numeric_limits<Scalar>::lowest();
        Scalar globalMin = std::numeric_limits<Scalar>::max();

        for (const auto& element : elements(gridView, Dune::Partitions::interior))
        {
            const auto globalIdxI = gridGeometry_->elementMapper().index(element);

            const auto geometry = element.geometry();
            const auto elemSol = elementSolution(element, sol, *gridGeometry_);
            const Scalar satI = evalSolution(element, geometry, *gridGeometry_, elemSol, geometry.center())[saturationIdx];

            using std::min;
            using std::max;
            globalMin = min(satI, globalMin);
            globalMax = max(satI, globalMax);

            for (const auto& intersection : intersections(gridView, element))
            {
                if (intersection.neighbor())
                {
                    const auto outside = intersection.outside();
                    const auto globalIdxJ = gridGeometry_->elementMapper().index(outside);

                    if (element.level() > outside.level() || (element.level() == outside.level() && globalIdxI < globalIdxJ))
                    {
                        const auto outsideGeometry = outside.geometry();
                        const auto elemSolJ = elementSolution(outside, sol, *gridGeometry_);
                        const Scalar satJ = evalSolution(outside, outsideGeometry, *gridGeometry_, elemSolJ, outsideGeometry.center())[saturationIdx];
 
                        using std::abs;
                        Scalar localdelta = abs(satI - satJ);
                        maxSaturationDelta_[globalIdxI] = max(maxSaturationDelta_[globalIdxI], localdelta);
                        maxSaturationDelta_[globalIdxJ] = max(maxSaturationDelta_[globalIdxJ], localdelta);
                    }
                }
            }
        }

        const auto globalDelta = globalMax - globalMin;

        refineBound_ = refineTol*globalDelta;
        coarsenBound_ = coarsenTol*globalDelta;
 
#if HAVE_MPI
        // communicate updated values
        using DataHandle = VectorCommDataHandleMax<typename GridGeometry::ElementMapper, std::vector<Scalar>, 0>;
        DataHandle dataHandle(gridGeometry_->elementMapper(), maxSaturationDelta_);
        if constexpr (Dune::Capabilities::canCommunicate<typename GridGeometry::GridView::Traits::Grid, 0>::v)
            gridView.template communicate<DataHandle>(
                dataHandle, Dune::InteriorBorder_All_Interface, Dune::ForwardCommunication
            );
        else
            DUNE_THROW(
                Dune::InvalidStateException,
                "Cannot call calculate indicator on grid that cannot communicate codim-0-entities."
            );
 
        using std::max;
        refineBound_ = gridView.comm().max(refineBound_);
        coarsenBound_ = gridView.comm().max(coarsenBound_);
#endif

        for (const auto& element : elements(gridView))
            if (this->operator()(element) > 0)
                checkNeighborsRefine_(element);
    }

    int operator() (const Element& element) const
    {
        if (element.hasFather()
            && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] < coarsenBound_)
        {
            return -1;
        }
        else if (element.level() < maxLevel_
                 && maxSaturationDelta_[gridGeometry_->elementMapper().index(element)] > refineBound_)
        {
            return 1;
        }
        else
            return 0;
    }
 
private:
    bool checkNeighborsRefine_(const Element &element, std::size_t level = 1)
    {
        for(const auto& intersection : intersections(gridGeometry_->gridView(), element))
        {
            if(!intersection.neighbor())
                continue;
 
            // obtain outside element
            const auto outside = intersection.outside();
 
            // only mark non-ghost elements
            if (outside.partitionType() == Dune::GhostEntity)
                continue;
 
            if (outside.level() < maxLevel_ && outside.level() < element.level())
            {
                // ensure refinement for outside element
                assert(outside.isLeaf());
                maxSaturationDelta_[gridGeometry_->elementMapper().index(outside)]
                    = std::numeric_limits<Scalar>::max();
 
                if (level < maxLevel_)
                    checkNeighborsRefine_(outside, ++level);
            }
        }
 
        return true;
    }
 
    std::shared_ptr<const GridGeometry> gridGeometry_;
 
    Scalar refineBound_;
    Scalar coarsenBound_;
    std::vector< Scalar > maxSaturationDelta_;
    std::size_t minLevel_;
    std::size_t maxLevel_;
};
 
} // end namespace Dumux
 
#endif