27 #ifndef OPM_IMMISCIBLE_FLASH_HPP 28 #define OPM_IMMISCIBLE_FLASH_HPP 34 #include <opm/common/Valgrind.hpp> 36 #include <opm/common/ErrorMacros.hpp> 37 #include <opm/common/Exceptions.hpp> 39 #include <dune/common/fvector.hh> 40 #include <dune/common/fmatrix.hh> 73 template <
class Scalar,
class Flu
idSystem>
76 static const int numPhases = FluidSystem::numPhases;
77 static const int numComponents = FluidSystem::numComponents;
78 static_assert(numPhases == numComponents,
79 "Immiscibility assumes that the number of phases" 80 " is equal to the number of components");
87 static const int numEq = numPhases;
93 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
95 const Dune::FieldVector<Evaluation, numComponents>& )
97 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
99 fluidState.setPressure(phaseIdx, 1e5);
102 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
112 template <
class MaterialLaw,
class Flu
idState>
113 static void solve(FluidState& fluidState,
114 const typename MaterialLaw::Params& matParams,
115 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
116 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
117 Scalar tolerance = -1)
119 typedef typename FluidState::Scalar InputEval;
124 bool allIncompressible =
true;
125 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
126 if (FluidSystem::isCompressible(phaseIdx)) {
127 allIncompressible =
false;
132 if (allIncompressible) {
136 paramCache.updateAll(fluidState);
137 solveAllIncompressible_(fluidState, paramCache, globalMolarities);
141 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
142 typedef Dune::FieldVector<InputEval, numEq> Vector;
145 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
148 Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
151 tolerance = std::min<Scalar>(1e-5,
152 1e8*std::numeric_limits<Scalar>::epsilon());
154 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
155 flashParamCache.assignPersistentData(paramCache);
168 Valgrind::SetUndefined(J);
169 Valgrind::SetUndefined(deltaX);
170 Valgrind::SetUndefined(b);
172 FlashFluidState flashFluidState;
173 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
178 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
179 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
180 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
182 FlashDefectVector defect;
183 const unsigned nMax = 50;
184 for (
unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
186 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
187 Valgrind::CheckDefined(defect);
191 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
192 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
193 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
195 b[eqIdx] = defect[eqIdx].value();
197 Valgrind::CheckDefined(J);
198 Valgrind::CheckDefined(b);
203 try { J.solve(deltaX, b); }
204 catch (Dune::FMatrixError e) {
205 throw Opm::NumericalProblem(e.what());
207 Valgrind::CheckDefined(deltaX);
210 Scalar relError = update_<MaterialLaw>(flashFluidState, flashParamCache, matParams, deltaX);
212 if (relError < tolerance) {
213 assignOutputFluidState_(flashFluidState, fluidState);
218 OPM_THROW(Opm::NumericalProblem,
219 "ImmiscibleFlash solver failed: " 220 "{c_alpha^kappa} = {" << globalMolarities <<
"}, " 221 <<
"T = " << fluidState.temperature(0));
226 template <
class Flu
idState>
227 static void printFluidState_(
const FluidState& fs)
229 std::cout <<
"saturations: ";
230 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
231 std::cout << fs.saturation(phaseIdx) <<
" ";
234 std::cout <<
"pressures: ";
235 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
236 std::cout << fs.pressure(phaseIdx) <<
" ";
239 std::cout <<
"densities: ";
240 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
241 std::cout << fs.density(phaseIdx) <<
" ";
244 std::cout <<
"global component molarities: ";
245 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
247 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
248 sum += fs.saturation(phaseIdx)*fs.molarity(phaseIdx, compIdx);
250 std::cout << sum <<
" ";
255 template <
class MaterialLaw,
class InputFlu
idState,
class FlashFlu
idState>
256 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
257 FlashFluidState& flashFluidState,
258 const typename MaterialLaw::Params& matParams,
259 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
261 typedef typename FlashFluidState::Scalar FlashEval;
268 flashFluidState.setTemperature(inputFluidState.temperature(0));
272 FlashEval Slast = 1.0;
273 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
274 FlashEval S = inputFluidState.saturation(phaseIdx);
275 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
279 flashFluidState.setSaturation(phaseIdx, S);
281 flashFluidState.setSaturation(numPhases - 1, Slast);
285 FlashEval p0 = inputFluidState.pressure(0);
286 p0.setDerivative(p0PvIdx, 1.0);
288 std::array<FlashEval, numPhases> pc;
289 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
290 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
291 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
293 flashParamCache.updateAll(flashFluidState);
296 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
297 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
298 flashFluidState.setDensity(phaseIdx, rho);
302 template <
class FlashFlu
idState,
class OutputFlu
idState>
303 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
304 OutputFluidState& outputFluidState)
306 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
309 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
310 const auto& S = flashFluidState.saturation(phaseIdx).value();
311 outputFluidState.setSaturation(phaseIdx, S);
313 const auto& p = flashFluidState.pressure(phaseIdx).value();
314 outputFluidState.setPressure(phaseIdx, p);
316 const auto& rho = flashFluidState.density(phaseIdx).value();
317 outputFluidState.setDensity(phaseIdx, rho);
321 template <
class Flu
idState>
322 static void solveAllIncompressible_(FluidState& fluidState,
323 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
324 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities)
326 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
327 Scalar rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
328 fluidState.setDensity(phaseIdx, rho);
331 globalMolarities[phaseIdx]
332 / fluidState.molarDensity(phaseIdx);
333 fluidState.setSaturation(phaseIdx, saturation);
337 template <
class Flu
idState,
class FlashDefectVector,
class FlashComponentVector>
338 static void evalDefect_(FlashDefectVector& b,
339 const FluidState& fluidState,
340 const FlashComponentVector& globalMolarities)
343 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
345 fluidState.saturation(phaseIdx)
346 * fluidState.molarity(phaseIdx, phaseIdx);
347 b[phaseIdx] -= globalMolarities[phaseIdx];
351 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
352 static Scalar update_(FlashFluidState& fluidState,
353 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
354 const typename MaterialLaw::Params& matParams,
355 const EvalVector& deltaX)
358 typedef typename FlashFluidState::Scalar FlashEval;
361 typedef typename FlashEvalToolbox::ValueType InnerEval;
365 assert(deltaX.dimension == numEq);
366 for (
unsigned i = 0; i < numEq; ++i)
367 assert(std::isfinite(Opm::scalarValue(deltaX[i])));
371 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
372 FlashEval tmp = getQuantity_(fluidState, pvIdx);
373 InnerEval delta = deltaX[pvIdx];
375 relError = std::max(relError,
376 std::abs(Opm::scalarValue(delta))
377 * quantityWeight_(fluidState, pvIdx));
379 if (isSaturationIdx_(pvIdx)) {
382 delta = Opm::min(0.25, Opm::max(-0.25, delta));
384 else if (isPressureIdx_(pvIdx)) {
387 delta = Opm::min(0.5*fluidState.pressure(0).value(),
388 Opm::max(-0.50*fluidState.pressure(0).value(), delta));
392 setQuantity_(fluidState, pvIdx, tmp);
395 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
400 template <
class MaterialLaw,
class FlashFlu
idState>
401 static void completeFluidState_(FlashFluidState& flashFluidState,
402 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
403 const typename MaterialLaw::Params& matParams)
405 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
407 typedef typename FlashFluidState::Scalar FlashEval;
411 FlashEval sumSat = 0.0;
412 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
413 sumSat += flashFluidState.saturation(phaseIdx);
414 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
419 Dune::FieldVector<FlashEval, numPhases> pC;
420 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
421 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
422 flashFluidState.setPressure(phaseIdx,
423 flashFluidState.pressure(0)
424 + (pC[phaseIdx] - pC[0]));
427 paramCache.updateAll(flashFluidState, ParamCache::Temperature|ParamCache::Composition);
430 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
431 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
432 flashFluidState.setDensity(phaseIdx, rho);
436 static bool isPressureIdx_(
unsigned pvIdx)
437 {
return pvIdx == 0; }
439 static bool isSaturationIdx_(
unsigned pvIdx)
440 {
return 1 <= pvIdx && pvIdx < numPhases; }
443 template <
class Flu
idState>
444 static const typename FluidState::Scalar& getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
446 assert(pvIdx < numEq);
450 unsigned phaseIdx = 0;
451 return fluidState.pressure(phaseIdx);
455 assert(pvIdx < numPhases);
456 unsigned phaseIdx = pvIdx - 1;
457 return fluidState.saturation(phaseIdx);
462 template <
class Flu
idState>
463 static void setQuantity_(FluidState& fluidState,
465 const typename FluidState::Scalar& value)
467 assert(pvIdx < numEq);
471 unsigned phaseIdx = 0;
472 fluidState.setPressure(phaseIdx, value);
476 assert(pvIdx < numPhases);
477 unsigned phaseIdx = pvIdx - 1;
478 fluidState.setSaturation(phaseIdx, value);
482 template <
class Flu
idState>
483 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
490 assert(pvIdx < numPhases);
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: ImmiscibleFluidState.hpp:49
Definition: Air_Mesitylene.hpp:33
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
static void guessInitial(FluidState &fluidState, const Dune::FieldVector< Evaluation, numComponents > &)
Guess initial values for all quantities.
Definition: ImmiscibleFlash.hpp:94
Determines the pressures and saturations of all fluid phases given the total mass of all components...
Definition: ImmiscibleFlash.hpp:74
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:57
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:113
Representation of an evaluation of a function and its derivatives w.r.t.