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>
72 template <
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);
136 solveAllIncompressible_(fluidState, paramCache, globalMolarities);
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);
169 Valgrind::SetUndefined(J);
170 Valgrind::SetUndefined(deltaX);
171 Valgrind::SetUndefined(b);
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);
188 Valgrind::CheckDefined(defect);
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();
198 Valgrind::CheckDefined(J);
199 Valgrind::CheckDefined(b);
204 try { J.solve(deltaX, b); }
205 catch (
const Dune::FMatrixError& e) {
208 Valgrind::CheckDefined(deltaX);
211 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
213 if (relError < tolerance) {
214 assignOutputFluidState_(flashFluidState, fluidState);
219 std::ostringstream oss;
220 oss <<
"ImmiscibleFlash solver failed:"
221 <<
" {c_alpha^kappa} = {" << globalMolarities <<
"},"
222 <<
" T = " << 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 InputFlu
idState,
class FlashFlu
idState>
258 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
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>
305 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
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>
324 static void solveAllIncompressible_(FluidState& fluidState,
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>
340 static void evalDefect_(FlashDefectVector& b,
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;
361 typedef MathToolbox<FlashEval> FlashEvalToolbox;
363 typedef typename FlashEvalToolbox::ValueType InnerEval;
367 assert(deltaX.dimension == numEq);
368 for (
unsigned i = 0; i < numEq; ++i)
369 assert(std::isfinite(scalarValue(deltaX[i])));
373 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
374 FlashEval tmp = getQuantity_(fluidState, pvIdx);
375 InnerEval delta = deltaX[pvIdx];
377 relError = std::max(relError,
378 std::abs(scalarValue(delta))
379 * quantityWeight_(fluidState, pvIdx));
381 if (isSaturationIdx_(pvIdx)) {
384 delta = min(0.25, max(-0.25, delta));
386 else if (isPressureIdx_(pvIdx)) {
389 delta = min(0.5*fluidState.pressure(0).value(),
390 max(-0.50*fluidState.pressure(0).value(), delta));
394 setQuantity_(fluidState, pvIdx, tmp);
397 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
402 template <
class MaterialLaw,
class FlashFlu
idState>
403 static void completeFluidState_(FlashFluidState& flashFluidState,
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);
438 static bool isPressureIdx_(
unsigned pvIdx)
439 {
return pvIdx == 0; }
441 static bool isSaturationIdx_(
unsigned pvIdx)
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>
465 static void setQuantity_(FluidState& fluidState,
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>
485 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
492 assert(pvIdx < numPhases);
Representation of an evaluation of a function and its derivatives w.r.t.
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.
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 void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &)
Guess initial values for all quantities.
Definition: ImmiscibleFlash.hpp:93
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
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: ImmiscibleFluidState.hpp:49
Definition: Exceptions.hpp:46