27#ifndef OPM_IMMISCIBLE_FLASH_HPP
28#define OPM_IMMISCIBLE_FLASH_HPP
37#include <dune/common/fvector.hh>
38#include <dune/common/fmatrix.hh>
39#include <dune/common/version.hh>
72template <
class Scalar,
class Flu
idSystem>
75 static const int numPhases = FluidSystem::numPhases;
76 static const int numComponents = FluidSystem::numComponents;
77 static_assert(numPhases == numComponents,
78 "Immiscibility assumes that the number of phases"
79 " is equal to the number of components");
86 static const int numEq = numPhases;
92 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
94 const Dune::FieldVector<Evaluation, numComponents>& )
96 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
98 fluidState.setPressure(phaseIdx, 1e5);
101 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
111 template <
class MaterialLaw,
class Flu
idState>
112 static void solve(FluidState& fluidState,
113 const typename MaterialLaw::Params& matParams,
114 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
115 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
116 Scalar tolerance = -1)
118 typedef typename FluidState::Scalar InputEval;
123 bool allIncompressible =
true;
124 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
125 if (FluidSystem::isCompressible(phaseIdx)) {
126 allIncompressible =
false;
131 if (allIncompressible) {
135 paramCache.updateAll(fluidState);
140 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
141 typedef Dune::FieldVector<InputEval, numEq> Vector;
144 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
147#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
148 Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
152 tolerance = std::min<Scalar>(1e-5,
153 1e8*std::numeric_limits<Scalar>::epsilon());
155 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
156 flashParamCache.assignPersistentData(paramCache);
173 FlashFluidState flashFluidState;
174 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
179 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
180 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
181 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
183 FlashDefectVector defect;
184 const unsigned nMax = 50;
185 for (
unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
187 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
192 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
193 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
194 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
196 b[eqIdx] = defect[eqIdx].value();
204 try { J.solve(deltaX, b); }
205 catch (
const Dune::FMatrixError& e) {
211 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
213 if (relError < tolerance) {
219 std::ostringstream oss;
220 oss <<
"ImmiscibleFlash solver failed:"
221 <<
" {c_alpha^kappa} = {" << globalMolarities <<
"},"
222 <<
" T = " << fluidState.temperature(0);
228 template <
class Flu
idState>
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 InputFlu
idState,
class FlashFlu
idState>
259 FlashFluidState& flashFluidState,
260 const typename MaterialLaw::Params& matParams,
261 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
263 typedef typename FlashFluidState::Scalar FlashEval;
270 flashFluidState.setTemperature(inputFluidState.temperature(0));
274 FlashEval Slast = 1.0;
275 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
276 FlashEval S = inputFluidState.saturation(phaseIdx);
277 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
281 flashFluidState.setSaturation(phaseIdx, S);
283 flashFluidState.setSaturation(numPhases - 1, Slast);
287 FlashEval p0 = inputFluidState.pressure(0);
288 p0.setDerivative(p0PvIdx, 1.0);
290 std::array<FlashEval, numPhases> pc;
291 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
292 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
293 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
295 flashParamCache.updateAll(flashFluidState);
298 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
299 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
300 flashFluidState.setDensity(phaseIdx, rho);
304 template <
class FlashFlu
idState,
class OutputFlu
idState>
306 OutputFluidState& outputFluidState)
308 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
311 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
312 const auto& S = flashFluidState.saturation(phaseIdx).value();
313 outputFluidState.setSaturation(phaseIdx, S);
315 const auto& p = flashFluidState.pressure(phaseIdx).value();
316 outputFluidState.setPressure(phaseIdx, p);
318 const auto& rho = flashFluidState.density(phaseIdx).value();
319 outputFluidState.setDensity(phaseIdx, rho);
323 template <
class Flu
idState>
325 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
326 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
328 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
329 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
330 fluidState.setDensity(phaseIdx, rho);
333 globalMolarities[phaseIdx]
334 / fluidState.molarDensity(phaseIdx);
335 fluidState.setSaturation(phaseIdx, saturation);
339 template <
class Flu
idState,
class FlashDefectVector,
class FlashComponentVector>
341 const FluidState& fluidState,
342 const FlashComponentVector& globalMolarities)
345 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
347 fluidState.saturation(phaseIdx)
348 * fluidState.molarity(phaseIdx, phaseIdx);
349 b[phaseIdx] -= globalMolarities[phaseIdx];
353 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
354 static Scalar
update_(FlashFluidState& fluidState,
355 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
356 const typename MaterialLaw::Params& matParams,
357 const EvalVector& deltaX)
360 typedef typename FlashFluidState::Scalar FlashEval;
363 typedef typename FlashEvalToolbox::ValueType InnerEval;
367 assert(deltaX.dimension == numEq);
368 for (
unsigned i = 0; i < numEq; ++i)
373 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
375 InnerEval delta = deltaX[pvIdx];
384 delta =
min(0.25,
max(-0.25, delta));
389 delta =
min(0.5*fluidState.pressure(0).value(),
390 max(-0.50*fluidState.pressure(0).value(), delta));
397 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
402 template <
class MaterialLaw,
class FlashFlu
idState>
404 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
405 const typename MaterialLaw::Params& matParams)
407 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
409 typedef typename FlashFluidState::Scalar FlashEval;
413 FlashEval sumSat = 0.0;
414 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
415 sumSat += flashFluidState.saturation(phaseIdx);
416 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
421 std::array<FlashEval, numPhases> pC;
422 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
423 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
424 flashFluidState.setPressure(phaseIdx,
425 flashFluidState.pressure(0)
426 + (pC[phaseIdx] - pC[0]));
429 paramCache.updateAll(flashFluidState, ParamCache::Temperature|ParamCache::Composition);
432 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
433 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
434 flashFluidState.setDensity(phaseIdx, rho);
439 {
return pvIdx == 0; }
442 {
return 1 <= pvIdx && pvIdx < numPhases; }
445 template <
class Flu
idState>
446 static const typename FluidState::Scalar&
getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
448 assert(pvIdx < numEq);
452 unsigned phaseIdx = 0;
453 return fluidState.pressure(phaseIdx);
457 assert(pvIdx < numPhases);
458 unsigned phaseIdx = pvIdx - 1;
459 return fluidState.saturation(phaseIdx);
464 template <
class Flu
idState>
467 const typename FluidState::Scalar& value)
469 assert(pvIdx < numEq);
473 unsigned phaseIdx = 0;
474 fluidState.setPressure(phaseIdx, value);
478 assert(pvIdx < numPhases);
479 unsigned phaseIdx = pvIdx - 1;
480 fluidState.setSaturation(phaseIdx, value);
484 template <
class Flu
idState>
492 assert(pvIdx < numPhases);
Representation of an evaluation of a function and its derivatives w.r.t. a set of variables in the lo...
Provides the opm-material specific exception classes.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Some templates to wrap the valgrind client request macros.
Represents a function evaluation and its derivatives w.r.t. a fixed set of variables.
Definition: Evaluation.hpp:59
Determines the pressures and saturations of all fluid phases given the total mass of all components.
Definition: ImmiscibleFlash.hpp:74
static bool isPressureIdx_(unsigned pvIdx)
Definition: ImmiscibleFlash.hpp:438
static void completeFluidState_(FlashFluidState &flashFluidState, typename FluidSystem::template ParameterCache< typename FlashFluidState::Scalar > ¶mCache, const typename MaterialLaw::Params &matParams)
Definition: ImmiscibleFlash.hpp:403
static Scalar quantityWeight_(const FluidState &, unsigned pvIdx)
Definition: ImmiscibleFlash.hpp:485
static void setQuantity_(FluidState &fluidState, unsigned pvIdx, const typename FluidState::Scalar &value)
Definition: ImmiscibleFlash.hpp:465
static void solveAllIncompressible_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &globalMolarities)
Definition: ImmiscibleFlash.hpp:324
static void printFluidState_(const FluidState &fs)
Definition: ImmiscibleFlash.hpp:229
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &)
Guess initial values for all quantities.
Definition: ImmiscibleFlash.hpp:93
static const FluidState::Scalar & getQuantity_(const FluidState &fluidState, unsigned pvIdx)
Definition: ImmiscibleFlash.hpp:446
static bool isSaturationIdx_(unsigned pvIdx)
Definition: ImmiscibleFlash.hpp:441
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)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: ImmiscibleFlash.hpp:112
static Scalar update_(FlashFluidState &fluidState, typename FluidSystem::template ParameterCache< typename FlashFluidState::Scalar > ¶mCache, const typename MaterialLaw::Params &matParams, const EvalVector &deltaX)
Definition: ImmiscibleFlash.hpp:354
static void evalDefect_(FlashDefectVector &b, const FluidState &fluidState, const FlashComponentVector &globalMolarities)
Definition: ImmiscibleFlash.hpp:340
static void assignOutputFluidState_(const FlashFluidState &flashFluidState, OutputFluidState &outputFluidState)
Definition: ImmiscibleFlash.hpp:305
static void assignFlashFluidState_(const InputFluidState &inputFluidState, FlashFluidState &flashFluidState, const typename MaterialLaw::Params &matParams, typename FluidSystem::template ParameterCache< typename FlashFluidState::Scalar > &flashParamCache)
Definition: ImmiscibleFlash.hpp:258
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: ImmiscibleFluidState.hpp:49
Definition: Exceptions.hpp:46
bool CheckDefined(const T &value)
Make valgrind complain if any of the memory occupied by an object is undefined.
Definition: Valgrind.hpp:74
void SetUndefined(const T &value)
Make the memory on which an object resides undefined in valgrind runs.
Definition: Valgrind.hpp:171
Definition: Air_Mesitylene.hpp:34
bool isfinite(const Evaluation &value)
Definition: MathToolbox.hpp:420
ReturnEval_< Evaluation1, Evaluation2 >::type min(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:346
ReturnEval_< Evaluation1, Evaluation2 >::type max(const Evaluation1 &arg1, const Evaluation2 &arg2)
Definition: MathToolbox.hpp:341
auto scalarValue(const Evaluation &val) -> decltype(MathToolbox< Evaluation >::scalarValue(val))
Definition: MathToolbox.hpp:335
Evaluation abs(const Evaluation &value)
Definition: MathToolbox.hpp:350