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 std::string twoPhaseMethod,
85 Scalar tolerance = -1.,
89 using InputEval =
typename FluidState::Scalar;
90 using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
93 tolerance = std::min<Scalar>(1e-3, 1e8 * std::numeric_limits<Scalar>::epsilon());
98 for(
int compIdx = 0; compIdx < numComponents; ++compIdx) {
99 K[compIdx] = fluid_state.K(compIdx);
106 if (verbosity >= 1) {
107 std::cout <<
"********" << std::endl;
108 std::cout <<
"Inputs are K = [" << K <<
"], L = [" << L <<
"], z = [" << z <<
"], P = " << fluid_state.pressure(0) <<
", and T = " << fluid_state.temperature(0) << std::endl;
111 using ScalarVector = Dune::FieldVector<Scalar, numComponents>;
112 Scalar L_scalar = Opm::getValue(L);
113 ScalarVector z_scalar;
114 ScalarVector K_scalar;
115 for (
unsigned i = 0; i < numComponents; ++i) {
116 z_scalar[i] = Opm::getValue(z[i]);
117 K_scalar[i] = Opm::getValue(K[i]);
120 ScalarFluidState fluid_state_scalar;
122 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
123 fluid_state_scalar.setMoleFraction(oilPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx)));
124 fluid_state_scalar.setMoleFraction(gasPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx)));
125 fluid_state_scalar.setKvalue(compIdx, Opm::getValue(fluid_state.K(compIdx)));
128 fluid_state_scalar.setLvalue(L_scalar);
130 fluid_state_scalar.setPressure(FluidSystem::oilPhaseIdx,
131 Opm::getValue(fluid_state.pressure(FluidSystem::oilPhaseIdx)));
132 fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
133 Opm::getValue(fluid_state.pressure(FluidSystem::gasPhaseIdx)));
135 fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
138 bool is_stable =
false;
139 if ( L <= 0 || L == 1 ) {
140 if (verbosity >= 1) {
141 std::cout <<
"Perform stability test (L <= 0 or L == 1)!" << std::endl;
143 phaseStabilityTest_(is_stable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
145 if (verbosity >= 1) {
146 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;
150 const bool is_single_phase = is_stable;
153 if ( !is_single_phase ) {
155 L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
156 flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
159 L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity);
161 fluid_state_scalar.setLvalue(L_scalar);
164 if (verbosity >= 1) {
165 std::cout <<
"********" << std::endl;
170 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
171 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
172 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x_i);
173 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, compIdx);
174 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y_i);
177 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
178 fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
179 fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
181 fluid_state.setLvalue(L_scalar);
183 updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
193 template <
class Flu
idState,
class ComponentVector>
194 static void solve(FluidState& fluid_state,
195 const ComponentVector& globalMolarities,
196 Scalar tolerance = 0.0)
200 using MaterialLawParams =
typename MaterialLaw::Params;
202 MaterialLawParams matParams;
203 solve<MaterialLaw>(fluid_state, matParams, globalMolarities, tolerance);
206 template <
class Vector>
207 static typename Vector::field_type solveRachfordRice_g_(
const Vector& K,
const Vector& z,
int verbosity)
211 using field_type =
typename Vector::field_type;
212 constexpr field_type tol = 1e-12;
213 constexpr int itmax = 10000;
214 field_type Kmin = K[0];
215 field_type Kmax = K[0];
216 for (
int compIdx = 1; compIdx < numComponents; ++compIdx){
217 if (K[compIdx] < Kmin)
219 else if (K[compIdx] >= Kmax)
223 auto Vmin = 1 / (1 - Kmax);
224 auto Vmax = 1 / (1 - Kmin);
226 auto V = (Vmin + Vmax)/2;
228 if (verbosity == 3 || verbosity == 4) {
229 std::cout <<
"Initial guess " << numComponents <<
"c : V = " << V <<
" and [Vmin, Vmax] = [" << Vmin <<
", " << Vmax <<
"]" << std::endl;
230 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"abs(step)" << std::setw(16) <<
"V" << std::endl;
233 for (
int iteration = 1; iteration < itmax; ++iteration) {
235 field_type denum = 0.0;
237 for (
int compIdx = 0; compIdx < numComponents; ++compIdx){
238 auto dK = K[compIdx] - 1.0;
239 auto a = z[compIdx] * dK;
240 auto b = (1 + V * dK);
242 denum += z[compIdx] * (dK*dK) / (b*b);
244 auto delta = r / denum;
248 if (V < Vmin || V > Vmax)
251 if (verbosity == 3 || verbosity == 4) {
252 std::cout <<
"V = " << V <<
" is not within the the range [Vmin, Vmax], solve using Bisection method!" << std::endl;
256 auto Lmin = 1 - Vmax;
257 auto Lmax = 1 - Vmin;
264 auto L = bisection_g_(K, Lmin, Lmax, z, verbosity);
267 if (verbosity >= 1) {
268 std::cout <<
"Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
274 if (verbosity == 3 || verbosity == 4) {
275 std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << V << std::endl;
278 if ( Opm::abs(r) < tol ) {
283 if (verbosity >= 1) {
284 std::cout <<
"Rachford-Rice converged to final solution L = " << L << std::endl;
291 throw std::runtime_error(
" Rachford-Rice did not converge within maximum number of iterations" );
294 template <
class Vector>
295 static typename Vector::field_type bisection_g_(
const Vector& K,
typename Vector::field_type Lmin,
296 typename Vector::field_type Lmax,
const Vector& z,
int verbosity)
299 typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
302 if (verbosity >= 3) {
303 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"g(Lmid)" << std::setw(16) <<
"L" << std::endl;
306 constexpr int max_it = 10000;
308 auto closeLmaxLmin = [](
double max_v,
double min_v) {
309 return Opm::abs(max_v - min_v) / 2. < 1e-10;
314 if (closeLmaxLmin(Lmax, Lmin) ){
315 throw std::runtime_error(fmt::format(
"Strange bisection with Lmax {} and Lmin {}?", Lmax, Lmin));
317 for (
int iteration = 0; iteration < max_it; ++iteration){
319 auto L = (Lmin + Lmax) / 2;
320 auto gMid = rachfordRice_g_(K, L, z);
322 if (verbosity == 3 || verbosity == 4) {
323 std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl;
327 if (Opm::abs(gMid) < 1e-16 || closeLmaxLmin(Lmax, Lmin)){
331 else if (Dune::sign(gMid) != Dune::sign(gLmin)) {
341 throw std::runtime_error(
" Rachford-Rice bisection failed with " + std::to_string(max_it) +
" iterations!");
344 template <
class Vector,
class FlashFlu
idState>
345 static typename Vector::field_type li_single_phase_label_(
const FlashFluidState& fluid_state,
const Vector& z,
int verbosity)
348 typename Vector::field_type sumVz = 0.0;
349 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
351 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
354 sumVz += (V_crit * z[compIdx]);
358 typename Vector::field_type Tc_est = 0.0;
359 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
361 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
362 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
365 Tc_est += (V_crit * T_crit * z[compIdx] / sumVz);
369 const auto& T = fluid_state.temperature(0);
372 typename Vector::field_type L;
378 if (verbosity >= 1) {
379 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;
387 if (verbosity >= 1) {
388 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;
396 template <
class FlashFlu
idState,
class ComponentVector>
397 static void phaseStabilityTest_(
bool& isStable, ComponentVector& K, FlashFluidState& fluid_state,
const ComponentVector& z,
int verbosity)
400 bool isTrivialL, isTrivialV;
401 ComponentVector x, y;
402 typename FlashFluidState::Scalar S_l, S_v;
403 ComponentVector K0 = K;
404 ComponentVector K1 = K;
407 if (verbosity == 3 || verbosity == 4) {
408 std::cout <<
"Stability test for vapor phase:" << std::endl;
410 checkStability_(fluid_state, isTrivialV, K0, y, S_v, z,
true, verbosity);
411 bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
414 if (verbosity == 3 || verbosity == 4) {
415 std::cout <<
"Stability test for liquid phase:" << std::endl;
417 checkStability_(fluid_state, isTrivialL, K1, x, S_l, z,
false, verbosity);
418 bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
421 isStable = L_stable && V_unstable;
425 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
426 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
427 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
432 for (
int compIdx = 0; compIdx<numComponents; ++compIdx) {
433 K[compIdx] = y[compIdx] / x[compIdx];
440 template <
class FlashFlu
idState>
441 static typename FlashFluidState::Scalar wilsonK_(
const FlashFluidState& fluid_state,
int compIdx)
443 const auto& acf = FluidSystem::acentricFactor(compIdx);
444 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
445 const auto& T = fluid_state.temperature(0);
446 const auto& p_crit = FluidSystem::criticalPressure(compIdx);
447 const auto& p = fluid_state.pressure(0);
449 const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
453 template <
class Vector>
454 static typename Vector::field_type rachfordRice_g_(
const Vector& K,
typename Vector::field_type L,
const Vector& z)
456 typename Vector::field_type g=0;
457 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
458 g += (z[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
463 template <
class Vector>
464 static typename Vector::field_type rachfordRice_dg_dL_(
const Vector& K,
const typename Vector::field_type L,
const Vector& z)
466 typename Vector::field_type dg=0;
467 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
468 dg += (z[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
473 template <
class FlashFlu
idState,
class ComponentVector>
474 static void checkStability_(
const FlashFluidState& fluid_state,
bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
475 typename FlashFluidState::Scalar& S_loc,
const ComponentVector& z,
bool isGas,
int verbosity)
477 using FlashEval =
typename FlashFluidState::Scalar;
481 FlashFluidState fluid_state_fake = fluid_state;
482 FlashFluidState fluid_state_global = fluid_state;
485 if (verbosity >= 3 || verbosity >= 4) {
486 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"K-Norm" << std::setw(16) <<
"R-Norm" << std::endl;
491 for (
int i = 0; i < 20000; ++i) {
494 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
495 xy_loc[compIdx] = K[compIdx] * z[compIdx];
496 S_loc += xy_loc[compIdx];
498 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
499 xy_loc[compIdx] /= S_loc;
500 fluid_state_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
504 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
505 xy_loc[compIdx] = z[compIdx]/K[compIdx];
506 S_loc += xy_loc[compIdx];
508 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
509 xy_loc[compIdx] /= S_loc;
510 fluid_state_fake.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
514 int phaseIdx = (isGas ?
static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
515 int phaseIdx2 = (isGas ?
static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
517 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
518 fluid_state_global.setMoleFraction(phaseIdx2, compIdx, z[compIdx]);
521 typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
522 paramCache_fake.updatePhase(fluid_state_fake, phaseIdx);
524 typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
525 paramCache_global.updatePhase(fluid_state_global, phaseIdx2);
528 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
532 fluid_state_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
533 fluid_state_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
538 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
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_global / fug_fake;
543 R[compIdx] = fug_ratio / S_loc;
546 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
547 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
548 auto fug_ratio = fug_fake / fug_global;
549 R[compIdx] = fug_ratio * S_loc;
553 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
554 K[compIdx] *= R[compIdx];
558 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
559 auto a = Opm::getValue(R[compIdx]) - 1.0;
560 auto b = Opm::log(Opm::getValue(K[compIdx]));
566 if (verbosity >= 3) {
567 std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
571 isTrivial = (K_norm < 1e-5);
572 if (isTrivial || R_norm < 1e-10)
577 throw std::runtime_error(
" Stability test did not converge");
581 template <
class FlashFlu
idState,
class ComponentVector>
582 static void computeLiquidVapor_(FlashFluidState& fluid_state,
typename FlashFluidState::Scalar& L, ComponentVector& K,
const ComponentVector& z)
587 typename FlashFluidState::Scalar sumx=0;
588 typename FlashFluidState::Scalar sumy=0;
589 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
590 x[compIdx] = z[compIdx]/(L + (1-L)*K[compIdx]);
592 y[compIdx] = (K[compIdx]*z[compIdx])/(L + (1-L)*K[compIdx]);
598 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
599 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
600 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
604 template <
class Flu
idState,
class ComponentVector>
605 static void flash_2ph(
const ComponentVector& z_scalar,
606 const std::string& flash_2p_method,
607 ComponentVector& K_scalar,
608 typename FluidState::Scalar& L_scalar,
609 FluidState& fluid_state_scalar,
611 if (verbosity >= 1) {
612 std::cout <<
"Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar <<
"]" << std::endl;
617 if (flash_2p_method ==
"newton") {
618 if (verbosity >= 1) {
619 std::cout <<
"Calculate composition using Newton." << std::endl;
621 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
622 }
else if (flash_2p_method ==
"ssi") {
624 if (verbosity >= 1) {
625 std::cout <<
"Calculate composition using Succcessive Substitution." << std::endl;
627 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
false, verbosity);
628 }
else if (flash_2p_method ==
"ssi+newton") {
629 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
true, verbosity);
630 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
632 throw std::runtime_error(
"unknown two phase flash method " + flash_2p_method +
" is specified");
636 template <
class FlashFlu
idState,
class ComponentVector>
637 static void newtonComposition_(ComponentVector& K,
typename FlashFluidState::Scalar& L,
638 FlashFluidState& fluid_state,
const ComponentVector& z,
643 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
644 constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
645 using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
646 using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
647 constexpr Scalar tolerance = 1.e-8;
649 NewtonVector soln(0.);
650 NewtonVector res(0.);
651 NewtonMatrix jac (0.);
654 computeLiquidVapor_(fluid_state, L, K, z);
655 if (verbosity >= 1) {
656 std::cout <<
" the current L is " << Opm::getValue(L) << std::endl;
660 if (verbosity >= 1) {
661 std::cout <<
"Initial guess: x = [";
662 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
663 if (compIdx < numComponents - 1)
664 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
666 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
668 std::cout <<
"], y = [";
669 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
670 if (compIdx < numComponents - 1)
671 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
673 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
675 std::cout <<
"], and " <<
"L = " << L << std::endl;
679 if (verbosity == 2 || verbosity == 4) {
680 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"Norm2(step)" << std::setw(16) <<
"Norm2(Residual)" << std::endl;
684 using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
686 std::vector<Eval> x(numComponents), y(numComponents);
690 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
691 x[compIdx] = Eval(fluid_state.moleFraction(oilPhaseIdx, compIdx), compIdx);
692 const unsigned idx = compIdx + numComponents;
693 y[compIdx] = Eval(fluid_state.moleFraction(gasPhaseIdx, compIdx), idx);
695 l = Eval(L, num_primary_variables - 1);
698 CompositionalFluidState<Eval, FluidSystem> flash_fluid_state;
699 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
700 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
701 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
703 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
705 flash_fluid_state.setLvalue(l);
707 flash_fluid_state.setPressure(FluidSystem::oilPhaseIdx,
708 fluid_state.pressure(FluidSystem::oilPhaseIdx));
709 flash_fluid_state.setPressure(FluidSystem::gasPhaseIdx,
710 fluid_state.pressure(FluidSystem::gasPhaseIdx));
713 flash_fluid_state.setSaturation(FluidSystem::gasPhaseIdx,
714 fluid_state.saturation(FluidSystem::gasPhaseIdx));
715 flash_fluid_state.setSaturation(FluidSystem::oilPhaseIdx,
716 fluid_state.saturation(FluidSystem::oilPhaseIdx));
717 flash_fluid_state.setTemperature(fluid_state.temperature(0));
719 using ParamCache =
typename FluidSystem::template ParameterCache<typename CompositionalFluidState<Eval, FluidSystem>::Scalar>;
720 ParamCache paramCache;
722 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
723 paramCache.updatePhase(flash_fluid_state, phaseIdx);
724 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
726 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
727 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
730 bool converged =
false;
732 constexpr unsigned max_iter = 1000;
733 while (iter < max_iter) {
734 assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
735 (flash_fluid_state, z, jac, res);
736 if (verbosity >= 1) {
737 std::cout <<
" newton residual is " << res.two_norm() << std::endl;
739 converged = res.two_norm() < tolerance;
744 jac.solve(soln, res);
745 constexpr Scalar damping_factor = 1.0;
747 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
748 x[compIdx] -= soln[compIdx] * damping_factor;
749 y[compIdx] -= soln[compIdx + numComponents] * damping_factor;
751 l -= soln[num_equations - 1] * damping_factor;
753 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
754 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
755 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
757 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
759 flash_fluid_state.setLvalue(l);
761 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
762 paramCache.updatePhase(flash_fluid_state, phaseIdx);
763 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
764 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
765 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);
798 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
799 static void assembleNewton_(
const FlashFluidState& fluid_state,
800 const ComponentVector& global_composition,
801 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
802 Dune::FieldVector<double, num_equation>& res)
804 using Eval = DenseAd::Evaluation<double, num_primary>;
805 std::vector<Eval> x(numComponents), y(numComponents);
806 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
807 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
808 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
810 const Eval& l = fluid_state.L();
814 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
817 auto local_res = -global_composition[compIdx] + l * x[compIdx] + (1 - l) * y[compIdx];
818 res[compIdx] = Opm::getValue(local_res);
819 for (
unsigned i = 0; i < num_primary; ++i) {
820 jac[compIdx][i] = local_res.derivative(i);
826 auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
827 fluid_state.fugacity(gasPhaseIdx, compIdx));
828 res[compIdx + numComponents] = Opm::getValue(local_res);
829 for (
unsigned i = 0; i < num_primary; ++i) {
830 jac[compIdx + numComponents][i] = local_res.derivative(i);
837 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
841 auto local_res = sumx - sumy;
842 res[num_equation - 1] = Opm::getValue(local_res);
843 for (
unsigned i = 0; i < num_primary; ++i) {
844 jac[num_equation - 1][i] = local_res.derivative(i);
849 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
850 static void assembleNewtonSingle_(
const FlashFluidState& fluid_state,
851 const ComponentVector& global_composition,
852 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
853 Dune::FieldVector<double, num_equation>& res)
855 using Eval = DenseAd::Evaluation<double, num_primary>;
856 std::vector<Eval> x(numComponents), y(numComponents);
857 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
858 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
859 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
861 const Eval& l = fluid_state.L();
866 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
869 auto local_res = -global_composition[compIdx] + x[compIdx];
870 res[compIdx] = Opm::getValue(local_res);
871 for (
unsigned i = 0; i < num_primary; ++i) {
872 jac[compIdx][i] = local_res.derivative(i);
878 auto local_res = -global_composition[compIdx] + y[compIdx];
879 res[compIdx + numComponents] = Opm::getValue(local_res);
880 for (
unsigned i = 0; i < num_primary; ++i) {
881 jac[compIdx + numComponents][i] = local_res.derivative(i);
887 const bool isGas = Opm::abs(l - 1.0) > std::numeric_limits<double>::epsilon();
895 res[num_equation - 1] = Opm::getValue(local_res);
896 for (
unsigned i = 0; i < num_primary; ++i) {
897 jac[num_equation - 1][i] = local_res.derivative(i);
901 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
902 static void updateDerivatives_(
const FlashFluidStateScalar& fluid_state_scalar,
903 const ComponentVector& z,
904 FluidState& fluid_state,
905 bool is_single_phase)
908 updateDerivativesTwoPhase_(fluid_state_scalar, z, fluid_state);
910 updateDerivativesSinglePhase_(fluid_state_scalar, z, fluid_state);
914 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
915 static void updateDerivativesTwoPhase_(
const FlashFluidStateScalar& fluid_state_scalar,
916 const ComponentVector& z,
917 FluidState& fluid_state)
920 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
921 constexpr size_t secondary_num_pv = numComponents + 1;
923 using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
926 SecondaryFlashFluidState secondary_fluid_state;
927 SecondaryComponentVector secondary_z;
930 const SecondaryEval sec_p = SecondaryEval(fluid_state_scalar.pressure(FluidSystem::oilPhaseIdx), 0);
931 secondary_fluid_state.setPressure(FluidSystem::oilPhaseIdx, sec_p);
932 secondary_fluid_state.setPressure(FluidSystem::gasPhaseIdx, sec_p);
935 secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));
937 for (
unsigned idx = 0; idx < numComponents; ++idx) {
938 secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
941 for (
unsigned idx = 0; idx < numComponents; ++idx) {
943 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, idx);
944 secondary_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
945 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, idx);
946 secondary_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
948 const auto K_i = fluid_state_scalar.K(idx);
949 secondary_fluid_state.setKvalue(idx, K_i);
951 const auto L = fluid_state_scalar.L();
952 secondary_fluid_state.setLvalue(L);
956 using SecondaryParamCache =
typename FluidSystem::template ParameterCache<typename SecondaryFlashFluidState::Scalar>;
957 SecondaryParamCache secondary_param_cache;
958 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
959 secondary_param_cache.updatePhase(secondary_fluid_state, phase_idx);
960 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
961 SecondaryEval phi = FluidSystem::fugacityCoefficient(secondary_fluid_state, secondary_param_cache, phase_idx, comp_idx);
962 secondary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
966 using SecondaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
967 using SecondaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv>;
968 SecondaryNewtonMatrix sec_jac;
969 SecondaryNewtonVector sec_res;
973 assembleNewton_<SecondaryFlashFluidState, SecondaryComponentVector, secondary_num_pv, num_equations>
974 (secondary_fluid_state, secondary_z, sec_jac, sec_res);
979 constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
981 using PrimaryComponentVector = Dune::FieldVector<double, numComponents>;
984 PrimaryFlashFluidState primary_fluid_state;
986 PrimaryComponentVector primary_z;
987 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
988 primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
990 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
991 const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
992 primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
993 const unsigned idx = comp_idx + numComponents;
994 const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
995 primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
996 primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
999 l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
1000 primary_fluid_state.setLvalue(l);
1001 primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
1002 primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
1003 primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
1006 using PrimaryParamCache =
typename FluidSystem::template ParameterCache<typename PrimaryFlashFluidState::Scalar>;
1007 PrimaryParamCache primary_param_cache;
1008 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
1009 primary_param_cache.updatePhase(primary_fluid_state, phase_idx);
1010 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
1011 PrimaryEval phi = FluidSystem::fugacityCoefficient(primary_fluid_state, primary_param_cache, phase_idx, comp_idx);
1012 primary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
1016 using PrimaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
1017 using PrimaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, primary_num_pv>;
1018 PrimaryNewtonVector pri_res;
1019 PrimaryNewtonMatrix pri_jac;
1022 assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
1023 (primary_fluid_state, primary_z, pri_jac, pri_res);
1029 sec_jac.template leftmultiply(pri_jac);
1031 using InputEval =
typename FluidState::Scalar;
1032 using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
1034 ComponentVectorMoleFraction x(numComponents), y(numComponents);
1035 InputEval L_eval = L;
1042 const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
1043 const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
1044 std::vector<double> K(numComponents);
1046 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1047 K[compIdx] = fluid_state_scalar.K(compIdx);
1048 x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);
1049 y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);
1056 constexpr size_t num_deri = numComponents;
1057 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1058 std::vector<double> deri(num_deri, 0.);
1060 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1061 deri[idx] = -sec_jac[compIdx][0] * p_l.derivative(idx);
1064 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1065 const double pz = -sec_jac[compIdx][cIdx + 1];
1066 const auto& zi = z[cIdx];
1067 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1068 deri[idx] += pz * zi.derivative(idx);
1071 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1072 x[compIdx].setDerivative(idx, deri[idx]);
1075 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1076 deri[idx] = -sec_jac[compIdx + numComponents][0] * p_v.derivative(idx);
1078 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1079 const double pz = -sec_jac[compIdx + numComponents][cIdx + 1];
1080 const auto& zi = z[cIdx];
1081 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1082 deri[idx] += pz * zi.derivative(idx);
1085 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1086 y[compIdx].setDerivative(idx, deri[idx]);
1090 std::vector<double> deriL(num_deri, 0.);
1091 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1092 deriL[idx] = -sec_jac[2 * numComponents][0] * p_v.derivative(idx);
1094 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1095 const double pz = -sec_jac[2 * numComponents][cIdx + 1];
1096 const auto& zi = z[cIdx];
1097 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1098 deriL[idx] += pz * zi.derivative(idx);
1102 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1103 L_eval.setDerivative(idx, deriL[idx]);
1108 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1109 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
1110 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
1112 fluid_state.setLvalue(L_eval);
1115 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
1116 static void updateDerivativesSinglePhase_(
const FlashFluidStateScalar& fluid_state_scalar,
1117 const ComponentVector& z,
1118 FluidState& fluid_state)
1120 using InputEval =
typename FluidState::Scalar;
1122 InputEval L_eval = fluid_state_scalar.L();;
1126 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1127 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
1128 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
1130 fluid_state.setLvalue(L_eval);
1135 template <
class FlashFlu
idState,
class ComponentVector>
1136 static void successiveSubstitutionComposition_(ComponentVector& K,
typename ComponentVector::field_type& L, FlashFluidState& fluid_state,
const ComponentVector& z,
1137 const bool newton_afterwards,
const int verbosity)
1140 const int maxIterations = newton_afterwards ? 3 : 100;
1143 std::ios_base::fmtflags f(std::cout.flags());
1147 std::cout <<
"Initial guess: K = [" << K <<
"] and L = " << L << std::endl;
1149 if (verbosity == 2 || verbosity == 4) {
1151 int fugWidth = (numComponents * 12)/2;
1152 int convWidth = fugWidth + 7;
1153 std::cout << std::setw(10) <<
"Iteration" << std::setw(fugWidth) <<
"fL/fV" << std::setw(convWidth) <<
"norm2(fL/fv-1)" << std::endl;
1158 for (
int i=0; i < maxIterations; ++i){
1160 computeLiquidVapor_(fluid_state, L, K, z);
1163 using ParamCache =
typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
1164 ParamCache paramCache;
1165 for (
int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
1166 paramCache.updatePhase(fluid_state, phaseIdx);
1167 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1168 auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, phaseIdx, compIdx);
1169 fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
1174 ComponentVector newFugRatio;
1175 ComponentVector convFugRatio;
1176 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1177 newFugRatio[compIdx] = fluid_state.fugacity(oilPhaseIdx, compIdx)/fluid_state.fugacity(gasPhaseIdx, compIdx);
1178 convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
1182 if (verbosity >= 2) {
1184 int fugWidth = (prec + 3);
1185 int convWidth = prec + 9;
1186 std::cout << std::defaultfloat;
1187 std::cout << std::fixed;
1188 std::cout << std::setw(5) << i;
1189 std::cout << std::setw(fugWidth);
1190 std::cout << std::setprecision(prec);
1191 std::cout << newFugRatio;
1192 std::cout << std::scientific;
1193 std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
1197 if (convFugRatio.two_norm() < 1e-6) {
1202 if (verbosity >= 1) {
1203 std::cout <<
"Solution converged to the following result :" << std::endl;
1204 std::cout <<
"x = [";
1205 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1206 if (compIdx < numComponents - 1)
1207 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
1209 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
1211 std::cout <<
"]" << std::endl;
1212 std::cout <<
"y = [";
1213 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1214 if (compIdx < numComponents - 1)
1215 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
1217 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
1219 std::cout <<
"]" << std::endl;
1220 std::cout <<
"K = [" << K <<
"]" << std::endl;
1221 std::cout <<
"L = " << L << std::endl;
1230 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1231 K[compIdx] *= newFugRatio[compIdx];
1235 L = solveRachfordRice_g_(K, z, 0);
1239 if (!newton_afterwards) {
1240 throw std::runtime_error(
1241 "Successive substitution composition update did not converge within maxIterations " + std::to_string(maxIterations));