My Project
PTFlash.hpp
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  Copyright 2022 NORCE.
5  Copyright 2022 SINTEF Digital, Mathematics and Cybernetics.
6 
7  This file is part of the Open Porous Media project (OPM).
8 
9  OPM is free software: you can redistribute it and/or modify
10  it under the terms of the GNU General Public License as published by
11  the Free Software Foundation, either version 2 of the License, or
12  (at your option) any later version.
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 
21  Consult the COPYING file in the top-level source directory of this
22  module for the precise wording of the license and the list of
23  copyright holders.
24 */
29 #ifndef OPM_CHI_FLASH_HPP
30 #define OPM_CHI_FLASH_HPP
31 
41 
43 
44 #include <dune/common/fvector.hh>
45 #include <dune/common/fmatrix.hh>
46 #include <dune/common/version.hh>
47 #include <dune/common/classname.hh>
48 
49 #include <limits>
50 #include <iostream>
51 #include <iomanip>
52 #include <type_traits>
53 
54 namespace Opm {
55 
61 template <class Scalar, class FluidSystem>
62 class PTFlash
63 {
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}; //oil, gas
70  enum {
71  numEq =
72  numMisciblePhases+
73  numMisciblePhases*numMiscibleComponents
74  };
75 
76 public:
81  template <class FluidState>
82  static void solve(FluidState& fluid_state,
83  const Dune::FieldVector<typename FluidState::Scalar, numComponents>& z,
84  int spatialIdx,
85  std::string twoPhaseMethod,
86  Scalar tolerance = -1.,
87  int verbosity = 0)
88  {
89 
90  using InputEval = typename FluidState::Scalar;
91  using ComponentVector = Dune::FieldVector<typename FluidState::Scalar, numComponents>;
92 
93 #if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
94  Dune::FMatrixPrecision<InputEval>::set_singular_limit(1e-35);
95 #endif
96  if (tolerance <= 0) {
97  tolerance = std::min<Scalar>(1e-3, 1e8 * std::numeric_limits<Scalar>::epsilon());
98  }
99 
100  // K and L from previous timestep (wilson and -1 initially)
101  ComponentVector K;
102  for(int compIdx = 0; compIdx < numComponents; ++compIdx) {
103  K[compIdx] = fluid_state.K(compIdx);
104  }
105  InputEval L;
106  // TODO: L has all the derivatives to be all ZEROs here.
107  L = fluid_state.L();
108 
109  // Print header
110  if (verbosity >= 1) {
111  std::cout << "********" << std::endl;
112  std::cout << "Flash calculations on Cell " << spatialIdx << std::endl;
113  std::cout << "Inputs are K = [" << K << "], L = [" << L << "], z = [" << z << "], P = " << fluid_state.pressure(0) << ", and T = " << fluid_state.temperature(0) << std::endl;
114  }
115 
116  using ScalarVector = Dune::FieldVector<Scalar, numComponents>;
117  Scalar L_scalar = Opm::getValue(L);
118  ScalarVector z_scalar;
119  ScalarVector K_scalar;
120  for (unsigned i = 0; i < numComponents; ++i) {
121  z_scalar[i] = Opm::getValue(z[i]);
122  K_scalar[i] = Opm::getValue(K[i]);
123  }
124  using ScalarFluidState = CompositionalFluidState<Scalar, FluidSystem>;
125  ScalarFluidState fluid_state_scalar;
126 
127  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
128  fluid_state_scalar.setMoleFraction(oilPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(oilPhaseIdx, compIdx)));
129  fluid_state_scalar.setMoleFraction(gasPhaseIdx, compIdx, Opm::getValue(fluid_state.moleFraction(gasPhaseIdx, compIdx)));
130  fluid_state_scalar.setKvalue(compIdx, Opm::getValue(fluid_state.K(compIdx)));
131  }
132 
133  fluid_state_scalar.setLvalue(L_scalar);
134  // other values need to be Scalar, but I guess the fluidstate does not support it yet.
135  fluid_state_scalar.setPressure(FluidSystem::oilPhaseIdx,
136  Opm::getValue(fluid_state.pressure(FluidSystem::oilPhaseIdx)));
137  fluid_state_scalar.setPressure(FluidSystem::gasPhaseIdx,
138  Opm::getValue(fluid_state.pressure(FluidSystem::gasPhaseIdx)));
139 
140  fluid_state_scalar.setTemperature(Opm::getValue(fluid_state.temperature(0)));
141 
142  // Do a stability test to check if cell is is_single_phase-phase (do for all cells the first time).
143  bool is_stable = false;
144  if ( L <= 0 || L == 1 ) {
145  if (verbosity >= 1) {
146  std::cout << "Perform stability test (L <= 0 or L == 1)!" << std::endl;
147  }
148  phaseStabilityTest_(is_stable, K_scalar, fluid_state_scalar, z_scalar, verbosity);
149  }
150  if (verbosity >= 1) {
151  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;
152  }
153  // TODO: we do not need two variables is_stable and is_single_hase, while lacking a good name
154  // TODO: from the later code, is good if we knows whether single_phase_gas or single_phase_oil here
155  const bool is_single_phase = is_stable;
156 
157  // Update the composition if cell is two-phase
158  if ( !is_single_phase ) {
159  // Rachford Rice equation to get initial L for composition solver
160  L_scalar = solveRachfordRice_g_(K_scalar, z_scalar, verbosity);
161  flash_2ph(z_scalar, twoPhaseMethod, K_scalar, L_scalar, fluid_state_scalar, verbosity);
162  } else {
163  // Cell is one-phase. Use Li's phase labeling method to see if it's liquid or vapor
164  L_scalar = li_single_phase_label_(fluid_state_scalar, z_scalar, verbosity);
165  }
166  fluid_state_scalar.setLvalue(L_scalar);
167 
168  // Print footer
169  if (verbosity >= 1) {
170  std::cout << "********" << std::endl;
171  }
172 
173  // the flash solution process were performed in scalar form, after the flash calculation finishes,
174  // ensure that things in fluid_state_scalar is transformed to fluid_state
175  for (int compIdx=0; compIdx<numComponents; ++compIdx){
176  const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, compIdx);
177  fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x_i);
178  const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, compIdx);
179  fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y_i);
180  }
181 
182  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
183  fluid_state.setKvalue(compIdx, K_scalar[compIdx]);
184  fluid_state_scalar.setKvalue(compIdx, K_scalar[compIdx]);
185  }
186  fluid_state.setLvalue(L_scalar);
187  // we update the derivatives in fluid_state
188  updateDerivatives_(fluid_state_scalar, z, fluid_state, is_single_phase);
189  }//end solve
190 
198  template <class FluidState, class ComponentVector>
199  static void solve(FluidState& fluid_state,
200  const ComponentVector& globalMolarities,
201  Scalar tolerance = 0.0)
202  {
203  using MaterialTraits = NullMaterialTraits<Scalar, numPhases>;
204  using MaterialLaw = NullMaterial<MaterialTraits>;
205  using MaterialLawParams = typename MaterialLaw::Params;
206 
207  MaterialLawParams matParams;
208  solve<MaterialLaw>(fluid_state, matParams, globalMolarities, tolerance);
209  }
210 
211 
212 protected:
213 
214  template <class FlashFluidState>
215  static typename FlashFluidState::Scalar wilsonK_(const FlashFluidState& fluid_state, int compIdx)
216  {
217  const auto& acf = FluidSystem::acentricFactor(compIdx);
218  const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
219  const auto& T = fluid_state.temperature(0);
220  const auto& p_crit = FluidSystem::criticalPressure(compIdx);
221  const auto& p = fluid_state.pressure(0); //for now assume no capillary pressure
222 
223  const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
224  return tmp;
225  }
226 
227  template <class Vector, class FlashFluidState>
228  static typename Vector::field_type li_single_phase_label_(const FlashFluidState& fluid_state, const Vector& z, int verbosity)
229  {
230  // Calculate intermediate sum
231  typename Vector::field_type sumVz = 0.0;
232  for (int compIdx=0; compIdx<numComponents; ++compIdx){
233  // Get component information
234  const auto& V_crit = FluidSystem::criticalVolume(compIdx);
235 
236  // Sum calculation
237  sumVz += (V_crit * z[compIdx]);
238  }
239 
240  // Calculate approximate (pseudo) critical temperature using Li's method
241  typename Vector::field_type Tc_est = 0.0;
242  for (int compIdx=0; compIdx<numComponents; ++compIdx){
243  // Get component information
244  const auto& V_crit = FluidSystem::criticalVolume(compIdx);
245  const auto& T_crit = FluidSystem::criticalTemperature(compIdx);
246 
247  // Sum calculation
248  Tc_est += (V_crit * T_crit * z[compIdx] / sumVz);
249  }
250 
251  // Get temperature
252  const auto& T = fluid_state.temperature(0);
253 
254  // If temperature is below estimated critical temperature --> phase = liquid; else vapor
255  typename Vector::field_type L;
256  if (T < Tc_est) {
257  // Liquid
258  L = 1.0;
259 
260  // Print
261  if (verbosity >= 1) {
262  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;
263  }
264  }
265  else {
266  // Vapor
267  L = 0.0;
268 
269  // Print
270  if (verbosity >= 1) {
271  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;
272  }
273  }
274 
275 
276  return L;
277  }
278 
279  template <class Vector>
280  static typename Vector::field_type rachfordRice_g_(const Vector& K, typename Vector::field_type L, const Vector& z)
281  {
282  typename Vector::field_type g=0;
283  for (int compIdx=0; compIdx<numComponents; ++compIdx){
284  g += (z[compIdx]*(K[compIdx]-1))/(K[compIdx]-L*(K[compIdx]-1));
285  }
286  return g;
287  }
288 
289  template <class Vector>
290  static typename Vector::field_type rachfordRice_dg_dL_(const Vector& K, const typename Vector::field_type L, const Vector& z)
291  {
292  typename Vector::field_type dg=0;
293  for (int compIdx=0; compIdx<numComponents; ++compIdx){
294  dg += (z[compIdx]*(K[compIdx]-1)*(K[compIdx]-1))/((K[compIdx]-L*(K[compIdx]-1))*(K[compIdx]-L*(K[compIdx]-1)));
295  }
296  return dg;
297  }
298 
299  template <class Vector>
300  static typename Vector::field_type solveRachfordRice_g_(const Vector& K, const Vector& z, int verbosity)
301  {
302  // Find min and max K. Have to do a laborious for loop to avoid water component (where K=0)
303  // TODO: Replace loop with Dune::min_value() and Dune::max_value() when water component is properly handled
304  typename Vector::field_type Kmin = K[0];
305  typename Vector::field_type Kmax = K[0];
306  for (int compIdx=1; compIdx<numComponents; ++compIdx){
307  if (K[compIdx] < Kmin)
308  Kmin = K[compIdx];
309  else if (K[compIdx] >= Kmax)
310  Kmax = K[compIdx];
311  }
312 
313  // Lower and upper bound for solution
314  auto Lmin = (Kmin / (Kmin - 1));
315  auto Lmax = Kmax / (Kmax - 1);
316 
317  // Check if Lmin < Lmax, and switch if not
318  if (Lmin > Lmax)
319  {
320  auto Ltmp = Lmin;
321  Lmin = Lmax;
322  Lmax = Ltmp;
323  }
324 
325  // Initial guess
326  auto L = (Lmin + Lmax)/2;
327 
328  // Print initial guess and header
329  if (verbosity == 3 || verbosity == 4) {
330  std::cout << "Initial guess: L = " << L << " and [Lmin, Lmax] = [" << Lmin << ", " << Lmax << "]" << std::endl;
331  std::cout << std::setw(10) << "Iteration" << std::setw(16) << "abs(step)" << std::setw(16) << "L" << std::endl;
332  }
333 
334  // Newton-Raphson loop
335  for (int iteration=1; iteration<100; ++iteration){
336  // Calculate function and derivative values
337  auto g = rachfordRice_g_(K, L, z);
338  auto dg_dL = rachfordRice_dg_dL_(K, L, z);
339 
340  // Lnew = Lold - g/dg;
341  auto delta = g/dg_dL;
342  L -= delta;
343 
344  // Check if L is within the bounds, and if not, we apply bisection method
345  if (L < Lmin || L > Lmax)
346  {
347  // Print info
348  if (verbosity == 3 || verbosity == 4) {
349  std::cout << "L is not within the the range [Lmin, Lmax], solve using Bisection method!" << std::endl;
350  }
351 
352  // Run bisection
353  L = bisection_g_(K, Lmin, Lmax, z, verbosity);
354 
355  // Ensure that L is in the range (0, 1)
356  L = Opm::min(Opm::max(L, 0.0), 1.0);
357 
358  // Print final result
359  if (verbosity >= 1) {
360  std::cout << "Rachford-Rice (Bisection) converged to final solution L = " << L << std::endl;
361  }
362  return L;
363  }
364 
365  // Print iteration info
366  if (verbosity == 3 || verbosity == 4) {
367  std::cout << std::setw(10) << iteration << std::setw(16) << Opm::abs(delta) << std::setw(16) << L << std::endl;
368  }
369  // Check for convergence
370  if ( Opm::abs(delta) < 1e-10 ) {
371  // Ensure that L is in the range (0, 1)
372  L = Opm::min(Opm::max(L, 0.0), 1.0);
373 
374  // Print final result
375  if (verbosity >= 1) {
376  std::cout << "Rachford-Rice converged to final solution L = " << L << std::endl;
377  }
378  return L;
379  }
380  }
381  // Throw error if Rachford-Rice fails
382  throw std::runtime_error(" Rachford-Rice did not converge within maximum number of iterations" );
383  }
384 
385  template <class Vector>
386  static typename Vector::field_type bisection_g_(const Vector& K, typename Vector::field_type Lmin,
387  typename Vector::field_type Lmax, const Vector& z, int verbosity)
388  {
389  // Calculate for g(Lmin) for first comparison with gMid = g(L)
390  typename Vector::field_type gLmin = rachfordRice_g_(K, Lmin, z);
391 
392  // Print new header
393  if (verbosity >= 3) {
394  std::cout << std::setw(10) << "Iteration" << std::setw(16) << "g(Lmid)" << std::setw(16) << "L" << std::endl;
395  }
396 
397  constexpr int max_it = 100;
398  // Bisection loop
399  for (int iteration = 0; iteration < max_it; ++iteration){
400  // New midpoint
401  auto L = (Lmin + Lmax) / 2;
402  auto gMid = rachfordRice_g_(K, L, z);
403  if (verbosity == 3 || verbosity == 4) {
404  std::cout << std::setw(10) << iteration << std::setw(16) << gMid << std::setw(16) << L << std::endl;
405  }
406 
407  // Check if midpoint fulfills g=0 or L - Lmin is sufficiently small
408  if (Opm::abs(gMid) < 1e-10 || Opm::abs((Lmax - Lmin) / 2) < 1e-10){
409  return L;
410  }
411 
412  // Else we repeat with midpoint being either Lmin og Lmax (depending on the signs).
413  else if (Dune::sign(gMid) != Dune::sign(gLmin)) {
414  // gMid has same sign as gLmax, so we set L as the new Lmax
415  Lmax = L;
416  }
417  else {
418  // gMid and gLmin have same sign so we set L as the new Lmin
419  Lmin = L;
420  gLmin = gMid;
421  }
422  }
423  throw std::runtime_error(" Rachford-Rice with bisection failed with " + std::to_string(max_it) + " iterations!");
424  }
425 
426  template <class FlashFluidState, class ComponentVector>
427  static void phaseStabilityTest_(bool& isStable, ComponentVector& K, FlashFluidState& fluid_state, const ComponentVector& z, int verbosity)
428  {
429  // Declarations
430  bool isTrivialL, isTrivialV;
431  ComponentVector x, y;
432  typename FlashFluidState::Scalar S_l, S_v;
433  ComponentVector K0 = K;
434  ComponentVector K1 = K;
435 
436  // Check for vapour instable phase
437  if (verbosity == 3 || verbosity == 4) {
438  std::cout << "Stability test for vapor phase:" << std::endl;
439  }
440  checkStability_(fluid_state, isTrivialV, K0, y, S_v, z, /*isGas=*/true, verbosity);
441  bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
442 
443  // Check for liquids stable phase
444  if (verbosity == 3 || verbosity == 4) {
445  std::cout << "Stability test for liquid phase:" << std::endl;
446  }
447  checkStability_(fluid_state, isTrivialL, K1, x, S_l, z, /*isGas=*/false, verbosity);
448  bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
449 
450  // L-stable means success in making liquid, V-unstable means no success in making vapour
451  isStable = L_stable && V_unstable;
452  if (isStable) {
453  // Single phase, i.e. phase composition is equivalent to the global composition
454  // Update fluid_state with mole fraction
455  for (int compIdx=0; compIdx<numComponents; ++compIdx){
456  fluid_state.setMoleFraction(gasPhaseIdx, compIdx, z[compIdx]);
457  fluid_state.setMoleFraction(oilPhaseIdx, compIdx, z[compIdx]);
458  }
459  }
460  // If not stable: use the mole fractions from Michelsen's test to update K
461  else {
462  for (int compIdx = 0; compIdx<numComponents; ++compIdx) {
463  K[compIdx] = y[compIdx] / x[compIdx];
464  }
465  }
466  }
467 
468  template <class FlashFluidState, class ComponentVector>
469  static void checkStability_(const FlashFluidState& fluid_state, bool& isTrivial, ComponentVector& K, ComponentVector& xy_loc,
470  typename FlashFluidState::Scalar& S_loc, const ComponentVector& z, bool isGas, int verbosity)
471  {
472  using FlashEval = typename FlashFluidState::Scalar;
473  using PengRobinsonMixture = typename Opm::PengRobinsonMixture<Scalar, FluidSystem>;
474 
475  // Declarations
476  FlashFluidState fluid_state_fake = fluid_state;
477  FlashFluidState fluid_state_global = fluid_state;
478 
479  // Setup output
480  if (verbosity >= 3 || verbosity >= 4) {
481  std::cout << std::setw(10) << "Iteration" << std::setw(16) << "K-Norm" << std::setw(16) << "R-Norm" << std::endl;
482  }
483 
484  // Michelsens stability test.
485  // Make two fake phases "inside" one phase and check for positive volume
486  for (int i = 0; i < 20000; ++i) {
487  S_loc = 0.0;
488  if (isGas) {
489  for (int compIdx=0; compIdx<numComponents; ++compIdx){
490  xy_loc[compIdx] = K[compIdx] * z[compIdx];
491  S_loc += xy_loc[compIdx];
492  }
493  for (int compIdx=0; compIdx<numComponents; ++compIdx){
494  xy_loc[compIdx] /= S_loc;
495  fluid_state_fake.setMoleFraction(gasPhaseIdx, compIdx, xy_loc[compIdx]);
496  }
497  }
498  else {
499  for (int compIdx=0; compIdx<numComponents; ++compIdx){
500  xy_loc[compIdx] = z[compIdx]/K[compIdx];
501  S_loc += xy_loc[compIdx];
502  }
503  for (int compIdx=0; compIdx<numComponents; ++compIdx){
504  xy_loc[compIdx] /= S_loc;
505  fluid_state_fake.setMoleFraction(oilPhaseIdx, compIdx, xy_loc[compIdx]);
506  }
507  }
508 
509  int phaseIdx = (isGas ? static_cast<int>(gasPhaseIdx) : static_cast<int>(oilPhaseIdx));
510  int phaseIdx2 = (isGas ? static_cast<int>(oilPhaseIdx) : static_cast<int>(gasPhaseIdx));
511  // TODO: not sure the following makes sense
512  for (int compIdx=0; compIdx<numComponents; ++compIdx){
513  fluid_state_global.setMoleFraction(phaseIdx2, compIdx, z[compIdx]);
514  }
515 
516  typename FluidSystem::template ParameterCache<FlashEval> paramCache_fake;
517  paramCache_fake.updatePhase(fluid_state_fake, phaseIdx);
518 
519  typename FluidSystem::template ParameterCache<FlashEval> paramCache_global;
520  paramCache_global.updatePhase(fluid_state_global, phaseIdx2);
521 
522  //fugacity for fake phases each component
523  for (int compIdx=0; compIdx<numComponents; ++compIdx){
524  auto phiFake = PengRobinsonMixture::computeFugacityCoefficient(fluid_state_fake, paramCache_fake, phaseIdx, compIdx);
525  auto phiGlobal = PengRobinsonMixture::computeFugacityCoefficient(fluid_state_global, paramCache_global, phaseIdx2, compIdx);
526 
527  fluid_state_fake.setFugacityCoefficient(phaseIdx, compIdx, phiFake);
528  fluid_state_global.setFugacityCoefficient(phaseIdx2, compIdx, phiGlobal);
529  }
530 
531 
532  ComponentVector R;
533  for (int compIdx=0; compIdx<numComponents; ++compIdx){
534  if (isGas){
535  auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
536  auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
537  auto fug_ratio = fug_global / fug_fake;
538  R[compIdx] = fug_ratio / S_loc;
539  }
540  else{
541  auto fug_fake = fluid_state_fake.fugacity(phaseIdx, compIdx);
542  auto fug_global = fluid_state_global.fugacity(phaseIdx2, compIdx);
543  auto fug_ratio = fug_fake / fug_global;
544  R[compIdx] = fug_ratio * S_loc;
545  }
546  }
547 
548  for (int compIdx=0; compIdx<numComponents; ++compIdx){
549  K[compIdx] *= R[compIdx];
550  }
551  Scalar R_norm = 0.0;
552  Scalar K_norm = 0.0;
553  for (int compIdx=0; compIdx<numComponents; ++compIdx){
554  auto a = Opm::getValue(R[compIdx]) - 1.0;
555  auto b = Opm::log(Opm::getValue(K[compIdx]));
556  R_norm += a*a;
557  K_norm += b*b;
558  }
559 
560  // Print iteration info
561  if (verbosity >= 3) {
562  std::cout << std::setw(10) << i << std::setw(16) << K_norm << std::setw(16) << R_norm << std::endl;
563  }
564 
565  // Check convergence
566  isTrivial = (K_norm < 1e-5);
567  if (isTrivial || R_norm < 1e-10)
568  return;
569  //todo: make sure that no mole fraction is smaller than 1e-8 ?
570  //todo: take care of water!
571  }
572  throw std::runtime_error(" Stability test did not converge");
573  }//end checkStability
574 
575  // TODO: basically FlashFluidState and ComponentVector are both depending on the one Scalar type
576  template <class FlashFluidState, class ComponentVector>
577  static void computeLiquidVapor_(FlashFluidState& fluid_state, typename FlashFluidState::Scalar& L, ComponentVector& K, const ComponentVector& z)
578  {
579  // Calculate x and y, and normalize
580  ComponentVector x;
581  ComponentVector y;
582  typename FlashFluidState::Scalar sumx=0;
583  typename FlashFluidState::Scalar sumy=0;
584  for (int compIdx=0; compIdx<numComponents; ++compIdx){
585  x[compIdx] = z[compIdx]/(L + (1-L)*K[compIdx]);
586  sumx += x[compIdx];
587  y[compIdx] = (K[compIdx]*z[compIdx])/(L + (1-L)*K[compIdx]);
588  sumy += y[compIdx];
589  }
590  x /= sumx;
591  y /= sumy;
592 
593  for (int compIdx=0; compIdx<numComponents; ++compIdx){
594  fluid_state.setMoleFraction(oilPhaseIdx, compIdx, x[compIdx]);
595  fluid_state.setMoleFraction(gasPhaseIdx, compIdx, y[compIdx]);
596  }
597  }
598 
599  template <class FluidState, class ComponentVector>
600  static void flash_2ph(const ComponentVector& z_scalar,
601  const std::string& flash_2p_method,
602  ComponentVector& K_scalar,
603  typename FluidState::Scalar& L_scalar,
604  FluidState& fluid_state_scalar,
605  int verbosity = 0) {
606  if (verbosity >= 1) {
607  std::cout << "Cell is two-phase! Solve Rachford-Rice with initial K = [" << K_scalar << "]" << std::endl;
608  }
609 
610  // Calculate composition using nonlinear solver
611  // Newton
612  if (flash_2p_method == "newton") {
613  if (verbosity >= 1) {
614  std::cout << "Calculate composition using Newton." << std::endl;
615  }
616  newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
617  } else if (flash_2p_method == "ssi") {
618  // Successive substitution
619  if (verbosity >= 1) {
620  std::cout << "Calculate composition using Succcessive Substitution." << std::endl;
621  }
622  successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, false, verbosity);
623  } else if (flash_2p_method == "ssi+newton") {
624  successiveSubstitutionComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, true, verbosity);
625  newtonComposition_(K_scalar, L_scalar, fluid_state_scalar, z_scalar, verbosity);
626  } else {
627  throw std::runtime_error("unknown two phase flash method " + flash_2p_method + " is specified");
628  }
629  }
630 
631  template <class FlashFluidState, class ComponentVector>
632  static void newtonComposition_(ComponentVector& K, typename FlashFluidState::Scalar& L,
633  FlashFluidState& fluid_state, const ComponentVector& z,
634  int verbosity)
635  {
636  // Note: due to the need for inverse flash update for derivatives, the following two can be different
637  // Looking for a good way to organize them
638  constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
639  constexpr size_t num_primary_variables = numMisciblePhases * numMiscibleComponents + 1;
640  using NewtonVector = Dune::FieldVector<Scalar, num_equations>;
641  using NewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, num_primary_variables>;
642  constexpr Scalar tolerance = 1.e-8;
643 
644  NewtonVector soln(0.);
645  NewtonVector res(0.);
646  NewtonMatrix jac (0.);
647 
648  // Compute x and y from K, L and Z
649  computeLiquidVapor_(fluid_state, L, K, z);
650  if (verbosity >= 1) {
651  std::cout << " the current L is " << Opm::getValue(L) << std::endl;
652  }
653 
654  // Print initial condition
655  if (verbosity >= 1) {
656  std::cout << "Initial guess: x = [";
657  for (int compIdx=0; compIdx<numComponents; ++compIdx){
658  if (compIdx < numComponents - 1)
659  std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
660  else
661  std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
662  }
663  std::cout << "], y = [";
664  for (int compIdx=0; compIdx<numComponents; ++compIdx){
665  if (compIdx < numComponents - 1)
666  std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
667  else
668  std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
669  }
670  std::cout << "], and " << "L = " << L << std::endl;
671  }
672 
673  // Print header
674  if (verbosity == 2 || verbosity == 4) {
675  std::cout << std::setw(10) << "Iteration" << std::setw(16) << "Norm2(step)" << std::setw(16) << "Norm2(Residual)" << std::endl;
676  }
677 
678  // AD type
679  using Eval = DenseAd::Evaluation<Scalar, num_primary_variables>;
680  // TODO: we might need to use numMiscibleComponents here
681  std::vector<Eval> x(numComponents), y(numComponents);
682  Eval l;
683 
684  // TODO: I might not need to set soln anything here.
685  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
686  x[compIdx] = Eval(fluid_state.moleFraction(oilPhaseIdx, compIdx), compIdx);
687  const unsigned idx = compIdx + numComponents;
688  y[compIdx] = Eval(fluid_state.moleFraction(gasPhaseIdx, compIdx), idx);
689  }
690  l = Eval(L, num_primary_variables - 1);
691 
692  // it is created for the AD calculation for the flash calculation
693  CompositionalFluidState<Eval, FluidSystem> flash_fluid_state;
694  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
695  flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
696  flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
697  // TODO: should we use wilsonK_ here?
698  flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
699  }
700  flash_fluid_state.setLvalue(l);
701  // other values need to be Scalar, but I guess the fluid_state does not support it yet.
702  flash_fluid_state.setPressure(FluidSystem::oilPhaseIdx,
703  fluid_state.pressure(FluidSystem::oilPhaseIdx));
704  flash_fluid_state.setPressure(FluidSystem::gasPhaseIdx,
705  fluid_state.pressure(FluidSystem::gasPhaseIdx));
706 
707  // TODO: not sure whether we need to set the saturations
708  flash_fluid_state.setSaturation(FluidSystem::gasPhaseIdx,
709  fluid_state.saturation(FluidSystem::gasPhaseIdx));
710  flash_fluid_state.setSaturation(FluidSystem::oilPhaseIdx,
711  fluid_state.saturation(FluidSystem::oilPhaseIdx));
712  flash_fluid_state.setTemperature(fluid_state.temperature(0));
713 
714  using ParamCache = typename FluidSystem::template ParameterCache<typename CompositionalFluidState<Eval, FluidSystem>::Scalar>;
715  ParamCache paramCache;
716 
717  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
718  paramCache.updatePhase(flash_fluid_state, phaseIdx);
719  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
720  // TODO: will phi here carry the correct derivatives?
721  Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
722  flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
723  }
724  }
725  bool converged = false;
726  unsigned iter = 0;
727  constexpr unsigned max_iter = 1000;
728  while (iter < max_iter) {
729  assembleNewton_<CompositionalFluidState<Eval, FluidSystem>, ComponentVector, num_primary_variables, num_equations>
730  (flash_fluid_state, z, jac, res);
731  if (verbosity >= 1) {
732  std::cout << " newton residual is " << res.two_norm() << std::endl;
733  }
734  converged = res.two_norm() < tolerance;
735  if (converged) {
736  break;
737  }
738 
739  jac.solve(soln, res);
740  constexpr Scalar damping_factor = 1.0;
741  // updating x and y
742  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
743  x[compIdx] -= soln[compIdx] * damping_factor;
744  y[compIdx] -= soln[compIdx + numComponents] * damping_factor;
745  }
746  l -= soln[num_equations - 1] * damping_factor;
747 
748  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
749  flash_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
750  flash_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
751  // TODO: should we use wilsonK_ here?
752  flash_fluid_state.setKvalue(compIdx, y[compIdx] / x[compIdx]);
753  }
754  flash_fluid_state.setLvalue(l);
755 
756  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
757  paramCache.updatePhase(flash_fluid_state, phaseIdx);
758  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
759  Eval phi = FluidSystem::fugacityCoefficient(flash_fluid_state, paramCache, phaseIdx, compIdx);
760  flash_fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
761  }
762  }
763  }
764  if (verbosity >= 1) {
765  for (unsigned i = 0; i < num_equations; ++i) {
766  for (unsigned j = 0; j < num_primary_variables; ++j) {
767  std::cout << " " << jac[i][j];
768  }
769  std::cout << std::endl;
770  }
771  std::cout << std::endl;
772  }
773  if (!converged) {
774  throw std::runtime_error(" Newton composition update did not converge within maxIterations " + std::to_string(max_iter));
775  }
776 
777  // fluid_state is scalar, we need to update all the values for fluid_state here
778  for (unsigned idx = 0; idx < numComponents; ++idx) {
779  const auto x_i = Opm::getValue(flash_fluid_state.moleFraction(oilPhaseIdx, idx));
780  fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
781  const auto y_i = Opm::getValue(flash_fluid_state.moleFraction(gasPhaseIdx, idx));
782  fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
783  const auto K_i = Opm::getValue(flash_fluid_state.K(idx));
784  fluid_state.setKvalue(idx, K_i);
785  // TODO: not sure we need K and L here, because they are in the flash_fluid_state anyway.
786  K[idx] = K_i;
787  }
788  L = Opm::getValue(l);
789  fluid_state.setLvalue(L);
790  }
791 
792  // TODO: the interface will need to refactor for later usage
793  template<typename FlashFluidState, typename ComponentVector, size_t num_primary, size_t num_equation >
794  static void assembleNewton_(const FlashFluidState& fluid_state,
795  const ComponentVector& global_composition,
796  Dune::FieldMatrix<double, num_equation, num_primary>& jac,
797  Dune::FieldVector<double, num_equation>& res)
798  {
799  using Eval = DenseAd::Evaluation<double, num_primary>;
800  std::vector<Eval> x(numComponents), y(numComponents);
801  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
802  x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
803  y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
804  }
805  const Eval& l = fluid_state.L();
806  // TODO: clearing zero whether necessary?
807  jac = 0.;
808  res = 0.;
809  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
810  {
811  // z - L*x - (1-L) * y
812  auto local_res = -global_composition[compIdx] + l * x[compIdx] + (1 - l) * y[compIdx];
813  res[compIdx] = Opm::getValue(local_res);
814  for (unsigned i = 0; i < num_primary; ++i) {
815  jac[compIdx][i] = local_res.derivative(i);
816  }
817  }
818 
819  {
820  // f_liquid - f_vapor = 0
821  auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
822  fluid_state.fugacity(gasPhaseIdx, compIdx));
823  res[compIdx + numComponents] = Opm::getValue(local_res);
824  for (unsigned i = 0; i < num_primary; ++i) {
825  jac[compIdx + numComponents][i] = local_res.derivative(i);
826  }
827  }
828  }
829  // sum(x) - sum(y) = 0
830  Eval sumx = 0.;
831  Eval sumy = 0.;
832  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
833  sumx += x[compIdx];
834  sumy += y[compIdx];
835  }
836  auto local_res = sumx - sumy;
837  res[num_equation - 1] = Opm::getValue(local_res);
838  for (unsigned i = 0; i < num_primary; ++i) {
839  jac[num_equation - 1][i] = local_res.derivative(i);
840  }
841  }
842 
843  // TODO: the interface will need to refactor for later usage
844  template<typename FlashFluidState, typename ComponentVector, size_t num_primary, size_t num_equation >
845  static void assembleNewtonSingle_(const FlashFluidState& fluid_state,
846  const ComponentVector& global_composition,
847  Dune::FieldMatrix<double, num_equation, num_primary>& jac,
848  Dune::FieldVector<double, num_equation>& res)
849  {
850  using Eval = DenseAd::Evaluation<double, num_primary>;
851  std::vector<Eval> x(numComponents), y(numComponents);
852  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
853  x[compIdx] = fluid_state.moleFraction(oilPhaseIdx, compIdx);
854  y[compIdx] = fluid_state.moleFraction(gasPhaseIdx, compIdx);
855  }
856  const Eval& l = fluid_state.L();
857 
858  // TODO: clearing zero whether necessary?
859  jac = 0.;
860  res = 0.;
861  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
862  {
863  // z - L*x - (1-L) * y ---> z - x;
864  auto local_res = -global_composition[compIdx] + x[compIdx];
865  res[compIdx] = Opm::getValue(local_res);
866  for (unsigned i = 0; i < num_primary; ++i) {
867  jac[compIdx][i] = local_res.derivative(i);
868  }
869  }
870 
871  {
872  // f_liquid - f_vapor = 0 -->z - y;
873  auto local_res = -global_composition[compIdx] + y[compIdx];
874  res[compIdx + numComponents] = Opm::getValue(local_res);
875  for (unsigned i = 0; i < num_primary; ++i) {
876  jac[compIdx + numComponents][i] = local_res.derivative(i);
877  }
878  }
879  }
880 
881  // TODO: better we have isGas or isLiquid here
882  const bool isGas = Opm::abs(l - 1.0) > std::numeric_limits<double>::epsilon();
883 
884  // sum(x) - sum(y) = 0
885  auto local_res = l;
886  if(isGas) {
887  local_res = l-1;
888  }
889 
890  res[num_equation - 1] = Opm::getValue(local_res);
891  for (unsigned i = 0; i < num_primary; ++i) {
892  jac[num_equation - 1][i] = local_res.derivative(i);
893  }
894  }
895 
896  template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
897  static void updateDerivatives_(const FlashFluidStateScalar& fluid_state_scalar,
898  const ComponentVector& z,
899  FluidState& fluid_state,
900  bool is_single_phase)
901  {
902  if(!is_single_phase)
903  updateDerivativesTwoPhase_(fluid_state_scalar, z, fluid_state);
904  else
905  updateDerivativesSinglePhase_(fluid_state_scalar, z, fluid_state);
906 
907  }
908 
909  template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
910  static void updateDerivativesTwoPhase_(const FlashFluidStateScalar& fluid_state_scalar,
911  const ComponentVector& z,
912  FluidState& fluid_state)
913  {
914  // getting the secondary Jocobian matrix
915  constexpr size_t num_equations = numMisciblePhases * numMiscibleComponents + 1;
916  constexpr size_t secondary_num_pv = numComponents + 1; // pressure, z for all the components
917  using SecondaryEval = Opm::DenseAd::Evaluation<double, secondary_num_pv>; // three z and one pressure
918  using SecondaryComponentVector = Dune::FieldVector<SecondaryEval, numComponents>;
919  using SecondaryFlashFluidState = Opm::CompositionalFluidState<SecondaryEval, FluidSystem>;
920 
921  SecondaryFlashFluidState secondary_fluid_state;
922  SecondaryComponentVector secondary_z;
923  // p and z are the primary variables here
924  // pressure
925  const SecondaryEval sec_p = SecondaryEval(fluid_state_scalar.pressure(FluidSystem::oilPhaseIdx), 0);
926  secondary_fluid_state.setPressure(FluidSystem::oilPhaseIdx, sec_p);
927  secondary_fluid_state.setPressure(FluidSystem::gasPhaseIdx, sec_p);
928 
929  // set the temperature // TODO: currently we are not considering the temperature derivatives
930  secondary_fluid_state.setTemperature(Opm::getValue(fluid_state_scalar.temperature(0)));
931 
932  for (unsigned idx = 0; idx < numComponents; ++idx) {
933  secondary_z[idx] = SecondaryEval(Opm::getValue(z[idx]), idx + 1);
934  }
935  // set up the mole fractions
936  for (unsigned idx = 0; idx < num_equations; ++idx) {
937  // TODO: double checking that fluid_state_scalar returns a scalar here
938  const auto x_i = fluid_state_scalar.moleFraction(oilPhaseIdx, idx);
939  secondary_fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, idx, x_i);
940  const auto y_i = fluid_state_scalar.moleFraction(gasPhaseIdx, idx);
941  secondary_fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, idx, y_i);
942  // TODO: double checking make sure those are consistent
943  const auto K_i = fluid_state_scalar.K(idx);
944  secondary_fluid_state.setKvalue(idx, K_i);
945  }
946  const auto L = fluid_state_scalar.L();
947  secondary_fluid_state.setLvalue(L);
948  // TODO: Do we need to update the saturations?
949  // compositions
950  // TODO: we probably can simplify SecondaryFlashFluidState::Scalar
951  using SecondaryParamCache = typename FluidSystem::template ParameterCache<typename SecondaryFlashFluidState::Scalar>;
952  SecondaryParamCache secondary_param_cache;
953  for (unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
954  secondary_param_cache.updatePhase(secondary_fluid_state, phase_idx);
955  for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
956  SecondaryEval phi = FluidSystem::fugacityCoefficient(secondary_fluid_state, secondary_param_cache, phase_idx, comp_idx);
957  secondary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
958  }
959  }
960 
961  using SecondaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
962  using SecondaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, secondary_num_pv>;
963  SecondaryNewtonMatrix sec_jac;
964  SecondaryNewtonVector sec_res;
965 
966 
967  //use the regular equations
968  assembleNewton_<SecondaryFlashFluidState, SecondaryComponentVector, secondary_num_pv, num_equations>
969  (secondary_fluid_state, secondary_z, sec_jac, sec_res);
970 
971 
972  // assembly the major matrix here
973  // primary variables are x, y and L
974  constexpr size_t primary_num_pv = numMisciblePhases * numMiscibleComponents + 1;
976  using PrimaryComponentVector = Dune::FieldVector<double, numComponents>;
977  using PrimaryFlashFluidState = Opm::CompositionalFluidState<PrimaryEval, FluidSystem>;
978 
979  PrimaryFlashFluidState primary_fluid_state;
980  // primary_z is not needed, because we use z will be okay here
981  PrimaryComponentVector primary_z;
982  for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
983  primary_z[comp_idx] = Opm::getValue(z[comp_idx]);
984  }
985  for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
986  const auto x_ii = PrimaryEval(fluid_state_scalar.moleFraction(oilPhaseIdx, comp_idx), comp_idx);
987  primary_fluid_state.setMoleFraction(oilPhaseIdx, comp_idx, x_ii);
988  const unsigned idx = comp_idx + numComponents;
989  const auto y_ii = PrimaryEval(fluid_state_scalar.moleFraction(gasPhaseIdx, comp_idx), idx);
990  primary_fluid_state.setMoleFraction(gasPhaseIdx, comp_idx, y_ii);
991  primary_fluid_state.setKvalue(comp_idx, y_ii / x_ii);
992  }
993  PrimaryEval l;
994  l = PrimaryEval(fluid_state_scalar.L(), primary_num_pv - 1);
995  primary_fluid_state.setLvalue(l);
996  primary_fluid_state.setPressure(oilPhaseIdx, fluid_state_scalar.pressure(oilPhaseIdx));
997  primary_fluid_state.setPressure(gasPhaseIdx, fluid_state_scalar.pressure(gasPhaseIdx));
998  primary_fluid_state.setTemperature(fluid_state_scalar.temperature(0));
999 
1000  // TODO: is PrimaryFlashFluidState::Scalar> PrimaryEval here?
1001  using PrimaryParamCache = typename FluidSystem::template ParameterCache<typename PrimaryFlashFluidState::Scalar>;
1002  PrimaryParamCache primary_param_cache;
1003  for (unsigned phase_idx = 0; phase_idx < numPhases; ++phase_idx) {
1004  primary_param_cache.updatePhase(primary_fluid_state, phase_idx);
1005  for (unsigned comp_idx = 0; comp_idx < numComponents; ++comp_idx) {
1006  PrimaryEval phi = FluidSystem::fugacityCoefficient(primary_fluid_state, primary_param_cache, phase_idx, comp_idx);
1007  primary_fluid_state.setFugacityCoefficient(phase_idx, comp_idx, phi);
1008  }
1009  }
1010 
1011  using PrimaryNewtonVector = Dune::FieldVector<Scalar, num_equations>;
1012  using PrimaryNewtonMatrix = Dune::FieldMatrix<Scalar, num_equations, primary_num_pv>;
1013  PrimaryNewtonVector pri_res;
1014  PrimaryNewtonMatrix pri_jac;
1015 
1016  //use the regular equations
1017  assembleNewton_<PrimaryFlashFluidState, PrimaryComponentVector, primary_num_pv, num_equations>
1018  (primary_fluid_state, primary_z, pri_jac, pri_res);
1019 
1020  // the following code does not compile with DUNE2.6
1021  // SecondaryNewtonMatrix xx;
1022  // pri_jac.solve(xx, sec_jac);
1023  pri_jac.invert();
1024  sec_jac.template leftmultiply(pri_jac);
1025 
1026  using InputEval = typename FluidState::Scalar;
1027  using ComponentVectorMoleFraction = Dune::FieldVector<InputEval, numComponents>;
1028 
1029  ComponentVectorMoleFraction x(numComponents), y(numComponents);
1030  InputEval L_eval = L;
1031 
1032  // use the chainrule (and using partial instead of total
1033  // derivatives, DF / Dp = dF / dp + dF / ds * ds/dp.
1034  // where p is the primary variables and s the secondary variables. We then obtain
1035  // ds / dp = -inv(dF / ds)*(DF / Dp)
1036 
1037  const auto p_l = fluid_state.pressure(FluidSystem::oilPhaseIdx);
1038  const auto p_v = fluid_state.pressure(FluidSystem::gasPhaseIdx);
1039  std::vector<double> K(numComponents);
1040 
1041  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1042  K[compIdx] = fluid_state_scalar.K(compIdx);
1043  x[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::oilPhaseIdx,compIdx);//;z[compIdx] * 1. / (L + (1 - L) * K[compIdx]);
1044  y[compIdx] = fluid_state_scalar.moleFraction(FluidSystem::gasPhaseIdx,compIdx);//;x[compIdx] * K[compIdx];
1045  }
1046 
1047  // then we try to set the derivatives for x, y and L against P and x.
1048  // p_l and p_v are the same here, in the future, there might be slightly more complicated scenarios when capillary
1049  // pressure joins
1050 
1051  constexpr size_t num_deri = numComponents;
1052  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1053  std::vector<double> deri(num_deri, 0.);
1054  // derivatives from P
1055  for (unsigned idx = 0; idx < num_deri; ++idx) {
1056  deri[idx] = -sec_jac[compIdx][0] * p_l.derivative(idx);
1057  }
1058 
1059  for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1060  const double pz = -sec_jac[compIdx][cIdx + 1];
1061  const auto& zi = z[cIdx];
1062  for (unsigned idx = 0; idx < num_deri; ++idx) {
1063  deri[idx] += pz * zi.derivative(idx);
1064  }
1065  }
1066  for (unsigned idx = 0; idx < num_deri; ++idx) {
1067  x[compIdx].setDerivative(idx, deri[idx]);
1068  }
1069  // handling y
1070  for (unsigned idx = 0; idx < num_deri; ++idx) {
1071  deri[idx] = -sec_jac[compIdx + numComponents][0] * p_v.derivative(idx);
1072  }
1073  for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1074  const double pz = -sec_jac[compIdx + numComponents][cIdx + 1];
1075  const auto& zi = z[cIdx];
1076  for (unsigned idx = 0; idx < num_deri; ++idx) {
1077  deri[idx] += pz * zi.derivative(idx);
1078  }
1079  }
1080  for (unsigned idx = 0; idx < num_deri; ++idx) {
1081  y[compIdx].setDerivative(idx, deri[idx]);
1082  }
1083 
1084  // handling derivatives of L
1085  std::vector<double> deriL(num_deri, 0.);
1086  for (unsigned idx = 0; idx < num_deri; ++idx) {
1087  deriL[idx] = -sec_jac[2 * numComponents][0] * p_v.derivative(idx);
1088  }
1089  for (unsigned cIdx = 0; cIdx < numComponents; ++cIdx) {
1090  const double pz = -sec_jac[2 * numComponents][cIdx + 1];
1091  const auto& zi = z[cIdx];
1092  for (unsigned idx = 0; idx < num_deri; ++idx) {
1093  deriL[idx] += pz * zi.derivative(idx);
1094  }
1095  }
1096 
1097  for (unsigned idx = 0; idx < num_deri; ++idx) {
1098  L_eval.setDerivative(idx, deriL[idx]);
1099  }
1100  }
1101 
1102  // set up the mole fractions
1103  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1104  fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, x[compIdx]);
1105  fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, y[compIdx]);
1106  }
1107  fluid_state.setLvalue(L_eval);
1108  } //end updateDerivativesTwoPhase
1109 
1110  template <typename FlashFluidStateScalar, typename FluidState, typename ComponentVector>
1111  static void updateDerivativesSinglePhase_(const FlashFluidStateScalar& fluid_state_scalar,
1112  const ComponentVector& z,
1113  FluidState& fluid_state)
1114  {
1115  using InputEval = typename FluidState::Scalar;
1116  // L_eval is converted from a scalar, so all derivatives are zero at this point
1117  InputEval L_eval = fluid_state_scalar.L();;
1118 
1119  // for single phase situation, x = y = z;
1120  // and L_eval have all zero derivatives
1121  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
1122  fluid_state.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, z[compIdx]);
1123  fluid_state.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, z[compIdx]);
1124  }
1125  fluid_state.setLvalue(L_eval);
1126  } //end updateDerivativesSinglePhase
1127 
1128 
1129  // TODO: or use typename FlashFluidState::Scalar
1130  template <class FlashFluidState, class ComponentVector>
1131  static void successiveSubstitutionComposition_(ComponentVector& K, typename ComponentVector::field_type& L, FlashFluidState& fluid_state, const ComponentVector& z,
1132  const bool newton_afterwards, const int verbosity)
1133  {
1134  // Determine max. iterations based on if it will be used as a standalone flash or as a pre-process to Newton (or other) method.
1135  const int maxIterations = newton_afterwards ? 3 : 10;
1136 
1137  // Store cout format before manipulation
1138  std::ios_base::fmtflags f(std::cout.flags());
1139 
1140  // Print initial guess
1141  if (verbosity >= 1)
1142  std::cout << "Initial guess: K = [" << K << "] and L = " << L << std::endl;
1143 
1144  if (verbosity == 2 || verbosity == 4) {
1145  // Print header
1146  int fugWidth = (numComponents * 12)/2;
1147  int convWidth = fugWidth + 7;
1148  std::cout << std::setw(10) << "Iteration" << std::setw(fugWidth) << "fL/fV" << std::setw(convWidth) << "norm2(fL/fv-1)" << std::endl;
1149  }
1150  //
1151  // Successive substitution loop
1152  //
1153  for (int i=0; i < maxIterations; ++i){
1154  // Compute (normalized) liquid and vapor mole fractions
1155  computeLiquidVapor_(fluid_state, L, K, z);
1156 
1157  // Calculate fugacity coefficient
1158  using ParamCache = typename FluidSystem::template ParameterCache<typename FlashFluidState::Scalar>;
1159  ParamCache paramCache;
1160  for (int phaseIdx=0; phaseIdx<numPhases; ++phaseIdx){
1161  paramCache.updatePhase(fluid_state, phaseIdx);
1162  for (int compIdx=0; compIdx<numComponents; ++compIdx){
1163  auto phi = FluidSystem::fugacityCoefficient(fluid_state, paramCache, phaseIdx, compIdx);
1164  fluid_state.setFugacityCoefficient(phaseIdx, compIdx, phi);
1165  }
1166  }
1167 
1168  // Calculate fugacity ratio
1169  ComponentVector newFugRatio;
1170  ComponentVector convFugRatio;
1171  for (int compIdx=0; compIdx<numComponents; ++compIdx){
1172  newFugRatio[compIdx] = fluid_state.fugacity(oilPhaseIdx, compIdx)/fluid_state.fugacity(gasPhaseIdx, compIdx);
1173  convFugRatio[compIdx] = newFugRatio[compIdx] - 1.0;
1174  }
1175 
1176  // Print iteration info
1177  if (verbosity == 2 || verbosity == 4) {
1178  int prec = 5;
1179  int fugWidth = (prec + 3);
1180  int convWidth = prec + 9;
1181  std::cout << std::defaultfloat;
1182  std::cout << std::fixed;
1183  std::cout << std::setw(5) << i;
1184  std::cout << std::setw(fugWidth);
1185  std::cout << std::setprecision(prec);
1186  std::cout << newFugRatio;
1187  std::cout << std::scientific;
1188  std::cout << std::setw(convWidth) << convFugRatio.two_norm() << std::endl;
1189  }
1190 
1191  // Check convergence
1192  if (convFugRatio.two_norm() < 1e-6){
1193  // Restore cout format
1194  std::cout.flags(f);
1195 
1196  // Print info
1197  if (verbosity >= 1) {
1198  std::cout << "Solution converged to the following result :" << std::endl;
1199  std::cout << "x = [";
1200  for (int compIdx=0; compIdx<numComponents; ++compIdx){
1201  if (compIdx < numComponents - 1)
1202  std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx) << " ";
1203  else
1204  std::cout << fluid_state.moleFraction(oilPhaseIdx, compIdx);
1205  }
1206  std::cout << "]" << std::endl;
1207  std::cout << "y = [";
1208  for (int compIdx=0; compIdx<numComponents; ++compIdx){
1209  if (compIdx < numComponents - 1)
1210  std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx) << " ";
1211  else
1212  std::cout << fluid_state.moleFraction(gasPhaseIdx, compIdx);
1213  }
1214  std::cout << "]" << std::endl;
1215  std::cout << "K = [" << K << "]" << std::endl;
1216  std::cout << "L = " << L << std::endl;
1217  }
1218  // Restore cout format format
1219  return;
1220  }
1221 
1222  // If convergence is not met, K is updated in a successive substitution manner
1223  else{
1224  // Update K
1225  for (int compIdx=0; compIdx<numComponents; ++compIdx){
1226  K[compIdx] *= newFugRatio[compIdx];
1227  }
1228 
1229  // Solve Rachford-Rice to get L from updated K
1230  L = solveRachfordRice_g_(K, z, 0);
1231  }
1232 
1233  }
1234  if (!newton_afterwards) {
1235  throw std::runtime_error(
1236  "Successive substitution composition update did not converge within maxIterations");
1237  }
1238  }
1239 
1240 };//end PTFlash
1241 
1242 } // namespace Opm
1243 
1244 #endif
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
A central place for various physical constants occuring in some equations.
Representation of an evaluation of a function and its derivatives w.r.t.
Provides the opm-material specific exception classes.
This file contains helper classes for the material laws.
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
A number of commonly used algebraic functions for the localized OPM automatic differentiation (AD) fr...
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Implements the Peng-Robinson equation of state for a mixture.
Some templates to wrap the valgrind client request macros.
Represents all relevant thermodynamic quantities of a multi-phase, multi-component fluid system assum...
Definition: CompositionalFluidState.hpp:46
Represents a function evaluation and its derivatives w.r.t.
Definition: Evaluation.hpp:59
A generic traits class which does not provide any indices.
Definition: MaterialTraits.hpp:45
Implements a dummy linear saturation-capillary pressure relation which just disables capillary pressu...
Definition: NullMaterial.hpp:46
Determines the phase compositions, pressures and saturations given the total mass of all components f...
Definition: PTFlash.hpp:63
static void solve(FluidState &fluid_state, const Dune::FieldVector< typename FluidState::Scalar, numComponents > &z, int spatialIdx, std::string twoPhaseMethod, Scalar tolerance=-1., int verbosity=0)
Calculates the fluid state from the global mole fractions of the components and the phase pressures.
Definition: PTFlash.hpp:82
static void solve(FluidState &fluid_state, const ComponentVector &globalMolarities, Scalar tolerance=0.0)
Calculates the chemical equilibrium from the component fugacities in a phase.
Definition: PTFlash.hpp:199
Implements the Peng-Robinson equation of state for a mixture.
Definition: PengRobinsonMixture.hpp:41
static LhsEval computeFugacityCoefficient(const FluidState &fs, const Params &params, unsigned phaseIdx, unsigned compIdx)
Returns the fugacity coefficient of an individual component in the phase.
Definition: PengRobinsonMixture.hpp:89