25 #ifndef OPM_IMMISCIBLE_FLASH_HPP
26 #define OPM_IMMISCIBLE_FLASH_HPP
28 #include <opm/common/ErrorMacros.hpp>
29 #include <opm/common/Exceptions.hpp>
32 #include <dune/common/fvector.hh>
33 #include <dune/common/fmatrix.hh>
66 template <
class Scalar,
class Flu
idSystem>
69 static const int numPhases = FluidSystem::numPhases;
70 static const int numComponents = FluidSystem::numComponents;
71 static_assert(numPhases == numComponents,
72 "Immiscibility assumes that the number of phases"
73 " is equal to the number of components");
75 typedef typename FluidSystem::ParameterCache ParameterCache;
77 static const int numEq = numPhases;
79 typedef Dune::FieldMatrix<Scalar, numEq, numEq> Matrix;
80 typedef Dune::FieldVector<Scalar, numEq> Vector;
83 typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
88 template <
class Flu
idState>
89 static void guessInitial(FluidState &fluidState,
91 const ComponentVector &globalMolarities)
95 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
96 sumMoles += globalMolarities[compIdx];
98 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
100 fluidState.setPressure(phaseIdx, 2e5);
103 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
113 template <
class MaterialLaw,
class Flu
idState>
114 static void solve(FluidState &fluidState,
115 ParameterCache ¶mCache,
116 const typename MaterialLaw::Params &matParams,
117 const ComponentVector &globalMolarities)
119 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-25);
124 bool allIncompressible =
true;
125 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
126 if (FluidSystem::isCompressible(phaseIdx)) {
127 allIncompressible =
false;
132 if (allIncompressible) {
133 paramCache.updateAll(fluidState);
134 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
135 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
136 fluidState.setDensity(phaseIdx, rho);
139 globalMolarities[phaseIdx]
140 / fluidState.molarDensity(phaseIdx);
141 fluidState.setSaturation(phaseIdx, saturation);
160 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
163 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
165 linearize_<MaterialLaw>(J, b, fluidState, paramCache, matParams, globalMolarities);
172 try { J.solve(deltaX, b); }
173 catch (Dune::FMatrixError e)
180 throw Opm::NumericalProblem(e.what());
206 Scalar relError = update_<MaterialLaw>(fluidState, paramCache, matParams, deltaX);
220 OPM_THROW(Opm::NumericalProblem,
221 "Flash calculation failed."
222 " {c_alpha^kappa} = {" << globalMolarities <<
"}, T = "
223 << fluidState.temperature(0));
228 template <
class Flu
idState>
229 static void printFluidState_(
const FluidState &fs)
231 std::cout <<
"saturations: ";
232 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
233 std::cout << fs.saturation(phaseIdx) <<
" ";
236 std::cout <<
"pressures: ";
237 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
238 std::cout << fs.pressure(phaseIdx) <<
" ";
241 std::cout <<
"densities: ";
242 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
243 std::cout << fs.density(phaseIdx) <<
" ";
246 std::cout <<
"global component molarities: ";
247 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
249 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
250 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
252 std::cout << sum <<
" ";
257 template <
class MaterialLaw,
class Flu
idState>
258 static void linearize_(Matrix &J,
260 FluidState &fluidState,
261 ParameterCache ¶mCache,
262 const typename MaterialLaw::Params &matParams,
263 const ComponentVector &globalMolarities)
265 FluidState origFluidState(fluidState);
266 ParameterCache origParamCache(paramCache);
274 calculateDefect_(b, fluidState, fluidState, globalMolarities);
278 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
286 Scalar xI = getQuantity_(fluidState, pvIdx);
287 const Scalar eps = 1e-10/quantityWeight_(fluidState, pvIdx);
288 setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, xI + eps);
289 assert(
std::abs(getQuantity_(fluidState, pvIdx) - (xI + eps))
293 calculateDefect_(tmp, origFluidState, fluidState, globalMolarities);
297 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
298 J[eqIdx][pvIdx] = tmp[eqIdx];
301 fluidState = origFluidState;
302 paramCache = origParamCache;
309 template <
class Flu
idState>
310 static void calculateDefect_(Vector &b,
312 const FluidState &fluidState,
313 const ComponentVector &globalMolarities)
316 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
318 fluidState.saturation(phaseIdx)
319 * fluidState.molarity(phaseIdx, phaseIdx);
321 b[phaseIdx] -= globalMolarities[phaseIdx];
325 template <
class MaterialLaw,
class Flu
idState>
326 static Scalar update_(FluidState &fluidState,
327 ParameterCache ¶mCache,
328 const typename MaterialLaw::Params &matParams,
329 const Vector &deltaX)
332 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
333 Scalar tmp = getQuantity_(fluidState, pvIdx);
334 Scalar delta = deltaX[pvIdx];
336 relError =
std::max(relError,
std::abs(delta)*quantityWeight_(fluidState, pvIdx));
338 if (isSaturationIdx_(pvIdx)) {
343 else if (isPressureIdx_(pvIdx)) {
346 delta =
std::min(0.30*fluidState.pressure(0),
std::max(-0.30*fluidState.pressure(0), delta));
349 setQuantityRaw_(fluidState, pvIdx, tmp - delta);
352 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
357 template <
class MaterialLaw,
class Flu
idState>
358 static void completeFluidState_(FluidState &fluidState,
359 ParameterCache ¶mCache,
360 const typename MaterialLaw::Params &matParams)
365 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
366 sumSat += fluidState.saturation(phaseIdx);
371 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
373 Scalar S = fluidState.saturation(phaseIdx);
374 fluidState.setSaturation(phaseIdx, S/sumSat);
380 fluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
386 MaterialLaw::capillaryPressures(pC, matParams, fluidState);
387 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
388 fluidState.setPressure(phaseIdx,
389 fluidState.pressure(0)
390 + (pC[phaseIdx] - pC[0]));
393 paramCache.updateAll(fluidState, ParameterCache::Temperature|ParameterCache::Composition);
396 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
397 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
398 fluidState.setDensity(phaseIdx, rho);
402 static bool isPressureIdx_(
unsigned pvIdx)
403 {
return pvIdx == 0; }
405 static bool isSaturationIdx_(
unsigned pvIdx)
406 {
return 1 <= pvIdx && pvIdx < numPhases; }
409 template <
class Flu
idState>
410 static Scalar getQuantity_(
const FluidState &fs,
unsigned pvIdx)
412 assert(pvIdx < numEq);
416 unsigned phaseIdx = 0;
417 return fs.pressure(phaseIdx);
421 assert(pvIdx < numPhases);
422 unsigned phaseIdx = pvIdx - 1;
423 return fs.saturation(phaseIdx);
428 template <
class MaterialLaw,
class Flu
idState>
429 static void setQuantity_(FluidState &fs,
430 ParameterCache ¶mCache,
431 const typename MaterialLaw::Params &matParams,
435 assert(0 <= pvIdx && pvIdx < numEq);
439 Scalar delta = value - fs.pressure(0);
443 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
444 fs.setPressure(phaseIdx, fs.pressure(phaseIdx) + delta);
445 paramCache.updateAllPressures(fs);
448 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
449 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
450 fs.setDensity(phaseIdx, rho);
454 assert(pvIdx < numPhases);
457 Scalar delta = value - fs.saturation(pvIdx - 1);
458 fs.setSaturation(pvIdx - 1, value);
462 fs.setSaturation(numPhases - 1,
463 fs.saturation(numPhases - 1) - delta);
468 MaterialLaw::capillaryPressures(pC, matParams, fs);
469 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
470 fs.setPressure(phaseIdx,
472 + (pC[phaseIdx] - pC[0]));
473 paramCache.updateAllPressures(fs);
476 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
477 Scalar rho = FluidSystem::density(fs, paramCache, phaseIdx);
478 fs.setDensity(phaseIdx, rho);
484 template <
class Flu
idState>
485 static void setQuantityRaw_(FluidState &fs,
unsigned pvIdx, Scalar value)
487 assert(pvIdx < numEq);
491 unsigned phaseIdx = 0;
492 fs.setPressure(phaseIdx, value);
496 assert(pvIdx < numPhases);
497 unsigned phaseIdx = pvIdx - 1;
502 fs.setSaturation(phaseIdx, value);
506 template <
class Flu
idState>
507 static Scalar quantityWeight_(
const FluidState &,
unsigned pvIdx)
514 assert(pvIdx < numPhases);
void SetUndefined(const T &value OPM_UNUSED)
Make the memory on which an object resides undefined in valgrind runs.
Definition: Valgrind.hpp:138
bool CheckDefined(const T &value OPM_UNUSED)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
Definition: Air_Mesitylene.hpp:31
Some templates to wrap the valgrind client request macros.
Evaluation< Scalar, VarSetTag, numVars > max(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:114
Determines the pressures and saturations of all fluid phases given the total mass of all components...
Definition: ImmiscibleFlash.hpp:67
Evaluation< Scalar, VarSetTag, numVars > abs(const Evaluation< Scalar, VarSetTag, numVars > &)
Definition: Math.hpp:41
Evaluation< Scalar, VarSetTag, numVars > min(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:61