P to dynes/cm2. Python install

.github/workflows/test.yaml

@@ -19,6 +19,11 @@ jobs:
         compiler: gcc
         version: 11
 
+    - name: Setup Python
+      uses: actions/setup-python@v3
+      with:
+        python-version: 3.9
+
     - name: Install dependencies
       run: |
         sudo apt-get install valgrind
@@ -32,4 +37,11 @@ jobs:
     - name: test
       working-directory: ${{github.workspace}}/build
       run: |
-        valgrind --error-exitcode=1 --leak-check=full ./test/test_equilibrate
+        valgrind --error-exitcode=1 --leak-check=full ./test/test_equilibrate
+
+    - name: build python
+      run: python -m pip install . -v
+
+    - name: test python
+      working-directory: ${{github.workspace}}/test
+      run: python test_equilibrate.py

CMakeLists.txt

@@ -1,6 +1,6 @@
 cmake_minimum_required(VERSION "3.14")
 
-project(EQUILIBRATE LANGUAGES Fortran C VERSION "0.1.0")
+project(EQUILIBRATE LANGUAGES Fortran C VERSION "0.1.1")
 
 set(CMAKE_Fortran_MODULE_DIRECTORY "${CMAKE_BINARY_DIR}/modules")
 include(cmake/CPM.cmake)

README.md

@@ -0,0 +1,8 @@
+# Equilibrate
+
+`Equilibrate` is a chemical equilibrium solver which is based on [EasyChem](https://gitlab.com/EliseLei/easychem), which itself is a clone of the NASA CEA tool (described in [Gordon and McBride 1994](https://ntrs.nasa.gov/api/citations/19950013764/downloads/19950013764.pdf)).
+
+
+
+
+

equilibrate/cython/_equilibrate.pyx

@@ -95,7 +95,7 @@ cdef class ChemEquiAnalysis:
     Parameters
     ----------
     P : double
-        Pressure in bars
+        Pressure in dynes/cm^2
     T : double
         Temperature in Kelvin
     molfracs_atoms : ndarray[double,ndim=1], optional

pyproject.toml

@@ -0,0 +1,2 @@
+[build-system]
+requires = ["setuptools>=42", "wheel", "scikit-build", "cmake>=3.12", "ninja", "numpy", "cython"]

setup.py

@@ -0,0 +1,33 @@
+from skbuild import setup
+from os import path
+
+this_directory = path.abspath(path.dirname(__file__))
+with open(path.join(this_directory, 'README.md'), encoding='utf-8') as f:
+    long_description = f.read()
+
+with open("CMakeLists.txt",'r') as fil:
+    lines = fil.readlines()
+    for line in lines:
+        if line.startswith("project(EQUILIBRATE"):
+            version = line.split('"')[1]
+            break
+
+setup(
+    name="equilibrate",
+    packages=['equilibrate'],
+    python_requires='>=3.6',
+    version=version,
+    license="GNU General Public License v3.0",
+    install_requires=['numpy'], 
+    author='Nicholas Wogan',
+    author_email = '[email protected]',
+    description = "A chemical equilibrium solver.",
+    long_description=long_description,
+    long_description_content_type='text/markdown',
+    url = "https://github.com/Nicholaswogan/Equilibrate",
+    cmake_args=['-DCMAKE_BUILD_TYPE=Release',\
+                '-DBUILD_PYTHON_EQUILIBRATE=ON',\
+                '-DBUILD_EXECUTABLES=OFF']
+)
+
+

src/equilibrate.f90

@@ -163,7 +163,7 @@ function create_ChemEquiAnalysis(thermopath, atoms, species, err) result(cea)
   !> `err` is allocated with an error message.
   function solve(self, P, T, molfracs_atoms, molfracs_species, err) result(converged)
     class(ChemEquiAnalysis), intent(inout) :: self
-    real(dp), intent(in) :: P !! Pressure in bars
+    real(dp), intent(in) :: P !! Pressure in dynes/cm^2
     real(dp), intent(in) :: T !! Temperature in Kelvin
     !> Atom mole fractions in the same order and length as self%atoms_names.
     real(dp), optional, intent(in) :: molfracs_atoms(:)
@@ -173,9 +173,11 @@ function solve(self, P, T, molfracs_atoms, molfracs_species, err) result(converg
     logical :: converged
 
     real(dp), allocatable :: molfracs_atoms_(:)
+    real(dp) :: P_bars
     integer :: i, j, jj
 
     converged = .false.
+    P_bars = P/1.0e6_dp
 
     if (present(molfracs_atoms) .and. present(molfracs_species)) then
       err = 'Both "molfracs_atoms" and "molfracs_species" are inputs, but only one is allowed.'
@@ -224,7 +226,7 @@ function solve(self, P, T, molfracs_atoms, molfracs_species, err) result(converg
                            mass_tol=self%mass_tol, &
                            molfracs_reactants=self%molfracs_species, &
                            massfracs_reactants=self%massfracs_species, &
-                           temp=T, press=P, &
+                           temp=T, press=P_bars, &
                            nabla_ad=self%nabla_ad, gamma2=self%gamma2, MMW=self%MMW, rho=self%rho, c_pe=self%c_pe)
     if (self%dat%error) then
       err = trim(self%dat%err_msg)

test/test_equilibrate.f90

@@ -264,14 +264,14 @@ subroutine test()
       stop 1
     endif
 
-    converged = cea%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
+    converged = cea%solve(1.0e6_dp, 1000.0_dp, molfracs_atoms=X, err=err)
     if (allocated(err)) then
       print*,err
       stop 1
     endif
     if (.not.converged) stop 1
 
-    converged = cea2%solve(1.0_dp, 1000.0_dp, molfracs_atoms=X, err=err)
+    converged = cea2%solve(1.0e6_dp, 1000.0_dp, molfracs_atoms=X, err=err)
     if (allocated(err)) then
       print*,err
       stop 1
@@ -294,7 +294,7 @@ subroutine test()
       endif
     enddo
 
-    converged = cea%solve(1.0_dp, 1000.0_dp, molfracs_species=cea%molfracs_species, err=err)
+    converged = cea%solve(1.0e6_dp, 1000.0_dp, molfracs_species=cea%molfracs_species, err=err)
     if (allocated(err)) then
       print*,err
       stop 1

test/test_equilibrate.py

@@ -0,0 +1,55 @@
+import numpy as np
+from equilibrate import ChemEquiAnalysis
+
+def main():
+    atoms = [
+        'H ',
+        'He',
+        'C ',
+        'N ',
+        'O ',
+        'Na',
+        'Mg',
+        'Al',
+        'Si',
+        'P ',
+        'S ',
+        'Cl',
+        'K ',
+        'Ca',
+        'Ti',
+        'V ',
+        'Fe',
+        'Ni'
+    ]
+    X = np.array([
+        9.207539305000000e-01,
+        7.836886940000000e-02,
+        2.478241000000000e-04,
+        6.225060569498810e-05,
+        4.509658000000000e-04,
+        1.600086943532050e-06,
+        3.665587420553620e-05,
+        2.595000000000000e-06,
+        2.979500000000000e-05,
+        2.366702019976680e-07,
+        1.213790073460400e-05,
+        2.911679584995890e-07,
+        9.866056119256769e-08,
+        2.014390114292550e-06,
+        8.206228043663590e-08,
+        7.836886940899920e-09,
+        2.911679584995890e-05,
+        1.528071168062810e-06
+    ])
+
+    # 
+    cea = ChemEquiAnalysis('thermo_easy_chem_simp_own.yaml', atoms=atoms)
+    converged = cea.solve(1.0e6, 1000.0, molfracs_atoms=X)
+
+    # Check the solution
+    ind = cea.species_names.index('H2')
+    assert np.isclose(cea.molfracs_species[ind],8.531731613887943e-01,rtol=1e-4)
+
+if __name__ == '__main__':
+    main()