74 using Toolbox = MathToolbox<Evaluation>;
76 using TabulatedFunction =
typename BlackOilFoamParams<Scalar>::TabulatedFunction;
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;
86 static constexpr unsigned numPhases = FluidSystem::numPhases;
110 if constexpr (enableFoam) {
111 if (Parameters::Get<Parameters::EnableVtkOutput>()) {
112 OpmLog::warning(
"VTK output requested, currently unsupported by the foam module.");
118 static bool primaryVarApplies(
unsigned pvIdx)
120 if constexpr (enableFoam)
121 return pvIdx == foamConcentrationIdx;
126 static std::string primaryVarName([[maybe_unused]]
unsigned pvIdx)
128 assert(primaryVarApplies(pvIdx));
129 return "foam_concentration";
132 static Scalar primaryVarWeight([[maybe_unused]]
unsigned pvIdx)
134 assert(primaryVarApplies(pvIdx));
137 return static_cast<Scalar
>(1.0);
140 static bool eqApplies(
unsigned eqIdx)
142 if constexpr (enableFoam)
143 return eqIdx == contiFoamEqIdx;
149 static std::string eqName([[maybe_unused]]
unsigned eqIdx)
151 assert(eqApplies(eqIdx));
156 static Scalar eqWeight([[maybe_unused]]
unsigned eqIdx)
158 assert(eqApplies(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;
206 static void computeFlux([[maybe_unused]] RateVector& flux,
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();
219 switch (transportPhase()) {
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>
291 static void serializeEntity([[maybe_unused]]
const Model& model,
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>
303 static void deserializeEntity([[maybe_unused]] Model& model,
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];
319 static const Scalar foamRockDensity(
const ElementContext& elemCtx,
323 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
324 return params_.foamRockDensity_[satnumRegionIdx];
327 static bool foamAllowDesorption(
const ElementContext& elemCtx,
331 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
332 return params_.foamAllowDesorption_[satnumRegionIdx];
335 static const TabulatedFunction& adsorbedFoamTable(
const ElementContext& elemCtx,
339 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
340 return params_.adsorbedFoamTable_[satnumRegionIdx];
343 static const TabulatedFunction& gasMobilityMultiplierTable(
const ElementContext& elemCtx,
347 unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
348 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
351 static const typename BlackOilFoamParams<Scalar>::FoamCoefficients&
352 foamCoefficients(
const ElementContext& elemCtx,
353 const unsigned scvIdx,
354 const unsigned timeIdx)
356 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
357 return params_.foamCoefficients_[satnumRegionIdx];
360 static Phase transportPhase() {
361 return params_.transport_phase_;
365 static BlackOilFoamParams<Scalar> params_;
413 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
414 foamConcentration_ = priVars.makeEvaluation(foamConcentrationIdx, timeIdx);
415 const auto& fs = asImp_().fluidState_;
418 Evaluation mobilityReductionFactor = 1.0;
423 const auto& foamCoefficients = FoamModule::foamCoefficients(elemCtx, dofIdx, timeIdx);
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;
440 const Evaluation C_surf = foamConcentration_;
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);
455 const auto& gasMobilityMultiplier = FoamModule::gasMobilityMultiplierTable(elemCtx, dofIdx, timeIdx);
456 mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_,
true);
460 switch (FoamModule::transportPhase()) {
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.");
482 foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
484 const auto& adsorbedFoamTable = FoamModule::adsorbedFoamTable(elemCtx, dofIdx, timeIdx);
485 foamAdsorbed_ = adsorbedFoamTable.eval(foamConcentration_,
true);
486 if (!FoamModule::foamAllowDesorption(elemCtx, dofIdx, timeIdx)) {
487 throw std::runtime_error(
"Foam module does not support the 'no desorption' option.");