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>
56template <class TypeTag, bool enableFoamV = getPropValue<TypeTag, Properties::EnableFoam>()>
72 using Toolbox = MathToolbox<Evaluation>;
76 static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
77 static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
78 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
79 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
81 static constexpr unsigned enableFoam = enableFoamV;
83 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
86 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
107 if constexpr (enableFoam) {
108 if (Parameters::Get<Parameters::EnableVtkOutput>()) {
109 OpmLog::warning(
"VTK output requested, currently unsupported by the foam module.");
117 if constexpr (enableFoam) {
118 return pvIdx == foamConcentrationIdx;
128 return "foam_concentration";
136 return static_cast<Scalar
>(1.0);
141 if constexpr (enableFoam) {
142 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 const LhsEval accumulationFoam = freeFoam + adsorbedFoam;
202 storage[contiFoamEqIdx] += accumulationFoam;
207 [[maybe_unused]]
const ElementContext& elemCtx,
208 [[maybe_unused]]
unsigned scvfIdx,
209 [[maybe_unused]]
unsigned timeIdx)
211 if constexpr (enableFoam) {
212 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
213 const unsigned inIdx = extQuants.interiorIndex();
220 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
221 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
222 if (upIdx == inIdx) {
223 flux[contiFoamEqIdx] =
224 extQuants.volumeFlux(waterPhaseIdx) *
225 up.fluidState().invB(waterPhaseIdx) *
226 up.foamConcentration();
228 flux[contiFoamEqIdx] =
229 extQuants.volumeFlux(waterPhaseIdx) *
230 decay<Scalar>(up.fluidState().invB(waterPhaseIdx)) *
231 decay<Scalar>(up.foamConcentration());
236 const unsigned upIdx = extQuants.upstreamIndex(gasPhaseIdx);
237 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
238 if (upIdx == inIdx) {
239 flux[contiFoamEqIdx] =
240 extQuants.volumeFlux(gasPhaseIdx) *
241 up.fluidState().invB(gasPhaseIdx) *
242 up.foamConcentration();
244 flux[contiFoamEqIdx] =
245 extQuants.volumeFlux(gasPhaseIdx) *
246 decay<Scalar>(up.fluidState().invB(gasPhaseIdx)) *
247 decay<Scalar>(up.foamConcentration());
252 if constexpr (enableSolvent) {
253 const unsigned upIdx = extQuants.solventUpstreamIndex();
254 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
255 if (upIdx == inIdx) {
256 flux[contiFoamEqIdx] =
257 extQuants.solventVolumeFlux() *
258 up.solventInverseFormationVolumeFactor() *
259 up.foamConcentration();
261 flux[contiFoamEqIdx] =
262 extQuants.solventVolumeFlux() *
263 decay<Scalar>(up.solventInverseFormationVolumeFactor()) *
264 decay<Scalar>(up.foamConcentration());
267 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
271 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
284 return static_cast<Scalar
>(0.0);
287 template <
class DofEntity>
289 [[maybe_unused]] std::ostream& outstream,
290 [[maybe_unused]]
const DofEntity& dof)
292 if constexpr (enableFoam) {
293 const unsigned dofIdx = model.dofMapper().index(dof);
294 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
295 outstream << priVars[foamConcentrationIdx];
299 template <
class DofEntity>
301 [[maybe_unused]] std::istream& instream,
302 [[maybe_unused]]
const DofEntity& dof)
304 if constexpr (enableFoam) {
305 const unsigned dofIdx = model.dofMapper().index(dof);
306 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
307 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
309 instream >> priVars0[foamConcentrationIdx];
312 priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
320 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
321 return params_.foamRockDensity_[satnumRegionIdx];
328 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
329 return params_.foamAllowDesorption_[satnumRegionIdx];
336 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
337 return params_.adsorbedFoamTable_[satnumRegionIdx];
344 const unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
345 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
350 const unsigned scvIdx,
351 const unsigned timeIdx)
353 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
354 return params_.foamCoefficients_[satnumRegionIdx];
358 {
return params_.transport_phase_; }
364template <
class TypeTag,
bool enableFoam>
365BlackOilFoamParams<typename BlackOilFoamModule<TypeTag, enableFoam>::Scalar>
366BlackOilFoamModule<TypeTag, enableFoam>::params_;
368template <
class TypeTag,
bool enableFoam>
378template <
class TypeTag>
393 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
394 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
396 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
397 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
398 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
399 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
411 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
412 foamConcentration_ = priVars.makeEvaluation(foamConcentrationIdx, timeIdx);
413 const auto& fs = asImp_().fluidState_;
416 Evaluation mobilityReductionFactor = 1.0;
417 if constexpr (
false) {
421 const auto& foamCoefficients = FoamModule::foamCoefficients(elemCtx, dofIdx, timeIdx);
423 const Scalar fm_mob = foamCoefficients.fm_mob;
425 const Scalar fm_surf = foamCoefficients.fm_surf;
426 const Scalar ep_surf = foamCoefficients.ep_surf;
428 const Scalar fm_oil = foamCoefficients.fm_oil;
429 const Scalar fl_oil = foamCoefficients.fl_oil;
430 const Scalar ep_oil = foamCoefficients.ep_oil;
432 const Scalar fm_dry = foamCoefficients.fm_dry;
433 const Scalar ep_dry = foamCoefficients.ep_dry;
435 const Scalar fm_cap = foamCoefficients.fm_cap;
436 const Scalar ep_cap = foamCoefficients.ep_cap;
438 const Evaluation C_surf = foamConcentration_;
439 const Evaluation Ca = 1e10;
440 const Evaluation S_o = fs.saturation(oilPhaseIdx);
441 const Evaluation S_w = fs.saturation(waterPhaseIdx);
443 const Evaluation F1 = pow(C_surf / fm_surf, ep_surf);
444 const Evaluation F2 = pow((fm_oil - S_o) / (fm_oil - fl_oil), ep_oil);
445 const Evaluation F3 = pow(fm_cap / Ca, ep_cap);
446 const Evaluation F7 = 0.5 + atan(ep_dry * (S_w - fm_dry)) / M_PI;
448 mobilityReductionFactor = 1. / (1. + fm_mob * F1 * F2 * F3 * F7);
453 const auto& gasMobilityMultiplier = FoamModule::gasMobilityMultiplierTable(elemCtx, dofIdx, timeIdx);
454 mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_,
true);
458 switch (FoamModule::transportPhase()) {
460 asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
463 asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
466 if constexpr (enableSolvent) {
467 asImp_().solventMobility_ *= mobilityReductionFactor;
469 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
473 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
476 foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
478 const auto& adsorbedFoamTable = FoamModule::adsorbedFoamTable(elemCtx, dofIdx, timeIdx);
479 foamAdsorbed_ = adsorbedFoamTable.eval(foamConcentration_,
true);
480 if (!FoamModule::foamAllowDesorption(elemCtx, dofIdx, timeIdx)) {
481 throw std::runtime_error(
"Foam module does not support the 'no desorption' option.");
486 {
return foamConcentration_; }
489 {
return foamRockDensity_; }
492 {
return foamAdsorbed_; }
496 {
return *
static_cast<Implementation*
>(
this); }
503template <
class TypeTag>
517 {
throw std::runtime_error(
"foamConcentration() called but foam is disabled"); }
520 {
throw std::runtime_error(
"foamRockDensity() called but foam is disabled"); }
523 {
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:519
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:516
void foamPropertiesUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilfoammodules.hh:511
Scalar foamAdsorbed() const
Definition: blackoilfoammodules.hh:522
Evaluation foamConcentration_
Definition: blackoilfoammodules.hh:498
void foamPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilfoammodules.hh:407
const Evaluation & foamAdsorbed() const
Definition: blackoilfoammodules.hh:491
Scalar foamRockDensity() const
Definition: blackoilfoammodules.hh:488
Evaluation foamAdsorbed_
Definition: blackoilfoammodules.hh:500
Implementation & asImp_()
Definition: blackoilfoammodules.hh:495
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:485
Scalar foamRockDensity_
Definition: blackoilfoammodules.hh:499
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilfoammodules.hh:369
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition: blackoilfoammodules.hh:58
static bool eqApplies(unsigned eqIdx)
Definition: blackoilfoammodules.hh:139
static void registerOutputModules(Model &, Simulator &)
Register all foam specific VTK and ECL output modules.
Definition: blackoilfoammodules.hh:104
static void registerParameters()
Register all run-time parameters for the black-oil foam module.
Definition: blackoilfoammodules.hh:98
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilfoammodules.hh:125
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilfoammodules.hh:156
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:300
static void setParams(BlackOilFoamParams< Scalar > &¶ms)
Set parameters.
Definition: blackoilfoammodules.hh:90
static std::string eqName(unsigned eqIdx)
Definition: blackoilfoammodules.hh:149
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilfoammodules.hh:115
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:324
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilfoammodules.hh:131
static const BlackOilFoamParams< Scalar >::FoamCoefficients & foamCoefficients(const ElementContext &elemCtx, const unsigned scvIdx, const unsigned timeIdx)
Definition: blackoilfoammodules.hh:349
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:279
static const Scalar foamRockDensity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:316
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:288
static const TabulatedFunction & adsorbedFoamTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:332
static Phase transportPhase()
Definition: blackoilfoammodules.hh:357
static const TabulatedFunction & gasMobilityMultiplierTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:340
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:39
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:233
Definition: blackoilfoamparams.hpp:64
Struct holding the parameters for the BlackoilFoamModule class.
Definition: blackoilfoamparams.hpp:46
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilfoamparams.hpp:47