blackoilfoammodules.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef EWOMS_BLACK_OIL_FOAM_MODULE_HH
29#define EWOMS_BLACK_OIL_FOAM_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
33#include <opm/common/OpmLog/OpmLog.hpp>
34
35#include <opm/input/eclipse/EclipseState/Phase.hpp>
36
39
42
43#include <cassert>
44#include <istream>
45#include <ostream>
46#include <stdexcept>
47#include <string>
48
49namespace Opm {
50
56template <class TypeTag, bool enableFoamV = getPropValue<TypeTag, Properties::EnableFoam>()>
58{
71
72 using Toolbox = MathToolbox<Evaluation>;
73
74 using TabulatedFunction = typename BlackOilFoamParams<Scalar>::TabulatedFunction;
75
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;
80
81 static constexpr unsigned enableFoam = enableFoamV;
82
83 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
84 static constexpr unsigned numPhases = FluidSystem::numPhases;
85
86 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
87
88public:
91 {
92 params_ = params;
93 }
94
98 static void registerParameters()
99 {}
100
104 static void registerOutputModules(Model&,
105 Simulator&)
106 {
107 if constexpr (enableFoam) {
108 if (Parameters::Get<Parameters::EnableVtkOutput>()) {
109 OpmLog::warning("VTK output requested, currently unsupported by the foam module.");
110 }
111 }
112 //model.addOutputModule(new VtkBlackOilFoamModule<TypeTag>(simulator));
113 }
114
115 static bool primaryVarApplies(unsigned pvIdx)
116 {
117 if constexpr (enableFoam) {
118 return pvIdx == foamConcentrationIdx;
119 }
120 else {
121 return false;
122 }
123 }
124
125 static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
126 {
127 assert(primaryVarApplies(pvIdx));
128 return "foam_concentration";
129 }
130
131 static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
132 {
133 assert(primaryVarApplies(pvIdx));
134
135 // TODO: it may be beneficial to chose this differently.
136 return static_cast<Scalar>(1.0);
137 }
138
139 static bool eqApplies(unsigned eqIdx)
140 {
141 if constexpr (enableFoam) {
142 return eqIdx == contiFoamEqIdx;
143 }
144 else {
145 return false;
146 }
147 }
148
149 static std::string eqName([[maybe_unused]] unsigned eqIdx)
150 {
151 assert(eqApplies(eqIdx));
152
153 return "conti^foam";
154 }
155
156 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
157 {
158 assert(eqApplies(eqIdx));
159
160 // TODO: it may be beneficial to chose this differently.
161 return static_cast<Scalar>(1.0);
162 }
163
164 // must be called after water storage is computed
165 template <class LhsEval>
166 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
167 const IntensiveQuantities& intQuants)
168 {
169 if constexpr (enableFoam) {
170 const auto& fs = intQuants.fluidState();
171
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());
183 }
184 } else {
185 throw std::runtime_error("Transport phase is GAS/WATER/SOLVENT");
186 }
187
188 // Avoid singular matrix if no gas is present.
189 surfaceVolume = max(surfaceVolume, 1e-10);
190
191 // Foam/surfactant in free phase.
192 const LhsEval freeFoam = surfaceVolume *
193 Toolbox::template decay<LhsEval>(intQuants.foamConcentration());
194
195 // Adsorbed foam/surfactant.
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());
200
201 const LhsEval accumulationFoam = freeFoam + adsorbedFoam;
202 storage[contiFoamEqIdx] += accumulationFoam;
203 }
204 }
205
206 static void computeFlux([[maybe_unused]] RateVector& flux,
207 [[maybe_unused]] const ElementContext& elemCtx,
208 [[maybe_unused]] unsigned scvfIdx,
209 [[maybe_unused]] unsigned timeIdx)
210 {
211 if constexpr (enableFoam) {
212 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
213 const unsigned inIdx = extQuants.interiorIndex();
214
215 // The effect of the mobility reduction factor is
216 // incorporated in the mobility for the relevant phase,
217 // so fluxes do not need modification here.
218 switch (transportPhase()) {
219 case Phase::WATER: {
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();
227 } else {
228 flux[contiFoamEqIdx] =
229 extQuants.volumeFlux(waterPhaseIdx) *
230 decay<Scalar>(up.fluidState().invB(waterPhaseIdx)) *
231 decay<Scalar>(up.foamConcentration());
232 }
233 break;
234 }
235 case Phase::GAS: {
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();
243 } else {
244 flux[contiFoamEqIdx] =
245 extQuants.volumeFlux(gasPhaseIdx) *
246 decay<Scalar>(up.fluidState().invB(gasPhaseIdx)) *
247 decay<Scalar>(up.foamConcentration());
248 }
249 break;
250 }
251 case Phase::SOLVENT:
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();
260 } else {
261 flux[contiFoamEqIdx] =
262 extQuants.solventVolumeFlux() *
263 decay<Scalar>(up.solventInverseFormationVolumeFactor()) *
264 decay<Scalar>(up.foamConcentration());
265 }
266 } else {
267 throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
268 }
269 break;
270 default:
271 throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
272 }
273 }
274 }
275
279 static Scalar computeUpdateError(const PrimaryVariables&,
280 const EqVector&)
281 {
282 // do not consider the change of foam primary variables for convergence
283 // TODO: maybe this should be changed
284 return static_cast<Scalar>(0.0);
285 }
286
287 template <class DofEntity>
288 static void serializeEntity([[maybe_unused]] const Model& model,
289 [[maybe_unused]] std::ostream& outstream,
290 [[maybe_unused]] const DofEntity& dof)
291 {
292 if constexpr (enableFoam) {
293 const unsigned dofIdx = model.dofMapper().index(dof);
294 const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
295 outstream << priVars[foamConcentrationIdx];
296 }
297 }
298
299 template <class DofEntity>
300 static void deserializeEntity([[maybe_unused]] Model& model,
301 [[maybe_unused]] std::istream& instream,
302 [[maybe_unused]] const DofEntity& dof)
303 {
304 if constexpr (enableFoam) {
305 const unsigned dofIdx = model.dofMapper().index(dof);
306 PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
307 PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
308
309 instream >> priVars0[foamConcentrationIdx];
310
311 // set the primary variables for the beginning of the current time step.
312 priVars1[foamConcentrationIdx] = priVars0[foamConcentrationIdx];
313 }
314 }
315
316 static const Scalar foamRockDensity(const ElementContext& elemCtx,
317 unsigned scvIdx,
318 unsigned timeIdx)
319 {
320 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
321 return params_.foamRockDensity_[satnumRegionIdx];
322 }
323
324 static bool foamAllowDesorption(const ElementContext& elemCtx,
325 unsigned scvIdx,
326 unsigned timeIdx)
327 {
328 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
329 return params_.foamAllowDesorption_[satnumRegionIdx];
330 }
331
332 static const TabulatedFunction& adsorbedFoamTable(const ElementContext& elemCtx,
333 unsigned scvIdx,
334 unsigned timeIdx)
335 {
336 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
337 return params_.adsorbedFoamTable_[satnumRegionIdx];
338 }
339
340 static const TabulatedFunction& gasMobilityMultiplierTable(const ElementContext& elemCtx,
341 unsigned scvIdx,
342 unsigned timeIdx)
343 {
344 const unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
345 return params_.gasMobilityMultiplierTable_[pvtnumRegionIdx];
346 }
347
349 foamCoefficients(const ElementContext& elemCtx,
350 const unsigned scvIdx,
351 const unsigned timeIdx)
352 {
353 const unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
354 return params_.foamCoefficients_[satnumRegionIdx];
355 }
356
357 static Phase transportPhase()
358 { return params_.transport_phase_; }
359
360private:
361 static BlackOilFoamParams<Scalar> params_;
362};
363
364template <class TypeTag, bool enableFoam>
365BlackOilFoamParams<typename BlackOilFoamModule<TypeTag, enableFoam>::Scalar>
366BlackOilFoamModule<TypeTag, enableFoam>::params_;
367
368template <class TypeTag, bool enableFoam>
370
378template <class TypeTag>
379class BlackOilFoamIntensiveQuantities<TypeTag, /*enableFoam=*/true>
380{
382
390
392
393 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
394 enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
395
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;
400
401public:
407 void foamPropertiesUpdate_(const ElementContext& elemCtx,
408 unsigned dofIdx,
409 unsigned timeIdx)
410 {
411 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
412 foamConcentration_ = priVars.makeEvaluation(foamConcentrationIdx, timeIdx);
413 const auto& fs = asImp_().fluidState_;
414
415 // Compute gas mobility reduction factor
416 Evaluation mobilityReductionFactor = 1.0;
417 if constexpr (false) {
418 // The functional model is used.
419 // TODO: allow this model.
420 // In order to do this we must allow transport to be in the water phase, not just the gas phase.
421 const auto& foamCoefficients = FoamModule::foamCoefficients(elemCtx, dofIdx, timeIdx);
422
423 const Scalar fm_mob = foamCoefficients.fm_mob;
424
425 const Scalar fm_surf = foamCoefficients.fm_surf;
426 const Scalar ep_surf = foamCoefficients.ep_surf;
427
428 const Scalar fm_oil = foamCoefficients.fm_oil;
429 const Scalar fl_oil = foamCoefficients.fl_oil;
430 const Scalar ep_oil = foamCoefficients.ep_oil;
431
432 const Scalar fm_dry = foamCoefficients.fm_dry;
433 const Scalar ep_dry = foamCoefficients.ep_dry;
434
435 const Scalar fm_cap = foamCoefficients.fm_cap;
436 const Scalar ep_cap = foamCoefficients.ep_cap;
437
438 const Evaluation C_surf = foamConcentration_;
439 const Evaluation Ca = 1e10; // TODO: replace with proper capillary number.
440 const Evaluation S_o = fs.saturation(oilPhaseIdx);
441 const Evaluation S_w = fs.saturation(waterPhaseIdx);
442
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;
447
448 mobilityReductionFactor = 1. / (1. + fm_mob * F1 * F2 * F3 * F7);
449 } else {
450 // The tabular model is used.
451 // Note that the current implementation only includes the effect of foam concentration (FOAMMOB),
452 // and not the optional pressure dependence (FOAMMOBP) or shear dependence (FOAMMOBS).
453 const auto& gasMobilityMultiplier = FoamModule::gasMobilityMultiplierTable(elemCtx, dofIdx, timeIdx);
454 mobilityReductionFactor = gasMobilityMultiplier.eval(foamConcentration_, /* extrapolate = */ true);
455 }
456
457 // adjust mobility
458 switch (FoamModule::transportPhase()) {
459 case Phase::WATER:
460 asImp_().mobility_[waterPhaseIdx] *= mobilityReductionFactor;
461 break;
462 case Phase::GAS:
463 asImp_().mobility_[gasPhaseIdx] *= mobilityReductionFactor;
464 break;
465 case Phase::SOLVENT:
466 if constexpr (enableSolvent) {
467 asImp_().solventMobility_ *= mobilityReductionFactor;
468 } else {
469 throw std::runtime_error("Foam transport phase is SOLVENT but SOLVENT is not activated.");
470 }
471 break;
472 default:
473 throw std::runtime_error("Foam transport phase must be GAS/WATER/SOLVENT.");
474 }
475
476 foamRockDensity_ = FoamModule::foamRockDensity(elemCtx, dofIdx, timeIdx);
477
478 const auto& adsorbedFoamTable = FoamModule::adsorbedFoamTable(elemCtx, dofIdx, timeIdx);
479 foamAdsorbed_ = adsorbedFoamTable.eval(foamConcentration_, /*extrapolate=*/true);
480 if (!FoamModule::foamAllowDesorption(elemCtx, dofIdx, timeIdx)) {
481 throw std::runtime_error("Foam module does not support the 'no desorption' option.");
482 }
483 }
484
485 const Evaluation& foamConcentration() const
486 { return foamConcentration_; }
487
488 Scalar foamRockDensity() const
489 { return foamRockDensity_; }
490
491 const Evaluation& foamAdsorbed() const
492 { return foamAdsorbed_; }
493
494protected:
495 Implementation& asImp_()
496 { return *static_cast<Implementation*>(this); }
497
500 Evaluation foamAdsorbed_;
501};
502
503template <class TypeTag>
505{
509
510public:
511 void foamPropertiesUpdate_(const ElementContext&,
512 unsigned,
513 unsigned)
514 {}
515
516 const Evaluation& foamConcentration() const
517 { throw std::runtime_error("foamConcentration() called but foam is disabled"); }
518
519 Scalar foamRockDensity() const
520 { throw std::runtime_error("foamRockDensity() called but foam is disabled"); }
521
522 Scalar foamAdsorbed() const
523 { throw std::runtime_error("foamAdsorbed() called but foam is disabled"); }
524};
525
526} // namespace Opm
527
528#endif
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 > &&params)
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