From 76ad99e4581c2932cd77817785351ba4b265b359 Mon Sep 17 00:00:00 2001
From: Kai Bao <kai.bao@sintef.no>
Date: Tue, 22 Mar 2022 14:26:47 +0100
Subject: [PATCH] using different forms of flash equations have better
 convergence

it is the fugacity coefficient equations. It begins with bigger
residual, while converges faster.

origionally,
 newton residual is 24.7536
 newton residual is 11.0045
 newton residual is 5.2655
 newton residual is 2.4048
 newton residual is 0.995777
 newton residual is 0.331732
 newton residual is 0.0663179
 newton residual is 0.00393044
 newton residual is 1.84752e-05
 newton residual is 1.55499e-08
 newton residual is 1.67839e-11

new form,
 newton residual is 694053
 newton residual is 5950.39
 newton residual is 52.1553
 newton residual is 0.470124
 newton residual is 0.00430846
 newton residual is 3.91461e-05
 newton residual is 3.55269e-07
 newton residual is 4.00957e-09
---
 opm/material/constraintsolvers/ChiFlash.hpp                | 7 +++++--
 .../chifluid/juliathreecomponentfluidsystem.hh             | 2 ++
 2 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/opm/material/constraintsolvers/ChiFlash.hpp b/opm/material/constraintsolvers/ChiFlash.hpp
index 8cc954537..77d691cd6 100644
--- a/opm/material/constraintsolvers/ChiFlash.hpp
+++ b/opm/material/constraintsolvers/ChiFlash.hpp
@@ -931,8 +931,11 @@ protected:
 
             {
                 // (f_liquid/f_vapor) - 1 = 0
-                auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) /
-                                  fluid_state.fugacity(gasPhaseIdx, compIdx)) - 1.0;
+                /* auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) /
+                                  fluid_state.fugacity(gasPhaseIdx, compIdx)) - 1.0; */
+                // TODO: it looks this formulation converges quicker while begin with bigger residual
+                auto local_res = (fluid_state.fugacity(oilPhaseIdx, compIdx) -
+                                  fluid_state.fugacity(gasPhaseIdx, compIdx));
                 res[compIdx + numComponents] = Opm::getValue(local_res);
                 for (unsigned i = 0; i < num_equation; ++i) {
                     jac[compIdx + numComponents][i] = local_res.derivative(i);
diff --git a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh
index 32b625972..fd9a04474 100644
--- a/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh
+++ b/opm/material/fluidsystems/chifluid/juliathreecomponentfluidsystem.hh
@@ -158,6 +158,8 @@ namespace Opm {
             assert(0 <= phaseIdx && phaseIdx < numPhases);
             assert(0 <= compIdx && compIdx < numComponents);
 
+            // TODO: here the derivatives for the phi are dropped. Should we keep the derivatives against the pressure
+            // and temperature?
             Scalar phi = Opm::getValue(
                     PengRobinsonMixture::computeFugacityCoefficient(fluidState, paramCache, phaseIdx, compIdx));
             return phi;