Implementation of minimization of tpd by using a second order newton

c_interface/yaeos_c.f90

@@ -339,12 +339,13 @@ subroutine set_qmrtd(ar_id, kij_0, kij_inf, t_star, lij)
       end associate
    end subroutine set_qmrtd
 
-   subroutine set_mhv(ar_id, ge_id, q)
+   subroutine set_mhv(ar_id, ge_id, q, lij)
       !! Michelsen's Modified Huron-Vidal 1 with constant `q_1` parameter
       use yaeos, only: MHV, CubicEoS
       integer(c_int), intent(in) :: ar_id
       integer(c_int), intent(in) :: ge_id
       real(c_double), intent(in) :: q
+      real(c_double), intent(in) :: lij(:, :)
 
       type(MHV) :: mixrule
 
@@ -353,7 +354,7 @@ subroutine set_mhv(ar_id, ge_id, q)
 
       select type(ar_model)
        class is(CubicEoS)
-         mixrule = MHV(ge=ge_model, b=ar_model%b, q=q)
+         mixrule = MHV(ge=ge_model, b=ar_model%b, q=q, lij=lij)
          deallocate(ar_model%mixrule)
          ar_model%mixrule = mixrule
       end select
@@ -1447,6 +1448,8 @@ subroutine generalized_isopleth(&
       ws_stab, &
       x_ls, ws, betas, Ts, Ps, w_more_stable, found_unstability, max_points &
       )
+      !! Calculates a generalized isopleth line for a given model and conditions.
+      !! This subroutine computes the generalized isopleth line for a given model and conditions,
       use yaeos, only: GeneralizedIsoZLine, create_generalized_isoz_line
       integer(c_int), intent(in) :: id
       integer(c_int), intent(in) :: nc
@@ -1538,11 +1541,12 @@ elemental subroutine convert_kind(numeric, string)
    end subroutine convert_kind
 
    subroutine convert_gerg_components(string, numeric)
+      !! Converts a string representation of a component to its numeric identifier
       use yaeos, only: G2008Components
       use iso_fortran_env, only: error_unit
 
-      character(len=*), intent(in) :: string
-      integer(c_int), intent(out) :: numeric
+      character(len=*), intent(in) :: string !! String representation of the component
+      integer(c_int), intent(out) :: numeric !! Numeric identifier for the component
 
       select case (string)
        case ("methane", "C1", "CH4", "nC1")

fpm.toml

@@ -27,7 +27,7 @@ json-fortran = {git="https://github.com/jacobwilliams/json-fortran"}
 forsus = {git="https://github.com/ipqa-research/forsus"}
 fortime = { git = "https://github.com/gha3mi/fortime.git" }
 PowellOpt = { git = "https://github.com/jacobwilliams/PowellOpt/" }
-#minpack = {git = "https://github.com/fortran-lang/minpack" }
+minpack = {git = "https://github.com/fortran-lang/minpack" }
 finterp = { git = "https://github.com/jacobwilliams/finterp.git" }
 
 [dev-dependencies]

python/docs/source/tutorial/instantiate_model.ipynb

@@ -20,17 +20,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "\u001b[0;31mType:\u001b[0m        module\n",
-      "\u001b[0;31mString form:\u001b[0m <module 'yaeos.models' from '/home/salvador/code/yaeos/python/yaeos/models/__init__.py'>\n",
-      "\u001b[0;31mFile:\u001b[0m        ~/code/yaeos/python/yaeos/models/__init__.py\n",
-      "\u001b[0;31mDocstring:\u001b[0m  \n",
+      "\u001b[31mType:\u001b[39m        module\n",
+      "\u001b[31mString form:\u001b[39m <module 'yaeos.models' from '/home/ruther/docs/programming/python/virtualenvs/thermo/lib/python3.13/site-packages/yaeos/models/__init__.py'>\n",
+      "\u001b[31mFile:\u001b[39m        ~/docs/programming/python/virtualenvs/thermo/lib/python3.13/site-packages/yaeos/models/__init__.py\n",
+      "\u001b[31mDocstring:\u001b[39m  \n",
       "Models module.\n",
       "\n",
       "Yaeos models module. This module provides the following submodules:\n",
@@ -48,14 +48,18 @@
       "        - RKPR: RKPR model\n",
       "    - Mixing rules: mixing rules for cubic EoS\n",
       "        - QMR: cuadratic mixing rule\n",
-      "        - MHV: modified Huron-Vidal mixing rule"
+      "        - MHV: modified Huron-Vidal mixing rule\n",
+      "        - HV: Huron-Vidal mixing rule\n",
+      "        - HVNRTL: Huron-Vidal NRTL mixing rule\n",
+      "    - Multifluid EoS:\n",
+      "        - GERG2008: GERG2008 Residual contribution"
      ]
     }
    ],
    "source": [
     "import yaeos\n",
     "\n",
-    "?yaeos.models"
+    "yaeos.models?"
    ]
   },
   {
@@ -911,7 +915,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "yaeos",
+   "display_name": "thermo",
    "language": "python",
    "name": "python3"
   },
@@ -925,7 +929,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.12"
+   "version": "3.13.5"
   }
  },
  "nbformat": 4,

python/meson.build

@@ -35,6 +35,7 @@ incdir_f2py = run_command(py,
 
 
 # fpm
+fflags = '-g -fPIC -funroll-loops -fstack-arrays -Ofast -frepack-arrays -faggressive-function-elimination -fopenmp'
 if get_option('debug')
   fflags = '-g -fPIC -Og'
 else

python/yaeos/core.py

@@ -4,7 +4,7 @@
 the models' thermoprops methods.
 """
 
-from abc import ABC
+from abc import ABC, abstractmethod
 from typing import Union
 from warnings import warn
 
@@ -44,22 +44,33 @@ def adjust_root_kind(number_of_phases, kinds_x=None, kind_w=None):
         kind_w_out : str
             Kind of the test phase
     """
-    if kinds_x:
-        kinds_x_out = [root_kinds[kind] for kind in kinds_x]
-    else:
+    if kinds_x is None:
         kinds_x_out = [root_kinds["stable"] for _ in range(number_of_phases)]
-
-    if kind_w:
-        kind_w_out = root_kinds[kind_w]
     else:
+        kinds_x_out = [root_kinds[kind] for kind in kinds_x]
+
+    if kind_w is None:
         kind_w_out = root_kinds["stable"]
+    else:
+        kind_w_out = root_kinds[kind_w]
 
     return kinds_x_out, kind_w_out
 
 
 class GeModel(ABC):
     """Excess Gibbs (Ge) model abstract class."""
 
+    @abstractmethod
+    def size(self) -> int:
+        """Get number of components in the model.
+
+        Raises
+        ------
+        NotImplementedError
+            Abstract error, this method must be implemented in the subclass
+        """
+        raise NotImplementedError
+
     def ln_gamma(
         self, moles, temperature: float, dt: bool = False, dn: bool = False
     ) -> Union[np.ndarray, tuple[np.ndarray, dict]]:
@@ -421,6 +432,17 @@ def __del__(self) -> None:
 class ArModel(ABC):
     """Residual Helmholtz (Ar) model abstract class."""
 
+    @abstractmethod
+    def size(self) -> int:
+        """Get number of components in the model.
+
+        Raises
+        ------
+        NotImplementedError
+            Abstract error, this method must be implemented in the subclass
+        """
+        raise NotImplementedError
+
     def lnphi_vt(
         self,
         moles,
@@ -2655,6 +2677,53 @@ def generalized_isopleth(
             "found_unstability": found_unstability,
         }
 
+    def precipitation_line_from_env(
+            self, env: PTEnvelope, spec_value, spec="P", ds0=1e-5,
+            ws_stability=None, max_points=100
+            ) -> dict:
+        """Calculate precipitation line from a PTEnvelope.
+        """
+        phases = env.number_of_phases
+        z = env.global_composition
+
+        x_l0 = env["x"][-1]
+        w0 = env["w"][-1]
+
+        betas0 = [*env.main_phases_molar_fractions[-1], 1e-15]
+
+        if spec == "P":
+            spec_variable = phases * len(z) + phases + 1 + 1
+            p0 = spec_value
+            spec_variable_value = np.log(p0)
+            loc = np.argmin(np.abs(env["P"] - p0))
+            t0 = env["T"][loc]
+        elif spec == "T":
+            spec_variable = phases * len(z) + phases + 1 + 2
+            t0 = spec_value
+            spec_variable_value = np.log(t0)
+            loc = np.argmin(np.abs(env["T"] - t0))
+            p0 = env["P"][loc]
+        else:
+            raise ValueError(
+                "spec must be either 'P' or 'T', got: {}".format(spec)
+            )
+
+        ns0 = phases * len(z) + phases + 1
+        s0 = betas0[-1]
+        kinds_x = env.main_phases_kinds[loc, :]
+        kind_w = env.reference_phase_kinds[loc]
+
+        prec = self.generalized_isopleth(
+            z, kinds_x=kinds_x, kind_w=kind_w,
+            spec_variable=spec_variable,
+            spec_variable_value=spec_variable_value,
+            x_l0=x_l0, w0=w0, betas0=betas0, t0=t0, p0=p0,
+            ns0=ns0, s0=s0, ds0=ds0, max_points=max_points,
+            ws_stability=ws_stability
+        )
+
+        return prec
+
     # =========================================================================
     # Stability analysis
     # -------------------------------------------------------------------------

python/yaeos/envelopes.py

@@ -398,6 +398,8 @@ def __getitem__(self, key):
                 alphas=self.alphas[key],
                 iterations=self.iterations[key],
                 specified_variable=self.specified_variable[key],
+                main_phases_kinds=self.main_phases_kinds[key],
+                reference_phase_kinds=self.reference_phase_kinds[key],
                 critical_pressures=self.critical_pressures,
                 critical_alphas=self.critical_alphas,
             )

python/yaeos/models/excess_gibbs/dortmund.py

@@ -37,9 +37,21 @@ class UNIFACDortmund(GeModel):
 
     def __init__(self, molecules: List[dict]) -> None:
 
+        self.molecules = molecules
+
         (number_of_groups, groups_ids, groups_ammounts) = groups_from_dicts(
             molecules
         )
         self.id = yaeos_c.unifac_dortmund(
             ngs=number_of_groups, g_ids=groups_ids, g_v=groups_ammounts
         )
+
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return len(self.molecules)

python/yaeos/models/excess_gibbs/nrtl.py

@@ -1,5 +1,7 @@
 """Non-random two-liquid model (NRTL) module."""
 
+import numpy as np
+
 from yaeos.core import GeModel
 from yaeos.lib import yaeos_c
 
@@ -43,7 +45,17 @@ class NRTL(GeModel):
     """
 
     def __init__(self, a, b, c) -> None:
-        self.a = a
-        self.b = b
-        self.c = c
+        self.a = np.array(a, order="F")
+        self.b = np.array(b, order="F")
+        self.c = np.array(c, order="F")
         self.id = yaeos_c.nrtl(a, b, c)
+
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return self.a.shape[0]

python/yaeos/models/excess_gibbs/psrk_unifac.py

@@ -37,9 +37,21 @@ class UNIFACPSRK(GeModel):
 
     def __init__(self, molecules: List[dict]) -> None:
 
+        self.molecules = molecules
+
         (number_of_groups, groups_ids, groups_ammounts) = groups_from_dicts(
             molecules
         )
         self.id = yaeos_c.unifac_psrk(
             ngs=number_of_groups, g_ids=groups_ids, g_v=groups_ammounts
         )
+
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return len(self.molecules)

python/yaeos/models/excess_gibbs/unifac.py

@@ -37,9 +37,21 @@ class UNIFACVLE(GeModel):
 
     def __init__(self, molecules: List[dict]) -> None:
 
+        self.molecules = molecules
+
         (number_of_groups, groups_ids, groups_ammounts) = groups_from_dicts(
             molecules
         )
         self.id = yaeos_c.unifac_vle(
             ngs=number_of_groups, g_ids=groups_ids, g_v=groups_ammounts
         )
+
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return len(self.molecules)

python/yaeos/models/excess_gibbs/uniquac.py

@@ -115,3 +115,14 @@ def __init__(
         self.id = yaeos_c.uniquac(
             qs, rs, self.aij, self.bij, self.cij, self.dij, self.eij
         )
+        self.nc = nc
+
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return self.nc

python/yaeos/models/residual_helmholtz/cubic_eos/cubic_eos.py

@@ -42,6 +42,16 @@ def __init__(
         self.pc = critical_pressures
         self.w = acentric_factors
 
+    def size(self) -> int:
+        """Get the number of components.
+
+        Returns
+        -------
+        int
+            Number of components
+        """
+        return self.nc
+
     def set_mixrule(self, mixrule: CubicMixRule) -> None:
         """Set the mixing rule for the EoS.