62 enum { numPhases = FluidSystem::numPhases };
63 enum { numComponents = FluidSystem::numComponents };
64 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx};
65 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx};
66 enum { numMiscibleComponents = FluidSystem::numMiscibleComponents};
67 enum { numMisciblePhases = FluidSystem::numMisciblePhases};
71 numMisciblePhases*numMiscibleComponents
79 template <
class Flu
idState>
80 static void solve(FluidState& fluid_state,
81 const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
83 std::string twoPhaseMethod,
84 Scalar tolerance = -1.,
88 using InputEval =
typename FluidState::Scalar;
89 using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
92 tolerance = std::min<Scalar>(1e-3, 1e8 * std::numeric_limits<Scalar>::epsilon());
97 for(
int compIdx = 0; compIdx < numComponents; ++compIdx) {
98 K[compIdx] = fluid_state.K(compIdx);
105 if (verbosity >= 1) {
106 std::cout <<
"********" << std::endl;
107 std::cout <<
"Flash calculations on Cell " << spatialIdx << 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);
209 template <
class FlashFlu
idState>
210 static typename FlashFluidState::Scalar wilsonK_(
const FlashFluidState& fluid_state,
int compIdx)
212 const auto& acf = FluidSystem::acentricFactor(compIdx);
213 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
214 const auto& T = fluid_state.temperature(0);
215 const auto& p_crit = FluidSystem::criticalPressure(compIdx);
216 const auto& p = fluid_state.pressure(0);
218 const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
222 template <
class Vector,
class FlashFlu
idState>
223 static typename Vector::field_type li_single_phase_label_(
const FlashFluidState& fluid_state,
const Vector& z,
int verbosity)
226 typename Vector::field_type sumVz = 0.0;
227 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
229 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
232 sumVz += (V_crit * z[compIdx]);
236 typename Vector::field_type Tc_est = 0.0;
237 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
239 const auto& V_crit = FluidSystem::criticalVolume(compIdx);
240 const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
243 Tc_est += (V_crit * T_crit * z[compIdx] / sumVz);
247 const auto& T = fluid_state.temperature(0);
250 typename Vector::field_type L;
256 if (verbosity >= 1) {
257 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;
265 if (verbosity >= 1) {
266 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;
274 template <
class Vector>
275 static typename Vector::field_type rachfordRice_g_(
const Vector& K,
typename Vector::field_type L,
const Vector& z)
277 typename Vector::field_type g=0;
278 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
279 g += (z[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
284 template <
class Vector>
285 static typename Vector::field_type rachfordRice_dg_dL_(
const Vector& K,
const typename Vector::field_type L,
const Vector& z)
287 typename Vector::field_type dg=0;
288 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
289 dg += (z[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
294 template <
class Vector>
295 static typename Vector::field_type solveRachfordRice_g_(
const Vector& K,
const Vector& z,
int verbosity)
299 typename Vector::field_type Kmin = K[0];
300 typename Vector::field_type Kmax = K[0];
301 for (
int compIdx=1; compIdx<numComponents; ++compIdx){
302 if (K[compIdx] < Kmin)
304 else if (K[compIdx] >= Kmax)
309 auto Lmin = (Kmin / (Kmin - 1));
310 auto Lmax = Kmax / (Kmax - 1);
321 auto L = (Lmin + Lmax)/2;
324 if (verbosity == 3 || verbosity == 4) {
325 std::cout <<
"Initial guess: L = " << L <<
" and [Lmin, Lmax] = [" << Lmin <<
", " << Lmax <<
"]" << std::endl;
326 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"abs(step)" << std::setw(16) <<
"L" << std::endl;
330 for (
int iteration=1; iteration<100; ++iteration){
332 auto g = rachfordRice_g_(K, L, z);
333 auto dg_dL = rachfordRice_dg_dL_(K, L, z);
336 auto delta = g/dg_dL;
340 if (L < Lmin || L > Lmax)
343 if (verbosity == 3 || verbosity == 4) {
344 std::cout <<
"L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl;
348 L = bisection_g_(K, Lmin, Lmax, z, verbosity);
351 L = Opm::min(Opm::max(L, 0.0), 1.0);
354 if (verbosity >= 1) {
355 std::cout <<
"Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
361 if (verbosity == 3 || verbosity == 4) {
362 std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << L << std::endl;
365 if ( Opm::abs(delta) < 1e-10 ) {
367 L = Opm::min(Opm::max(L, 0.0), 1.0);
370 if (verbosity >= 1) {
371 std::cout <<
"Rachford-Rice converged to final solution L = " << L << std::endl;
377 throw std::runtime_error(
" Rachford-Rice did not converge within maximum number of iterations" );
380 template <
class Vector>
381 static typename Vector::field_type bisection_g_(
const Vector& K,
typename Vector::field_type Lmin,
382 typename Vector::field_type Lmax,
const Vector& z,
int verbosity)
385 typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
388 if (verbosity >= 3) {
389 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"g(Lmid)" << std::setw(16) <<
"L" << std::endl;
392 constexpr int max_it = 100;
394 for (
int iteration = 0; iteration < max_it; ++iteration){
396 auto L = (Lmin + Lmax) / 2;
397 auto gMid = rachfordRice_g_(K, L, z);
398 if (verbosity == 3 || verbosity == 4) {
399 std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl;
403 if (Opm::abs(gMid) < 1e-10 || Opm::abs((Lmax - Lmin) / 2) < 1e-10){
408 else if (Dune::sign(gMid) != Dune::sign(gLmin)) {
418 throw std::runtime_error(
" Rachford-Rice with bisection failed with " + std::to_string(max_it) +
" iterations!");
421 template <
class FlashFlu
idState,
class ComponentVector>
422 static void phaseStabilityTest_(
bool& isStable, ComponentVector& K, FlashFluidState& fluid_state,
const ComponentVector& z,
int verbosity)
425 bool isTrivialL, isTrivialV;
426 ComponentVector x, y;
427 typename FlashFluidState::Scalar S_l, S_v;
428 ComponentVector K0 = K;
429 ComponentVector K1 = K;
432 if (verbosity == 3 || verbosity == 4) {
433 std::cout <<
"Stability test for vapor phase:" << std::endl;
435 checkStability_(fluid_state, isTrivialV, K0, y, S_v, z,
true, verbosity);
436 bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
439 if (verbosity == 3 || verbosity == 4) {
440 std::cout <<
"Stability test for liquid phase:" << std::endl;
442 checkStability_(fluid_state, isTrivialL, K1, x, S_l, z,
false, verbosity);
443 bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
446 isStable = L_stable && V_unstable;
450 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
451 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
452 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
457 for (
int compIdx = 0; compIdx<numComponents; ++compIdx) {
458 K[compIdx] = y[compIdx] / x[compIdx];
463 template <
class FlashFlu
idState,
class ComponentVector>
464 static void checkStability_(
const FlashFluidState& fluid_state,
bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
465 typename FlashFluidState::Scalar& S_loc,
const ComponentVector& z,
bool isGas,
int verbosity)
467 using FlashEval =
typename FlashFluidState::Scalar;
471 FlashFluidState fluid_state_fake = fluid_state;
472 FlashFluidState fluid_state_global = fluid_state;
475 if (verbosity >= 3 || verbosity >= 4) {
476 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"K-Norm" << std::setw(16) <<
"R-Norm" << std::endl;
481 for (
int i = 0; i < 20000; ++i) {
484 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
485 xy_loc[compIdx] = K[compIdx] * z[compIdx];
486 S_loc += xy_loc[compIdx];
488 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
489 xy_loc[compIdx] /= S_loc;
490 fluid_state_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
494 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
495 xy_loc[compIdx] = z[compIdx]/K[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(oilPhaseIdx, compIdx, xy_loc[compIdx]);
504 int phaseIdx = (isGas ?
static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
505 int phaseIdx2 = (isGas ?
static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
507 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
508 fluid_state_global.setMoleFraction(phaseIdx2, compIdx, z[compIdx]);
511 typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
512 paramCache_fake.updatePhase(fluid_state_fake, phaseIdx);
514 typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
515 paramCache_global.updatePhase(fluid_state_global, phaseIdx2);
518 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
522 fluid_state_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
523 fluid_state_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
528 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
530 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
531 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
532 auto fug_ratio = fug_global / fug_fake;
533 R[compIdx] = fug_ratio / S_loc;
536 auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
537 auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
538 auto fug_ratio = fug_fake / fug_global;
539 R[compIdx] = fug_ratio * S_loc;
543 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
544 K[compIdx] *= R[compIdx];
548 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
549 auto a = Opm::getValue(R[compIdx]) - 1.0;
550 auto b = Opm::log(Opm::getValue(K[compIdx]));
556 if (verbosity >= 3) {
557 std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
561 isTrivial = (K_norm < 1e-5);
562 if (isTrivial || R_norm < 1e-10)
567 throw std::runtime_error(
" Stability test did not converge");
571 template <
class FlashFlu
idState,
class ComponentVector>
572 static void computeLiquidVapor_(FlashFluidState& fluid_state,
typename FlashFluidState::Scalar& L, ComponentVector& K,
const ComponentVector& z)
577 typename FlashFluidState::Scalar sumx=0;
578 typename FlashFluidState::Scalar sumy=0;
579 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
580 x[compIdx] = z[compIdx]/(L + (1-L)*K[compIdx]);
582 y[compIdx] = (K[compIdx]*z[compIdx])/(L + (1-L)*K[compIdx]);
588 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
589 fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
590 fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
594 template <
class Flu
idState,
class ComponentVector>
595 static void flash_2ph(
const ComponentVector& z_scalar,
596 const std::string& flash_2p_method,
597 ComponentVector& K_scalar,
598 typename FluidState::Scalar& L_scalar,
599 FluidState& fluid_state_scalar,
601 if (verbosity >= 1) {
602 std::cout <<
"Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar <<
"]" << std::endl;
607 if (flash_2p_method ==
"newton") {
608 if (verbosity >= 1) {
609 std::cout <<
"Calculate composition using Newton." << std::endl;
611 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
612 }
else if (flash_2p_method ==
"ssi") {
614 if (verbosity >= 1) {
615 std::cout <<
"Calculate composition using Succcessive Substitution." << std::endl;
617 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
false, verbosity);
618 }
else if (flash_2p_method ==
"ssi+newton") {
619 successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar,
true, verbosity);
620 newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
622 throw std::runtime_error(
"unknown two phase flash method " + flash_2p_method +
" is specified");
626 template <
class FlashFlu
idState,
class ComponentVector>
627 static void newtonComposition_(ComponentVector& K,
typename FlashFluidState::Scalar& L,
628 FlashFluidState& fluid_state,
const ComponentVector& z,
633 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
634 constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
635 using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
636 using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
637 constexpr Scalar tolerance = 1.e-8;
639 NewtonVector soln(0.);
640 NewtonVector res(0.);
641 NewtonMatrix jac (0.);
644 computeLiquidVapor_(fluid_state, L, K, z);
645 if (verbosity >= 1) {
646 std::cout <<
" the current L is " << Opm::getValue(L) << std::endl;
650 if (verbosity >= 1) {
651 std::cout <<
"Initial guess: x = [";
652 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
653 if (compIdx < numComponents - 1)
654 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
656 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
658 std::cout <<
"], y = [";
659 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
660 if (compIdx < numComponents - 1)
661 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
663 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
665 std::cout <<
"], and " <<
"L = " << L << std::endl;
669 if (verbosity == 2 || verbosity == 4) {
670 std::cout << std::setw(10) <<
"Iteration" << std::setw(16) <<
"Norm2(step)" << std::setw(16) <<
"Norm2(Residual)" << std::endl;
674 using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
676 std::vector<Eval> x(numComponents), y(numComponents);
680 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
681 x[compIdx] = Eval(fluid_state.moleFraction(oilPhaseIdx, compIdx), compIdx);
682 const unsigned idx = compIdx + numComponents;
683 y[compIdx] = Eval(fluid_state.moleFraction(gasPhaseIdx, compIdx), idx);
685 l = Eval(L, num_primary_variables - 1);
688 CompositionalFluidState<Eval, FluidSystem> flash_fluid_state;
689 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
690 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
691 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
693 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
695 flash_fluid_state.setLvalue(l);
697 flash_fluid_state.setPressure(FluidSystem::oilPhaseIdx,
698 fluid_state.pressure(FluidSystem::oilPhaseIdx));
699 flash_fluid_state.setPressure(FluidSystem::gasPhaseIdx,
700 fluid_state.pressure(FluidSystem::gasPhaseIdx));
703 flash_fluid_state.setSaturation(FluidSystem::gasPhaseIdx,
704 fluid_state.saturation(FluidSystem::gasPhaseIdx));
705 flash_fluid_state.setSaturation(FluidSystem::oilPhaseIdx,
706 fluid_state.saturation(FluidSystem::oilPhaseIdx));
707 flash_fluid_state.setTemperature(fluid_state.temperature(0));
709 using ParamCache =
typename FluidSystem::template ParameterCache<typename CompositionalFluidState<Eval, FluidSystem>::Scalar>;
710 ParamCache paramCache;
712 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
713 paramCache.updatePhase(flash_fluid_state, phaseIdx);
714 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
716 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
717 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
720 bool converged =
false;
722 constexpr unsigned max_iter = 1000;
723 while (iter < max_iter) {
724 assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
725 (flash_fluid_state, z, jac, res);
726 if (verbosity >= 1) {
727 std::cout <<
" newton residual is " << res.two_norm() << std::endl;
729 converged = res.two_norm() < tolerance;
734 jac.solve(soln, res);
735 constexpr Scalar damping_factor = 1.0;
737 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
738 x[compIdx] -= soln[compIdx] * damping_factor;
739 y[compIdx] -= soln[compIdx + numComponents] * damping_factor;
741 l -= soln[num_equations - 1] * damping_factor;
743 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
744 flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
745 flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
747 flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
749 flash_fluid_state.setLvalue(l);
751 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
752 paramCache.updatePhase(flash_fluid_state, phaseIdx);
753 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
754 Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
755 flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
759 if (verbosity >= 1) {
760 for (
unsigned i = 0; i < num_equations; ++i) {
761 for (
unsigned j = 0; j < num_primary_variables; ++j) {
762 std::cout <<
" " << jac[i][j];
764 std::cout << std::endl;
766 std::cout << std::endl;
769 throw std::runtime_error(
" Newton composition update did not converge within maxIterations " + std::to_string(max_iter));
773 for (
unsigned idx = 0; idx < numComponents; ++idx) {
774 const auto x_i = Opm::getValue(flash_fluid_state.moleFraction(oilPhaseIdx, idx));
775 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
776 const auto y_i = Opm::getValue(flash_fluid_state.moleFraction(gasPhaseIdx, idx));
777 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
778 const auto K_i = Opm::getValue(flash_fluid_state.K(idx));
779 fluid_state.setKvalue(idx, K_i);
783 L = Opm::getValue(l);
784 fluid_state.setLvalue(L);
788 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
789 static void assembleNewton_(
const FlashFluidState& fluid_state,
790 const ComponentVector& global_composition,
791 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
792 Dune::FieldVector<double, num_equation>& res)
794 using Eval = DenseAd::Evaluation<double, num_primary>;
795 std::vector<Eval> x(numComponents), y(numComponents);
796 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
797 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
798 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
800 const Eval& l = fluid_state.L();
804 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
807 auto local_res = -global_composition[compIdx] + l * x[compIdx] + (1 - l) * y[compIdx];
808 res[compIdx] = Opm::getValue(local_res);
809 for (
unsigned i = 0; i < num_primary; ++i) {
810 jac[compIdx][i] = local_res.derivative(i);
816 auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
817 fluid_state.fugacity(gasPhaseIdx, compIdx));
818 res[compIdx + numComponents] = Opm::getValue(local_res);
819 for (
unsigned i = 0; i < num_primary; ++i) {
820 jac[compIdx + numComponents][i] = local_res.derivative(i);
827 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
831 auto local_res = sumx - sumy;
832 res[num_equation - 1] = Opm::getValue(local_res);
833 for (
unsigned i = 0; i < num_primary; ++i) {
834 jac[num_equation - 1][i] = local_res.derivative(i);
839 template<
typename FlashFlu
idState,
typename ComponentVector,
size_t num_primary,
size_t num_equation >
840 static void assembleNewtonSingle_(
const FlashFluidState& fluid_state,
841 const ComponentVector& global_composition,
842 Dune::FieldMatrix<double, num_equation, num_primary>& jac,
843 Dune::FieldVector<double, num_equation>& res)
845 using Eval = DenseAd::Evaluation<double, num_primary>;
846 std::vector<Eval> x(numComponents), y(numComponents);
847 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
848 x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
849 y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
851 const Eval& l = fluid_state.L();
856 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
859 auto local_res = -global_composition[compIdx] + x[compIdx];
860 res[compIdx] = Opm::getValue(local_res);
861 for (
unsigned i = 0; i < num_primary; ++i) {
862 jac[compIdx][i] = local_res.derivative(i);
868 auto local_res = -global_composition[compIdx] + y[compIdx];
869 res[compIdx + numComponents] = Opm::getValue(local_res);
870 for (
unsigned i = 0; i < num_primary; ++i) {
871 jac[compIdx + numComponents][i] = local_res.derivative(i);
877 const bool isGas = Opm::abs(l - 1.0) > std::numeric_limits<double>::epsilon();
885 res[num_equation - 1] = Opm::getValue(local_res);
886 for (
unsigned i = 0; i < num_primary; ++i) {
887 jac[num_equation - 1][i] = local_res.derivative(i);
891 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
892 static void updateDerivatives_(
const FlashFluidStateScalar& fluid_state_scalar,
893 const ComponentVector& z,
894 FluidState& fluid_state,
895 bool is_single_phase)
898 updateDerivativesTwoPhase_(fluid_state_scalar, z, fluid_state);
900 updateDerivativesSinglePhase_(fluid_state_scalar, z, fluid_state);
904 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
905 static void updateDerivativesTwoPhase_(
const FlashFluidStateScalar& fluid_state_scalar,
906 const ComponentVector& z,
907 FluidState& fluid_state)
910 constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
911 constexpr size_t secondary_num_pv = numComponents + 1;
913 using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
916 SecondaryFlashFluidState secondary_fluid_state;
917 SecondaryComponentVector secondary_z;
920 const SecondaryEval sec_p = SecondaryEval(fluid_state_scalar.pressure(FluidSystem::oilPhaseIdx), 0);
921 secondary_fluid_state.setPressure(FluidSystem::oilPhaseIdx, sec_p);
922 secondary_fluid_state.setPressure(FluidSystem::gasPhaseIdx, sec_p);
925 secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));
927 for (
unsigned idx = 0; idx < numComponents; ++idx) {
928 secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
931 for (
unsigned idx = 0; idx < numComponents; ++idx) {
933 const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, idx);
934 secondary_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
935 const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, idx);
936 secondary_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
938 const auto K_i = fluid_state_scalar.K(idx);
939 secondary_fluid_state.setKvalue(idx, K_i);
941 const auto L = fluid_state_scalar.L();
942 secondary_fluid_state.setLvalue(L);
946 using SecondaryParamCache =
typename FluidSystem::template ParameterCache<typename SecondaryFlashFluidState::Scalar>;
947 SecondaryParamCache secondary_param_cache;
948 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
949 secondary_param_cache.updatePhase(secondary_fluid_state, phase_idx);
950 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
951 SecondaryEval phi = FluidSystem::fugacityCoefficient(secondary_fluid_state, secondary_param_cache, phase_idx, comp_idx);
952 secondary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
956 using SecondaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
957 using SecondaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv>;
958 SecondaryNewtonMatrix sec_jac;
959 SecondaryNewtonVector sec_res;
963 assembleNewton_<SecondaryFlashFluidState, SecondaryComponentVector, secondary_num_pv, num_equations>
964 (secondary_fluid_state, secondary_z, sec_jac, sec_res);
969 constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
971 using PrimaryComponentVector = Dune::FieldVector<double, numComponents>;
974 PrimaryFlashFluidState primary_fluid_state;
976 PrimaryComponentVector primary_z;
977 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
978 primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
980 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
981 const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
982 primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
983 const unsigned idx = comp_idx + numComponents;
984 const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
985 primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
986 primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
989 l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
990 primary_fluid_state.setLvalue(l);
991 primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
992 primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
993 primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
996 using PrimaryParamCache =
typename FluidSystem::template ParameterCache<typename PrimaryFlashFluidState::Scalar>;
997 PrimaryParamCache primary_param_cache;
998 for (
unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
999 primary_param_cache.updatePhase(primary_fluid_state, phase_idx);
1000 for (
unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
1001 PrimaryEval phi = FluidSystem::fugacityCoefficient(primary_fluid_state, primary_param_cache, phase_idx, comp_idx);
1002 primary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
1006 using PrimaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
1007 using PrimaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, primary_num_pv>;
1008 PrimaryNewtonVector pri_res;
1009 PrimaryNewtonMatrix pri_jac;
1012 assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
1013 (primary_fluid_state, primary_z, pri_jac, pri_res);
1019 sec_jac.template leftmultiply(pri_jac);
1021 using InputEval =
typename FluidState::Scalar;
1022 using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
1024 ComponentVectorMoleFraction x(numComponents), y(numComponents);
1025 InputEval L_eval = L;
1032 const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
1033 const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
1034 std::vector<double> K(numComponents);
1036 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1037 K[compIdx] = fluid_state_scalar.K(compIdx);
1038 x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);
1039 y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);
1046 constexpr size_t num_deri = numComponents;
1047 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1048 std::vector<double> deri(num_deri, 0.);
1050 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1051 deri[idx] = -sec_jac[compIdx][0] * p_l.derivative(idx);
1054 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1055 const double pz = -sec_jac[compIdx][cIdx + 1];
1056 const auto& zi = z[cIdx];
1057 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1058 deri[idx] += pz * zi.derivative(idx);
1061 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1062 x[compIdx].setDerivative(idx, deri[idx]);
1065 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1066 deri[idx] = -sec_jac[compIdx + numComponents][0] * p_v.derivative(idx);
1068 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1069 const double pz = -sec_jac[compIdx + numComponents][cIdx + 1];
1070 const auto& zi = z[cIdx];
1071 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1072 deri[idx] += pz * zi.derivative(idx);
1075 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1076 y[compIdx].setDerivative(idx, deri[idx]);
1080 std::vector<double> deriL(num_deri, 0.);
1081 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1082 deriL[idx] = -sec_jac[2 * numComponents][0] * p_v.derivative(idx);
1084 for (
unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1085 const double pz = -sec_jac[2 * numComponents][cIdx + 1];
1086 const auto& zi = z[cIdx];
1087 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1088 deriL[idx] += pz * zi.derivative(idx);
1092 for (
unsigned idx = 0; idx < num_deri; ++idx) {
1093 L_eval.setDerivative(idx, deriL[idx]);
1098 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1099 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
1100 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
1102 fluid_state.setLvalue(L_eval);
1105 template <
typename FlashFlu
idStateScalar,
typename Flu
idState,
typename ComponentVector>
1106 static void updateDerivativesSinglePhase_(
const FlashFluidStateScalar& fluid_state_scalar,
1107 const ComponentVector& z,
1108 FluidState& fluid_state)
1110 using InputEval =
typename FluidState::Scalar;
1112 InputEval L_eval = fluid_state_scalar.L();;
1116 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1117 fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
1118 fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
1120 fluid_state.setLvalue(L_eval);
1125 template <
class FlashFlu
idState,
class ComponentVector>
1126 static void successiveSubstitutionComposition_(ComponentVector& K,
typename ComponentVector::field_type& L, FlashFluidState& fluid_state,
const ComponentVector& z,
1127 const bool newton_afterwards,
const int verbosity)
1130 const int maxIterations = newton_afterwards ? 3 : 10;
1133 std::ios_base::fmtflags f(std::cout.flags());
1137 std::cout <<
"Initial guess: K = [" << K <<
"] and L = " << L << std::endl;
1139 if (verbosity == 2 || verbosity == 4) {
1141 int fugWidth = (numComponents * 12)/2;
1142 int convWidth = fugWidth + 7;
1143 std::cout << std::setw(10) <<
"Iteration" << std::setw(fugWidth) <<
"fL/fV" << std::setw(convWidth) <<
"norm2(fL/fv-1)" << std::endl;
1148 for (
int i=0; i < maxIterations; ++i){
1150 computeLiquidVapor_(fluid_state, L, K, z);
1153 using ParamCache =
typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
1154 ParamCache paramCache;
1155 for (
int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
1156 paramCache.updatePhase(fluid_state, phaseIdx);
1157 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1158 auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, phaseIdx, compIdx);
1159 fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
1164 ComponentVector newFugRatio;
1165 ComponentVector convFugRatio;
1166 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1167 newFugRatio[compIdx] = fluid_state.fugacity(oilPhaseIdx, compIdx)/fluid_state.fugacity(gasPhaseIdx, compIdx);
1168 convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
1172 if (verbosity == 2 || verbosity == 4) {
1174 int fugWidth = (prec + 3);
1175 int convWidth = prec + 9;
1176 std::cout << std::defaultfloat;
1177 std::cout << std::fixed;
1178 std::cout << std::setw(5) << i;
1179 std::cout << std::setw(fugWidth);
1180 std::cout << std::setprecision(prec);
1181 std::cout << newFugRatio;
1182 std::cout << std::scientific;
1183 std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
1187 if (convFugRatio.two_norm() < 1e-6){
1192 if (verbosity >= 1) {
1193 std::cout <<
"Solution converged to the following result :" << std::endl;
1194 std::cout <<
"x = [";
1195 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1196 if (compIdx < numComponents - 1)
1197 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) <<
" ";
1199 std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
1201 std::cout <<
"]" << std::endl;
1202 std::cout <<
"y = [";
1203 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1204 if (compIdx < numComponents - 1)
1205 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) <<
" ";
1207 std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
1209 std::cout <<
"]" << std::endl;
1210 std::cout <<
"K = [" << K <<
"]" << std::endl;
1211 std::cout <<
"L = " << L << std::endl;
1220 for (
int compIdx=0; compIdx<numComponents; ++compIdx){
1221 K[compIdx] *= newFugRatio[compIdx];
1225 L = solveRachfordRice_g_(K, z, 0);
1229 if (!newton_afterwards) {
1230 throw std::runtime_error(
1231 "Successive substitution composition update did not converge within maxIterations");