27 #ifndef OPM_NCP_FLASH_HPP
28 #define OPM_NCP_FLASH_HPP
40 #include <dune/common/fvector.hh>
41 #include <dune/common/fmatrix.hh>
42 #include <dune/common/version.hh>
88 template <
class Scalar,
class Flu
idSystem>
91 enum { numPhases = FluidSystem::numPhases };
92 enum { numComponents = FluidSystem::numComponents };
97 x00PvIdx = S0PvIdx + numPhases - 1
100 static const int numEq = numPhases*(numComponents + 1);
106 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
108 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
111 Evaluation sumMoles = 0;
112 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
113 sumMoles += globalMolarities[compIdx];
115 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
117 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
118 fluidState.setMoleFraction(phaseIdx,
120 globalMolarities[compIdx]/sumMoles);
123 fluidState.setPressure(phaseIdx, 1.0135e5);
126 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
130 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
131 paramCache.updateAll(fluidState);
132 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
133 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
134 const typename FluidState::Scalar phi =
135 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
136 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
147 template <
class MaterialLaw,
class Flu
idState>
148 static void solve(FluidState& fluidState,
149 const typename MaterialLaw::Params& matParams,
150 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
151 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
152 Scalar tolerance = -1.0)
154 typedef typename FluidState::Scalar InputEval;
156 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
157 typedef Dune::FieldVector<InputEval, numEq> Vector;
162 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
165 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
166 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) {
225 Valgrind::CheckDefined(deltaX);
228 Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
230 if (relError < tolerance) {
231 assignOutputFluidState_(flashFluidState, fluidState);
236 std::ostringstream oss;
237 oss <<
"NcpFlash solver failed:"
238 <<
" {c_alpha^kappa} = {" << globalMolarities <<
"}, "
239 <<
" T = " << fluidState.temperature(0);
250 template <
class Flu
idState,
class ComponentVector>
251 static void solve(FluidState& fluidState,
252 const ComponentVector& globalMolarities,
253 Scalar tolerance = 0.0)
257 typedef typename MaterialLaw::Params MaterialLawParams;
259 MaterialLawParams matParams;
260 solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
265 template <
class Flu
idState>
266 static void printFluidState_(
const FluidState& fluidState)
268 typedef typename FluidState::Scalar FsScalar;
270 std::cout <<
"saturations: ";
271 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
272 std::cout << fluidState.saturation(phaseIdx) <<
" ";
275 std::cout <<
"pressures: ";
276 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
277 std::cout << fluidState.pressure(phaseIdx) <<
" ";
280 std::cout <<
"densities: ";
281 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
282 std::cout << fluidState.density(phaseIdx) <<
" ";
285 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
286 std::cout <<
"composition " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
287 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
288 std::cout << fluidState.moleFraction(phaseIdx, compIdx) <<
" ";
293 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
294 std::cout <<
"fugacities " << FluidSystem::phaseName(phaseIdx) <<
"Phase: ";
295 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
296 std::cout << fluidState.fugacity(phaseIdx, compIdx) <<
" ";
301 std::cout <<
"global component molarities: ";
302 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
304 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
305 sum += fluidState.saturation(phaseIdx)*fluidState.molarity(phaseIdx, compIdx);
307 std::cout << sum <<
" ";
312 template <
class MaterialLaw,
class InputFlu
idState,
class FlashFlu
idState>
313 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
314 FlashFluidState& flashFluidState,
315 const typename MaterialLaw::Params& matParams,
316 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
318 typedef typename FlashFluidState::Scalar FlashEval;
325 flashFluidState.setTemperature(inputFluidState.temperature(0));
329 FlashEval Slast = 1.0;
330 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
331 FlashEval S = inputFluidState.saturation(phaseIdx);
332 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
336 flashFluidState.setSaturation(phaseIdx, S);
338 flashFluidState.setSaturation(numPhases - 1, Slast);
342 FlashEval p0 = inputFluidState.pressure(0);
343 p0.setDerivative(p0PvIdx, 1.0);
345 std::array<FlashEval, numPhases> pc;
346 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
347 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
348 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
351 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
352 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
353 FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
354 x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
355 flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
359 flashParamCache.updateAll(flashFluidState);
363 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
364 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
365 flashFluidState.setDensity(phaseIdx, rho);
367 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
368 const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
369 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
374 template <
class FlashFlu
idState,
class OutputFlu
idState>
375 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
376 OutputFluidState& outputFluidState)
378 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
381 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
382 const auto& S = flashFluidState.saturation(phaseIdx).value();
383 outputFluidState.setSaturation(phaseIdx, S);
385 const auto& p = flashFluidState.pressure(phaseIdx).value();
386 outputFluidState.setPressure(phaseIdx, p);
388 const auto& rho = flashFluidState.density(phaseIdx).value();
389 outputFluidState.setDensity(phaseIdx, rho);
393 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
394 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
395 const auto& moleFrac =
396 flashFluidState.moleFraction(phaseIdx, compIdx).value();
397 outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
399 const auto& fugCoeff =
400 flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
401 outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
406 template <
class FlashFlu
idState,
class FlashDefectVector,
class FlashComponentVector>
407 static void evalDefect_(FlashDefectVector& b,
408 const FlashFluidState& fluidState,
409 const FlashComponentVector& globalMolarities)
411 typedef typename FlashFluidState::Scalar FlashEval;
416 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
417 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
419 fluidState.fugacity(0, compIdx)
420 - fluidState.fugacity(phaseIdx, compIdx);
424 assert(eqIdx == numComponents*(numPhases - 1));
430 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
432 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
434 fluidState.saturation(phaseIdx)
435 * fluidState.molarity(phaseIdx, compIdx);
438 b[eqIdx] -= globalMolarities[compIdx];
443 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
444 FlashEval oneMinusSumMoleFrac = 1.0;
445 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
446 oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
448 if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
449 b[eqIdx] = fluidState.saturation(phaseIdx);
451 b[eqIdx] = oneMinusSumMoleFrac;
457 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
458 static Scalar update_(FlashFluidState& fluidState,
459 const typename MaterialLaw::Params& matParams,
460 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
461 const EvalVector& deltaX)
464 typedef typename FlashFluidState::Scalar FlashEval;
465 typedef typename FlashEval::ValueType InnerEval;
469 assert(deltaX.dimension == numEq);
470 for (
unsigned i = 0; i < numEq; ++i)
471 assert(std::isfinite(scalarValue(deltaX[i])));
475 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
476 FlashEval tmp = getQuantity_(fluidState, pvIdx);
477 InnerEval delta = deltaX[pvIdx];
479 relError = std::max(relError,
480 std::abs(scalarValue(delta))
481 * quantityWeight_(fluidState, pvIdx));
483 if (isSaturationIdx_(pvIdx)) {
485 delta = min(0.25, max(-0.25, delta));
487 else if (isMoleFracIdx_(pvIdx)) {
489 delta = min(0.20, max(-0.20, delta));
491 else if (isPressureIdx_(pvIdx)) {
493 delta = min(0.5*fluidState.pressure(0).value(),
494 max(-0.5*fluidState.pressure(0).value(),
499 setQuantity_(fluidState, pvIdx, tmp);
502 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
507 template <
class MaterialLaw,
class FlashFlu
idState>
508 static void completeFluidState_(FlashFluidState& flashFluidState,
509 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
510 const typename MaterialLaw::Params& matParams)
512 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
514 typedef typename FlashFluidState::Scalar FlashEval;
518 FlashEval sumSat = 0.0;
519 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
520 sumSat += flashFluidState.saturation(phaseIdx);
521 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
526 std::array<FlashEval, numPhases> pC;
527 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
528 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
529 flashFluidState.setPressure(phaseIdx,
530 flashFluidState.pressure(0)
531 + (pC[phaseIdx] - pC[0]));
534 paramCache.updateAll(flashFluidState, ParamCache::Temperature);
537 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
538 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
539 flashFluidState.setDensity(phaseIdx, rho);
541 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
542 const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
543 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
548 static bool isPressureIdx_(
unsigned pvIdx)
549 {
return pvIdx == 0; }
551 static bool isSaturationIdx_(
unsigned pvIdx)
552 {
return 1 <= pvIdx && pvIdx < numPhases; }
554 static bool isMoleFracIdx_(
unsigned pvIdx)
555 {
return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
558 template <
class Flu
idState>
559 static const typename FluidState::Scalar& getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
561 assert(pvIdx < numEq);
565 unsigned phaseIdx = 0;
566 return fluidState.pressure(phaseIdx);
569 else if (pvIdx < numPhases) {
570 unsigned phaseIdx = pvIdx - 1;
571 return fluidState.saturation(phaseIdx);
576 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
577 unsigned compIdx = (pvIdx - numPhases)%numComponents;
578 return fluidState.moleFraction(phaseIdx, compIdx);
583 template <
class Flu
idState>
584 static void setQuantity_(FluidState& fluidState,
586 const typename FluidState::Scalar& value)
588 assert(pvIdx < numEq);
590 Valgrind::CheckDefined(value);
593 unsigned phaseIdx = 0;
594 fluidState.setPressure(phaseIdx, value);
597 else if (pvIdx < numPhases) {
598 unsigned phaseIdx = pvIdx - 1;
599 fluidState.setSaturation(phaseIdx, value);
603 assert(pvIdx < numPhases + numPhases*numComponents);
604 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
605 unsigned compIdx = (pvIdx - numPhases)%numComponents;
606 fluidState.setMoleFraction(phaseIdx, compIdx, value);
610 template <
class Flu
idState>
611 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
617 else if (pvIdx < numPhases)
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Representation of an evaluation of a function and its derivatives w.r.t.
Provides the opm-material specific exception classes.
This file contains helper classes for the material laws.
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Some templates to wrap the valgrind client request macros.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: CompositionalFluidState.hpp:46
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:59
Determines the phase compositions, pressures and saturations given the total mass of all components.
Definition: NcpFlash.hpp:90
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &globalMolarities)
Guess initial values for all quantities.
Definition: NcpFlash.hpp:107
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:148
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:251
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:45
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:46
Definition: Exceptions.hpp:46