28#ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH
29#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
33#include <opm/common/OpmLog/OpmLog.hpp>
38#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
39#include <opm/input/eclipse/EclipseState/Tables/FoamadsTable.hpp>
40#include <opm/input/eclipse/EclipseState/Tables/FoammobTable.hpp>
43#include <dune/common/fvector.hh>
55template <class TypeTag, bool enableFoamV = getPropValue<TypeTag, Properties::EnableFoam>()>
71 using Toolbox = MathToolbox<Evaluation>;
75 static constexpr unsigned foamConcentrationIdx = Indices::foamConcentrationIdx;
76 static constexpr unsigned contiFoamEqIdx = Indices::contiFoamEqIdx;
77 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
78 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
80 static constexpr unsigned enableFoam = enableFoamV;
81 static constexpr bool enableVtkOutput = getPropValue<TypeTag, Properties::EnableVtkOutput>();
83 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
86 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
93 static void initFromState(
const EclipseState& eclState)
97 if (enableFoam && !eclState.runspec().phases().active(Phase::FOAM)) {
98 throw std::runtime_error(
"Non-trivial foam treatment requested at compile time, but "
99 "the deck does not contain the FOAM keyword");
101 else if (!enableFoam && eclState.runspec().phases().active(Phase::FOAM)) {
102 throw std::runtime_error(
"Foam treatment disabled at compile time, but the deck "
103 "contains the FOAM keyword");
106 if (!eclState.runspec().phases().active(Phase::FOAM)) {
110 params_.transport_phase_ = eclState.getInitConfig().getFoamConfig().getTransportPhase();
112 if (eclState.getInitConfig().getFoamConfig().getMobilityModel() != FoamConfig::MobilityModel::TAB) {
113 throw std::runtime_error(
"In FOAMOPTS, only TAB is allowed for the gas mobility factor reduction model.");
116 const auto& tableManager = eclState.getTableManager();
117 const unsigned int numSatRegions = tableManager.getTabdims().getNumSatTables();
118 params_.setNumSatRegions(numSatRegions);
119 const unsigned int numPvtRegions = tableManager.getTabdims().getNumPVTTables();
120 params_.gasMobilityMultiplierTable_.resize(numPvtRegions);
123 const FoamConfig& foamConf = eclState.getInitConfig().getFoamConfig();
124 if (numSatRegions != foamConf.size()) {
125 throw std::runtime_error(
"Inconsistent sizes, number of saturation regions differ from the number of elements "
126 "in FoamConfig, which typically corresponds to the number of records in FOAMROCK.");
130 const auto& foamadsTables = tableManager.getFoamadsTables();
131 if (foamadsTables.empty()) {
132 throw std::runtime_error(
"FOAMADS must be specified in FOAM runs");
134 if (numSatRegions != foamadsTables.size()) {
135 throw std::runtime_error(
"Inconsistent sizes, number of saturation regions differ from the "
136 "number of FOAMADS tables.");
140 for (std::size_t satReg = 0; satReg < numSatRegions; ++satReg) {
141 const auto& rec = foamConf.getRecord(satReg);
143 params_.foamCoefficients_[satReg].fm_min = rec.minimumSurfactantConcentration();
144 params_.foamCoefficients_[satReg].fm_surf = rec.referenceSurfactantConcentration();
145 params_.foamCoefficients_[satReg].ep_surf = rec.exponent();
146 params_.foamRockDensity_[satReg] = rec.rockDensity();
147 params_.foamAllowDesorption_[satReg] = rec.allowDesorption();
148 const auto& foamadsTable = foamadsTables.template getTable<FoamadsTable>(satReg);
149 const auto& conc = foamadsTable.getFoamConcentrationColumn();
150 const auto& ads = foamadsTable.getAdsorbedFoamColumn();
151 params_.adsorbedFoamTable_[satReg].setXYContainers(conc, ads);
155 const auto& foammobTables = tableManager.getFoammobTables();
156 if (foammobTables.empty()) {
160 throw std::runtime_error(
"FOAMMOB must be specified in FOAM runs");
162 if (numPvtRegions != foammobTables.size()) {
163 throw std::runtime_error(
"Inconsistent sizes, number of PVT regions differ from the "
164 "number of FOAMMOB tables.");
168 for (std::size_t pvtReg = 0; pvtReg < numPvtRegions; ++pvtReg) {
169 const auto& foammobTable = foammobTables.template getTable<FoammobTable>(pvtReg);
170 const auto& conc = foammobTable.getFoamConcentrationColumn();
171 const auto& mobMult = foammobTable.getMobilityMultiplierColumn();
172 params_.gasMobilityMultiplierTable_[pvtReg].setXYContainers(conc, mobMult);
190 if constexpr (enableFoam) {
191 if (enableVtkOutput) {
192 OpmLog::warning(
"VTK output requested, currently unsupported by the foam module.");
200 if constexpr (enableFoam)
201 return pvIdx == foamConcentrationIdx;
209 return "foam_concentration";
217 return static_cast<Scalar
>(1.0);
222 if constexpr (enableFoam)
223 return eqIdx == contiFoamEqIdx;
229 static std::string
eqName([[maybe_unused]]
unsigned eqIdx)
236 static Scalar
eqWeight([[maybe_unused]]
unsigned eqIdx)
241 return static_cast<Scalar
>(1.0);
245 template <
class LhsEval>
246 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
247 const IntensiveQuantities& intQuants)
249 if constexpr (enableFoam) {
250 const auto& fs = intQuants.fluidState();
252 LhsEval surfaceVolume = Toolbox::template decay<LhsEval>(intQuants.porosity());
253 if (params_.transport_phase_ == Phase::WATER) {
254 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
255 * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx)));
256 }
else if (params_.transport_phase_ == Phase::GAS) {
257 surfaceVolume *= (Toolbox::template decay<LhsEval>(fs.saturation(gasPhaseIdx))
258 * Toolbox::template decay<LhsEval>(fs.invB(gasPhaseIdx)));
259 }
else if (params_.transport_phase_ == Phase::SOLVENT) {
260 if constexpr (enableSolvent) {
261 surfaceVolume *= (Toolbox::template decay<LhsEval>( intQuants.solventSaturation())
262 * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor()));
265 throw std::runtime_error(
"Transport phase is GAS/WATER/SOLVENT");
269 surfaceVolume = max(surfaceVolume, 1e-10);
272 const LhsEval freeFoam = surfaceVolume
273 * Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
276 const LhsEval adsorbedFoam =
277 Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
278 * Toolbox::template decay<LhsEval>(intQuants.foamRockDensity())
279 * Toolbox::template decay<LhsEval>(intQuants.foamAdsorbed());
281 LhsEval accumulationFoam = freeFoam + adsorbedFoam;
282 storage[contiFoamEqIdx] += accumulationFoam;
287 [[maybe_unused]]
const ElementContext& elemCtx,
288 [[maybe_unused]]
unsigned scvfIdx,
289 [[maybe_unused]]
unsigned timeIdx)
292 if constexpr (enableFoam) {
293 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
294 const unsigned inIdx = extQuants.interiorIndex();
301 const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
302 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
303 if (upIdx == inIdx) {
304 flux[contiFoamEqIdx] =
305 extQuants.volumeFlux(waterPhaseIdx)
306 *up.fluidState().invB(waterPhaseIdx)
307 *up.foamConcentration();
309 flux[contiFoamEqIdx] =
310 extQuants.volumeFlux(waterPhaseIdx)
311 *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
312 *decay<Scalar>(up.foamConcentration());
317 const unsigned upIdx = extQuants.upstreamIndex(gasPhaseIdx);
318 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
319 if (upIdx == inIdx) {
320 flux[contiFoamEqIdx] =
321 extQuants.volumeFlux(gasPhaseIdx)
322 *up.fluidState().invB(gasPhaseIdx)
323 *up.foamConcentration();
325 flux[contiFoamEqIdx] =
326 extQuants.volumeFlux(gasPhaseIdx)
327 *decay<Scalar>(up.fluidState().invB(gasPhaseIdx))
328 *decay<Scalar>(up.foamConcentration());
332 case Phase::SOLVENT: {
333 if constexpr (enableSolvent) {
334 const unsigned upIdx = extQuants.solventUpstreamIndex();
335 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
336 if (upIdx == inIdx) {
337 flux[contiFoamEqIdx] =
338 extQuants.solventVolumeFlux()
339 *up.solventInverseFormationVolumeFactor()
340 *up.foamConcentration();
342 flux[contiFoamEqIdx] =
343 extQuants.solventVolumeFlux()
344 *decay<Scalar>(up.solventInverseFormationVolumeFactor())
345 *decay<Scalar>(up.foamConcentration());
348 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
353 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
367 return static_cast<Scalar
>(0.0);
370 template <
class DofEntity>
372 [[maybe_unused]] std::ostream& outstream,
373 [[maybe_unused]]
const DofEntity& dof)
375 if constexpr (enableFoam) {
376 unsigned dofIdx = model.dofMapper().index(dof);
377 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
378 outstream << priVars[foamConcentrationIdx];
382 template <
class DofEntity>
384 [[maybe_unused]] std::istream& instream,
385 [[maybe_unused]]
const DofEntity& dof)
387 if constexpr (enableFoam) {
388 unsigned dofIdx = model.dofMapper().index(dof);
389 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
390 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
392 instream >> priVars0[foamConcentrationIdx];
395 priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
403 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
404 return params_.foamRockDensity_[satnumRegionIdx];
411 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
412 return params_.foamAllowDesorption_[satnumRegionIdx];
419 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
420 return params_.adsorbedFoamTable_[satnumRegionIdx];
427 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
428 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
433 const unsigned scvIdx,
434 const unsigned timeIdx)
436 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
437 return params_.foamCoefficients_[satnumRegionIdx];
441 return params_.transport_phase_;
448template <
class TypeTag,
bool enableFoam>
449BlackOilFoamParams<typename BlackOilFoamModule<TypeTag, enableFoam>::Scalar>
450BlackOilFoamModule<TypeTag, enableFoam>::params_;
459template <class TypeTag, bool enableFoam = getPropValue<TypeTag, Properties::EnableFoam>()>
474 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
475 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
477 static constexpr int foamConcentrationIdx = Indices::foamConcentrationIdx;
478 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
479 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
480 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
493 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
495 const auto& fs =
asImp_().fluidState_;
498 Evaluation mobilityReductionFactor = 1.0;
505 const Scalar fm_mob = foamCoefficients.fm_mob;
507 const Scalar fm_surf = foamCoefficients.fm_surf;
508 const Scalar ep_surf = foamCoefficients.ep_surf;
510 const Scalar fm_oil = foamCoefficients.fm_oil;
511 const Scalar fl_oil = foamCoefficients.fl_oil;
512 const Scalar ep_oil = foamCoefficients.ep_oil;
514 const Scalar fm_dry = foamCoefficients.fm_dry;
515 const Scalar ep_dry = foamCoefficients.ep_dry;
517 const Scalar fm_cap = foamCoefficients.fm_cap;
518 const Scalar ep_cap = foamCoefficients.ep_cap;
521 const Evaluation Ca = 1e10;
522 const Evaluation S_o = fs.saturation(oilPhaseIdx);
523 const Evaluation S_w = fs.saturation(waterPhaseIdx);
525 Evaluation F1 = pow(C_surf/fm_surf, ep_surf);
526 Evaluation F2 = pow((fm_oil-S_o)/(fm_oil-fl_oil), ep_oil);
527 Evaluation F3 = pow(fm_cap/Ca, ep_cap);
528 Evaluation F7 = 0.5 + atan(ep_dry*(S_w-fm_dry))/M_PI;
530 mobilityReductionFactor = 1./(1. + fm_mob*F1*F2*F3*F7);
542 asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
546 asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
549 case Phase::SOLVENT: {
550 if constexpr (enableSolvent) {
551 asImp_().solventMobility_ *= mobilityReductionFactor;
553 throw std::runtime_error(
"Foam transport phase is SOLVENT but SOLVENT is not activated.");
558 throw std::runtime_error(
"Foam transport phase must be GAS/WATER/SOLVENT.");
567 throw std::runtime_error(
"Foam module does not support the 'no desorption' option.");
582 {
return *
static_cast<Implementation*
>(
this); }
589template <
class TypeTag>
604 {
throw std::runtime_error(
"foamConcentration() called but foam is disabled"); }
607 {
throw std::runtime_error(
"foamRockDensity() called but foam is disabled"); }
610 {
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:606
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:603
void foamPropertiesUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilfoammodules.hh:597
Scalar foamAdsorbed() const
Definition: blackoilfoammodules.hh:609
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilfoammodules.hh:461
const Evaluation & foamConcentration() const
Definition: blackoilfoammodules.hh:571
Implementation & asImp_()
Definition: blackoilfoammodules.hh:581
Scalar foamRockDensity_
Definition: blackoilfoammodules.hh:585
const Evaluation & foamAdsorbed() const
Definition: blackoilfoammodules.hh:577
Evaluation foamAdsorbed_
Definition: blackoilfoammodules.hh:586
Evaluation foamConcentration_
Definition: blackoilfoammodules.hh:584
Scalar foamRockDensity() const
Definition: blackoilfoammodules.hh:574
void foamPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilfoammodules.hh:489
Contains the high level supplements required to extend the black oil model to include the effects of ...
Definition: blackoilfoammodules.hh:57
static bool eqApplies(unsigned eqIdx)
Definition: blackoilfoammodules.hh:220
static void registerOutputModules(Model &, Simulator &)
Register all foam specific VTK and ECL output modules.
Definition: blackoilfoammodules.hh:187
static void registerParameters()
Register all run-time parameters for the black-oil foam module.
Definition: blackoilfoammodules.hh:180
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilfoammodules.hh:206
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilfoammodules.hh:236
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:383
static std::string eqName(unsigned eqIdx)
Definition: blackoilfoammodules.hh:229
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilfoammodules.hh:198
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilfoammodules.hh:246
static bool foamAllowDesorption(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:407
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilfoammodules.hh:212
static const BlackOilFoamParams< Scalar >::FoamCoefficients & foamCoefficients(const ElementContext &elemCtx, const unsigned scvIdx, const unsigned timeIdx)
Definition: blackoilfoammodules.hh:432
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:286
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilfoammodules.hh:362
static const Scalar foamRockDensity(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:399
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilfoammodules.hh:371
static const TabulatedFunction & adsorbedFoamTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:415
static Phase transportPhase()
Definition: blackoilfoammodules.hh:440
static const TabulatedFunction & gasMobilityMultiplierTable(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Definition: blackoilfoammodules.hh:423
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:242
Definition: blackoilfoamparams.hh:59
Struct holding the parameters for the BlackoilFoamModule class.
Definition: blackoilfoamparams.hh:40
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilfoamparams.hh:41