25 #ifndef OPM_NCP_FLASH_HPP
26 #define OPM_NCP_FLASH_HPP
28 #include <dune/common/fvector.hh>
29 #include <dune/common/fmatrix.hh>
34 #include <opm/common/ErrorMacros.hpp>
35 #include <opm/common/Exceptions.hpp>
83 template <
class Scalar,
class Flu
idSystem>
86 enum { numPhases = FluidSystem::numPhases };
87 enum { numComponents = FluidSystem::numComponents };
89 typedef typename FluidSystem::ParameterCache ParameterCache;
91 static const int numEq = numPhases*(numComponents + 1);
97 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
99 ParameterCache ¶mCache,
100 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
103 Evaluation sumMoles = 0;
104 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
105 sumMoles += globalMolarities[compIdx];
107 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
109 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
110 fluidState.setMoleFraction(phaseIdx,
112 globalMolarities[compIdx]/sumMoles);
115 fluidState.setPressure(phaseIdx, 1.0135e5);
118 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
122 paramCache.updateAll(fluidState);
123 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
124 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
125 const typename FluidState::Scalar phi =
126 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
127 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
138 template <
class MaterialLaw,
class Flu
idState>
139 static void solve(FluidState &fluidState,
140 ParameterCache ¶mCache,
141 const typename MaterialLaw::Params &matParams,
142 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
143 Scalar tolerance = 0.0)
145 typedef typename FluidState::Scalar Evaluation;
146 typedef Dune::FieldMatrix<Evaluation, numEq, numEq> Matrix;
147 typedef Dune::FieldVector<Evaluation, numEq> Vector;
149 Dune::FMatrixPrecision<Scalar>::set_singular_limit(1e-35);
151 if (tolerance <= 0.0) {
152 tolerance = std::min<Scalar>(1e-10,
154 std::numeric_limits<Scalar>::epsilon()));
173 completeFluidState_<MaterialLaw>(fluidState,
185 for (
int nIdx = 0; nIdx < nMax; ++nIdx) {
187 linearize_<MaterialLaw>(J,
199 try { J.solve(deltaX, b); }
200 catch (Dune::FMatrixError e)
209 throw Opm::NumericalProblem(e.what());
236 Scalar relError = update_<MaterialLaw>(fluidState, paramCache, matParams, deltaX);
250 OPM_THROW(NumericalProblem,
251 "Flash calculation failed."
252 " {c_alpha^kappa} = {" << globalMolarities <<
"}, T = "
253 << fluidState.temperature(0));
263 template <
class Flu
idState,
class ComponentVector>
264 static void solve(FluidState &fluidState,
265 const ComponentVector &globalMolarities,
266 Scalar tolerance = 0.0)
268 ParameterCache paramCache;
269 paramCache.updateAll(fluidState);
273 typedef typename MaterialLaw::Params MaterialLawParams;
275 MaterialLawParams matParams;
276 solve<MaterialLaw>(fluidState, paramCache, matParams, globalMolarities, tolerance);
281 template <
class Flu
idState>
284 std::cout <<
"saturations: ";
285 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
286 std::cout << fluidState.saturation(phaseIdx) <<
" ";
289 std::cout <<
"pressures: ";
290 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
291 std::cout << fluidState.pressure(phaseIdx) <<
" ";
294 std::cout <<
"densities: ";
295 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
296 std::cout << fluidState.density(phaseIdx) <<
" ";
299 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
300 std::cout <<
"composition " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
301 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
302 std::cout << fluidState.moleFraction(phaseIdx, compIdx) <<
" ";
307 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
308 std::cout <<
"fugacities " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
309 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
310 std::cout << fluidState.fugacity(phaseIdx, compIdx) <<
" ";
315 std::cout <<
"global component molarities: ";
316 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
318 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
319 sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
321 std::cout << sum <<
" ";
326 template <
class MaterialLaw,
330 class ComponentVector>
333 FluidState &fluidState,
334 ParameterCache ¶mCache,
335 const typename MaterialLaw::Params &matParams,
336 const ComponentVector &globalMolarities)
338 typedef typename FluidState::Scalar Evaluation;
340 FluidState origFluidState(fluidState);
341 ParameterCache origParamCache(paramCache);
355 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
363 const Evaluation& x_i =
getQuantity_(fluidState, pvIdx);
364 const Scalar eps = std::numeric_limits<Scalar>::epsilon()*1e7/(
quantityWeight_(fluidState, pvIdx));
366 setQuantity_<MaterialLaw>(fluidState, paramCache, matParams, pvIdx, x_i + eps);
371 for (
unsigned i = 0; i < numEq; ++ i)
375 for (
unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx)
376 J[eqIdx][pvIdx] = tmp[eqIdx];
379 fluidState = origFluidState;
380 paramCache = origParamCache;
387 template <
class Flu
idState,
class Vector,
class ComponentVector>
389 const FluidState &fluidStateEval,
390 const FluidState &fluidState,
391 const ComponentVector &globalMolarities)
393 typedef typename FluidState::Scalar Evaluation;
398 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
399 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
401 fluidState.fugacity(0, compIdx)
402 - fluidState.fugacity(phaseIdx, compIdx);
407 assert(eqIdx == numComponents*(numPhases - 1));
413 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
415 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
417 fluidState.saturation(phaseIdx)
418 * fluidState.molarity(phaseIdx, compIdx);
421 b[eqIdx] -= globalMolarities[compIdx];
426 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
427 Evaluation sumMoleFracEval = 0.0;
428 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
429 sumMoleFracEval += fluidStateEval.moleFraction(phaseIdx, compIdx);
431 if (1.0 - sumMoleFracEval > fluidStateEval.saturation(phaseIdx)) {
432 b[eqIdx] = fluidState.saturation(phaseIdx);
435 Evaluation sumMoleFrac = 0.0;
436 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
437 sumMoleFrac += fluidState.moleFraction(phaseIdx, compIdx);
438 b[eqIdx] = 1.0 - sumMoleFrac;
445 template <
class MaterialLaw,
class Flu
idState,
class Vector>
447 ParameterCache ¶mCache,
448 const typename MaterialLaw::Params &matParams,
449 const Vector &deltaX)
451 typedef typename FluidState::Scalar Evaluation;
456 assert(deltaX.dimension == numEq);
457 for (
unsigned i = 0; i < numEq; ++i)
458 assert(std::isfinite(Toolbox::value(deltaX[i])));
462 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
463 const Evaluation& tmp =
getQuantity_(fluidState, pvIdx);
464 Evaluation delta = deltaX[pvIdx];
466 relError = std::max<Scalar>(relError,
518 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
523 template <
class MaterialLaw,
class Flu
idState>
525 ParameterCache ¶mCache,
526 const typename MaterialLaw::Params &matParams)
528 typedef typename FluidState::Scalar Evaluation;
532 Evaluation sumSat = 0.0;
533 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
534 sumSat += fluidState.saturation(phaseIdx);
535 fluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
540 Dune::FieldVector<Scalar, numPhases> pC;
541 MaterialLaw::capillaryPressures(pC, matParams, fluidState);
542 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
543 fluidState.setPressure(phaseIdx,
544 fluidState.pressure(0)
545 + (pC[phaseIdx] - pC[0]));
548 paramCache.updateAll(fluidState, ParameterCache::Temperature);
551 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
552 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
553 fluidState.setDensity(phaseIdx, rho);
555 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
556 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
557 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
563 {
return pvIdx == 0; }
566 {
return 1 <= pvIdx && pvIdx < numPhases; }
569 {
return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
572 template <
class Flu
idState>
573 static const typename FluidState::Scalar&
getQuantity_(
const FluidState &fluidState,
unsigned pvIdx)
575 assert(pvIdx < numEq);
579 unsigned phaseIdx = 0;
580 return fluidState.pressure(phaseIdx);
583 else if (pvIdx < numPhases) {
584 unsigned phaseIdx = pvIdx - 1;
585 return fluidState.saturation(phaseIdx);
590 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
591 unsigned compIdx = (pvIdx - numPhases)%numComponents;
592 return fluidState.moleFraction(phaseIdx, compIdx);
597 template <
class MaterialLaw,
class Flu
idState>
599 ParameterCache ¶mCache,
600 const typename MaterialLaw::Params &matParams,
602 const typename FluidState::Scalar& value)
604 typedef typename FluidState::Scalar Evaluation;
606 assert(0 <= pvIdx && pvIdx < numEq);
609 Evaluation delta = value - fluidState.pressure(0);
613 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
614 fluidState.setPressure(phaseIdx, fluidState.pressure(phaseIdx) + delta);
615 paramCache.updateAllPressures(fluidState);
618 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
619 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
621 fluidState.setDensity(phaseIdx, rho);
623 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
624 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
626 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
630 else if (pvIdx < numPhases) {
631 const Evaluation& delta = value - fluidState.saturation(pvIdx - 1);
633 fluidState.setSaturation(pvIdx - 1, value);
637 fluidState.setSaturation(numPhases - 1,
638 fluidState.saturation(numPhases - 1) - delta);
641 Dune::FieldVector<Evaluation, numPhases> pC;
642 MaterialLaw::capillaryPressures(pC, matParams, fluidState);
644 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
645 fluidState.setPressure(phaseIdx,
646 fluidState.pressure(0)
647 + (pC[phaseIdx] - pC[0]));
648 paramCache.updateAllPressures(fluidState);
651 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
652 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
654 fluidState.setDensity(phaseIdx, rho);
656 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
657 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
659 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
663 else if (pvIdx < numPhases + numPhases*numComponents)
665 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
666 unsigned compIdx = (pvIdx - numPhases)%numComponents;
669 fluidState.setMoleFraction(phaseIdx, compIdx, value);
670 paramCache.updateSingleMoleFraction(fluidState, phaseIdx, compIdx);
673 const Evaluation& rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
675 fluidState.setDensity(phaseIdx, rho);
679 if (!FluidSystem::isIdealMixture(phaseIdx)) {
680 for (
unsigned fugCompIdx = 0; fugCompIdx < numComponents; ++fugCompIdx) {
681 const Evaluation& phi = FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, fugCompIdx);
683 fluidState.setFugacityCoefficient(phaseIdx, fugCompIdx, phi);
693 template <
class Flu
idState>
696 const typename FluidState::Scalar& value)
698 assert(pvIdx < numEq);
703 unsigned phaseIdx = 0;
704 fluidState.setPressure(phaseIdx, value);
707 else if (pvIdx < numPhases) {
708 unsigned phaseIdx = pvIdx - 1;
709 fluidState.setSaturation(phaseIdx, value);
714 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
715 unsigned compIdx = (pvIdx - numPhases)%numComponents;
716 fluidState.setMoleFraction(phaseIdx, compIdx, value);
720 template <
class Flu
idState>
727 else if (pvIdx < numPhases)
static void setQuantity_(FluidState &fluidState, ParameterCache ¶mCache, const typename MaterialLaw::Params &matParams, unsigned pvIdx, const typename FluidState::Scalar &value)
Definition: NcpFlash.hpp:598
static bool isMoleFracIdx_(unsigned pvIdx)
Definition: NcpFlash.hpp:568
void SetUndefined(const T &value OPM_UNUSED)
Make the memory on which an object resides undefined in valgrind runs.
Definition: Valgrind.hpp:138
static Scalar update_(FluidState &fluidState, ParameterCache ¶mCache, const typename MaterialLaw::Params &matParams, const Vector &deltaX)
Definition: NcpFlash.hpp:446
static void calculateDefect_(Vector &b, const FluidState &fluidStateEval, const FluidState &fluidState, const ComponentVector &globalMolarities)
Definition: NcpFlash.hpp:388
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
static bool isPressureIdx_(unsigned pvIdx)
Definition: NcpFlash.hpp:562
Definition: Air_Mesitylene.hpp:31
static void setQuantityRaw_(FluidState &fluidState, unsigned pvIdx, const typename FluidState::Scalar &value)
Definition: NcpFlash.hpp:694
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:41
Implements some common averages.
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
static void solve(FluidState &fluidState, const ComponentVector &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: NcpFlash.hpp:264
static void completeFluidState_(FluidState &fluidState, ParameterCache ¶mCache, const typename MaterialLaw::Params &matParams)
Definition: NcpFlash.hpp:524
Scalar geometricMean(Scalar x, Scalar y)
Computes the geometric average of two values.
Definition: Means.hpp:55
static void guessInitial(FluidState &fluidState, ParameterCache ¶mCache, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition: NcpFlash.hpp:98
static void solve(FluidState &fluidState, ParameterCache ¶mCache, const typename MaterialLaw::Params &matParams, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: NcpFlash.hpp:139
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:45
static void linearize_(Matrix &J, Vector &b, FluidState &fluidState, ParameterCache ¶mCache, const typename MaterialLaw::Params &matParams, const ComponentVector &globalMolarities)
Definition: NcpFlash.hpp:331
Evaluation< Scalar, VarSetTag, numVars > abs(const Evaluation< Scalar, VarSetTag, numVars > &)
Definition: Math.hpp:41
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Evaluation< Scalar, VarSetTag, numVars > min(const Evaluation< Scalar, VarSetTag, numVars > &x1, const Evaluation< Scalar, VarSetTag, numVars > &x2)
Definition: Math.hpp:61
This file contains helper classes for the material laws.
Determines the phase compositions, pressures and saturations given the total mass of all components...
Definition: NcpFlash.hpp:84
static void printFluidState_(const FluidState &fluidState)
Definition: NcpFlash.hpp:282
static bool isSaturationIdx_(unsigned pvIdx)
Definition: NcpFlash.hpp:565
static const FluidState::Scalar & getQuantity_(const FluidState &fluidState, unsigned pvIdx)
Definition: NcpFlash.hpp:573
static Scalar quantityWeight_(const FluidState &, unsigned pvIdx)
Definition: NcpFlash.hpp:721