NewTranFluxModule.hpp

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 OPM_NEWTRAN_FLUX_MODULE_HPP
#define OPM_NEWTRAN_FLUX_MODULE_HPP

#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>

#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
#include <opm/common/utility/gpuDecorators.hpp>

#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>

#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/thermal/EnergyModuleType.hpp>

#include <opm/models/blackoil/blackoilmoduleparams.hh>
#include <opm/models/blackoil/blackoilmodules.hpp>
#include <opm/models/blackoil/blackoilconvectivemixingmoduleparam.hpp>
#include <opm/models/blackoil/blackoilproperties.hh>

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

#include <opm/models/utils/signum.hh>

#include <array>

namespace Opm {

template <class TypeTag>
class NewTranIntensiveQuantities;

template <class TypeTag>
class NewTranExtensiveQuantities;

template <class TypeTag>
class NewTranBaseProblem;

template <class TypeTag>
struct NewTranFluxModule
{
    using FluxIntensiveQuantities = NewTranIntensiveQuantities<TypeTag>;
    using FluxExtensiveQuantities = NewTranExtensiveQuantities<TypeTag>;
    using FluxBaseProblem = NewTranBaseProblem<TypeTag>;

    static void registerParameters()
    { }
};

template <class TypeTag>
class NewTranBaseProblem
{ };

template <class TypeTag>
class NewTranIntensiveQuantities
{
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
protected:
    void update_(const ElementContext&, unsigned, unsigned)
    { }
};

template <class TypeTag>
class NewTranExtensiveQuantities
{
    using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;

    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using GridView = GetPropType<TypeTag, Properties::GridView>;
    using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;

    enum { dimWorld = GridView::dimensionworld };
    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
    enum { numPhases = FluidSystem::numPhases };

    static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
    static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool enableEnergy =
        getPropValue<TypeTag, Properties::EnergyModuleType>() == EnergyModules::FullyImplicitThermal;
    static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();

    using Toolbox = MathToolbox<Evaluation>;
    using DimVector = Dune::FieldVector<Scalar, dimWorld>;
    using EvalDimVector = Dune::FieldVector<Evaluation, dimWorld>;
    using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;

    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
    using ModuleParams = BlackoilModuleParams<ConvectiveMixingModuleParam<Scalar>>;

public:
    OPM_HOST_DEVICE const DimMatrix& intrinsicPermeability() const
    {
        OPM_THROW(std::invalid_argument, "The ECL transmissibility module does not provide an explicit intrinsic permeability");
    }

    OPM_HOST_DEVICE const EvalDimVector& potentialGrad(unsigned) const
    {
        OPM_THROW(std::invalid_argument, "The ECL transmissibility module does not provide explicit potential gradients");
    }

    OPM_HOST_DEVICE const Evaluation& pressureDifference(unsigned phaseIdx) const
    { return pressureDifference_[phaseIdx]; }

    OPM_HOST_DEVICE const EvalDimVector& filterVelocity(unsigned) const
    {
        OPM_THROW(std::invalid_argument, "The ECL transmissibility module does not provide explicit filter velocities");
    }

    OPM_HOST_DEVICE const Evaluation& volumeFlux(unsigned phaseIdx) const
    { return volumeFlux_[phaseIdx]; }

protected:
    OPM_HOST_DEVICE unsigned upstreamIndex_(unsigned phaseIdx) const
    {
        assert(phaseIdx < numPhases);

        return upIdx_[phaseIdx];
    }

    OPM_HOST_DEVICE unsigned downstreamIndex_(unsigned phaseIdx) const
    {
        assert(phaseIdx < numPhases);

        return dnIdx_[phaseIdx];
    }

    OPM_HOST_DEVICE void updateSolvent(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
    {
        if constexpr (enableSolvent) {
            asImp_().updateVolumeFluxTrans(elemCtx, scvfIdx, timeIdx);
        }
    }

    OPM_HOST_DEVICE void updatePolymer(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
    {
        if constexpr (enablePolymer) {
            asImp_().updateShearMultipliers(elemCtx, scvfIdx, timeIdx);
        }
    }

public:

    OPM_HOST_DEVICE static void volumeAndPhasePressureDifferences(std::array<short, numPhases>& upIdx,
                                                                  std::array<short, numPhases>& dnIdx,
                                                                  Evaluation (&volumeFlux)[numPhases],
                                                                  Evaluation (&pressureDifferences)[numPhases],
                                                                  const ElementContext& elemCtx,
                                                                  unsigned scvfIdx,
                                                                  unsigned timeIdx)
    {
        const auto& problem = elemCtx.problem();
        const auto& stencil = elemCtx.stencil(timeIdx);
        const auto& scvf = stencil.interiorFace(scvfIdx);
        unsigned interiorDofIdx = scvf.interiorIndex();
        unsigned exteriorDofIdx = scvf.exteriorIndex();

        assert(interiorDofIdx != exteriorDofIdx);

        unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
        unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
        Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
        Scalar faceArea = scvf.area();
        Scalar thpresInToEx = problem.thresholdPressure(I, J);
        Scalar thpresExToIn = problem.thresholdPressure(J, I);

        // estimate the gravity correction: for performance reasons we use a simplified
        // approach for this flux module that assumes that gravity is constant and always
        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
        Scalar g = elemCtx.problem().gravity()[dimWorld - 1];

        const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
        const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);

        // this is quite hacky because the dune grid interface does not provide a
        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
        // solution would be to take the Z coordinate of the element centroids, but since
        // ECL seems to like to be inconsistent on that front, it needs to be done like
        // here...
        Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
        Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);

        // the distances from the DOF's depths. (i.e., the additional depth of the
        // exterior DOF)
        Scalar distZ = zIn - zEx;

        Scalar Vin = elemCtx.dofVolume(interiorDofIdx, /*timeIdx=*/0);
        Scalar Vex = elemCtx.dofVolume(exteriorDofIdx, /*timeIdx=*/0);

        for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
            if (!FluidSystem::phaseIsActive(phaseIdx))
                continue;
            calculatePhasePressureDiff_(upIdx[phaseIdx],
                                        dnIdx[phaseIdx],
                                        pressureDifferences[phaseIdx],
                                        intQuantsIn,
                                        intQuantsEx,
                                        phaseIdx,//input
                                        interiorDofIdx,//input
                                        exteriorDofIdx,//input
                                        Vin,
                                        Vex,
                                        I,
                                        J,
                                        distZ*g,
                                        thpresInToEx,
                                        thpresExToIn,
                                        problem.moduleParams());

            const bool upwindIsInterior = (static_cast<unsigned>(upIdx[phaseIdx]) == interiorDofIdx);
            const IntensiveQuantities& up = upwindIsInterior ? intQuantsIn : intQuantsEx;
            // Use arithmetic average (more accurate with harmonic, but that requires recomputing the transmissbility)
            const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;

            const auto& materialLawManager = problem.materialLawManager();
            FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
            if (materialLawManager->hasDirectionalRelperms()) {
                facedir = scvf.faceDirFromDirId();  // direction (X, Y, or Z) of the face
            }
            if (upwindIsInterior)
                volumeFlux[phaseIdx] =
                    pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
            else
                volumeFlux[phaseIdx] =
                    pressureDifferences[phaseIdx]*
                        (Toolbox::value(up.mobility(phaseIdx, facedir))*transMult*(-trans/faceArea));
        }
    }

    template<class EvalType, class ModuleParamsT = ModuleParams>
    OPM_HOST_DEVICE static void calculatePhasePressureDiff_(short& upIdx,
                                                            short& dnIdx,
                                                            EvalType& pressureDifference,
                                                            const IntensiveQuantities& intQuantsIn,
                                                            const IntensiveQuantities& intQuantsEx,
                                                            const unsigned phaseIdx,
                                                            const unsigned interiorDofIdx,
                                                            const unsigned exteriorDofIdx,
                                                            const Scalar Vin,
                                                            const Scalar Vex,
                                                            const unsigned globalIndexIn,
                                                            const unsigned globalIndexEx,
                                                            const Scalar distZg,
                                                            const Scalar thpresInToEx,
                                                            const Scalar thpresExToIn,
                                                            const ModuleParamsT& moduleParams)
    {

        // check shortcut: if the mobility of the phase is zero in the interior as
        // well as the exterior DOF, we can skip looking at the phase.
        if (intQuantsIn.mobility(phaseIdx) <= 0.0 &&
            intQuantsEx.mobility(phaseIdx) <= 0.0)
        {
            upIdx = interiorDofIdx;
            dnIdx = exteriorDofIdx;
            pressureDifference = 0.0;
            return;
        }

        // do the gravity correction: compute the hydrostatic pressure for the
        // external at the depth of the internal one
        const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
        Scalar rhoEx = Toolbox::value(intQuantsEx.fluidState().density(phaseIdx));
        Evaluation rhoAvg = (rhoIn + rhoEx)/2;

        if constexpr (enableConvectiveMixing) {
            ConvectiveMixingModule::modifyAvgDensity(rhoAvg, intQuantsIn, intQuantsEx, phaseIdx, moduleParams.convectiveMixingModuleParam);
        }

        const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
        Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(phaseIdx));
        if (enableExtbo) // added stability; particulary useful for solvent migrating in pure water
                         // where the solvent fraction displays a 0/1 behaviour ...
            pressureExterior += Toolbox::value(rhoAvg)*(distZg);
        else
            pressureExterior += rhoAvg*(distZg);

        pressureDifference = pressureExterior - pressureInterior;

        // decide the upstream index for the phase. for this we make sure that the
        // degree of freedom which is regarded upstream if both pressures are equal
        // is always the same: if the pressure is equal, the DOF with the lower
        // global index is regarded to be the upstream one.
        if (pressureDifference > 0.0) {
            upIdx = exteriorDofIdx;
            dnIdx = interiorDofIdx;
        }
        else if (pressureDifference < 0.0) {
            upIdx = interiorDofIdx;
            dnIdx = exteriorDofIdx;
        }
        else {
            // if the pressure difference is zero, we chose the DOF which has the
            // larger volume associated to it as upstream DOF
            if (Vin > Vex) {
                upIdx = interiorDofIdx;
                dnIdx = exteriorDofIdx;
            }
            else if (Vin < Vex) {
                upIdx = exteriorDofIdx;
                dnIdx = interiorDofIdx;
            }
            else {
                assert(Vin == Vex);
                // if the volumes are also equal, we pick the DOF which exhibits the
                // smaller global index
                if (globalIndexIn < globalIndexEx) {
                    upIdx = interiorDofIdx;
                    dnIdx = exteriorDofIdx;
                }
                else {
                    upIdx  = exteriorDofIdx;
                    dnIdx = interiorDofIdx;
                }
            }
        }

        // apply the threshold pressure for the intersection. note that the concept
        // of threshold pressure is a quite big hack that only makes sense for ECL
        // datasets. (and even there, its physical justification is quite
        // questionable IMO.)
        const Scalar thpres = pressureDifference < 0.0 ? thpresInToEx : thpresExToIn;
        if (thpres > 0.0) {
            if (std::abs(Toolbox::value(pressureDifference)) > thpres) {
                if (pressureDifference < 0.0)
                    pressureDifference += thpres;
                else
                    pressureDifference -= thpres;
            }
            else {
                pressureDifference = 0.0;
            }
        }
    }

protected:
    OPM_HOST_DEVICE void calculateGradients_(const ElementContext& elemCtx, unsigned scvfIdx, unsigned timeIdx)
    {
        Valgrind::SetUndefined(*this);

        volumeAndPhasePressureDifferences(upIdx_ , dnIdx_, volumeFlux_, pressureDifference_, elemCtx, scvfIdx, timeIdx);
    }

    template <class FluidState>
    OPM_HOST_DEVICE void calculateBoundaryGradients_(const ElementContext& elemCtx,
                                                     unsigned scvfIdx,
                                                     unsigned timeIdx,
                                                     const FluidState& exFluidState)
    {
        const auto& scvf = elemCtx.stencil(timeIdx).boundaryFace(scvfIdx);
        const Scalar faceArea = scvf.area();
        const Scalar zEx = scvf.integrationPos()[dimWorld - 1];
        const auto& problem = elemCtx.problem();
        const unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(0, timeIdx);
        const auto& intQuantsIn = elemCtx.intensiveQuantities(0, timeIdx);

        calculateBoundaryGradients_(problem,
                                    globalSpaceIdx,
                                    intQuantsIn,
                                    scvfIdx,
                                    faceArea,
                                    zEx,
                                    exFluidState,
                                    upIdx_,
                                    dnIdx_,
                                    volumeFlux_,
                                    pressureDifference_);

        // Treating solvent here and not in the static method, since that would require more
        // extensive refactoring. It means that the TpfaLinearizer will not support bcs for solvent until this is
        // addressed.
        if constexpr (enableSolvent) {
            if (upIdx_[gasPhaseIdx] == 0) {
                const Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, scvfIdx);
                const Scalar transModified = trans * Toolbox::value(intQuantsIn.rockCompTransMultiplier());
                const auto solventFlux = pressureDifference_[gasPhaseIdx] * intQuantsIn.mobility(gasPhaseIdx) * (-transModified/faceArea);
                asImp_().setSolventVolumeFlux(solventFlux);
            } else {
                asImp_().setSolventVolumeFlux(0.0);
            }
        }
    }

public:
    template <class Problem, class FluidState, class EvaluationContainer>
    OPM_HOST_DEVICE static void calculateBoundaryGradients_(const Problem& problem,
                                                            const unsigned globalSpaceIdx,
                                                            const IntensiveQuantities& intQuantsIn,
                                                            const unsigned bfIdx,
                                                            const double faceArea,
                                                            const double zEx,
                                                            const FluidState& exFluidState,
                                                            std::array<short, numPhases>& upIdx,
                                                            std::array<short, numPhases>& dnIdx,
                                                            EvaluationContainer& volumeFlux,
                                                            EvaluationContainer& pressureDifference)
    {

        bool enableBoundaryMassFlux = problem.nonTrivialBoundaryConditions();
        if (!enableBoundaryMassFlux)
            return;

        Scalar trans = problem.transmissibilityBoundary(globalSpaceIdx, bfIdx);

        // estimate the gravity correction: for performance reasons we use a simplified
        // approach for this flux module that assumes that gravity is constant and always
        // acts into the downwards direction. (i.e., no centrifuge experiments, sorry.)
        Scalar g = problem.gravity()[dimWorld - 1];

        // this is quite hacky because the dune grid interface does not provide a
        // cellCenterDepth() method (so we ask the problem to provide it). The "good"
        // solution would be to take the Z coordinate of the element centroids, but since
        // ECL seems to like to be inconsistent on that front, it needs to be done like
        // here...
        Scalar zIn = problem.dofCenterDepth(globalSpaceIdx);

        // the distances from the DOF's depths. (i.e., the additional depth of the
        // exterior DOF)
        Scalar distZ = zIn - zEx;

        for (unsigned phaseIdx=0; phaseIdx < numPhases; phaseIdx++) {
            if (!FluidSystem::phaseIsActive(phaseIdx))
                continue;

            // do the gravity correction: compute the hydrostatic pressure for the
            // integration position
            const Evaluation& rhoIn = intQuantsIn.fluidState().density(phaseIdx);
            const auto& rhoEx = exFluidState.density(phaseIdx);
            Evaluation rhoAvg = (rhoIn + rhoEx)/2;

            const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(phaseIdx);
            Evaluation pressureExterior = exFluidState.pressure(phaseIdx);
            pressureExterior += rhoAvg*(distZ*g);

            pressureDifference[phaseIdx] = pressureExterior - pressureInterior;

            // decide the upstream index for the phase. for this we make sure that the
            // degree of freedom which is regarded upstream if both pressures are equal
            // is always the same: if the pressure is equal, the DOF with the lower
            // global index is regarded to be the upstream one.
            const unsigned interiorDofIdx = 0; // Valid only for cell-centered FV.
            if (pressureDifference[phaseIdx] > 0.0) {
                upIdx[phaseIdx] = -1;
                dnIdx[phaseIdx] = interiorDofIdx;
            }
            else {
                upIdx[phaseIdx] = interiorDofIdx;
                dnIdx[phaseIdx] = -1;
            }

            Evaluation transModified = trans;

            if (upIdx[phaseIdx] == interiorDofIdx) {

                // this is slightly hacky because in the automatic differentiation case, it
                // only works for the element centered finite volume method.
                const auto& up = intQuantsIn;

                // deal with water induced rock compaction
                const Scalar transMult = Toolbox::value(up.rockCompTransMultiplier());
                transModified *= transMult;

                volumeFlux[phaseIdx] =
                    pressureDifference[phaseIdx]*up.mobility(phaseIdx)*(-transModified/faceArea);
            }
            else {
                // compute the phase mobility using the material law parameters of the
                // interior element. TODO: this could probably be done more efficiently
                const auto& matParams = problem.materialLawParams(globalSpaceIdx);
                std::array<typename FluidState::ValueType,numPhases> kr;
                MaterialLaw::relativePermeabilities(kr, matParams, exFluidState);

                const auto& mob = kr[phaseIdx]/exFluidState.viscosity(phaseIdx);
                volumeFlux[phaseIdx] =
                    pressureDifference[phaseIdx]*mob*(-transModified/faceArea);
            }
        }
    }

protected:

    OPM_HOST_DEVICE void calculateFluxes_(const ElementContext&, unsigned, unsigned)
    { }

    OPM_HOST_DEVICE void calculateBoundaryFluxes_(const ElementContext&, unsigned, unsigned)
    {}

private:
    Implementation& asImp_()
    { return *static_cast<Implementation*>(this); }

    const Implementation& asImp_() const
    { return *static_cast<const Implementation*>(this); }

    // the volumetric flux of all phases [m^3/s]
    Evaluation volumeFlux_[numPhases];

    // the difference in effective pressure between the exterior and the interior degree
    // of freedom [Pa]
    Evaluation pressureDifference_[numPhases];

    // the local indices of the interior and exterior degrees of freedom
    std::array<short, numPhases> upIdx_;
    std::array<short, numPhases> dnIdx_;
 };

} // namespace Opm

#endif // OPM_NEWTRAN_FLUX_MODULE_HPP