blackoillocalresidual.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_BLACK_OIL_LOCAL_RESIDUAL_HH
#define EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH

#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/fluidstates/BlackOilFluidState.hpp>

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

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

#include <cassert>
#include <limits>

namespace Opm {

template <class TypeTag>
class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalResidual>
{
    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
    using Indices = GetPropType<TypeTag, Properties::Indices>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;

    enum { conti0EqIdx = Indices::conti0EqIdx };
    enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
    enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };

    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };

    enum { gasCompIdx = FluidSystem::gasCompIdx };
    enum { oilCompIdx = FluidSystem::oilCompIdx };
    enum { waterCompIdx = FluidSystem::waterCompIdx };
    static constexpr unsigned compositionSwitchIdx = Indices::compositionSwitchIdx;

    static constexpr bool waterEnabled = Indices::waterEnabled;
    static constexpr bool gasEnabled = Indices::gasEnabled;
    static constexpr bool oilEnabled = Indices::oilEnabled;
    static constexpr bool compositionSwitchEnabled =
        compositionSwitchIdx != std::numeric_limits<unsigned>::max();

    static constexpr bool blackoilConserveSurfaceVolume = 
        getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
    static constexpr bool enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>();
    static constexpr bool enableBrine = getPropValue<TypeTag, Properties::EnableBrine>();
    static constexpr bool enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>();
    static constexpr bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
    static constexpr bool enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>();
    static constexpr bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>();
    static constexpr EnergyModules energyModuleType = getPropValue<TypeTag, Properties::EnergyModuleType>();
    static constexpr bool enableFullyImplicitThermal = energyModuleType == EnergyModules::FullyImplicitThermal;
    static constexpr bool enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>();
    static constexpr bool enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>();

    using Toolbox = MathToolbox<Evaluation>;

    using BioeffectsModule = BlackOilBioeffectsModule<TypeTag, enableBioeffects>;
    using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
    using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
    using EnergyModule = BlackOilEnergyModule<TypeTag, energyModuleType>;
    using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
    using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
    using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
    using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;

public:
    template <class LhsEval>
    void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                        const ElementContext& elemCtx,
                        unsigned dofIdx,
                        unsigned timeIdx) const
    {
        // retrieve the intensive quantities for the SCV at the specified point in time
        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
        const auto& fs = intQuants.fluidState();

        storage = 0.0;

        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                if (Indices::numPhases == 3) { // add trivial equation for the pseudo phase
                    const unsigned activeCompIdx =
                        FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
                    if (timeIdx == 0) {
                        storage[conti0EqIdx + activeCompIdx] = variable<LhsEval>(0.0, conti0EqIdx + activeCompIdx);
                    }
                    else {
                        storage[conti0EqIdx + activeCompIdx] = 0.0;
                    }
                }
                continue;
            }

            const unsigned activeCompIdx =
                FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));
            const LhsEval surfaceVolume =
                Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx)) *
                Toolbox::template decay<LhsEval>(fs.invB(phaseIdx)) *
                Toolbox::template decay<LhsEval>(intQuants.porosity());

            storage[conti0EqIdx + activeCompIdx] += surfaceVolume;

            // account for dissolved gas
            if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) {
                const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
                storage[conti0EqIdx + activeGasCompIdx] +=
                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs()) *
                    surfaceVolume;
            }

            // account for dissolved gas in water phase
            if (phaseIdx == waterPhaseIdx && FluidSystem::enableDissolvedGasInWater()) {
                const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
                storage[conti0EqIdx + activeGasCompIdx] +=
                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rsw()) *
                    surfaceVolume;
            }

            // account for vaporized oil
            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) {
                const unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
                storage[conti0EqIdx + activeOilCompIdx] +=
                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rv()) *
                    surfaceVolume;
            }

            // account for vaporized water
            if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) {
                const unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
                storage[conti0EqIdx + activeWaterCompIdx] +=
                    Toolbox::template decay<LhsEval>(intQuants.fluidState().Rvw()) *
                    surfaceVolume;
            }
        }

        adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex());

        // deal with solvents (if present)
        if constexpr (enableSolvent) {
            SolventModule::addStorage(storage, intQuants);
        }

        // deal with zFracton (if present)
        if constexpr (enableExtbo) {
            ExtboModule::addStorage(storage, intQuants);
        }

        // deal with polymer (if present)
        if constexpr (enablePolymer) {
            PolymerModule::addStorage(storage, intQuants);
        }

        // deal with energy (if present)
        if constexpr (enableFullyImplicitThermal) {
            EnergyModule::addStorage(storage, intQuants);
        }

        // deal with foam (if present)
        if constexpr (enableFoam) {
            FoamModule::addStorage(storage, intQuants);
        }

        // deal with salt (if present)
        if constexpr (enableBrine) {
            BrineModule::addStorage(storage, intQuants);
        }

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addStorage(storage, intQuants);
        }
    }

    void computeFlux(RateVector& flux,
                     const ElementContext& elemCtx,
                     unsigned scvfIdx,
                     unsigned timeIdx) const
    {
        assert(timeIdx == 0);

        flux = 0.0;

        const ExtensiveQuantities& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
        unsigned focusDofIdx = elemCtx.focusDofIndex();
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if (!FluidSystem::phaseIsActive(phaseIdx)) {
                continue;
            }

            const unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(phaseIdx));
            const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
            const unsigned pvtRegionIdx = up.pvtRegionIndex();
            if (upIdx == focusDofIdx) {
                evalPhaseFluxes_<Evaluation>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
            }
            else {
                evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
            }
        }

        // deal with solvents (if present)
        if constexpr (enableSolvent) {
            SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with zFracton (if present)
        if constexpr (enableExtbo) {
            ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with polymer (if present)
        if constexpr (enablePolymer) {
            PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with energy (if present)
        if constexpr (enableFullyImplicitThermal) {
            EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with foam (if present)
        if constexpr (enableFoam) {
            FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with salt (if present)
        if constexpr (enableBrine) {
            BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with bioeffects (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with diffusion (if present)
        if constexpr (enableDiffusion) {
            DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
        }

        // deal with convective mixing (if enabled)
        if constexpr (enableConvectiveMixing) {
            ConvectiveMixingModule::addConvectiveMixingFlux(flux, elemCtx, scvfIdx, timeIdx);
        }
    }

    void computeSource(RateVector& source,
                       const ElementContext& elemCtx,
                       unsigned dofIdx,
                       unsigned timeIdx) const
    {
        // retrieve the source term intrinsic to the problem
        elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);

        // deal with MICP (if present)
        if constexpr (enableBioeffects) {
            BioeffectsModule::addSource(source, elemCtx, dofIdx, timeIdx);
        }

        // scale the source term of the energy equation
        if constexpr (enableFullyImplicitThermal) {
            source[Indices::contiEnergyEqIdx] *=
                getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
        }
    }

    template <class UpEval, class FluidState>
    static void evalPhaseFluxes_(RateVector& flux,
                                 unsigned phaseIdx,
                                 unsigned pvtRegionIdx,
                                 const ExtensiveQuantities& extQuants,
                                 const FluidState& upFs)
    {
        const auto& invB = getInvB_<FluidSystem, FluidState, UpEval>(upFs, phaseIdx, pvtRegionIdx);
        const auto& surfaceVolumeFlux = invB*extQuants.volumeFlux(phaseIdx);
        const unsigned activeCompIdx =
            FluidSystem::canonicalToActiveCompIdx(FluidSystem::solventComponentIndex(phaseIdx));

        if constexpr (blackoilConserveSurfaceVolume) {
            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux;
        }
        else {
            flux[conti0EqIdx + activeCompIdx] += surfaceVolumeFlux *
                                                 FluidSystem::referenceDensity(phaseIdx, pvtRegionIdx);
        }

        if (phaseIdx == oilPhaseIdx) {
            // dissolved gas (in the oil phase).
            if (FluidSystem::enableDissolvedGas()) {
                const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeGasCompIdx] += Rs * surfaceVolumeFlux;
                }
                else {
                    flux[conti0EqIdx + activeGasCompIdx] +=
                        Rs * surfaceVolumeFlux *
                        FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
                }
            }
        } else if (phaseIdx == waterPhaseIdx) {
            // dissolved gas (in the water phase).
            if (FluidSystem::enableDissolvedGasInWater()) {
                const auto& Rsw = BlackOil::getRsw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeGasCompIdx] += Rsw * surfaceVolumeFlux;
                }
                else {
                    flux[conti0EqIdx + activeGasCompIdx] +=
                        Rsw * surfaceVolumeFlux *
                        FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
                }
            }
        }
        else if (phaseIdx == gasPhaseIdx) {
            // vaporized oil (in the gas phase).
            if (FluidSystem::enableVaporizedOil()) {
                const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeOilCompIdx] += Rv * surfaceVolumeFlux;
                }
                else {
                    flux[conti0EqIdx + activeOilCompIdx] +=
                        Rv * surfaceVolumeFlux *
                        FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
                }
            }
             // vaporized water (in the gas phase).
            if (FluidSystem::enableVaporizedWater()) {
                const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, UpEval>(upFs, pvtRegionIdx);

                const unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
                if constexpr (blackoilConserveSurfaceVolume) {
                    flux[conti0EqIdx + activeWaterCompIdx] += Rvw * surfaceVolumeFlux;
                }
                else {
                    flux[conti0EqIdx + activeWaterCompIdx] +=
                        Rvw * surfaceVolumeFlux *
                        FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
                }
            }
        }
    }

    template <class Scalar>
    static void adaptMassConservationQuantities_(Dune::FieldVector<Scalar, numEq>& container,
                                                 unsigned pvtRegionIdx)
    {
        if constexpr (!blackoilConserveSurfaceVolume) {
            // convert "surface volume" to mass. this is complicated a bit by the fact that
            // not all phases are necessarily enabled. (we here assume that if a fluid phase
            // is disabled, its respective "main" component is not considered as well.)

            if constexpr (waterEnabled) {
                const unsigned activeWaterCompIdx = FluidSystem::canonicalToActiveCompIdx(waterCompIdx);
                container[conti0EqIdx + activeWaterCompIdx] *=
                    FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
            }

            if constexpr (gasEnabled) {
                const unsigned activeGasCompIdx = FluidSystem::canonicalToActiveCompIdx(gasCompIdx);
                container[conti0EqIdx + activeGasCompIdx] *=
                    FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
            }

            if constexpr (oilEnabled) {
                const unsigned activeOilCompIdx = FluidSystem::canonicalToActiveCompIdx(oilCompIdx);
                container[conti0EqIdx + activeOilCompIdx] *=
                    FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
            }
        }
    }
};

} // namespace Opm

#endif