Add phase equilibrium

.github/workflows/release.yml

@@ -34,7 +34,7 @@ jobs:
   build-macos:
     name: Build macOS Wheels
     if: ${{ github.event.inputs.build_os == 'Both' || github.event.inputs.build_os == 'MacOS' }}
-    runs-on: macos-latest
+    runs-on: macos-15
     strategy:
       fail-fast: true
       matrix:

.gitignore

@@ -81,3 +81,6 @@ IMPLICIT_INTEGRATION_GAP_ANALYSIS.md
 ION_CHEMISTRY_SUPPORT_PLAN.md
 KINETICS_BASE_TITAN_GAP_STATUS.md
 python/kinetics_base/titan/ARCHITECTURE.md
+
+# pdfs
+*.pdf

README.md

@@ -23,6 +23,34 @@ KINTERA provides efficient implementations of:
 - Phase equilibrium computations
 - Atmospheric chemistry models
 
+### Multiphase Equilibrium
+
+`EquilibriumTP` is a fixed-temperature, fixed-pressure constrained chemistry
+solver. The C++/CUDA core accepts component moles and precomputed logarithmic
+equilibrium constants; the module derives phase membership and stoichiometry
+from its options. Case-specific thermodynamics remains in Python under
+`kintera.equilibrium`.
+
+Equilibrium networks use the repository's top-level `phases`, `species`, and
+`reactions` YAML layout. Phase species determine component ordering, species
+compositions validate elemental balance, and reactions with `type: equilibrium`
+generate the module's stoichiometric buffer:
+
+```python
+from kintera import EquilibriumOptions, EquilibriumTP
+
+options = EquilibriumOptions.from_yaml("equilibrium.yaml")
+solver = EquilibriumTP(options)
+```
+
+`Nasa9LogK` evaluates ideal-gas equilibrium constants from the bundled NASA-9
+database. See `examples/equilibrium_nasa9.yaml` and
+`examples/equilibrium_nasa9.py` for a complete YAML-defined sample:
+
+```bash
+python examples/equilibrium_nasa9.py
+```
+
 The library is written in C++17 with Python bindings, leveraging PyTorch for tensor operations and providing GPU acceleration support via CUDA.
 
 ## Features

examples/equilibrium_nasa9.py

@@ -0,0 +1,31 @@
+"""Gas equilibrium sample using the standard NASA-9 database."""
+
+from pathlib import Path
+
+import torch
+from kintera import EquilibriumOptions, EquilibriumTP
+from kintera.equilibrium import Nasa9LogK
+
+
+def main() -> None:
+    config = Path(__file__).with_suffix(".yaml")
+    options = EquilibriumOptions.from_yaml(str(config))
+    species = options.components()
+    solver = EquilibriumTP(options)
+    log_k_model = Nasa9LogK(options)
+
+    temp = torch.tensor(4000.0, dtype=torch.float64)
+    pressure = torch.tensor(1.0e5, dtype=torch.float64)
+    initial_moles = torch.tensor([0.60, 0.20, 0.20], dtype=torch.float64)
+    log_k = log_k_model(temp, pressure)
+    moles, _, diagnostics = solver(temp, pressure, initial_moles, log_k)
+
+    print(f"T = {temp.item():.0f} K, P = {pressure.item():.0f} Pa")
+    print("log(K):", dict(zip(options.reactions(), log_k.tolist())))
+    print("initial moles:", dict(zip(species, initial_moles.tolist())))
+    print("equilibrium moles:", dict(zip(species, moles.tolist())))
+    print("diagnostics [status, iterations, error]:", diagnostics.tolist())
+
+
+if __name__ == "__main__":
+    main()

examples/equilibrium_nasa9.yaml

@@ -0,0 +1,24 @@
+phases:
+  - name: gas
+    thermo: ideal-gas
+    species: [H2, O2, H2O]
+
+species:
+  - name: H2
+    composition: {H: 2}
+
+  - name: O2
+    composition: {O: 2}
+
+  - name: H2O
+    composition: {H: 2, O: 1}
+
+reactions:
+  - equation: "H2 + 0.5 O2 <=> H2O"
+    type: equilibrium
+
+equilibrium:
+  standard-pressure: 1.0e5
+  max-iter: 100
+  ftol: 1.0e-8
+  mole-floor: 1.0e-30

pyproject.toml

@@ -54,6 +54,7 @@ packages = [
   "kintera.atm2d.conservation",
   "kintera.atm2d.newton",
   "kintera.atm2d.sources",
+  "kintera.equilibrium",
   "kintera.kinetics_base",
   "kintera.kinetics_base.titan",
   "kintera.kinetics_base_titan",
@@ -64,6 +65,7 @@ include-package-data = false
 "kintera" = [
   "kintera/**/*",
   "data/*.dat",
+  "equilibrium/*.yaml",
   "lib/*.so",
   "lib/*.dylib",
   "*.pyi",

python/__init__.py

@@ -18,6 +18,7 @@ def _add_packaged_resource_directory() -> None:
 _add_packaged_resource_directory()
 
 from .atm2d import *
+from .equilibrium import *
 from .kinetics_base_titan import *
 
 try:

python/csrc/kintera.cpp

@@ -11,6 +11,7 @@
 #include <kintera/photochem/kinetics_base_reader.hpp>
 #include <kintera/species.hpp>
 #include <kintera/utils/find_resource.hpp>
+#include <kintera/utils/molar_mass.hpp>
 
 // python
 #include "pyoptions.hpp"
@@ -27,24 +28,34 @@ extern std::vector<double> species_sref_R;
 
 }  // namespace kintera
 
-void bind_thermo(py::module& m);
-void bind_constants(py::module& m);
-void bind_kinetics(py::module& m);
-void bind_photolysis(py::module& m);
-void bind_diffusion(py::module& m);
-void bind_sparse_solver(py::module& m);
+void bind_thermo(py::module &m);
+void bind_constants(py::module &m);
+void bind_kinetics(py::module &m);
+void bind_photolysis(py::module &m);
+void bind_diffusion(py::module &m);
+void bind_sparse_solver(py::module &m);
+void bind_equilibrium(py::module &m);
 
 PYBIND11_MODULE(kintera, m) {
   m.attr("__name__") = "kintera";
   m.doc() = R"(Atmospheric Thermodynamics and Chemistry Library)";
 
+  m.def("atomic_mass", &kintera::atomic_mass, py::arg("element"));
+  m.def("molar_mass", &kintera::molar_mass, py::arg("composition"));
+  m.def("molar_masses", &kintera::molar_masses, py::arg("elements"),
+        py::arg("element_matrix"));
+  m.def("molar_masses_from_yaml", &kintera::molar_masses_from_yaml,
+        py::arg("filename"));
+  m.def("nasa9_gibbs_rt", &kintera::nasa9_gibbs_rt, py::arg("temp"),
+        py::arg("species"));
+
   auto pySpeciesThermo =
       py::class_<kintera::SpeciesThermoImpl, kintera::SpeciesThermo>(
           m, "SpeciesThermo");
 
   pySpeciesThermo.def(py::init<>())
       .def("__repr__",
-           [](const kintera::SpeciesThermo& self) {
+           [](const kintera::SpeciesThermo &self) {
              return fmt::format("SpeciesThermo({})", self);
            })
       .def("species", &kintera::SpeciesThermoImpl::species)
@@ -59,9 +70,9 @@ PYBIND11_MODULE(kintera, m) {
   auto pyReaction = py::class_<kintera::Reaction>(m, "Reaction");
 
   pyReaction.def(py::init<>())
-      .def(py::init<const std::string&>())
+      .def(py::init<const std::string &>())
       .def("__repr__",
-           [](const kintera::Reaction& self) {
+           [](const kintera::Reaction &self) {
              return fmt::format("Reaction({})", self);
            })
       .def("equation", &kintera::Reaction::equation)
@@ -209,7 +220,7 @@ PYBIND11_MODULE(kintera, m) {
                      &kintera::KBTitanReactionReport::unsupported_reaction_ids);
 
   m.def("classify_kinetics_base_titan_reactions",
-        py::overload_cast<std::string const&, std::string const&>(
+        py::overload_cast<std::string const &, std::string const &>(
             &kintera::classify_kinetics_base_titan_reactions),
         py::arg("pun_path"), py::arg("run_input_path"));
 
@@ -219,45 +230,49 @@ PYBIND11_MODULE(kintera, m) {
   bind_photolysis(m);
   bind_diffusion(m);
   bind_sparse_solver(m);
+  bind_equilibrium(m);
 
-  m.def("species_names", []() -> const std::vector<std::string>& {
+  m.def("species_names", []() -> const std::vector<std::string> & {
     return kintera::species_names;
   });
 
-  m.def("set_species_names", [](const std::vector<std::string>& names) {
+  m.def("set_species_names", [](const std::vector<std::string> &names) {
     kintera::species_names = names;
     return kintera::species_names;
   });
 
-  m.def("species_weights", []() -> const std::vector<double>& {
+  m.def("species_weights", []() -> const std::vector<double> & {
     return kintera::species_weights;
   });
 
-  m.def("set_species_weights", [](const std::vector<double>& weights) {
+  m.def("set_species_weights", [](const std::vector<double> &weights) {
     kintera::species_weights = weights;
     return kintera::species_weights;
   });
 
-  m.def("species_cref_R",
-        []() -> const std::vector<double>& { return kintera::species_cref_R; });
+  m.def("species_cref_R", []() -> const std::vector<double> & {
+    return kintera::species_cref_R;
+  });
 
-  m.def("set_species_cref_R", [](const std::vector<double>& cref_R) {
+  m.def("set_species_cref_R", [](const std::vector<double> &cref_R) {
     kintera::species_cref_R = cref_R;
     return kintera::species_cref_R;
   });
 
-  m.def("species_uref_R",
-        []() -> const std::vector<double>& { return kintera::species_uref_R; });
+  m.def("species_uref_R", []() -> const std::vector<double> & {
+    return kintera::species_uref_R;
+  });
 
-  m.def("set_species_uref_R", [](const std::vector<double>& uref_R) {
+  m.def("set_species_uref_R", [](const std::vector<double> &uref_R) {
     kintera::species_uref_R = uref_R;
     return kintera::species_uref_R;
   });
 
-  m.def("species_sref_R",
-        []() -> const std::vector<double>& { return kintera::species_sref_R; });
+  m.def("species_sref_R", []() -> const std::vector<double> & {
+    return kintera::species_sref_R;
+  });
 
-  m.def("set_species_sref_R", [](const std::vector<double>& sref_R) {
+  m.def("set_species_sref_R", [](const std::vector<double> &sref_R) {
     kintera::species_sref_R = sref_R;
     return kintera::species_sref_R;
   });

python/csrc/pyequilibrium.cpp

@@ -0,0 +1,44 @@
+#include <fmt/format.h>
+#include <pybind11/stl.h>
+#include <torch/extension.h>
+
+#include <kintera/equilibrium/equilibrium.hpp>
+#include <sstream>
+
+#include "pyoptions.hpp"
+
+namespace py = pybind11;
+
+void bind_equilibrium(py::module &m) {
+  auto pyOptions =
+      py::class_<kintera::EquilibriumOptionsImpl, kintera::EquilibriumOptions>(
+          m, "EquilibriumOptions");
+
+  pyOptions.def(py::init<>())
+      .def("__repr__",
+           [](kintera::EquilibriumOptions const &self) {
+             std::stringstream ss;
+             self->report(ss);
+             return ss.str();
+           })
+      .def_static("from_yaml", &kintera::EquilibriumOptionsImpl::from_yaml,
+                  py::arg("filename"), py::arg("verbose") = false)
+      .ADD_OPTION(std::vector<std::string>, kintera::EquilibriumOptionsImpl,
+                  components)
+      .ADD_OPTION(std::vector<std::string>, kintera::EquilibriumOptionsImpl,
+                  phases)
+      .ADD_OPTION(std::vector<std::string>, kintera::EquilibriumOptionsImpl,
+                  reactions)
+      .ADD_OPTION(std::vector<int>, kintera::EquilibriumOptionsImpl, phase_ids)
+      .ADD_OPTION(int, kintera::EquilibriumOptionsImpl, gas_phase)
+      .ADD_OPTION(double, kintera::EquilibriumOptionsImpl, standard_pressure)
+      .ADD_OPTION(int, kintera::EquilibriumOptionsImpl, max_iter)
+      .ADD_OPTION(double, kintera::EquilibriumOptionsImpl, ftol)
+      .ADD_OPTION(double, kintera::EquilibriumOptionsImpl, mole_floor)
+      .def("validate", &kintera::EquilibriumOptionsImpl::validate);
+
+  ADD_KINTERA_MODULE(EquilibriumTP, EquilibriumOptions,
+                     &kintera::EquilibriumTPImpl::forward, py::arg("temp"),
+                     py::arg("pres"), py::arg("moles"), py::arg("log_k"),
+                     py::arg("warm_start") = false);
+}

python/equilibrium/__init__.py

@@ -0,0 +1,3 @@
+from .nasa9 import Nasa9LogK
+
+__all__ = ["Nasa9LogK"]

python/equilibrium/nasa9.py

@@ -0,0 +1,28 @@
+"""NASA-9 equilibrium-constant models."""
+
+from collections.abc import Sequence
+
+import torch
+
+from ..kintera import EquilibriumOptions, EquilibriumTP, nasa9_gibbs_rt
+
+
+class Nasa9LogK:
+    """Compute reaction log(K) from bundled ideal-gas NASA-9 data."""
+
+    def __init__(
+        self,
+        options: EquilibriumOptions,
+        nasa9_species: Sequence[str] | None = None,
+    ) -> None:
+        self.species = list(nasa9_species or options.components())
+        if len(self.species) != len(options.components()):
+            raise ValueError("nasa9_species must contain one name per component")
+        self.stoich = EquilibriumTP(options).buffer("stoich")
+
+    def __call__(self, temp: torch.Tensor, pressure: torch.Tensor) -> torch.Tensor:
+        """Return natural-log equilibrium constants in option reaction order."""
+        del pressure
+        gibbs_rt = nasa9_gibbs_rt(temp, self.species)
+        stoich = self.stoich.to(device=temp.device, dtype=temp.dtype)
+        return -(gibbs_rt @ stoich)