27 #ifndef OPM_NCP_FLASH_HPP
28 #define OPM_NCP_FLASH_HPP
36 #include <opm/common/Valgrind.hpp>
38 #include <opm/common/ErrorMacros.hpp>
39 #include <opm/common/Exceptions.hpp>
41 #include <opm/common/utility/platform_dependent/disable_warnings.h>
42 #include <dune/common/fvector.hh>
43 #include <dune/common/fmatrix.hh>
44 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
90 template <
class Scalar,
class Flu
idSystem>
93 enum { numPhases = FluidSystem::numPhases };
94 enum { numComponents = FluidSystem::numComponents };
99 x00PvIdx = S0PvIdx + numPhases - 1
102 static const int numEq = numPhases*(numComponents + 1);
108 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
110 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
113 Evaluation sumMoles = 0;
114 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
115 sumMoles += globalMolarities[compIdx];
117 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
119 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
120 fluidState.setMoleFraction(phaseIdx,
122 globalMolarities[compIdx]/sumMoles);
125 fluidState.setPressure(phaseIdx, 1.0135e5);
128 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
132 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
133 paramCache.updateAll(fluidState);
134 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
135 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
136 const typename FluidState::Scalar phi =
137 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
138 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
149 template <
class MaterialLaw,
class Flu
idState>
150 static void solve(FluidState& fluidState,
151 const typename MaterialLaw::Params& matParams,
152 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
153 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
154 Scalar tolerance = -1.0)
156 typedef typename FluidState::Scalar InputEval;
158 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
159 typedef Dune::FieldVector<InputEval, numEq> Vector;
164 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
167 Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
170 tolerance = std::min<Scalar>(1e-3,
171 1e8*std::numeric_limits<Scalar>::epsilon());
173 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
174 flashParamCache.assignPersistentData(paramCache);
187 Valgrind::SetUndefined(J);
188 Valgrind::SetUndefined(deltaX);
189 Valgrind::SetUndefined(b);
191 FlashFluidState flashFluidState;
192 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
197 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
198 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
199 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
201 FlashDefectVector defect;
202 const unsigned nMax = 50;
203 for (
unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
205 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
206 Valgrind::CheckDefined(defect);
210 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
211 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
212 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
214 b[eqIdx] = defect[eqIdx].value();
216 Valgrind::CheckDefined(J);
217 Valgrind::CheckDefined(b);
221 try { J.solve(deltaX, b); }
222 catch (
const Dune::FMatrixError& e) {
223 throw Opm::NumericalProblem(e.what());
225 Valgrind::CheckDefined(deltaX);
228 Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
230 if (relError < tolerance) {
231 assignOutputFluidState_(flashFluidState, fluidState);
236 OPM_THROW(NumericalProblem,
237 "NcpFlash solver failed: "
238 "{c_alpha^kappa} = {" << globalMolarities <<
"}, "
239 <<
"T = " << fluidState.temperature(0));
249 template <
class Flu
idState,
class ComponentVector>
250 static void solve(FluidState& fluidState,
251 const ComponentVector& globalMolarities,
252 Scalar tolerance = 0.0)
256 typedef typename MaterialLaw::Params MaterialLawParams;
258 MaterialLawParams matParams;
259 solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
264 template <
class Flu
idState>
265 static void printFluidState_(
const FluidState& fluidState)
267 typedef typename FluidState::Scalar FsScalar;
269 std::cout <<
"saturations: ";
270 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
271 std::cout << fluidState.saturation(phaseIdx) <<
" ";
274 std::cout <<
"pressures: ";
275 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
276 std::cout << fluidState.pressure(phaseIdx) <<
" ";
279 std::cout <<
"densities: ";
280 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
281 std::cout << fluidState.density(phaseIdx) <<
" ";
284 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
285 std::cout <<
"composition " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
286 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
287 std::cout << fluidState.moleFraction(phaseIdx, compIdx) <<
" ";
292 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
293 std::cout <<
"fugacities " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
294 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
295 std::cout << fluidState.fugacity(phaseIdx, compIdx) <<
" ";
300 std::cout <<
"global component molarities: ";
301 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
303 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
304 sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
306 std::cout << sum <<
" ";
311 template <
class MaterialLaw,
class InputFlu
idState,
class FlashFlu
idState>
312 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
313 FlashFluidState& flashFluidState,
314 const typename MaterialLaw::Params& matParams,
315 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
317 typedef typename FlashFluidState::Scalar FlashEval;
324 flashFluidState.setTemperature(inputFluidState.temperature(0));
328 FlashEval Slast = 1.0;
329 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
330 FlashEval S = inputFluidState.saturation(phaseIdx);
331 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
335 flashFluidState.setSaturation(phaseIdx, S);
337 flashFluidState.setSaturation(numPhases - 1, Slast);
341 FlashEval p0 = inputFluidState.pressure(0);
342 p0.setDerivative(p0PvIdx, 1.0);
344 std::array<FlashEval, numPhases> pc;
345 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
346 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
347 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
350 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
351 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
352 FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
353 x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
354 flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
358 flashParamCache.updateAll(flashFluidState);
362 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
363 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
364 flashFluidState.setDensity(phaseIdx, rho);
366 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
367 const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
368 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
373 template <
class FlashFlu
idState,
class OutputFlu
idState>
374 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
375 OutputFluidState& outputFluidState)
377 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
380 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
381 const auto& S = flashFluidState.saturation(phaseIdx).value();
382 outputFluidState.setSaturation(phaseIdx, S);
384 const auto& p = flashFluidState.pressure(phaseIdx).value();
385 outputFluidState.setPressure(phaseIdx, p);
387 const auto& rho = flashFluidState.density(phaseIdx).value();
388 outputFluidState.setDensity(phaseIdx, rho);
392 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
393 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
394 const auto& moleFrac =
395 flashFluidState.moleFraction(phaseIdx, compIdx).value();
396 outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
398 const auto& fugCoeff =
399 flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
400 outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
405 template <
class FlashFlu
idState,
class FlashDefectVector,
class FlashComponentVector>
406 static void evalDefect_(FlashDefectVector& b,
407 const FlashFluidState& fluidState,
408 const FlashComponentVector& globalMolarities)
410 typedef typename FlashFluidState::Scalar FlashEval;
415 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
416 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
418 fluidState.fugacity(0, compIdx)
419 - fluidState.fugacity(phaseIdx, compIdx);
423 assert(eqIdx == numComponents*(numPhases - 1));
429 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
431 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
433 fluidState.saturation(phaseIdx)
434 * fluidState.molarity(phaseIdx, compIdx);
437 b[eqIdx] -= globalMolarities[compIdx];
442 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
443 FlashEval oneMinusSumMoleFrac = 1.0;
444 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
445 oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
447 if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
448 b[eqIdx] = fluidState.saturation(phaseIdx);
450 b[eqIdx] = oneMinusSumMoleFrac;
456 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
457 static Scalar update_(FlashFluidState& fluidState,
458 const typename MaterialLaw::Params& matParams,
459 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
460 const EvalVector& deltaX)
463 typedef typename FlashFluidState::Scalar FlashEval;
464 typedef typename FlashEval::ValueType InnerEval;
468 assert(deltaX.dimension == numEq);
469 for (
unsigned i = 0; i < numEq; ++i)
470 assert(std::isfinite(Opm::scalarValue(deltaX[i])));
474 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
475 FlashEval tmp = getQuantity_(fluidState, pvIdx);
476 InnerEval delta = deltaX[pvIdx];
478 relError = std::max(relError,
479 std::abs(Opm::scalarValue(delta))
480 * quantityWeight_(fluidState, pvIdx));
482 if (isSaturationIdx_(pvIdx)) {
484 delta = Opm::min(0.25, Opm::max(-0.25, delta));
486 else if (isMoleFracIdx_(pvIdx)) {
488 delta = Opm::min(0.20, Opm::max(-0.20, delta));
490 else if (isPressureIdx_(pvIdx)) {
492 delta = Opm::min(0.5*fluidState.pressure(0).value(),
493 Opm::max(-0.5*fluidState.pressure(0).value(),
498 setQuantity_(fluidState, pvIdx, tmp);
501 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
506 template <
class MaterialLaw,
class FlashFlu
idState>
507 static void completeFluidState_(FlashFluidState& flashFluidState,
508 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
509 const typename MaterialLaw::Params& matParams)
511 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
513 typedef typename FlashFluidState::Scalar FlashEval;
517 FlashEval sumSat = 0.0;
518 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
519 sumSat += flashFluidState.saturation(phaseIdx);
520 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
525 Dune::FieldVector<FlashEval, numPhases> pC;
526 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
527 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
528 flashFluidState.setPressure(phaseIdx,
529 flashFluidState.pressure(0)
530 + (pC[phaseIdx] - pC[0]));
533 paramCache.updateAll(flashFluidState, ParamCache::Temperature);
536 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
537 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
538 flashFluidState.setDensity(phaseIdx, rho);
540 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
541 const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
542 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
547 static bool isPressureIdx_(
unsigned pvIdx)
548 {
return pvIdx == 0; }
550 static bool isSaturationIdx_(
unsigned pvIdx)
551 {
return 1 <= pvIdx && pvIdx < numPhases; }
553 static bool isMoleFracIdx_(
unsigned pvIdx)
554 {
return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
557 template <
class Flu
idState>
558 static const typename FluidState::Scalar& getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
560 assert(pvIdx < numEq);
564 unsigned phaseIdx = 0;
565 return fluidState.pressure(phaseIdx);
568 else if (pvIdx < numPhases) {
569 unsigned phaseIdx = pvIdx - 1;
570 return fluidState.saturation(phaseIdx);
575 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
576 unsigned compIdx = (pvIdx - numPhases)%numComponents;
577 return fluidState.moleFraction(phaseIdx, compIdx);
582 template <
class Flu
idState>
583 static void setQuantity_(FluidState& fluidState,
585 const typename FluidState::Scalar& value)
587 assert(pvIdx < numEq);
589 Valgrind::CheckDefined(value);
592 unsigned phaseIdx = 0;
593 fluidState.setPressure(phaseIdx, value);
596 else if (pvIdx < numPhases) {
597 unsigned phaseIdx = pvIdx - 1;
598 fluidState.setSaturation(phaseIdx, value);
602 assert(pvIdx < numPhases + numPhases*numComponents);
603 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
604 unsigned compIdx = (pvIdx - numPhases)%numComponents;
605 fluidState.setMoleFraction(phaseIdx, compIdx, value);
609 template <
class Flu
idState>
610 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
616 else if (pvIdx < numPhases)
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition: NcpFlash.hpp:109
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
This file contains helper classes for the material laws.
static void solve(FluidState &fluidState, const typename MaterialLaw::Params &matParams, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &globalMolarities, Scalar tolerance=-1.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: NcpFlash.hpp:150
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:43
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:47
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
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:250
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
Determines the phase compositions, pressures and saturations given the total mass of all components...
Definition: NcpFlash.hpp:91
Representation of an evaluation of a function and its derivatives w.r.t.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: CompositionalFluidState.hpp:46