28#ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH
29#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/common/OpmLog/OpmLog.hpp>
35#include <opm/input/eclipse/EclipseState/Phase.hpp>
44#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
45#include <opm/input/eclipse/EclipseState/Tables/FoamadsTable.hpp>
46#include <opm/input/eclipse/EclipseState/Tables/FoammobTable.hpp>
58template <class TypeTag, bool enableFoamV = getPropValue<TypeTag, Properties::EnableFoam>()>
74 using Toolbox = MathToolbox<Evaluation>;
78 static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
79 static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
80 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
81 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
83 static constexpr unsigned enableFoam = enableFoamV;
85 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
86 static constexpr unsigned numPhases = FluidSystem::numPhases;
88 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
110 if constexpr (enableFoam) {
111 if (Parameters::Get<Parameters::EnableVtkOutput>()) {
112 OpmLog::warning(
"VTK output requested, currently unsupported by the foam module.");
120 if constexpr (enableFoam)
121 return pvIdx == foamConcentrationIdx;
129 return "foam_concentration";
137 return static_cast<Scalar
>(1.0);
142 if constexpr (enableFoam)
143 return eqIdx == contiFoamEqIdx;
149 static std::string
eqName([[maybe_unused]]
unsigned eqIdx)
156 static Scalar
eqWeight([[maybe_unused]]
unsigned eqIdx)
161 return static_cast<Scalar
>(1.0);
165 template <
class LhsEval>
166 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
167 const IntensiveQuantities& intQuants)
169 if constexpr (enableFoam) {
170 const auto& fs = intQuants.fluidState();
172 LhsEval surfaceVolume = Toolbox::template decay<LhsEval>(intQuants.porosity());
173 if (params_.transport_phase_ == Phase::WATER) {
174 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
175 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx)));
176 }
else if (params_.transport_phase_ == Phase::GAS) {
177 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
178 * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx)));
179 }
else if (params_.transport_phase_ == Phase::SOLVENT) {
180 if constexpr (enableSolvent) {
181 surfaceVolume *= (Toolbox::template decay<LhsEval>( intQuants.solventSaturation())
182 * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor()));
185 throw std::runtime_error(
"Transport phase is GAS/WATER/SOLVENT");
189 surfaceVolume = max(surfaceVolume, 1e-10);
192 const LhsEval freeFoam = surfaceVolume
193 * Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
196 const LhsEval adsorbedFoam =
197 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
198 * Toolbox::template decay<LhsEval>(intQuants.foamRockDensity())
199 * Toolbox::template decay<LhsEval>(intQuants.foamAdsorbed());
201 LhsEval accumulationFoam = freeFoam + adsorbedFoam;
202 storage[contiFoamEqIdx] += accumulationFoam;
207 [[maybe_unused]]
const ElementContext& elemCtx,
208 [[maybe_unused]]
unsigned scvfIdx,
209 [[maybe_unused]]
unsigned timeIdx)
212 if constexpr (enableFoam) {
213 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
214 const unsigned inIdx = extQuants.interiorIndex();
221 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
222 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
223 if (upIdx == inIdx) {
224 flux[contiFoamEqIdx] =
225 extQuants.volumeFlux(waterPhaseIdx)
226 *up.fluidState().invB(waterPhaseIdx)
227 *up.foamConcentration();
229 flux[contiFoamEqIdx] =
230 extQuants.volumeFlux(waterPhaseIdx)
231 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
232 *decay<Scalar>(up.foamConcentration());
237 const unsigned upIdx = extQuants.upstreamIndex(gasPhaseIdx);
238 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
239 if (upIdx == inIdx) {
240 flux[contiFoamEqIdx] =
241 extQuants.volumeFlux(gasPhaseIdx)
242 *up.fluidState().invB(gasPhaseIdx)
243 *up.foamConcentration();
245 flux[contiFoamEqIdx] =
246 extQuants.volumeFlux(gasPhaseIdx)
247 *decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
248 *decay<Scalar>(up.foamConcentration());
252 case Phase::SOLVENT: {
253 if constexpr (enableSolvent) {
254 const unsigned upIdx = extQuants.solventUpstreamIndex();
255 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
256 if (upIdx == inIdx) {
257 flux[contiFoamEqIdx] =
258 extQuants.solventVolumeFlux()
259 *up.solventInverseFormationVolumeFactor()
260 *up.foamConcentration();
262 flux[contiFoamEqIdx] =
263 extQuants.solventVolumeFlux()
264 *decay<Scalar>(up.solventInverseFormationVolumeFactor())
265 *decay<Scalar>(up.foamConcentration());
268 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
273 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
287 return static_cast<Scalar
>(0.0);
290 template <
class DofEntity>
292 [[maybe_unused]] std::ostream& outstream,
293 [[maybe_unused]]
const DofEntity& dof)
295 if constexpr (enableFoam) {
296 unsigned dofIdx = model.dofMapper().index(dof);
297 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
298 outstream << priVars[foamConcentrationIdx];
302 template <
class DofEntity>
304 [[maybe_unused]] std::istream& instream,
305 [[maybe_unused]]
const DofEntity& dof)
307 if constexpr (enableFoam) {
308 unsigned dofIdx = model.dofMapper().index(dof);
309 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
310 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
312 instream >> priVars0[foamConcentrationIdx];
315 priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
323 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
324 return params_.foamRockDensity_[satnumRegionIdx];
331 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
332 return params_.foamAllowDesorption_[satnumRegionIdx];
339 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
340 return params_.adsorbedFoamTable_[satnumRegionIdx];
347 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
348 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
353 const unsigned scvIdx,
354 const unsigned timeIdx)
356 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
357 return params_.foamCoefficients_[satnumRegionIdx];
361 return params_.transport_phase_;
368template <
class TypeTag,
bool enableFoam>
369BlackOilFoamParams<typename BlackOilFoamModule<TypeTag, enableFoam>::Scalar>
370BlackOilFoamModule<TypeTag, enableFoam>::params_;
379template <class TypeTag, bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>()>
394 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
395 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
397 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
398 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
399 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
400 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
413 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
415 const auto& fs =
asImp_().fluidState_;
418 Evaluation mobilityReductionFactor = 1.0;
425 const Scalar fm_mob = foamCoefficients.fm_mob;
427 const Scalar fm_surf = foamCoefficients.fm_surf;
428 const Scalar ep_surf = foamCoefficients.ep_surf;
430 const Scalar fm_oil = foamCoefficients.fm_oil;
431 const Scalar fl_oil = foamCoefficients.fl_oil;
432 const Scalar ep_oil = foamCoefficients.ep_oil;
434 const Scalar fm_dry = foamCoefficients.fm_dry;
435 const Scalar ep_dry = foamCoefficients.ep_dry;
437 const Scalar fm_cap = foamCoefficients.fm_cap;
438 const Scalar ep_cap = foamCoefficients.ep_cap;
441 const Evaluation Ca = 1e10;
442 const Evaluation S_o = fs.saturation(oilPhaseIdx);
443 const Evaluation S_w = fs.saturation(waterPhaseIdx);
445 Evaluation F1 = pow(C_surf/fm_surf, ep_surf);
446 Evaluation F2 = pow((fm_oil-S_o)/(fm_oil-fl_oil), ep_oil);
447 Evaluation F3 = pow(fm_cap/Ca, ep_cap);
448 Evaluation F7 = 0.5 + atan(ep_dry*(S_w-fm_dry))/M_PI;
450 mobilityReductionFactor = 1./(1. + fm_mob*F1*F2*F3*F7);
462 asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
466 asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
469 case Phase::SOLVENT: {
470 if constexpr (enableSolvent) {
471 asImp_().solventMobility_ *= mobilityReductionFactor;
473 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
478 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
487 throw std::runtime_error(
"Foam module does not support the 'no desorption' option.");
502 {
return *
static_cast<Implementation*
>(
this); }
509template <
class TypeTag>
524 {
throw std::runtime_error(
"foamConcentration() called but foam is disabled"); }
527 {
throw std::runtime_error(
"foamRockDensity() called but foam is disabled"); }
530 {
throw std::runtime_error(
"foamAdsorbed() called but foam is disabled"); }
Contains the parameters to extend the black-oil model to include the effects of foam.
Declares the properties required by the black oil model.
Scalar foamRockDensity() const
Definition: blackoilfoammodules.hh:526
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:523
void foamPropertiesUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilfoammodules.hh:517
Scalar foamAdsorbed() const
Definition: blackoilfoammodules.hh:529
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilfoammodules.hh:381
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:491
Implementation & asImp_()
Definition: blackoilfoammodules.hh:501
Scalar foamRockDensity_
Definition: blackoilfoammodules.hh:505
const Evaluation & foamAdsorbed() const
Definition: blackoilfoammodules.hh:497
Evaluation foamAdsorbed_
Definition: blackoilfoammodules.hh:506
Evaluation foamConcentration_
Definition: blackoilfoammodules.hh:504
Scalar foamRockDensity() const
Definition: blackoilfoammodules.hh:494
void foamPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilfoammodules.hh:409
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition: blackoilfoammodules.hh:60
static bool eqApplies(unsigned eqIdx)
Definition: blackoilfoammodules.hh:140
static void registerOutputModules(Model &, Simulator &)
Register all foam specific VTK and ECL output modules.
Definition: blackoilfoammodules.hh:107
static void registerParameters()
Register all run-time parameters for the black-oil foam module.
Definition: blackoilfoammodules.hh:100
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilfoammodules.hh:126
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilfoammodules.hh:156
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:303
static void setParams(BlackOilFoamParams< Scalar > &¶ms)
Set parameters.
Definition: blackoilfoammodules.hh:92
static std::string eqName(unsigned eqIdx)
Definition: blackoilfoammodules.hh:149
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilfoammodules.hh:118
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilfoammodules.hh:166
static bool foamAllowDesorption(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:327
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilfoammodules.hh:132
static const BlackOilFoamParams< Scalar >::FoamCoefficients & foamCoefficients(const ElementContext &elemCtx, const unsigned scvIdx, const unsigned timeIdx)
Definition: blackoilfoammodules.hh:352
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:206
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilfoammodules.hh:282
static const Scalar foamRockDensity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:319
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:291
static const TabulatedFunction & adsorbedFoamTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:335
static Phase transportPhase()
Definition: blackoilfoammodules.hh:360
static const TabulatedFunction & gasMobilityMultiplierTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:343
Declare the properties used by the infrastructure code of the finite volume discretizations.
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
Definition: blackoilfoamparams.hpp:64
Struct holding the parameters for the BlackoilFoamModule class.
Definition: blackoilfoamparams.hpp:46
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilfoamparams.hpp:47