fvbasegradientcalculator.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:
/*
  This file is part of the Open Porous Media project (OPM).

  OPM 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 2 of the License, or
  (at your option) any later version.

  OPM 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 OPM.  If not, see <http://www.gnu.org/licenses/>.

  Consult the COPYING file in the top-level source directory of this
  module for the precise wording of the license and the list of
  copyright holders.
*/
#ifndef EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH
#define EWOMS_FV_BASE_GRADIENT_CALCULATOR_HH

#include <dune/common/fvector.hh>

#include <opm/models/discretization/common/fvbaseproperties.hh>
#include <opm/models/utils/basicproperties.hh>

#include <cmath>
#include <type_traits>

namespace Opm {

template<class TypeTag>
class EcfvDiscretization;

template<class TypeTag>
class FvBaseGradientCalculator
{
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using Discretization = GetPropType<TypeTag, Properties::Discretization>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;

    enum { dim = GridView::dimension };
    enum { dimWorld = GridView::dimensionworld };

    // maximum number of flux approximation points. to calculate this,
    // we assume that the geometry with the most points is a cube.
    enum { maxFap = 2 << dim };

    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
    using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;

public:
    static void registerParameters()
    {}

    template <bool prepareValues = true, bool prepareGradients = true>
    void prepare(const ElementContext&, unsigned)
    { /* nothing to do */ }

    template <class QuantityCallback>
    auto calculateScalarValue(const ElementContext& elemCtx,
                              unsigned fapIdx,
                              const QuantityCallback& quantityCallback) const
        -> std::remove_reference_t<decltype(quantityCallback.operator()(0))>
    {
        using RawReturnType = decltype(quantityCallback.operator()(0));
        using ReturnType = std::remove_const_t<std::remove_reference_t<RawReturnType>>;

        Scalar interiorDistance;
        Scalar exteriorDistance;
        computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);

        const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
        const auto i = face.interiorIndex();
        const auto j = face.exteriorIndex();
        const auto focusDofIdx = elemCtx.focusDofIndex();

        // use the average weighted by distance...
        ReturnType value;
        if (i == focusDofIdx) {
            value = quantityCallback(i) * interiorDistance;
        } else {
            value = getValue(quantityCallback(i)) * interiorDistance;
        }

        if (j == focusDofIdx) {
            value += quantityCallback(j) * exteriorDistance;
        }
        else {
            value += getValue(quantityCallback(j)) * exteriorDistance;
        }

        value /= interiorDistance + exteriorDistance;

        return value;
    }

    template <class QuantityCallback>
    auto calculateVectorValue(const ElementContext& elemCtx,
                              unsigned fapIdx,
                              const QuantityCallback& quantityCallback) const
        -> std::remove_reference_t<decltype(quantityCallback.operator()(0))>
    {
        using RawReturnType = decltype(quantityCallback.operator()(0));
        using ReturnType = std::remove_const_t<std::remove_reference_t<RawReturnType>>;

        Scalar interiorDistance;
        Scalar exteriorDistance;
        computeDistances_(interiorDistance, exteriorDistance, elemCtx, fapIdx);

        const auto& face = elemCtx.stencil(/*timeIdx=*/0).interiorFace(fapIdx);
        const auto i = face.interiorIndex();
        const auto j = face.exteriorIndex();
        const auto focusDofIdx = elemCtx.focusDofIndex();

        // use the average weighted by distance...
        ReturnType value;
        if (i == focusDofIdx) {
            value = quantityCallback(i);
            for (int k = 0; k < value.size(); ++k) {
                value[k] *= interiorDistance;
            }
        }
        else {
            const auto& dofVal = getValue(quantityCallback(i));
            for (int k = 0; k < dofVal.size(); ++k) {
                value[k] = getValue(dofVal[k]) * interiorDistance;
            }
        }

        if (j == focusDofIdx) {
            const auto& dofVal = quantityCallback(j);
            for (int k = 0; k < dofVal.size(); ++k) {
                value[k] += dofVal[k] * exteriorDistance;
            }
        }
        else {
            const auto& dofVal = quantityCallback(j);
            for (int k = 0; k < dofVal.size(); ++k) {
                value[k] += getValue(dofVal[k]) * exteriorDistance;
            }
        }

        const Scalar totDistance = interiorDistance + exteriorDistance;
        for (int k = 0; k < value.size(); ++k) {
            value[k] /= totDistance;
        }

        return value;
    }

    template <class QuantityCallback>
    void calculateGradient(EvalDimVector& quantityGrad,
                           const ElementContext& elemCtx,
                           unsigned fapIdx,
                           const QuantityCallback& quantityCallback) const
    {
        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
        const auto& face = stencil.interiorFace(fapIdx);

        const auto i = face.interiorIndex();
        const auto j = face.exteriorIndex();
        const auto focusIdx = elemCtx.focusDofIndex();

        const auto& interiorPos = stencil.subControlVolume(i).globalPos();
        const auto& exteriorPos = stencil.subControlVolume(j).globalPos();

        Evaluation deltay;
        if (i == focusIdx) {
            deltay = getValue(quantityCallback(j)) - quantityCallback(i);
        }
        else if (j == focusIdx) {
            deltay = quantityCallback(j) - getValue(quantityCallback(i));
        }
        else {
            deltay = getValue(quantityCallback(j)) - getValue(quantityCallback(i));
        }

        Scalar distSquared = 0.0;
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            const Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
            distSquared += tmp * tmp;
        }

        // divide the gradient by the squared distance between the centers of the
        // sub-control volumes: the gradient is the normalized directional vector between
        // the two centers times the ratio of the difference of the values and their
        // distance, i.e., d/abs(d) * delta y / abs(d) = d*delta y / abs(d)^2.
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            const Scalar tmp = exteriorPos[dimIdx] - interiorPos[dimIdx];
            quantityGrad[dimIdx] = deltay *( tmp / distSquared);
        }
    }

    template <class QuantityCallback>
    auto calculateBoundaryValue(const ElementContext&,
                                unsigned,
                                const QuantityCallback& quantityCallback)
        -> decltype(quantityCallback.boundaryValue())
    { return quantityCallback.boundaryValue(); }

    template <class QuantityCallback>
    void calculateBoundaryGradient(EvalDimVector& quantityGrad,
                                   const ElementContext& elemCtx,
                                   unsigned faceIdx,
                                   const QuantityCallback& quantityCallback) const
    {
        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
        const auto& face = stencil.boundaryFace(faceIdx);

        Evaluation deltay;
        if (face.interiorIndex() == elemCtx.focusDofIndex()) {
            deltay = quantityCallback.boundaryValue() - quantityCallback(face.interiorIndex());
        }
        else {
            deltay = getValue(quantityCallback.boundaryValue()) -
                     getValue(quantityCallback(face.interiorIndex()));
        }

        const auto& boundaryFacePos = face.integrationPos();
        const auto& interiorPos = stencil.subControlVolume(face.interiorIndex()).center();

        Scalar distSquared = 0;
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            const Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
            distSquared += tmp * tmp;
        }

        // divide the gradient by the squared distance between the center of the
        // sub-control and the center of the boundary face: the gradient is the
        // normalized directional vector between the two centers times the ratio of the
        // difference of the values and their distance, i.e., d/abs(d) * deltay / abs(d)
        // = d*deltay / abs(d)^2.
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            const Scalar tmp = boundaryFacePos[dimIdx] - interiorPos[dimIdx];
            quantityGrad[dimIdx] = deltay * (tmp / distSquared);
        }
    }

private:
    void computeDistances_(Scalar& interiorDistance,
                           Scalar& exteriorDistance,
                           const ElementContext& elemCtx,
                           unsigned fapIdx) const
    {
        const auto& stencil = elemCtx.stencil(/*timeIdx=*/0);
        const auto& face = stencil.interiorFace(fapIdx);

        // calculate the distances of the position of the interior and of the exterior
        // finite volume to the position of the integration point.
        const auto& normal = face.normal();
        const auto i = face.interiorIndex();
        const auto j = face.exteriorIndex();
        const auto& interiorPos = stencil.subControlVolume(i).globalPos();
        const auto& exteriorPos = stencil.subControlVolume(j).globalPos();
        const auto& integrationPos = face.integrationPos();

        interiorDistance = 0.0;
        exteriorDistance = 0.0;
        for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
            interiorDistance += (interiorPos[dimIdx] - integrationPos[dimIdx]) * normal[dimIdx];
            exteriorDistance += (exteriorPos[dimIdx] - integrationPos[dimIdx]) * normal[dimIdx];
        }

        interiorDistance = std::sqrt(std::abs(interiorDistance));
        exteriorDistance = std::sqrt(std::abs(exteriorDistance));
    }
};

} // namespace Opm

#endif