64 enum { numPhases = FluidSystem::numPhases };
65 enum { numComponents = FluidSystem::numComponents };
66 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
67 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx};
68 enum { numMiscibleComponents = FluidSystem::numMiscibleComponents};
69 enum { numMisciblePhases = FluidSystem::numMisciblePhases};
73 numMisciblePhases*numMiscibleComponents
81 template <
class Flu
idState>
82 static void solve(FluidState& fluid_state,
83 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
84 const std::string& twoPhaseMethod,
89 using InputEval =
typename FluidState::Scalar;
90 using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
94 for(
int compIdx = 0; compIdx < numComponents; ++compIdx) {
95 K[compIdx] = fluid_state.K(compIdx);
102 if (verbosity >= 1) {
103 std::cout <<
"********" << std::endl;
104 std::cout <<
"Inputs are K = [" << K <<
"], L = [" << L <<
"], z = [" << z <<
"], P = " << fluid_state.pressure(0) <<
", and T = " << fluid_state.temperature(0) << std::endl;
107 using ScalarVector = Dune::FieldVector<Scalar, numComponents>;
108 Scalar L_scalar = Opm::getValue(L);
109 ScalarVector z_scalar;
110 ScalarVector K_scalar;
111 for (
unsigned i = 0; i < numComponents; ++i) {
112 z_scalar[i] = Opm::getValue(z[i]);
113 K_scalar[i] = Opm::getValue(K[i]);
116 ScalarFluidState fluid_state_scalar;
118 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
119 fluid_state_scalar.setMoleFraction(oilPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx)));
120 fluid_state_scalar.setMoleFraction(gasPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx)));
121 fluid_state_scalar.setKvalue(compIdx, Opm::getValue(fluid_state.K(compIdx)));
124 fluid_state_scalar.setLvalue(L_scalar);
126 fluid_state_scalar.setPressure(FluidSystem::oilPhaseIdx,
127 Opm::getValue(fluid_state.pressure(FluidSystem::oilPhaseIdx)));
128 fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
129 Opm::getValue(fluid_state.pressure(FluidSystem::gasPhaseIdx)));
131 fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
134 bool is_stable =
false;
135 if ( L <= 0 || L == 1 ) {
136 if (verbosity >= 1) {
137 std::cout <<
"Perform stability test (L <= 0 or L == 1)!" << std::endl;
139 phaseStabilityTest_(is_stable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
141 if (verbosity >= 1) {
142 std::cout <<
"Inputs after stability test are K = [" << K_scalar <<
"], L = [" << L_scalar <<
"], z = [" << z_scalar <<
"], P = " << fluid_state.pressure(0) <<
", and T = " << fluid_state.temperature(0) << std::endl;
146 const bool is_single_phase = is_stable;
149 if ( !is_single_phase ) {
151 L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
152 flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
155 L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity);
157 fluid_state_scalar.setLvalue(L_scalar);
160 if (verbosity >= 1) {
161 std::cout <<
"********" << std::endl;
166 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
167 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
168 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x_i);
169 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, compIdx);
170 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y_i);
173 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
174 fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
175 fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
177 fluid_state.setLvalue(L_scalar);
179 updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
189 template <
class Flu
idState,
class ComponentVector>
190 static void solve(FluidState& fluid_state,
191 const ComponentVector& globalMolarities,
192 Scalar tolerance = 0.0)
196 using MaterialLawParams =
typename MaterialLaw::Params;
198 MaterialLawParams matParams;
199 solve<MaterialLaw>(fluid_state, matParams, globalMolarities, tolerance);
202 template <
class Vector>
203 static typename Vector::field_type solveRachfordRice_g_(
const Vector& K,
const Vector& z,
int verbosity)
207 using field_type =
typename Vector::field_type;
208 constexpr field_type tol = 1e-12;
209 constexpr int itmax = 10000;
210 field_type Kmin = K[0];
211 field_type Kmax = K[0];
212 for (
int compIdx = 1; compIdx < numComponents; ++compIdx){
213 if (K[compIdx] < Kmin)
215 else if (K[compIdx] >= Kmax)
219 auto Vmin = 1 / (1 - Kmax);
220 auto Vmax = 1 / (1 - Kmin);
222 auto V = (Vmin + Vmax)/2;
224 if (verbosity == 3 || verbosity == 4) {
225 std::cout <<
"Initial guess " << numComponents <<
"c : V = " << V <<
" and [Vmin, Vmax] = [" << Vmin <<
", " << Vmax <<
"]" << std::endl;
226 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"abs(step)" << std::setw(16) <<
"V" << std::endl;
229 for (
int iteration = 1; iteration < itmax; ++iteration) {
231 field_type denum = 0.0;
233 for (
int compIdx = 0; compIdx < numComponents; ++compIdx){
234 auto dK = K[compIdx] - 1.0;
235 auto a = z[compIdx] * dK;
236 auto b = (1 + V * dK);
238 denum += z[compIdx] * (dK*dK) / (b*b);
240 auto delta = r / denum;
244 if (V < Vmin || V > Vmax)
247 if (verbosity == 3 || verbosity == 4) {
248 std::cout <<
"V = " << V <<
" is not within the the range [Vmin, Vmax], solve using Bisection method!" << std::endl;
256 decltype(Vmax) Lmin = 1.0;
257 decltype(Vmin) Lmax = 0.0;
258 auto L = bisection_g_(K, Lmin, Lmax, z, verbosity);
261 if (verbosity >= 1) {
262 std::cout <<
"Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
268 if (verbosity == 3 || verbosity == 4) {
269 std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << V << std::endl;
272 if ( Opm::abs(r) < tol ) {
277 if (verbosity >= 1) {
278 std::cout <<
"Rachford-Rice converged to final solution L = " << L << std::endl;
285 throw std::runtime_error(
" Rachford-Rice did not converge within maximum number of iterations" );
288 template <
class Vector>
289 static typename Vector::field_type bisection_g_(
const Vector& K,
typename Vector::field_type Lmin,
290 typename Vector::field_type Lmax,
const Vector& z,
int verbosity)
293 typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
296 if (verbosity >= 3) {
297 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"g(Lmid)" << std::setw(16) <<
"L" << std::endl;
300 constexpr int max_it = 10000;
302 auto closeLmaxLmin = [](
double max_v,
double min_v) {
303 return Opm::abs(max_v - min_v) / 2. < 1e-10;
308 if (closeLmaxLmin(Lmax, Lmin) ){
309 throw std::runtime_error(fmt::format(
"Strange bisection with Lmax {} and Lmin {}?", Lmax, Lmin));
311 for (
int iteration = 0; iteration < max_it; ++iteration){
313 auto L = (Lmin + Lmax) / 2;
314 auto gMid = rachfordRice_g_(K, L, z);
316 if (verbosity == 3 || verbosity == 4) {
317 std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl;
321 if (Opm::abs(gMid) < 1e-16 || closeLmaxLmin(Lmax, Lmin)){
325 else if (Dune::sign(gMid) != Dune::sign(gLmin)) {
335 throw std::runtime_error(
" Rachford-Rice bisection failed with " + std::to_string(max_it) +
" iterations!");
338 template <
class Vector,
class FlashFlu
idState>
339 static typename Vector::field_type li_single_phase_label_(
const FlashFluidState& fluid_state,
const Vector& z,
int verbosity)
342 typename Vector::field_type sumVz = 0.0;
343 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
345 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
348 sumVz += (V_crit * z[compIdx]);
352 typename Vector::field_type Tc_est = 0.0;
353 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
355 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
356 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
359 Tc_est += (V_crit * T_crit * z[compIdx] / sumVz);
363 const auto& T = fluid_state.temperature(0);
366 typename Vector::field_type L;
372 if (verbosity >= 1) {
373 std::cout <<
"Cell is single-phase, liquid (L = 1.0) due to Li's phase labeling method giving T < Tc_est (" << T <<
" < " << Tc_est <<
")!" << std::endl;
381 if (verbosity >= 1) {
382 std::cout <<
"Cell is single-phase, vapor (L = 0.0) due to Li's phase labeling method giving T >= Tc_est (" << T <<
" >= " << Tc_est <<
")!" << std::endl;
390 template <
class FlashFlu
idState,
class ComponentVector>
391 static void phaseStabilityTest_(
bool& isStable, ComponentVector& K, FlashFluidState& fluid_state,
const ComponentVector& z,
int verbosity)
394 bool isTrivialL, isTrivialV;
395 ComponentVector x, y;
396 typename FlashFluidState::Scalar S_l, S_v;
397 ComponentVector K0 = K;
398 ComponentVector K1 = K;
401 if (verbosity == 3 || verbosity == 4) {
402 std::cout <<
"Stability test for vapor phase:" << std::endl;
404 checkStability_(fluid_state, isTrivialV, K0, y, S_v, z,
true, verbosity);
405 bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
408 if (verbosity == 3 || verbosity == 4) {
409 std::cout <<
"Stability test for liquid phase:" << std::endl;
411 checkStability_(fluid_state, isTrivialL, K1, x, S_l, z,
false, verbosity);
412 bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
415 isStable = L_stable && V_unstable;
419 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
420 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
421 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
426 for (
int compIdx = 0; compIdx<numComponents; ++compIdx) {
427 K[compIdx] = y[compIdx] / x[compIdx];
434 template <
class FlashFlu
idState>
435 static typename FlashFluidState::Scalar wilsonK_(
const FlashFluidState& fluid_state,
int compIdx)
437 const auto& acf = FluidSystem::acentricFactor(compIdx);
438 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
439 const auto& T = fluid_state.temperature(0);
440 const auto& p_crit = FluidSystem::criticalPressure(compIdx);
441 const auto& p = fluid_state.pressure(0);
443 const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
447 template <
class Vector>
448 static typename Vector::field_type rachfordRice_g_(
const Vector& K,
typename Vector::field_type L,
const Vector& z)
450 typename Vector::field_type g=0;
451 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
452 g += (z[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
457 template <
class Vector>
458 static typename Vector::field_type rachfordRice_dg_dL_(
const Vector& K,
const typename Vector::field_type L,
const Vector& z)
460 typename Vector::field_type dg=0;
461 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
462 dg += (z[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
467 template <
class FlashFlu
idState,
class ComponentVector>
468 static void checkStability_(
const FlashFluidState& fluid_state,
bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
469 typename FlashFluidState::Scalar& S_loc,
const ComponentVector& z,
bool isGas,
int verbosity)
471 using FlashEval =
typename FlashFluidState::Scalar;
475 FlashFluidState fluid_state_fake = fluid_state;
476 FlashFluidState fluid_state_global = fluid_state;
479 if (verbosity >= 3) {
480 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"K-Norm" << std::setw(16) <<
"R-Norm" << std::endl;
485 for (
int i = 0; i < 20000; ++i) {
488 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
489 xy_loc[compIdx] = K[compIdx] * z[compIdx];
490 S_loc += xy_loc[compIdx];
492 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
493 xy_loc[compIdx] /= S_loc;
494 fluid_state_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
498 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
499 xy_loc[compIdx] = z[compIdx]/K[compIdx];
500 S_loc += xy_loc[compIdx];
502 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
503 xy_loc[compIdx] /= S_loc;
504 fluid_state_fake.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
508 int phaseIdx = (isGas ?
static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
509 int phaseIdx2 = (isGas ?
static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
511 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
512 fluid_state_global.setMoleFraction(phaseIdx2, compIdx, z[compIdx]);
515 typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
516 paramCache_fake.updatePhase(fluid_state_fake, phaseIdx);
518 typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
519 paramCache_global.updatePhase(fluid_state_global, phaseIdx2);
522 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
526 fluid_state_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
527 fluid_state_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
532 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
534 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
535 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
536 auto fug_ratio = fug_global / fug_fake;
537 R[compIdx] = fug_ratio / S_loc;
540 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
541 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
542 auto fug_ratio = fug_fake / fug_global;
543 R[compIdx] = fug_ratio * S_loc;
547 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
548 K[compIdx] *= R[compIdx];
552 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
553 auto a = Opm::getValue(R[compIdx]) - 1.0;
554 auto b = Opm::log(Opm::getValue(K[compIdx]));
560 if (verbosity >= 3) {
561 std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
565 isTrivial = (K_norm < 1e-5);
566 if (isTrivial || R_norm < 1e-10)
571 throw std::runtime_error(
" Stability test did not converge");
575 template <
class FlashFlu
idState,
class ComponentVector>
576 static void computeLiquidVapor_(FlashFluidState& fluid_state,
typename FlashFluidState::Scalar& L, ComponentVector& K,
const ComponentVector& z)
581 typename FlashFluidState::Scalar sumx=0;
582 typename FlashFluidState::Scalar sumy=0;
583 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
584 x[compIdx] = z[compIdx]/(L + (1-L)*K[compIdx]);
586 y[compIdx] = (K[compIdx]*z[compIdx])/(L + (1-L)*K[compIdx]);
592 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
593 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
594 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
598 template <
class Flu
idState,
class ComponentVector>
599 static void flash_2ph(
const ComponentVector& z_scalar,
600 const std::string& flash_2p_method,
601 ComponentVector& K_scalar,
602 typename FluidState::Scalar& L_scalar,
603 FluidState& fluid_state_scalar,
605 if (verbosity >= 1) {
606 std::cout <<
"Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar <<
"]" << std::endl;
611 bool converged =
false;
612 if (flash_2p_method ==
"newton") {
613 if (verbosity >= 1) {
614 std::cout <<
"Calculate composition using Newton." << std::endl;
616 converged = newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
617 }
else if (flash_2p_method ==
"ssi") {
619 if (verbosity >= 1) {
620 std::cout <<
"Calculate composition using Succcessive Substitution." << std::endl;
622 converged = successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
false, verbosity);
623 }
else if (flash_2p_method ==
"ssi+newton") {
624 converged = successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
true, verbosity);
626 converged = newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
629 throw std::logic_error(
"unknown two phase flash method " + flash_2p_method +
" is specified");
633 throw std::runtime_error(
"flash calculation did not get converged with " + flash_2p_method);
637 template <
class FlashFlu
idState,
class ComponentVector>
638 static bool newtonComposition_(ComponentVector& K,
typename FlashFluidState::Scalar& L,
639 FlashFluidState& fluid_state,
const ComponentVector& z,
644 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
645 constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
646 using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
647 using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
648 constexpr Scalar tolerance = 1.e-8;
650 NewtonVector soln(0.);
651 NewtonVector res(0.);
652 NewtonMatrix jac (0.);
655 computeLiquidVapor_(fluid_state, L, K, z);
658 if (verbosity >= 1) {
659 std::cout <<
" the current L is " << Opm::getValue(L) << std::endl;
660 std::cout <<
"Initial guess: x = [";
661 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
662 if (compIdx < numComponents - 1)
663 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
665 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
667 std::cout <<
"], y = [";
668 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
669 if (compIdx < numComponents - 1)
670 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
672 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
674 std::cout <<
"], and " <<
"L = " << L << std::endl;
678 if (verbosity == 2 || verbosity == 4) {
679 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"Norm2(step)" << std::setw(16) <<
"Norm2(Residual)" << std::endl;
683 using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
685 std::vector<Eval> x(numComponents), y(numComponents);
689 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
690 x[compIdx] = Eval(fluid_state.moleFraction(oilPhaseIdx, compIdx), compIdx);
691 const unsigned idx = compIdx + numComponents;
692 y[compIdx] = Eval(fluid_state.moleFraction(gasPhaseIdx, compIdx), idx);
694 l = Eval(L, num_primary_variables - 1);
697 CompositionalFluidState<Eval, FluidSystem> flash_fluid_state;
698 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
699 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
700 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
702 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
704 flash_fluid_state.setLvalue(l);
706 flash_fluid_state.setPressure(FluidSystem::oilPhaseIdx,
707 fluid_state.pressure(FluidSystem::oilPhaseIdx));
708 flash_fluid_state.setPressure(FluidSystem::gasPhaseIdx,
709 fluid_state.pressure(FluidSystem::gasPhaseIdx));
712 flash_fluid_state.setSaturation(FluidSystem::gasPhaseIdx,
713 fluid_state.saturation(FluidSystem::gasPhaseIdx));
714 flash_fluid_state.setSaturation(FluidSystem::oilPhaseIdx,
715 fluid_state.saturation(FluidSystem::oilPhaseIdx));
716 flash_fluid_state.setTemperature(fluid_state.temperature(0));
718 using ParamCache =
typename FluidSystem::template ParameterCache<typename CompositionalFluidState<Eval, FluidSystem>::Scalar>;
719 ParamCache paramCache;
721 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
722 paramCache.updatePhase(flash_fluid_state, phaseIdx);
723 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
725 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
726 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
729 bool converged =
false;
731 constexpr unsigned max_iter = 1000;
732 while (iter < max_iter) {
733 assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
734 (flash_fluid_state, z, jac, res);
735 if (verbosity >= 1) {
736 std::cout <<
" newton residual is " << res.two_norm() << std::endl;
738 converged = res.two_norm() < tolerance;
743 jac.solve(soln, res);
744 constexpr Scalar damping_factor = 1.0;
746 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
747 x[compIdx] -= soln[compIdx] * damping_factor;
748 y[compIdx] -= soln[compIdx + numComponents] * damping_factor;
750 l -= soln[num_equations - 1] * damping_factor;
752 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
753 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
754 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
756 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
758 flash_fluid_state.setLvalue(l);
760 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
761 paramCache.updatePhase(flash_fluid_state, phaseIdx);
762 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
763 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
764 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
769 if (verbosity >= 1) {
770 for (
unsigned i = 0; i < num_equations; ++i) {
771 for (
unsigned j = 0; j < num_primary_variables; ++j) {
772 std::cout <<
" " << jac[i][j];
774 std::cout << std::endl;
776 std::cout << std::endl;
779 throw std::runtime_error(
" Newton composition update did not converge within maxIterations " + std::to_string(max_iter));
783 for (
unsigned idx = 0; idx < numComponents; ++idx) {
784 const auto x_i = Opm::getValue(flash_fluid_state.moleFraction(oilPhaseIdx, idx));
785 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
786 const auto y_i = Opm::getValue(flash_fluid_state.moleFraction(gasPhaseIdx, idx));
787 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
788 const auto K_i = Opm::getValue(flash_fluid_state.K(idx));
789 fluid_state.setKvalue(idx, K_i);
793 L = Opm::getValue(l);
794 fluid_state.setLvalue(L);
799 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
800 static void assembleNewton_(
const FlashFluidState& fluid_state,
801 const ComponentVector& global_composition,
802 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
803 Dune::FieldVector<double, num_equation>& res)
805 using Eval = DenseAd::Evaluation<double, num_primary>;
806 std::vector<Eval> x(numComponents), y(numComponents);
807 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
808 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
809 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
811 const Eval& l = fluid_state.L();
815 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
818 auto local_res = -global_composition[compIdx] + l * x[compIdx] + (1 - l) * y[compIdx];
819 res[compIdx] = Opm::getValue(local_res);
820 for (
unsigned i = 0; i < num_primary; ++i) {
821 jac[compIdx][i] = local_res.derivative(i);
827 auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
828 fluid_state.fugacity(gasPhaseIdx, compIdx));
829 res[compIdx + numComponents] = Opm::getValue(local_res);
830 for (
unsigned i = 0; i < num_primary; ++i) {
831 jac[compIdx + numComponents][i] = local_res.derivative(i);
838 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
842 auto local_res = sumx - sumy;
843 res[num_equation - 1] = Opm::getValue(local_res);
844 for (
unsigned i = 0; i < num_primary; ++i) {
845 jac[num_equation - 1][i] = local_res.derivative(i);
850 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
851 static void assembleNewtonSingle_(
const FlashFluidState& fluid_state,
852 const ComponentVector& global_composition,
853 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
854 Dune::FieldVector<double, num_equation>& res)
856 using Eval = DenseAd::Evaluation<double, num_primary>;
857 std::vector<Eval> x(numComponents), y(numComponents);
858 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
859 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
860 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
862 const Eval& l = fluid_state.L();
867 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
870 auto local_res = -global_composition[compIdx] + x[compIdx];
871 res[compIdx] = Opm::getValue(local_res);
872 for (
unsigned i = 0; i < num_primary; ++i) {
873 jac[compIdx][i] = local_res.derivative(i);
879 auto local_res = -global_composition[compIdx] + y[compIdx];
880 res[compIdx + numComponents] = Opm::getValue(local_res);
881 for (
unsigned i = 0; i < num_primary; ++i) {
882 jac[compIdx + numComponents][i] = local_res.derivative(i);
888 const bool isGas = Opm::abs(l - 1.0) > std::numeric_limits<double>::epsilon();
896 res[num_equation - 1] = Opm::getValue(local_res);
897 for (
unsigned i = 0; i < num_primary; ++i) {
898 jac[num_equation - 1][i] = local_res.derivative(i);
902 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
903 static void updateDerivatives_(
const FlashFluidStateScalar& fluid_state_scalar,
904 const ComponentVector& z,
905 FluidState& fluid_state,
906 bool is_single_phase)
909 updateDerivativesTwoPhase_(fluid_state_scalar, z, fluid_state);
911 updateDerivativesSinglePhase_(fluid_state_scalar, z, fluid_state);
915 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
916 static void updateDerivativesTwoPhase_(
const FlashFluidStateScalar& fluid_state_scalar,
917 const ComponentVector& z,
918 FluidState& fluid_state)
921 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
922 constexpr size_t secondary_num_pv = numComponents + 1;
924 using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
927 SecondaryFlashFluidState secondary_fluid_state;
928 SecondaryComponentVector secondary_z;
931 const SecondaryEval sec_p = SecondaryEval(fluid_state_scalar.pressure(FluidSystem::oilPhaseIdx), 0);
932 secondary_fluid_state.setPressure(FluidSystem::oilPhaseIdx, sec_p);
933 secondary_fluid_state.setPressure(FluidSystem::gasPhaseIdx, sec_p);
936 secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));
938 for (
unsigned idx = 0; idx < numComponents; ++idx) {
939 secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
942 for (
unsigned idx = 0; idx < numComponents; ++idx) {
944 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, idx);
945 secondary_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
946 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, idx);
947 secondary_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
949 const auto K_i = fluid_state_scalar.K(idx);
950 secondary_fluid_state.setKvalue(idx, K_i);
952 const auto L = fluid_state_scalar.L();
953 secondary_fluid_state.setLvalue(L);
957 using SecondaryParamCache =
typename FluidSystem::template ParameterCache<typename SecondaryFlashFluidState::Scalar>;
958 SecondaryParamCache secondary_param_cache;
959 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
960 secondary_param_cache.updatePhase(secondary_fluid_state, phase_idx);
961 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
962 SecondaryEval phi = FluidSystem::fugacityCoefficient(secondary_fluid_state, secondary_param_cache, phase_idx, comp_idx);
963 secondary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
967 using SecondaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
968 using SecondaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv>;
969 SecondaryNewtonMatrix sec_jac;
970 SecondaryNewtonVector sec_res;
974 assembleNewton_<SecondaryFlashFluidState, SecondaryComponentVector, secondary_num_pv, num_equations>
975 (secondary_fluid_state, secondary_z, sec_jac, sec_res);
980 constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
982 using PrimaryComponentVector = Dune::FieldVector<double, numComponents>;
985 PrimaryFlashFluidState primary_fluid_state;
987 PrimaryComponentVector primary_z;
988 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
989 primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
991 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
992 const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
993 primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
994 const unsigned idx = comp_idx + numComponents;
995 const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
996 primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
997 primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
1000 l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
1001 primary_fluid_state.setLvalue(l);
1002 primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
1003 primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
1004 primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
1007 using PrimaryParamCache =
typename FluidSystem::template ParameterCache<typename PrimaryFlashFluidState::Scalar>;
1008 PrimaryParamCache primary_param_cache;
1009 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
1010 primary_param_cache.updatePhase(primary_fluid_state, phase_idx);
1011 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
1012 PrimaryEval phi = FluidSystem::fugacityCoefficient(primary_fluid_state, primary_param_cache, phase_idx, comp_idx);
1013 primary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
1017 using PrimaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
1018 using PrimaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, primary_num_pv>;
1019 PrimaryNewtonVector pri_res;
1020 PrimaryNewtonMatrix pri_jac;
1023 assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
1024 (primary_fluid_state, primary_z, pri_jac, pri_res);
1030 sec_jac.template leftmultiply(pri_jac);
1032 using InputEval =
typename FluidState::Scalar;
1033 using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
1035 ComponentVectorMoleFraction x(numComponents), y(numComponents);
1036 InputEval L_eval = L;
1043 const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
1044 const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
1046 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1047 x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);
1048 y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);
1055 constexpr size_t num_deri = numComponents;
1056 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1057 std::vector<double> deri(num_deri, 0.);
1059 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1060 deri[idx] = -sec_jac[compIdx][0] * p_l.derivative(idx);
1063 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1064 const double pz = -sec_jac[compIdx][cIdx + 1];
1065 const auto& zi = z[cIdx];
1066 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1067 deri[idx] += pz * zi.derivative(idx);
1070 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1071 x[compIdx].setDerivative(idx, deri[idx]);
1074 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1075 deri[idx] = -sec_jac[compIdx + numComponents][0] * p_v.derivative(idx);
1077 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1078 const double pz = -sec_jac[compIdx + numComponents][cIdx + 1];
1079 const auto& zi = z[cIdx];
1080 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1081 deri[idx] += pz * zi.derivative(idx);
1084 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1085 y[compIdx].setDerivative(idx, deri[idx]);
1089 std::vector<double> deriL(num_deri, 0.);
1090 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1091 deriL[idx] = -sec_jac[2 * numComponents][0] * p_v.derivative(idx);
1093 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1094 const double pz = -sec_jac[2 * numComponents][cIdx + 1];
1095 const auto& zi = z[cIdx];
1096 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1097 deriL[idx] += pz * zi.derivative(idx);
1101 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1102 L_eval.setDerivative(idx, deriL[idx]);
1107 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1108 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
1109 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
1111 fluid_state.setLvalue(L_eval);
1114 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
1115 static void updateDerivativesSinglePhase_(
const FlashFluidStateScalar& fluid_state_scalar,
1116 const ComponentVector& z,
1117 FluidState& fluid_state)
1119 using InputEval =
typename FluidState::Scalar;
1121 InputEval L_eval = fluid_state_scalar.L();;
1125 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1126 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
1127 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
1129 fluid_state.setLvalue(L_eval);
1134 template <
class FlashFlu
idState,
class ComponentVector>
1135 static bool successiveSubstitutionComposition_(ComponentVector& K,
typename ComponentVector::field_type& L, FlashFluidState& fluid_state,
const ComponentVector& z,
1136 const bool newton_afterwards,
const int verbosity)
1139 const int maxIterations = newton_afterwards ? 5 : 100;
1142 std::ios_base::fmtflags f(std::cout.flags());
1146 std::cout <<
"Initial guess: K = [" << K <<
"] and L = " << L << std::endl;
1148 if (verbosity == 2 || verbosity == 4) {
1150 int fugWidth = (numComponents * 12)/2;
1151 int convWidth = fugWidth + 7;
1152 std::cout << std::setw(10) <<
"Iteration" << std::setw(fugWidth) <<
"fL/fV" << std::setw(convWidth) <<
"norm2(fL/fv-1)" << std::endl;
1157 for (
int i=0; i < maxIterations; ++i){
1159 computeLiquidVapor_(fluid_state, L, K, z);
1162 using ParamCache =
typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
1163 ParamCache paramCache;
1164 for (
int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
1165 paramCache.updatePhase(fluid_state, phaseIdx);
1166 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1167 auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, phaseIdx, compIdx);
1168 fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
1173 ComponentVector newFugRatio;
1174 ComponentVector convFugRatio;
1175 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1176 newFugRatio[compIdx] = fluid_state.fugacity(oilPhaseIdx, compIdx)/fluid_state.fugacity(gasPhaseIdx, compIdx);
1177 convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
1181 if (verbosity >= 2) {
1183 int fugWidth = (prec + 3);
1184 int convWidth = prec + 9;
1185 std::cout << std::defaultfloat;
1186 std::cout << std::fixed;
1187 std::cout << std::setw(5) << i;
1188 std::cout << std::setw(fugWidth);
1189 std::cout << std::setprecision(prec);
1190 std::cout << newFugRatio;
1191 std::cout << std::scientific;
1192 std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
1196 if (convFugRatio.two_norm() < 1.e-6) {
1201 if (verbosity >= 1) {
1202 std::cout <<
"Solution converged to the following result :" << std::endl;
1203 std::cout <<
"x = [";
1204 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
1205 if (compIdx < numComponents - 1)
1206 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
1208 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
1210 std::cout <<
"]" << std::endl;
1211 std::cout <<
"y = [";
1212 for (
int compIdx = 0; compIdx < numComponents; ++compIdx) {
1213 if (compIdx < numComponents - 1)
1214 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
1216 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
1218 std::cout <<
"]" << std::endl;
1219 std::cout <<
"K = [" << K <<
"]" << std::endl;
1220 std::cout <<
"L = " << L << std::endl;
1228 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1229 K[compIdx] *= newFugRatio[compIdx];
1233 L = solveRachfordRice_g_(K, z, 0);
1238 const std::string msg = fmt::format(
"Successive substitution composition update did not converge within maxIterations {}.", maxIterations);
1239 if (!newton_afterwards) {
1240 throw std::runtime_error(msg);
1241 }
else if (verbosity > 0) {
1242 std::cout << msg << std::endl;