Remove penalization in fugacity computations and improve Gibbs energy calculation #5

volpatto · 2021-05-27T18:33:22Z

RES-166

This is an intermediate step required to define a BIP function wrapper, since the result has to provide k, kT and kTT as attributes

This default values allow the user to provide, for instance, only k values without the need to define kT and kTT when there is no T dependency

Asserts to check that kT and kTT are empty are now added

…ests

To compare such results, a custom database is defined containing the hydrocarbon data. See tests/data/hydrocarbons.xml

…ests

The pip option --force-reinstall is passed in order to reinstall reaktoro python bindings in every build

Such tests are failing due to the changes in the fugacity calculations. The workaround using a penalization is now removed. (RES-166)

volpatto · 2021-05-27T18:41:57Z

.dependencies/python.devenv.yml

+  - jupyter
+  - jupyter_contrib_nbextensions
+  - jupyter_nbextensions_configurator


This should be removed.

volpatto · 2021-05-27T18:43:28Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

        ChemicalScalar Z(nspecies);
+        // Z = selectCompressibilityFactorByGibbsEnergy(nspecies, Zs, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, epsilon, sigma);


Hmm.. I don't remember why I commented here. Needs investigation.

Great ideia,
Selecting the Z factor based on the GibbsEnergy seems to be much more consistent instead of getting min values if is liquid and max if it is gas.
But I realy don't know what would happen if we were computing properties of a liquid and selectCompressibilityFactorByGibbsEnergy return the max(Z1, Z2). Especialy considering that the output will say to the user that this phase is liquid but their properties will be computed as they were a gas 🤔

It is a great idea but, since it will affect the selection of Z even for problems that has only Water and Gas, the use of selectCompressibilityFactorByGibbsEnergy broke more than 80 test :/
That made me decided to keep this line commented and make the selection of Z based on min and max of the computed Z values.

I also got some great results of test_pedersen63.py by making PhaseIdentificationMethod = PhaseIdentificationMethod:None which makes the solver not penalize any phase

alexandrebbruno

Nicely done,
I only added some minor comments which could add more value to the work

Reaktoro/Thermodynamics/EOS/CubicEOS.hpp

Reaktoro/Thermodynamics/EOS/PhaseIdentification.cpp

...ibrium_liq_gas/test_equilibrium_CH4_H2S_liq_gas_temperature_equal_273_15_K_and_30_0_bar_.csv

...ibrium_liq_gas/test_equilibrium_CH4_H2S_liq_gas_temperature_equal_293_15_K_and_70_0_bar_.csv

alexandrebbruno · 2021-06-01T14:12:35Z

tests/regression/test_pedersen63.py

@@ -0,0 +1,380 @@
+import numpy as np


Great job,
I think that would be nice to add the reference used to make this test. Maybe as documentation in one of the tests

alexandrebbruno · 2021-06-01T14:12:48Z

tests/regression/test_whitson18.py

@@ -0,0 +1,375 @@
+import numpy as np


tests/test_oil_equilibrium.py

alexandrebbruno · 2021-06-01T14:39:17Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

        ChemicalScalar Z(nspecies);
+        // Z = selectCompressibilityFactorByGibbsEnergy(nspecies, Zs, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, epsilon, sigma);


Great ideia,
Selecting the Z factor based on the GibbsEnergy seems to be much more consistent instead of getting min values if is liquid and max if it is gas.
But I realy don't know what would happen if we were computing properties of a liquid and selectCompressibilityFactorByGibbsEnergy return the max(Z1, Z2). Especialy considering that the output will say to the user that this phase is liquid but their properties will be computed as they were a gas 🤔

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

Reaktoro/Thermodynamics/EOS/PhaseIdentification.cpp

(RES-166)

volpatto · 2021-07-01T17:13:09Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

-            Z.val = *std::max_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
-        else
-            Z.val = *std::min_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
+         Z = selectCompressibilityFactorByGibbsEnergy(nspecies, Zs, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, epsilon, sigma);


Please note that I have enabled this part, so more tests are now failing.

(RES-166)

volpatto · 2021-07-01T19:44:05Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+        // When using selectCompressibilityFactorByGibbsEnergy, the switch-case below
+        // is not necessary.
+        // TODO: remove switch-case if selectCompressibilityFactorByGibbsEnergy is used
+        auto identified_phase_type = input_phase_type;


Please have a look on this. The cubic EoS model is attributed to a phase. When calculating Z for such phase, a cubic equation is solved and the root is selected according to Gibbs Energy criterion, so there is no need to identify the phase here anymore. For thermodynamic consistency, the Z related to the minimum Gibbs energy is always selected despite of the phase.

volpatto · 2021-07-01T19:45:30Z

Nicely done,
I only added some minor comments which could add more value to the work

Thanks!

…:ESSS/reaktoro into fb-RES-166-improve-cubiceos-gibbs-energy

…mparing floats

OtavioPellicano · 2021-07-06T16:47:14Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

-            result.ln_fugacity_coefficients.fill(100.0);
-            return result;
-        }
+		if (identified_phase_type != input_phase_type) {


It seems there is more indentation than is needed.

It's because spaces were converted to tabs (probably because of editor configuration). It seems that these lines 524-538 could be reset.

OtavioPellicano · 2021-07-06T17:42:21Z

Reaktoro/Core/ChemicalProperties.cpp

@@ -321,6 +321,13 @@ auto ChemicalProperties::phaseSpecificVolumes() const -> ChemicalVector
    return phaseAmounts()/phaseMasses() % phaseMolarVolumes();
 }

+auto ChemicalProperties::phaseCompressibilityFactors() const -> ChemicalVector
+{
+    const auto& R = universalGasConstant;


I am curious. Why did you use const auto& since universalGasConstant is just a const double? Maybe const auto fit better here.

OtavioPellicano · 2021-07-06T17:45:15Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+{
+    const double R = universalGasConstant;
+
+	const double almost_zero = 1e-20;


Fix indentation

OtavioPellicano · 2021-07-06T17:53:34Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+
+    // Calculate the integration factor I
+    ChemicalScalar I;
+    if(std::abs(epsilon - sigma) > almost_zero) I = log((Z + sigma*beta)/(Z + epsilon*beta))/(sigma - epsilon);


IMHO, assignments like this would be more legible if you use a ternary operator, or even curly brackets in each block.

OtavioPellicano · 2021-07-06T17:55:08Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+
+    // Calculate the integration factor I for each component
+    ChemicalScalar Ii;
+    if(std::abs(epsilon - sigma) > almost_zero) Ii = I + ((Zi + sigma*betai)/(Z + sigma*beta) - (Zi + epsilon*betai)/(Z + epsilon*beta))/(sigma - epsilon);


OtavioPellicano · 2021-07-06T17:57:05Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+    {
+        auto xi = x[i];
+        auto fi = computeSpeciesFugacity(P, T, xi, abar[i], bbar[i], abarT[i], amix, amixT, bmix, A, B, C, Z, epsilon, sigma);
+		G_normalized += xi.val * log(fi).val;


Fix indentation

OtavioPellicano · 2021-07-06T18:01:14Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

            Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[0]));
        }
+        else if (cubic_size == 2)
+        {
+            Zs.push_back(ChemicalScalar(nspecies, cubicEOS_roots[0]));


Avoid temporary objects, use emplace_back instead of push_back

OtavioPellicano · 2021-07-06T18:03:46Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

        ChemicalScalar Z(nspecies);
+
+		// Selecting compressibility factor - Z_liq < Z_gas
+		//TODO: study the possibilty to use selectCompressibilityFactorByGibbsEnergy for Z selection


Normally TODO must be followed by the Task number

OtavioPellicano · 2021-07-06T18:05:12Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

        if (isvapor)
            Z.val = *std::max_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
        else
-            Z.val = *std::min_element(cubicEOS_roots.begin(), cubicEOS_roots.end());
+			Z.val = *std::min_element(cubicEOS_roots.begin(), cubicEOS_roots.end());


Fix indentation

OtavioPellicano · 2021-07-06T18:25:34Z

tests/regression/test_equilibrium_liq_gas.py


    solver.setOptions(options)

    state = ChemicalState(system)
+    state.setSpeciesAmounts(0.001)  # start will all having 0.001 moles
+    state.setSpeciesAmount("CH4(g)", overall_composition[0])  # overwrite amount of C1(g) (same below)


I didn't get C1(g). That means Carbon dissociated? If so, the dissociation should be C4+, ins't it?

tadeu · 2021-07-06T22:17:46Z

.dependencies/python.devenv.yml

@@ -3,6 +3,7 @@
 dependencies:
  - numpy
  - pandas
+  - matplotlib


Perhaps this one could also be removed?

tadeu · 2021-07-06T22:22:31Z

tests/regression/test_equilibrium_liq_gas.py

    solver.setOptions(options)

    state = ChemicalState(system)
+    state.setSpeciesAmounts(0.001)  # start will all having 0.001 moles


will -> with?

tadeu · 2021-07-06T22:25:07Z

Reaktoro/Thermodynamics/EOS/PhaseIdentification.cpp

+    auto const& T = temperature;
+
+    // Computing the values of residual Gibbs energy for all Zs
+    ChemicalScalar Gs;


Minor, but it doesn't seem necessary to have this uninitialized variable here, line 138 could be ChemicalScalar Gs = R * temperature*(Z - 1 - log(Z - beta) - q * I); or line 140 could be return R * temperature*(Z - 1 - log(Z - beta) - q * I).

tadeu · 2021-07-06T22:27:18Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

+    return G_normalized;
+}
+
+auto selectCompressibilityFactorByGibbsEnergy(


Is this function unused?

tadeu · 2021-07-06T22:30:22Z

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp

-            result.ln_fugacity_coefficients.fill(100.0);
-            return result;
-        }
+		if (identified_phase_type != input_phase_type) {


It's because spaces were converted to tabs (probably because of editor configuration). It seems that these lines 524-538 could be reset.

volpatto and others added 30 commits January 29, 2021 16:46

Separate CubicEOS python bindings in a proper file

afea63e

Remove redundant includes in PyCubicEos.cpp

bb7ed6b

Rename CubicEos -> CubicEOS

f64d2fb

Expose an interface to define BIPs values

6b6af9e

This is an intermediate step required to define a BIP function wrapper, since the result has to provide k, kT and kTT as attributes

Add some docs toInteractionParamsResult

5310496

Expose function that calculate BIPs through CubicEOSParams

a71d8ca

Use MatrixXd as default initialization in BinaryInteractionParams

7d5fbf8

This default values allow the user to provide, for instance, only k values without the need to define kT and kTT when there is no T dependency

Merge branch 'add-hydrocarbon-tests' into hydrocarbon-mixtures-tests

4cbf627

Improve BIPs asserts and add docs for the tests

5dd15d4

Add a sanity check for BIPs setup in CubicEOS

96cfee3

Add the BIPs sanity check and improve docs

60f9206

Improve test_bips_setup_without_derivatives

adff826

Asserts to check that kT and kTT are empty are now added

Remove leftovers in test_cubiceos.py

837d5e7

Merge branch 'refactor-cubiceos-bindings' into hydrocarbon-mixtures-t…

12ee1e8

…ests

Fix BIPs sanity check logic

b0fc3ae

Add test comparing results with and without BIPs setup

45c92b4

To compare such results, a custom database is defined containing the hydrocarbon data. See tests/data/hydrocarbons.xml

Fix database parses for Liquid type species

39a35ef

Add HC equilibrium test cases with and without BIPs

8215c84

Add forgotten assert test_equilibrium_CH4_CO2_pedersen

951b4e6

Merge branch 'refactor-cubiceos-bindings' into hydrocarbon-mixtures-t…

b464378

…ests

Add ex C1-CO2 from Pedersen, table 6.3

fba8e38

Update python deps to allow pandas to read xlsx

8466980

Fix activity indicies usage

8b22d0d

Add phase composition comparison for pedersen 63 ex

8e0e84e

Add plots improvements to the script of the cases

96da719

Remove unnecessary tracked file

5ad0328

Merge branch 'master' into hydrocarbon-mixtures-tests

fd4eb12

Merge branch 'master' into hydrocarbon-mixtures-tests

4919ff2

Update CubicEOS with BIPs sanitizer from master

2f18f46

Add a TODO for commented code in CubicEOS.cpp

b23efeb

Diego Tavares Volpatto added 4 commits March 9, 2021 22:05

Remove hanging "private" in CubicEOS Impl

0e29074

Improve Reaktoro python install while developing

4cc33af

The pip option --force-reinstall is passed in order to reinstall reaktoro python bindings in every build

Merge branch 'master' into fb-RES-166-improve-cubiceos-gibbs-energy

a08d878

Put xfail mark in tests with liq phase which are failing

669c193

Such tests are failing due to the changes in the fugacity calculations. The workaround using a penalization is now removed. (RES-166)

volpatto requested review from tadeu, alexandrebbruno, jaimeambrus and OtavioPellicano May 27, 2021 18:33

volpatto commented May 27, 2021

View reviewed changes

alexandrebbruno approved these changes Jun 1, 2021

View reviewed changes

volpatto marked this pull request as draft June 14, 2021 01:22

jaimeambrus reviewed Jun 14, 2021

View reviewed changes

Reaktoro/Thermodynamics/EOS/CubicEOS.cpp Outdated Show resolved Hide resolved

jaimeambrus reviewed Jun 14, 2021

View reviewed changes

Reaktoro/Thermodynamics/EOS/PhaseIdentification.cpp Outdated Show resolved Hide resolved

Enable Z selection by Gibbs energy

1f91ad0

(RES-166)

volpatto commented Jul 1, 2021

View reviewed changes

Add comment about phase identification switch-case

d3a73fd

(RES-166)

volpatto commented Jul 1, 2021

View reviewed changes

alexandrebbruno changed the title ~~[WIP] Remove penalization in fugacity computations and improve Gibbs energy calculation~~ Remove penalization in fugacity computations and improve Gibbs energy calculation Jul 5, 2021

alexandrebbruno added 7 commits July 5, 2021 18:12

remove unnecessary dependencies

dcb2e8e

removing xfail and update results

d64160e

selecting Z based on min and max and using penalty

35a34a7

removing phase_identification_method to set default value

ff05efb

Merge branch 'fb-RES-166-improve-cubiceos-gibbs-energy' of github.com…

f5db28c

…:ESSS/reaktoro into fb-RES-166-improve-cubiceos-gibbs-energy

remove unnecessary comment

ff3c7af

absolute difference is bigger than a threshold instead of directly co…

c7a595b

…mparing floats

OtavioPellicano approved these changes Jul 6, 2021

View reviewed changes

tadeu reviewed Jul 6, 2021

View reviewed changes

fix indentation

a2a6b0b

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Remove penalization in fugacity computations and improve Gibbs energy calculation #5

Remove penalization in fugacity computations and improve Gibbs energy calculation #5

volpatto commented May 27, 2021

volpatto May 27, 2021

alexandrebbruno Jul 5, 2021

volpatto May 27, 2021

alexandrebbruno Jun 1, 2021

alexandrebbruno Jul 5, 2021

alexandrebbruno left a comment

alexandrebbruno Jun 1, 2021

alexandrebbruno Jun 1, 2021

alexandrebbruno Jun 1, 2021

volpatto Jul 1, 2021

volpatto Jul 1, 2021

volpatto commented Jul 1, 2021

OtavioPellicano Jul 6, 2021

tadeu Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021 •

edited

Loading

OtavioPellicano Jul 6, 2021

OtavioPellicano Jul 6, 2021

tadeu Jul 6, 2021

tadeu Jul 6, 2021

tadeu Jul 6, 2021

tadeu Jul 6, 2021

tadeu Jul 6, 2021

		ChemicalScalar Z(nspecies);
		// Z = selectCompressibilityFactorByGibbsEnergy(nspecies, Zs, P, T, x, abar, bbar, abarT, amix, amixT, bmix, A, B, C, epsilon, sigma);

Remove penalization in fugacity computations and improve Gibbs energy calculation #5

Are you sure you want to change the base?

Remove penalization in fugacity computations and improve Gibbs energy calculation #5

Conversation

volpatto commented May 27, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

alexandrebbruno left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

volpatto commented Jul 1, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

OtavioPellicano Jul 6, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

OtavioPellicano Jul 6, 2021 •

edited

Loading