89 enum { numPhases = FluidSystem::numPhases };
90 enum { numComponents = FluidSystem::numComponents };
95 x00PvIdx = S0PvIdx + numPhases - 1
98 static const int numEq = numPhases*(numComponents + 1);
104 template <
class Flu
idState,
class Evaluation =
typename Flu
idState::Scalar>
106 const Dune::FieldVector<Evaluation, numComponents>& globalMolarities)
109 Evaluation sumMoles = 0;
110 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
111 sumMoles += globalMolarities[compIdx];
113 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
115 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
116 fluidState.setMoleFraction(phaseIdx,
118 globalMolarities[compIdx]/sumMoles);
121 fluidState.setPressure(phaseIdx, 1.0135e5);
124 fluidState.setSaturation(phaseIdx, 1.0/numPhases);
128 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
129 paramCache.updateAll(fluidState);
130 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
131 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
132 const typename FluidState::Scalar phi =
133 FluidSystem::fugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx);
134 fluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
145 template <
class MaterialLaw,
class Flu
idState>
146 static void solve(FluidState& fluidState,
147 const typename MaterialLaw::Params& matParams,
148 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
149 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& globalMolarities,
150 Scalar tolerance = -1.0)
152 typedef typename FluidState::Scalar InputEval;
154 typedef Dune::FieldMatrix<InputEval, numEq, numEq> Matrix;
155 typedef Dune::FieldVector<InputEval, numEq> Vector;
160 typedef Dune::FieldVector<FlashEval, numEq> FlashDefectVector;
164 tolerance = std::min<Scalar>(1e-3,
165 1e8*std::numeric_limits<Scalar>::epsilon());
167 typename FluidSystem::template ParameterCache<FlashEval> flashParamCache;
168 flashParamCache.assignPersistentData(paramCache);
181 Valgrind::SetUndefined(J);
182 Valgrind::SetUndefined(deltaX);
183 Valgrind::SetUndefined(b);
185 FlashFluidState flashFluidState;
186 assignFlashFluidState_<MaterialLaw>(fluidState, flashFluidState, matParams, flashParamCache);
191 Dune::FieldVector<FlashEval, numComponents> flashGlobalMolarities;
192 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
193 flashGlobalMolarities[compIdx] = globalMolarities[compIdx];
195 FlashDefectVector defect;
196 const unsigned nMax = 50;
197 for (
unsigned nIdx = 0; nIdx < nMax; ++nIdx) {
199 evalDefect_(defect, flashFluidState, flashGlobalMolarities);
200 Valgrind::CheckDefined(defect);
204 for (
unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
205 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx)
206 J[eqIdx][pvIdx] = defect[eqIdx].derivative(pvIdx);
208 b[eqIdx] = defect[eqIdx].value();
210 Valgrind::CheckDefined(J);
211 Valgrind::CheckDefined(b);
215 try { J.solve(deltaX, b); }
216 catch (
const Dune::FMatrixError& e) {
219 Valgrind::CheckDefined(deltaX);
222 Scalar relError = update_<MaterialLaw>(flashFluidState, matParams, flashParamCache, deltaX);
224 if (relError < tolerance) {
225 assignOutputFluidState_(flashFluidState, fluidState);
230 auto cast = [](
const auto d)
233 if constexpr (std::is_same_v<
decltype(d),
const quad>)
234 return static_cast<double>(d);
240 std::string msg =
"NcpFlash solver failed: "
241 "{c_alpha^kappa} = {";
242 for (
const auto& v : globalMolarities)
243 msg +=
" " + std::to_string(cast(getValue(v)));
245 msg += std::to_string(cast(getValue(fluidState.temperature(0))));
256 template <
class Flu
idState,
class ComponentVector>
257 static void solve(FluidState& fluidState,
258 const ComponentVector& globalMolarities,
259 Scalar tolerance = 0.0)
263 typedef typename MaterialLaw::Params MaterialLawParams;
265 MaterialLawParams matParams;
266 solve<MaterialLaw>(fluidState, matParams, globalMolarities, tolerance);
271 template <
class MaterialLaw,
class InputFlu
idState,
class FlashFlu
idState>
272 static void assignFlashFluidState_(
const InputFluidState& inputFluidState,
273 FlashFluidState& flashFluidState,
274 const typename MaterialLaw::Params& matParams,
275 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& flashParamCache)
277 typedef typename FlashFluidState::Scalar FlashEval;
284 flashFluidState.setTemperature(inputFluidState.temperature(0));
288 FlashEval Slast = 1.0;
289 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx) {
290 FlashEval S = inputFluidState.saturation(phaseIdx);
291 S.setDerivative(S0PvIdx + phaseIdx, 1.0);
295 flashFluidState.setSaturation(phaseIdx, S);
297 flashFluidState.setSaturation(numPhases - 1, Slast);
301 FlashEval p0 = inputFluidState.pressure(0);
302 p0.setDerivative(p0PvIdx, 1.0);
304 std::array<FlashEval, numPhases> pc;
305 MaterialLaw::capillaryPressures(pc, matParams, flashFluidState);
306 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
307 flashFluidState.setPressure(phaseIdx, p0 + (pc[phaseIdx] - pc[0]));
310 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
311 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
312 FlashEval x = inputFluidState.moleFraction(phaseIdx, compIdx);
313 x.setDerivative(x00PvIdx + phaseIdx*numComponents + compIdx, 1.0);
314 flashFluidState.setMoleFraction(phaseIdx, compIdx, x);
318 flashParamCache.updateAll(flashFluidState);
322 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
323 const FlashEval& rho = FluidSystem::density(flashFluidState, flashParamCache, phaseIdx);
324 flashFluidState.setDensity(phaseIdx, rho);
326 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
327 const FlashEval& fugCoeff = FluidSystem::fugacityCoefficient(flashFluidState, flashParamCache, phaseIdx, compIdx);
328 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
333 template <
class FlashFlu
idState,
class OutputFlu
idState>
334 static void assignOutputFluidState_(
const FlashFluidState& flashFluidState,
335 OutputFluidState& outputFluidState)
337 outputFluidState.setTemperature(flashFluidState.temperature(0).value());
340 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
341 const auto& S = flashFluidState.saturation(phaseIdx).value();
342 outputFluidState.setSaturation(phaseIdx, S);
344 const auto& p = flashFluidState.pressure(phaseIdx).value();
345 outputFluidState.setPressure(phaseIdx, p);
347 const auto& rho = flashFluidState.density(phaseIdx).value();
348 outputFluidState.setDensity(phaseIdx, rho);
352 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
353 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
354 const auto& moleFrac =
355 flashFluidState.moleFraction(phaseIdx, compIdx).value();
356 outputFluidState.setMoleFraction(phaseIdx, compIdx, moleFrac);
358 const auto& fugCoeff =
359 flashFluidState.fugacityCoefficient(phaseIdx, compIdx).value();
360 outputFluidState.setFugacityCoefficient(phaseIdx, compIdx, fugCoeff);
365 template <
class FlashFlu
idState,
class FlashDefectVector,
class FlashComponentVector>
366 static void evalDefect_(FlashDefectVector& b,
367 const FlashFluidState& fluidState,
368 const FlashComponentVector& globalMolarities)
370 typedef typename FlashFluidState::Scalar FlashEval;
375 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
376 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx) {
378 fluidState.fugacity(0, compIdx)
379 - fluidState.fugacity(phaseIdx, compIdx);
383 assert(eqIdx == numComponents*(numPhases - 1));
389 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
391 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
393 fluidState.saturation(phaseIdx)
394 * fluidState.molarity(phaseIdx, compIdx);
397 b[eqIdx] -= globalMolarities[compIdx];
402 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
403 FlashEval oneMinusSumMoleFrac = 1.0;
404 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
405 oneMinusSumMoleFrac -= fluidState.moleFraction(phaseIdx, compIdx);
407 if (oneMinusSumMoleFrac > fluidState.saturation(phaseIdx))
408 b[eqIdx] = fluidState.saturation(phaseIdx);
410 b[eqIdx] = oneMinusSumMoleFrac;
416 template <
class MaterialLaw,
class FlashFlu
idState,
class EvalVector>
417 static Scalar update_(FlashFluidState& fluidState,
418 const typename MaterialLaw::Params& matParams,
419 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
420 const EvalVector& deltaX)
423 typedef typename FlashFluidState::Scalar FlashEval;
424 typedef typename FlashEval::ValueType InnerEval;
428 assert(deltaX.dimension == numEq);
429 for (
unsigned i = 0; i < numEq; ++i)
430 assert(std::isfinite(scalarValue(deltaX[i])));
434 for (
unsigned pvIdx = 0; pvIdx < numEq; ++ pvIdx) {
435 FlashEval tmp = getQuantity_(fluidState, pvIdx);
436 InnerEval delta = deltaX[pvIdx];
438 relError = std::max(relError,
439 std::abs(scalarValue(delta))
440 * quantityWeight_(fluidState, pvIdx));
442 if (isSaturationIdx_(pvIdx)) {
444 delta = min(0.25, max(-0.25, delta));
446 else if (isMoleFracIdx_(pvIdx)) {
448 delta = min(0.20, max(-0.20, delta));
450 else if (isPressureIdx_(pvIdx)) {
452 delta = min(0.5*fluidState.pressure(0).value(),
453 max(-0.5*fluidState.pressure(0).value(),
458 setQuantity_(fluidState, pvIdx, tmp);
461 completeFluidState_<MaterialLaw>(fluidState, paramCache, matParams);
466 template <
class MaterialLaw,
class FlashFlu
idState>
467 static void completeFluidState_(FlashFluidState& flashFluidState,
468 typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>& paramCache,
469 const typename MaterialLaw::Params& matParams)
471 typedef typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar> ParamCache;
473 typedef typename FlashFluidState::Scalar FlashEval;
477 FlashEval sumSat = 0.0;
478 for (
unsigned phaseIdx = 0; phaseIdx < numPhases - 1; ++phaseIdx)
479 sumSat += flashFluidState.saturation(phaseIdx);
480 flashFluidState.setSaturation(numPhases - 1, 1.0 - sumSat);
485 std::array<FlashEval, numPhases> pC;
486 MaterialLaw::capillaryPressures(pC, matParams, flashFluidState);
487 for (
unsigned phaseIdx = 1; phaseIdx < numPhases; ++phaseIdx)
488 flashFluidState.setPressure(phaseIdx,
489 flashFluidState.pressure(0)
490 + (pC[phaseIdx] - pC[0]));
493 paramCache.updateAll(flashFluidState, ParamCache::Temperature);
496 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
497 const FlashEval& rho = FluidSystem::density(flashFluidState, paramCache, phaseIdx);
498 flashFluidState.setDensity(phaseIdx, rho);
500 for (
unsigned compIdx = 0; compIdx < numComponents; ++ compIdx) {
501 const FlashEval& phi = FluidSystem::fugacityCoefficient(flashFluidState, paramCache, phaseIdx, compIdx);
502 flashFluidState.setFugacityCoefficient(phaseIdx, compIdx, phi);
507 static bool isPressureIdx_(
unsigned pvIdx)
508 {
return pvIdx == 0; }
510 static bool isSaturationIdx_(
unsigned pvIdx)
511 {
return 1 <= pvIdx && pvIdx < numPhases; }
513 static bool isMoleFracIdx_(
unsigned pvIdx)
514 {
return numPhases <= pvIdx && pvIdx < numPhases + numPhases*numComponents; }
517 template <
class Flu
idState>
518 static const typename FluidState::Scalar& getQuantity_(
const FluidState& fluidState,
unsigned pvIdx)
520 assert(pvIdx < numEq);
524 unsigned phaseIdx = 0;
525 return fluidState.pressure(phaseIdx);
528 else if (pvIdx < numPhases) {
529 unsigned phaseIdx = pvIdx - 1;
530 return fluidState.saturation(phaseIdx);
535 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
536 unsigned compIdx = (pvIdx - numPhases)%numComponents;
537 return fluidState.moleFraction(phaseIdx, compIdx);
542 template <
class Flu
idState>
543 static void setQuantity_(FluidState& fluidState,
545 const typename FluidState::Scalar& value)
547 assert(pvIdx < numEq);
549 Valgrind::CheckDefined(value);
552 unsigned phaseIdx = 0;
553 fluidState.setPressure(phaseIdx, value);
556 else if (pvIdx < numPhases) {
557 unsigned phaseIdx = pvIdx - 1;
558 fluidState.setSaturation(phaseIdx, value);
562 assert(pvIdx < numPhases + numPhases*numComponents);
563 unsigned phaseIdx = (pvIdx - numPhases)/numComponents;
564 unsigned compIdx = (pvIdx - numPhases)%numComponents;
565 fluidState.setMoleFraction(phaseIdx, compIdx, value);
569 template <
class Flu
idState>
570 static Scalar quantityWeight_(
const FluidState& ,
unsigned pvIdx)
576 else if (pvIdx < numPhases)