From 3856a5a84eba36800af19020da3a53056b47e6a9 Mon Sep 17 00:00:00 2001 From: hnil Date: Tue, 7 Jun 2022 21:07:01 +0200 Subject: [PATCH] added tpfa variant --- .../blackoil/blackoillocalresidualtpfa.hh | 379 ++++++++++++++++++ 1 file changed, 379 insertions(+) create mode 100644 opm/models/blackoil/blackoillocalresidualtpfa.hh diff --git a/opm/models/blackoil/blackoillocalresidualtpfa.hh b/opm/models/blackoil/blackoillocalresidualtpfa.hh new file mode 100644 index 000000000..eda11d0bf --- /dev/null +++ b/opm/models/blackoil/blackoillocalresidualtpfa.hh @@ -0,0 +1,379 @@ +// -*- 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 . + + 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. +*/ +/*! + * \file + * + * \copydoc Opm::BlackOilLocalResidual + */ +#ifndef EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH +#define EWOMS_BLACK_OIL_LOCAL_RESIDUAL_HH + +#include "blackoilproperties.hh" +#include "blackoilsolventmodules.hh" +#include "blackoilextbomodules.hh" +#include "blackoilpolymermodules.hh" +#include "blackoilenergymodules.hh" +#include "blackoilfoammodules.hh" +#include "blackoilbrinemodules.hh" +#include "blackoildiffusionmodule.hh" +#include "blackoilmicpmodules.hh" +#include + +namespace Opm { +/*! + * \ingroup BlackOilModel + * + * \brief Calculates the local residual of the black oil model. + */ +template +class BlackOilLocalResidual : public GetPropType +{ + using IntensiveQuantities = GetPropType; + using ExtensiveQuantities = GetPropType; + using ElementContext = GetPropType; + using Indices = GetPropType; + using Scalar = GetPropType; + using Evaluation = GetPropType; + using EqVector = GetPropType; + using RateVector = GetPropType; + using FluidSystem = GetPropType; + + + enum { conti0EqIdx = Indices::conti0EqIdx }; + enum { numEq = getPropValue() }; + enum { numPhases = getPropValue() }; + enum { numComponents = getPropValue() }; + + enum { gasPhaseIdx = FluidSystem::gasPhaseIdx }; + enum { oilPhaseIdx = FluidSystem::oilPhaseIdx }; + enum { waterPhaseIdx = FluidSystem::waterPhaseIdx }; + + enum { gasCompIdx = FluidSystem::gasCompIdx }; + enum { oilCompIdx = FluidSystem::oilCompIdx }; + enum { waterCompIdx = FluidSystem::waterCompIdx }; + enum { compositionSwitchIdx = Indices::compositionSwitchIdx }; + + static const bool waterEnabled = Indices::waterEnabled; + static const bool gasEnabled = Indices::gasEnabled; + static const bool oilEnabled = Indices::oilEnabled; + static const bool compositionSwitchEnabled = (compositionSwitchIdx >= 0); + + static constexpr bool blackoilConserveSurfaceVolume = getPropValue(); + static constexpr bool enableEnergy = getPropValue(); + static constexpr bool enableDiffusion = getPropValue(); + + using Toolbox = MathToolbox; + using SolventModule = BlackOilSolventModule; + using ExtboModule = BlackOilExtboModule; + using PolymerModule = BlackOilPolymerModule; + using EnergyModule = BlackOilEnergyModule; + using FoamModule = BlackOilFoamModule; + using BrineModule = BlackOilBrineModule; + using DiffusionModule = BlackOilDiffusionModule; + using MICPModule = BlackOilMICPModule; + +public: + /*! + * \copydoc FvBaseLocalResidual::computeStorage + */ + template + void computeStorage(Dune::FieldVector& storage, + const ElementContext& elemCtx, + unsigned dofIdx, + unsigned timeIdx) const + { + const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx); + computeStorage(storage, + intQuants, + timeIdx); + } + template + void computeStorage(Dune::FieldVector& storage, + const IntensiveQuantities& intQuants, + unsigned timeIdx) const{ + // retrieve the intensive quantities for the SCV at the specified point in time + 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 + unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + if (timeIdx == 0) + storage[conti0EqIdx + activeCompIdx] = variable(0.0, conti0EqIdx + activeCompIdx); + else + storage[conti0EqIdx + activeCompIdx] = 0.0; + } + continue; + } + + unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + LhsEval surfaceVolume = + Toolbox::template decay(fs.saturation(phaseIdx)) + * Toolbox::template decay(fs.invB(phaseIdx)) + * Toolbox::template decay(intQuants.porosity()); + + storage[conti0EqIdx + activeCompIdx] += surfaceVolume; + + // account for dissolved gas + if (phaseIdx == oilPhaseIdx && FluidSystem::enableDissolvedGas()) { + unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx); + storage[conti0EqIdx + activeGasCompIdx] += + Toolbox::template decay(intQuants.fluidState().Rs()) + * surfaceVolume; + } + + // account for vaporized oil + if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedOil()) { + unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx); + storage[conti0EqIdx + activeOilCompIdx] += + Toolbox::template decay(intQuants.fluidState().Rv()) + * surfaceVolume; + } + + // account for vaporized water + if (phaseIdx == gasPhaseIdx && FluidSystem::enableVaporizedWater()) { + unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx); + storage[conti0EqIdx + activeWaterCompIdx] += + Toolbox::template decay(intQuants.fluidState().Rvw()) + * surfaceVolume; + } + } + + adaptMassConservationQuantities_(storage, intQuants.pvtRegionIndex()); + + // deal with solvents (if present) + SolventModule::addStorage(storage, intQuants); + + // deal with zFracton (if present) + ExtboModule::addStorage(storage, intQuants); + + // deal with polymer (if present) + PolymerModule::addStorage(storage, intQuants); + + // deal with energy (if present) + EnergyModule::addStorage(storage, intQuants); + + // deal with foam (if present) + FoamModule::addStorage(storage, intQuants); + + // deal with salt (if present) + BrineModule::addStorage(storage, intQuants); + + // deal with micp (if present) + MICPModule::addStorage(storage, intQuants); + } + + /*! + * \copydoc FvBaseLocalResidual::computeFlux + */ + 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 auto& darcyFlux = extQuants.volumeFlux(phaseIdx); + unsigned upIdx = static_cast(extQuants.upstreamIndex(phaseIdx)); + const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx); + unsigned pvtRegionIdx = up.pvtRegionIndex(); + using FluidState = typename IntensiveQuantities::FluidState; + if (upIdx == focusDofIdx){ + const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); + const auto& surfaceVolumeFlux = invB*darcyFlux; + evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); + }else{ + const auto& invB = getInvB_(up.fluidState(), phaseIdx, pvtRegionIdx); + const auto& surfaceVolumeFlux = invB*darcyFlux; + evalPhaseFluxes_(flux, phaseIdx, pvtRegionIdx, surfaceVolumeFlux, up.fluidState()); + } + } + + // deal with solvents (if present) + SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with zFracton (if present) + ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with polymer (if present) + PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with energy (if present) + EnergyModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with foam (if present) + FoamModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with salt (if present) + BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + // deal with micp (if present) + MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx); + + DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx); + } + + /*! + * \copydoc FvBaseLocalResidual::computeSource + */ + 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) + MICPModule::addSource(source, elemCtx, dofIdx, timeIdx); + + // scale the source term of the energy equation + if (enableEnergy) + source[Indices::contiEnergyEqIdx] *= getPropValue(); + } + + template + static void evalPhaseFluxes_(RateVector& flux, + unsigned phaseIdx, + unsigned pvtRegionIdx, + const ExtensiveQuantities& extQuants, + const FluidState& upFs){ + + const auto& invB = getInvB_(upFs, phaseIdx, pvtRegionIdx); + const auto& surfaceVolumeFlux = invB*extQuants.volumeFlux(phaseIdx); + evalPhaseFluxes_(flux, + phaseIdx, + pvtRegionIdx, + surfaceVolumeFlux, + upFs); + + } + /*! + * \brief Helper function to calculate the flux of mass in terms of conservation + * quantities via specific fluid phase over a face. + */ + template + static void evalPhaseFluxes_(RateVector& flux, + unsigned phaseIdx, + unsigned pvtRegionIdx, + const Eval& surfaceVolumeFlux,//value of RateVector + const FluidState& upFs) + { + //const auto& invB = getInvB_(upFs, phaseIdx, pvtRegionIdx); + //const auto& surfaceVolumeFlux = invB*extQuants.volumeFlux(phaseIdx); + unsigned activeCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::solventComponentIndex(phaseIdx)); + + if (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_(upFs, pvtRegionIdx); + + unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx); + if (blackoilConserveSurfaceVolume) + flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux; + else + flux[conti0EqIdx + activeGasCompIdx] += Rs*surfaceVolumeFlux*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); + } + } + else if (phaseIdx == gasPhaseIdx) { + // vaporized oil (in the gas phase). + if (FluidSystem::enableVaporizedOil()) { + const auto& Rv = BlackOil::getRv_(upFs, pvtRegionIdx); + + unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx); + if (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_(upFs, pvtRegionIdx); + + unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx); + if (blackoilConserveSurfaceVolume) + flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux; + else + flux[conti0EqIdx + activeWaterCompIdx] += Rvw*surfaceVolumeFlux*FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx); + } + } + } + + /*! + * \brief Helper function to convert the mass-related parts of a Dune::FieldVector + * that stores conservation quantities in terms of "surface-volume" to the + * conservation quantities used by the model. + * + * Depending on the value of the BlackoilConserveSurfaceVolume property, the model + * either conserves mass by means of "surface volume" of the components or mass + * directly. In the former case, this method is a no-op; in the latter, the values + * passed are multiplied by their respective pure component's density at surface + * conditions. + */ + template + static void adaptMassConservationQuantities_(Dune::FieldVector& container, unsigned pvtRegionIdx) + { + if (blackoilConserveSurfaceVolume) + return; + + // 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 (waterEnabled) { + unsigned activeWaterCompIdx = Indices::canonicalToActiveComponentIndex(waterCompIdx); + container[conti0EqIdx + activeWaterCompIdx] *= + FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx); + } + + if (gasEnabled) { + unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(gasCompIdx); + container[conti0EqIdx + activeGasCompIdx] *= + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx); + } + + if (oilEnabled) { + unsigned activeOilCompIdx = Indices::canonicalToActiveComponentIndex(oilCompIdx); + container[conti0EqIdx + activeOilCompIdx] *= + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx); + } + } +}; + +} // namespace Opm + +#endif