diff --git a/CHANGELOG.md b/CHANGELOG.md index 6053719..d968ce0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,10 +5,20 @@ All notable user-visible changes to this project are documented here. ## [Unreleased] ### Changed +- Python `Mixture` input validation now accepts `str` and `cea.Reactant` entries (including mixed lists) and no longer accepts raw `bytes` species names (`#53`). +- Added SI-focused custom-reactant handling at the Python API layer: `Reactant.temperature` is specified in K and `Reactant.enthalpy` in J/kg (converted internally for core input) (`#53`). +- Legacy input parsing now supports repeated `outp` dataset keywords (including multiline forms) by merging successive `outp` entries during dataset assembly (`#52`). ### Fixed +- Legacy CLI equilibrium/rocket/shock workflows now propagate `include_ions` into generated product mixtures so ionized products are retained when requested (`#52`). ### Added +- Added C and Python support for custom reactant data (including species not present in `thermo.lib`) in parity with the main interface workflow used by RP-1311 Example 5 (`#53`). +- Added new C-API constructors for generating product mixtures from input-reactant payloads: + - `cea_mixture_create_products_from_input_reactants` (`#53`). + - `cea_mixture_create_products_from_input_reactants_w_ions` (`#53`). +- Added a shared bindc parser path for `cea_reactant_input -> ReactantInput` conversion to reduce duplicated C-binding logic (`#53`). +- Added Python `cea.Reactant` and mixed-input `Mixture(...)` support in the Cython binding (`#53`). ## [3.1.0] - 2026-03-02 diff --git a/docs/source/examples/equilibrium/example5.rst b/docs/source/examples/equilibrium/example5.rst new file mode 100644 index 0000000..ece71f0 --- /dev/null +++ b/docs/source/examples/equilibrium/example5.rst @@ -0,0 +1,54 @@ +Example 5 from RP-1311 +====================== +.. note:: The python script for this example is available in the `source/bind/python/cea/samples` directory of the CEA repository. + +Here we describe how to run example 5 from RP-1311 [1]_ using the Python API. +This is an HP equilibrium problem for a solid propellant blend with five reactants, +including one custom reactant (``CHOS-Binder``) that is not present in ``thermo.lib``. + +First import the required libraries: + +.. code-block:: python + + import numpy as np + import cea + +Define the pressure schedule and reactant composition by weight fraction: + +.. code-block:: python + + pressures = np.array([34.473652, 17.236826, 8.618413, 3.447365, 0.344737]) + weights = np.array([0.7206, 0.1858, 0.09, 0.002, 0.0016], dtype=np.float64) + T_reac = np.array([298.15, 298.15, 298.15, 298.15, 298.15], dtype=np.float64) + +Create a :class:`~cea.Reactant` for ``CHOS-Binder`` using SI values. +In the Python API, custom reactant ``enthalpy`` is in J/kg and ``temperature`` is in K. +Use :mod:`cea.units` helpers for pre-conversion as needed. + +.. code-block:: python + + chos_binder_mw_kg_per_mol = 14.6652984484e-3 + chos_binder_h_si = cea.units.cal_to_joule(-2999.082) / chos_binder_mw_kg_per_mol + + reactants = [ + "NH4CLO4(I)", + cea.Reactant( + name="CHOS-Binder", + formula={"C": 1.0, "H": 1.86955, "O": 0.031256, "S": 0.008415}, + molecular_weight=14.6652984484, + enthalpy=chos_binder_h_si, + temperature=298.15, + ), + "AL(cr)", + "MgO(cr)", + "H2O(L)", + ] + +Then build reactant/product mixtures, solve the HP problem, and print the full table output: + +.. literalinclude:: ../../../../source/bind/python/cea/samples/rp1311/example5.py + :language: python + +.. rubric:: References + +.. [1] B. J. McBride and S. Gordon, *Computer Program for Calculation of Complex Chemical Equilibrium Compositions and Applications*, NASA RP-1311, 1996. diff --git a/docs/source/examples/equilibrium_examples.rst b/docs/source/examples/equilibrium_examples.rst index 6a6d13a..7968a70 100644 --- a/docs/source/examples/equilibrium_examples.rst +++ b/docs/source/examples/equilibrium_examples.rst @@ -10,4 +10,5 @@ This section describes how to solve equilibrium problems using the Python interf equilibrium/example2 equilibrium/example3 equilibrium/example4 - equilibrium/example14 \ No newline at end of file + equilibrium/example5 + equilibrium/example14 diff --git a/docs/source/interfaces/python_api.rst b/docs/source/interfaces/python_api.rst index 0b8fae5..442a58b 100644 --- a/docs/source/interfaces/python_api.rst +++ b/docs/source/interfaces/python_api.rst @@ -9,6 +9,11 @@ Mixture ------- The :class:`~cea.Mixture` class is used to define a mixture of product or reactant species. It allows the user to specify the composition of the mixture and provides methods to compute thermodynamic curve fit properties. The instances of this class are then passed as inputs to the available solver classes (e.g., :class:`~cea.EqSolver`, :class:`~cea.RocketSolver`, :class:`~cea.ShockSolver`, or :class:`~cea.DetonationSolver`). +Custom reactants can be provided through :class:`~cea.Reactant` objects (including mixed lists of strings and Reactant objects). +For :class:`~cea.Reactant`, ``enthalpy`` and ``temperature`` are SI-only in the Python API (J/kg and K, respectively); use :mod:`cea.units` helpers for pre-conversion. + +.. autoclass:: cea.Reactant + :members: .. autoclass:: cea.Mixture :members: diff --git a/source/bind/c/CMakeLists.txt b/source/bind/c/CMakeLists.txt index 3509795..5a0a981 100644 --- a/source/bind/c/CMakeLists.txt +++ b/source/bind/c/CMakeLists.txt @@ -56,6 +56,15 @@ if (CEA_BUILD_TESTING) COMMAND cea_bindc_rp1311_ex3 WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} ) + + add_executable(cea_bindc_rp1311_ex5 samples/rp1311_example5.c) + target_link_libraries(cea_bindc_rp1311_ex5 PRIVATE cea::bindc) + add_test( + NAME cea_bindc_rp1311_ex5 + COMMAND cea_bindc_rp1311_ex5 + WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} + ) + add_executable(cea_bindc_rp1311_ex6 samples/rp1311_example6.c) target_link_libraries(cea_bindc_rp1311_ex6 PRIVATE cea::bindc) add_test( @@ -89,4 +98,3 @@ if (CEA_BUILD_TESTING) ) endif() - diff --git a/source/bind/c/bindc.F90 b/source/bind/c/bindc.F90 index a3d1bd4..bcd62a0 100644 --- a/source/bind/c/bindc.F90 +++ b/source/bind/c/bindc.F90 @@ -574,227 +574,119 @@ function cea_mixture_create_from_input_reactants(mptr, nreac, creac) result(ierr type(cea_reactant_input), intent(in) :: creac(*) type(Mixture), pointer :: mix type(ReactantInput), allocatable :: input_reactants(:) - type(c_ptr), pointer :: c_elements(:) - real(c_double), pointer :: c_coeffs(:) - character(:), allocatable :: name, units, elem - integer :: i, j, ne ierr = CEA_SUCCESS - if (nreac <= 0) then - ierr = CEA_INVALID_SIZE - return - end if - - allocate(input_reactants(nreac)) - - do i = 1, nreac - if (.not. c_associated(creac(i)%name)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%name, name) - if (len_trim(name) == 0 .or. len(name) > snl) then - ierr = CEA_INVALID_SIZE - return - end if - input_reactants(i)%name = name - - ne = creac(i)%num_elements - if (ne < 0) then - ierr = CEA_INVALID_SIZE - return - end if + call parse_c_reactant_inputs(nreac, creac, input_reactants, ierr) + if (ierr /= CEA_SUCCESS) return - if (ne > 0) then - if (.not. c_associated(creac(i)%elements) .or. .not. c_associated(creac(i)%coefficients)) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%formula) - allocate(input_reactants(i)%formula%elements(ne)) - allocate(input_reactants(i)%formula%coefficients(ne)) - call c_f_pointer(creac(i)%elements, c_elements, [ne]) - call c_f_pointer(creac(i)%coefficients, c_coeffs, [ne]) - do j = 1, ne - if (.not. c_associated(c_elements(j))) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(c_elements(j), elem) - if (len_trim(elem) == 0 .or. len(elem) > enl) then - ierr = CEA_INVALID_SIZE - return - end if - input_reactants(i)%formula%elements(j) = elem - input_reactants(i)%formula%coefficients(j) = c_coeffs(j) - end do - else - if (.not. thermodb_has_species(name)) then - ierr = CEA_INVALID_SIZE - return - end if - end if + allocate(mix) + mix = Mixture(global_thermodb, input_reactants=input_reactants) + mptr = c_loc(mix) + call log_info('BINDC: Created Mixture object at '//to_str(mptr)) - if (creac(i)%has_molecular_weight) then - allocate(input_reactants(i)%molecular_weight) - input_reactants(i)%molecular_weight = creac(i)%molecular_weight - end if + end function - if (creac(i)%has_enthalpy) then - if (.not. c_associated(creac(i)%enthalpy_units)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%enthalpy_units, units) - if (len_trim(units) == 0) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%enthalpy) - input_reactants(i)%enthalpy%name = 'h' - input_reactants(i)%enthalpy%units = units - allocate(input_reactants(i)%enthalpy%values(1)) - input_reactants(i)%enthalpy%values(1) = creac(i)%enthalpy - end if + function cea_mixture_create_from_input_reactants_w_ions(mptr, nreac, creac) result(ierr) bind(c) + integer(c_int) :: ierr + type(c_ptr), intent(out) :: mptr + integer(c_int), value :: nreac + type(cea_reactant_input), intent(in) :: creac(*) + type(Mixture), pointer :: mix + type(ReactantInput), allocatable :: input_reactants(:) - if (creac(i)%has_temperature) then - if (.not. c_associated(creac(i)%temperature_units)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%temperature_units, units) - if (len_trim(units) == 0) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%temperature) - input_reactants(i)%temperature%name = 't' - input_reactants(i)%temperature%units = units - allocate(input_reactants(i)%temperature%values(1)) - input_reactants(i)%temperature%values(1) = creac(i)%temperature - end if - end do + ierr = CEA_SUCCESS + call parse_c_reactant_inputs(nreac, creac, input_reactants, ierr) + if (ierr /= CEA_SUCCESS) return allocate(mix) - mix = Mixture(global_thermodb, input_reactants=input_reactants) + mix = Mixture(global_thermodb, input_reactants=input_reactants, ions=.true.) mptr = c_loc(mix) call log_info('BINDC: Created Mixture object at '//to_str(mptr)) end function - function cea_mixture_create_from_input_reactants_w_ions(mptr, nreac, creac) result(ierr) bind(c) + function cea_mixture_create_products_from_input_reactants(mptr, nreac, creac, nomit, comit) result(ierr) bind(c) integer(c_int) :: ierr type(c_ptr), intent(out) :: mptr integer(c_int), value :: nreac type(cea_reactant_input), intent(in) :: creac(*) + integer(c_int), value :: nomit + type(c_ptr), intent(in) :: comit(*) type(Mixture), pointer :: mix + type(Mixture) :: reactants type(ReactantInput), allocatable :: input_reactants(:) - type(c_ptr), pointer :: c_elements(:) - real(c_double), pointer :: c_coeffs(:) - character(:), allocatable :: name, units, elem - integer :: i, j, ne + character(snl), allocatable :: product_names(:) + character(snl), allocatable :: omit(:) + character(:), allocatable :: name + integer :: n ierr = CEA_SUCCESS - if (nreac <= 0) then + if (nomit < 0) then ierr = CEA_INVALID_SIZE return end if + call parse_c_reactant_inputs(nreac, creac, input_reactants, ierr) + if (ierr /= CEA_SUCCESS) return - allocate(input_reactants(nreac)) - - do i = 1, nreac - if (.not. c_associated(creac(i)%name)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%name, name) - if (len_trim(name) == 0 .or. len(name) > snl) then + allocate(omit(nomit)) + do n = 1,nomit + call c_copy(comit(n), name) + if (len(name) > snl) then ierr = CEA_INVALID_SIZE return end if - input_reactants(i)%name = name + omit(n) = name + end do - ne = creac(i)%num_elements - if (ne < 0) then - ierr = CEA_INVALID_SIZE - return - end if + reactants = Mixture(global_thermodb, input_reactants=input_reactants) + product_names = reactants%get_products(global_thermodb, omit) - if (ne > 0) then - if (.not. c_associated(creac(i)%elements) .or. .not. c_associated(creac(i)%coefficients)) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%formula) - allocate(input_reactants(i)%formula%elements(ne)) - allocate(input_reactants(i)%formula%coefficients(ne)) - call c_f_pointer(creac(i)%elements, c_elements, [ne]) - call c_f_pointer(creac(i)%coefficients, c_coeffs, [ne]) - do j = 1, ne - if (.not. c_associated(c_elements(j))) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(c_elements(j), elem) - if (len_trim(elem) == 0 .or. len(elem) > enl) then - ierr = CEA_INVALID_SIZE - return - end if - input_reactants(i)%formula%elements(j) = elem - input_reactants(i)%formula%coefficients(j) = c_coeffs(j) - end do - else - if (.not. thermodb_has_species(name)) then - ierr = CEA_INVALID_SIZE - return - end if - end if + allocate(mix) + mix = Mixture(global_thermodb, product_names) + mptr = c_loc(mix) + call log_info('BINDC: Created product Mixture object at '//to_str(mptr)) + end function - if (creac(i)%has_molecular_weight) then - allocate(input_reactants(i)%molecular_weight) - input_reactants(i)%molecular_weight = creac(i)%molecular_weight - end if + function cea_mixture_create_products_from_input_reactants_w_ions(mptr, nreac, creac, nomit, comit) result(ierr) bind(c) + integer(c_int) :: ierr + type(c_ptr), intent(out) :: mptr + integer(c_int), value :: nreac + type(cea_reactant_input), intent(in) :: creac(*) + integer(c_int), value :: nomit + type(c_ptr), intent(in) :: comit(*) + type(Mixture), pointer :: mix + type(Mixture) :: reactants + type(ReactantInput), allocatable :: input_reactants(:) + character(snl), allocatable :: product_names(:) + character(snl), allocatable :: omit(:) + character(:), allocatable :: name + integer :: n - if (creac(i)%has_enthalpy) then - if (.not. c_associated(creac(i)%enthalpy_units)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%enthalpy_units, units) - if (len_trim(units) == 0) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%enthalpy) - input_reactants(i)%enthalpy%name = 'h' - input_reactants(i)%enthalpy%units = units - allocate(input_reactants(i)%enthalpy%values(1)) - input_reactants(i)%enthalpy%values(1) = creac(i)%enthalpy - end if + ierr = CEA_SUCCESS + if (nomit < 0) then + ierr = CEA_INVALID_SIZE + return + end if + call parse_c_reactant_inputs(nreac, creac, input_reactants, ierr) + if (ierr /= CEA_SUCCESS) return - if (creac(i)%has_temperature) then - if (.not. c_associated(creac(i)%temperature_units)) then - ierr = CEA_INVALID_SIZE - return - end if - call c_copy(creac(i)%temperature_units, units) - if (len_trim(units) == 0) then - ierr = CEA_INVALID_SIZE - return - end if - allocate(input_reactants(i)%temperature) - input_reactants(i)%temperature%name = 't' - input_reactants(i)%temperature%units = units - allocate(input_reactants(i)%temperature%values(1)) - input_reactants(i)%temperature%values(1) = creac(i)%temperature + allocate(omit(nomit)) + do n = 1,nomit + call c_copy(comit(n), name) + if (len(name) > snl) then + ierr = CEA_INVALID_SIZE + return end if + omit(n) = name end do + reactants = Mixture(global_thermodb, input_reactants=input_reactants, ions=.true.) + product_names = reactants%get_products(global_thermodb, omit) + allocate(mix) - mix = Mixture(global_thermodb, input_reactants=input_reactants, ions=.true.) + mix = Mixture(global_thermodb, product_names, ions=.true.) mptr = c_loc(mix) - call log_info('BINDC: Created Mixture object at '//to_str(mptr)) - + call log_info('BINDC: Created product Mixture object at '//to_str(mptr)) end function function cea_mixture_destroy(mptr) result(ierr) bind(c) @@ -4041,6 +3933,114 @@ function cea_detonation_solution_get_converged(slptr, converged) result(ierr) bi !----------------------------------------------------------------- ! Helper Functions !----------------------------------------------------------------- + subroutine parse_c_reactant_inputs(nreac, creac, input_reactants, ierr) + integer(c_int), intent(in), value :: nreac + type(cea_reactant_input), intent(in) :: creac(*) + type(ReactantInput), allocatable, intent(out) :: input_reactants(:) + integer(c_int), intent(out) :: ierr + + type(c_ptr), pointer :: c_elements(:) + real(c_double), pointer :: c_coeffs(:) + character(:), allocatable :: name, units, elem + integer :: i, j, ne + + ierr = CEA_SUCCESS + if (nreac <= 0) then + ierr = CEA_INVALID_SIZE + return + end if + + allocate(input_reactants(nreac)) + + do i = 1, nreac + if (.not. c_associated(creac(i)%name)) then + ierr = CEA_INVALID_SIZE + return + end if + call c_copy(creac(i)%name, name) + if (len_trim(name) == 0 .or. len(name) > snl) then + ierr = CEA_INVALID_SIZE + return + end if + input_reactants(i)%name = name + + ne = creac(i)%num_elements + if (ne < 0) then + ierr = CEA_INVALID_SIZE + return + end if + + if (ne > 0) then + if (.not. c_associated(creac(i)%elements) .or. .not. c_associated(creac(i)%coefficients)) then + ierr = CEA_INVALID_SIZE + return + end if + allocate(input_reactants(i)%formula) + allocate(input_reactants(i)%formula%elements(ne)) + allocate(input_reactants(i)%formula%coefficients(ne)) + call c_f_pointer(creac(i)%elements, c_elements, [ne]) + call c_f_pointer(creac(i)%coefficients, c_coeffs, [ne]) + do j = 1, ne + if (.not. c_associated(c_elements(j))) then + ierr = CEA_INVALID_SIZE + return + end if + call c_copy(c_elements(j), elem) + if (len_trim(elem) == 0 .or. len(elem) > enl) then + ierr = CEA_INVALID_SIZE + return + end if + input_reactants(i)%formula%elements(j) = elem + input_reactants(i)%formula%coefficients(j) = c_coeffs(j) + end do + else + if (.not. thermodb_has_species(name)) then + ierr = CEA_INVALID_SIZE + return + end if + end if + + if (creac(i)%has_molecular_weight) then + allocate(input_reactants(i)%molecular_weight) + input_reactants(i)%molecular_weight = creac(i)%molecular_weight + end if + + if (creac(i)%has_enthalpy) then + if (.not. c_associated(creac(i)%enthalpy_units)) then + ierr = CEA_INVALID_SIZE + return + end if + call c_copy(creac(i)%enthalpy_units, units) + if (len_trim(units) == 0) then + ierr = CEA_INVALID_SIZE + return + end if + allocate(input_reactants(i)%enthalpy) + input_reactants(i)%enthalpy%name = 'h' + input_reactants(i)%enthalpy%units = units + allocate(input_reactants(i)%enthalpy%values(1)) + input_reactants(i)%enthalpy%values(1) = creac(i)%enthalpy + end if + + if (creac(i)%has_temperature) then + if (.not. c_associated(creac(i)%temperature_units)) then + ierr = CEA_INVALID_SIZE + return + end if + call c_copy(creac(i)%temperature_units, units) + if (len_trim(units) == 0) then + ierr = CEA_INVALID_SIZE + return + end if + allocate(input_reactants(i)%temperature) + input_reactants(i)%temperature%name = 't' + input_reactants(i)%temperature%units = units + allocate(input_reactants(i)%temperature%values(1)) + input_reactants(i)%temperature%values(1) = creac(i)%temperature + end if + end do + end subroutine + logical function thermodb_has_species(name) result(found) character(*), intent(in) :: name integer :: i diff --git a/source/bind/c/cea.h b/source/bind/c/cea.h index 9dc518a..0be18ae 100644 --- a/source/bind/c/cea.h +++ b/source/bind/c/cea.h @@ -170,6 +170,20 @@ extern "C" const cea_int nreactants, const cea_reactant_input reactants[]); + cea_err cea_mixture_create_products_from_input_reactants( + cea_mixture *mix, + const cea_int nreactants, + const cea_reactant_input reactants[], + const cea_int nomit, + const cea_string omit[]); + + cea_err cea_mixture_create_products_from_input_reactants_w_ions( + cea_mixture *mix, + const cea_int nreactants, + const cea_reactant_input reactants[], + const cea_int nomit, + const cea_string omit[]); + cea_err cea_mixture_destroy( cea_mixture *mix); diff --git a/source/bind/c/samples/rp1311_example5.c b/source/bind/c/samples/rp1311_example5.c new file mode 100644 index 0000000..f32cf16 --- /dev/null +++ b/source/bind/c/samples/rp1311_example5.c @@ -0,0 +1,76 @@ +#include "stdio.h" +#include "cea.h" + +#define LEN(x) (sizeof(x) / sizeof((x)[0])) + +int main(void) +{ + const cea_real pressures[] = {34.473652, 17.236826, 8.618413, 3.447365, 0.344737}; + const cea_real reactant_weights[] = {0.7206, 0.1858, 0.09, 0.002, 0.0016}; + const cea_real reactant_temps[] = {298.15, 298.15, 298.15, 298.15, 298.15}; + const cea_string omitted_products[] = { + "CCN", "CNC", "C2N2", "C2O", "C3H4,allene", "C3H4,propyne", "C3H4,cyclo-", "C3", + "C3H5,allyl", "C3H6,propylene", "C3H6,cyclo-", "C3H3,propargyl", "C3H6O", "C3H7,n-propyl", + "C3H7,i-propyl", "Jet-A(g)", "C3O2", "C4", "C4H2", "C3H8O,2propanol", "C4H4,1,3-cyclo-", + "C4H6,butadiene", "C4H6,2-butyne", "C3H8O,1propanol", "C4H8,tr2-butene", "C4H8,isobutene", + "C4H8,cyclo-", "C4H6,cyclo-", "(CH3COOH)2", "C4H9,n-butyl", "C4H9,i-butyl", "C4H8,1-butene", + "C4H9,s-butyl", "C4H9,t-butyl", "C4H10,isobutane", "C4H8,cis2-buten", "C4H10,n-butane", + "C4N2", "C5", "C3H8", "Jet-A(L)", "C6H6(L)", "H2O(s)", "H2O(L)"}; + + const cea_string binder_elements[] = {"C", "H", "O", "S"}; + const cea_real binder_coeffs[] = {1.0, 1.86955, 0.031256, 0.008415}; + const cea_string unit_enthalpy = "j/mole"; + const cea_string unit_temp = "k"; + + cea_reactant_input reactants[] = { + {.name = "NH4CLO4(I)"}, + {.name = "CHOS-Binder", + .num_elements = 4, + .elements = binder_elements, + .coefficients = binder_coeffs, + .has_enthalpy = true, + .enthalpy = -12548.159088, + .enthalpy_units = unit_enthalpy, + .has_temperature = true, + .temperature = 298.15, + .temperature_units = unit_temp}, + {.name = "AL(cr)"}, + {.name = "MgO(cr)"}, + {.name = "H2O(L)"}}; + + cea_set_log_level(CEA_LOG_NONE); + cea_init(); + + cea_mixture reac, prod; + cea_mixture_create_from_input_reactants(&reac, LEN(reactants), reactants); + cea_mixture_create_products_from_input_reactants( + &prod, LEN(reactants), reactants, LEN(omitted_products), omitted_products); + + cea_eqsolver solver; + cea_eqsolver_create_with_reactants(&solver, prod, reac); + + cea_eqsolution soln; + cea_eqsolution_create(&soln, solver); + + cea_real h0; + cea_mixture_calc_property_multitemp( + reac, CEA_ENTHALPY, LEN(reactants), reactant_weights, LEN(reactants), reactant_temps, &h0); + + printf("%12s %12s %12s %12s\n", "P(bar)", "T(K)", "Gamma_s", "Cp(cal/g-K)"); + for (int i = 0; i < LEN(pressures); ++i) + { + cea_eqsolver_solve(solver, CEA_HP, h0 / 8314.51, pressures[i], reactant_weights, soln); + cea_real temperature, gamma_s, cp_eq; + cea_eqsolution_get_property(soln, CEA_TEMPERATURE, &temperature); + cea_eqsolution_get_property(soln, CEA_GAMMA_S, &gamma_s); + cea_eqsolution_get_property(soln, CEA_EQUILIBRIUM_CP, &cp_eq); + cp_eq = cp_eq / 4.184; + printf("%12.6f %12.3f %12.5f %12.5f\n", pressures[i], temperature, gamma_s, cp_eq); + } + + cea_eqsolution_destroy(&soln); + cea_eqsolver_destroy(&solver); + cea_mixture_destroy(&prod); + cea_mixture_destroy(&reac); + return 0; +} diff --git a/source/bind/python/CEA.pyx b/source/bind/python/CEA.pyx index 092bc5d..697e835 100644 --- a/source/bind/python/CEA.pyx +++ b/source/bind/python/CEA.pyx @@ -1,5 +1,6 @@ # Import the necessary Cython modules from libc.stdlib cimport malloc, free +from libc.stddef cimport size_t from cpython.bytes cimport PyBytes_AsString import cython import ctypes @@ -516,6 +517,74 @@ def is_initialized(): return bool(initialized) +cdef class Reactant: + """ + Reactant specification for custom or overridden reactant thermo data. + + Parameters + ---------- + name : str + Reactant name. + formula : dict[str, float], optional + Exploded formula mapping element symbols to coefficients. + molecular_weight : float, optional + Molecular weight in kg/kmol (numerically equivalent to g/mol). + enthalpy : float, optional + Reference enthalpy in SI units (J/kg). + temperature : float, optional + Reference temperature in SI units (K). + + Notes + ----- + Python Reactant inputs are SI-only. Use :mod:`cea.units` helpers to pre-convert + non-SI values before constructing a Reactant. + """ + cdef public object name + cdef public object formula + cdef public object molecular_weight + cdef public object enthalpy + cdef public object temperature + + def __init__(self, name, formula=None, molecular_weight=None, enthalpy=None, temperature=None): + if not isinstance(name, str) or len(name.strip()) == 0: + raise TypeError("Reactant name must be a non-empty str") + + if formula is not None: + if not isinstance(formula, dict): + raise TypeError("Reactant formula must be a dict[str, float] when provided") + if len(formula) == 0: + raise ValueError("Reactant formula cannot be empty when provided") + for key, value in formula.items(): + if not isinstance(key, str) or len(key.strip()) == 0: + raise TypeError("Reactant formula element names must be non-empty str") + try: + fval = float(value) + except (TypeError, ValueError): + raise TypeError("Reactant formula coefficients must be numeric") + if not np.isfinite(fval): + raise ValueError("Reactant formula coefficients must be finite") + + for field_name, field_val in (("molecular_weight", molecular_weight), + ("enthalpy", enthalpy), + ("temperature", temperature)): + if field_val is None: + continue + try: + fval = float(field_val) + except (TypeError, ValueError): + raise TypeError(f"Reactant {field_name} must be numeric") + if not np.isfinite(fval): + raise ValueError(f"Reactant {field_name} must be finite") + if enthalpy is not None and molecular_weight is None: + raise ValueError("Reactant molecular_weight is required when enthalpy is specified") + + self.name = name + self.formula = formula + self.molecular_weight = molecular_weight + self.enthalpy = enthalpy + self.temperature = temperature + + cdef class Mixture: """ Mixture base class for managing chemical species compositions. @@ -526,8 +595,8 @@ cdef class Mixture: Parameters ---------- - species : list of str - List of chemical species names + species : list of str | list of Reactant | mixed list + List of chemical species names and/or Reactant objects products_from_reactants : bool, default False If True, treat species as reactants and generate corresponding products omit : list of str, default [] @@ -539,6 +608,7 @@ cdef class Mixture: cdef public int num_species cdef object _keepalive_species cdef object _keepalive_omit + cdef object _keepalive_reactants def __init__(self, list species=None, bint products_from_reactants=False, list omit=[], bint ions=False): """ Create the Mixture object @@ -546,15 +616,49 @@ cdef class Mixture: returned object is the corresponding product Mixture """ - cdef cea_err ierr + cdef cea_err ierr = SUCCESS cdef cea_int num_products - self.num_species = len(species) + cdef int i, j + cdef int ne + cdef int nomit + cdef bint has_custom = False + cdef object entry + cdef object formula_dict + cdef object fkey + cdef object fval + cdef Reactant reactant + cdef _CString _encoded cdef cea_string* cea_species = NULL - cdef int nomit = len(omit) cdef cea_string* cea_omit = NULL + cdef cea_reactant_input* cea_reactants = NULL + cdef cea_string* elem_ptrs = NULL + cdef cea_real* coeff_ptrs = NULL cdef list _species_keepalive = [] cdef list _omit_keepalive = [] - cdef _CString _encoded + cdef list _reactant_keepalive = [] + cdef list _element_ptr_buffers = [] + cdef list _coeff_ptr_buffers = [] + cdef _CString _enthalpy_units + cdef _CString _temperature_units + + if species is None: + raise TypeError("Mixture species must be provided") + if omit is None: + omit = [] + if type(species) is not list: + raise TypeError("Mixture species must be a list") + if type(omit) is not list: + raise TypeError("Mixture omit must be a list") + _enthalpy_units = _CString("j/mole", "Reactant enthalpy units") + _temperature_units = _CString("k", "Reactant temperature units") + + self.num_species = len(species) + nomit = len(omit) + for entry in species: + if isinstance(entry, Reactant): + has_custom = True + elif not isinstance(entry, str): + raise TypeError("Mixture species entries must be str or Reactant") cea_species = malloc(sizeof(cea_string)*len(species)) if cea_species == NULL: @@ -565,8 +669,11 @@ cdef class Mixture: raise MemoryError("Failed to allocate omit pointer buffer") try: - for i, val in enumerate(species): - _encoded = _CString(val, "Mixture species entries") + for i, entry in enumerate(species): + if isinstance(entry, Reactant): + _encoded = _CString(entry.name, "Mixture species entries") + else: + _encoded = _CString(entry, "Mixture species entries") _species_keepalive.append(_encoded) cea_species[i] = _encoded.ptr for i, val in enumerate(omit): @@ -575,28 +682,120 @@ cdef class Mixture: cea_omit[i] = _encoded.ptr self._keepalive_species = _species_keepalive self._keepalive_omit = _omit_keepalive - - if ions: - if products_from_reactants: - # Create the products mixture from the provided reactants - ierr = cea_mixture_create_from_reactants_w_ions(&self.ptr, self.num_species, cea_species, nomit, cea_omit) - ierr = cea_mixture_get_num_species(self.ptr, &num_products) - self.num_species = num_products + self._keepalive_reactants = _reactant_keepalive + + if has_custom: + cea_reactants = malloc(sizeof(cea_reactant_input) * self.num_species) + if cea_reactants == NULL: + raise MemoryError("Failed to allocate reactant input buffer") + + for i, entry in enumerate(species): + cea_reactants[i].name = cea_species[i] + cea_reactants[i].num_elements = 0 + cea_reactants[i].elements = NULL + cea_reactants[i].coefficients = NULL + cea_reactants[i].has_molecular_weight = 0 + cea_reactants[i].molecular_weight = 0.0 + cea_reactants[i].has_enthalpy = 0 + cea_reactants[i].enthalpy = 0.0 + cea_reactants[i].enthalpy_units = NULL + cea_reactants[i].has_temperature = 0 + cea_reactants[i].temperature = 0.0 + cea_reactants[i].temperature_units = NULL + + if isinstance(entry, Reactant): + reactant = entry + else: + reactant = Reactant(entry) + + formula_dict = reactant.formula + if formula_dict is not None: + ne = len(formula_dict) + elem_ptrs = malloc(sizeof(cea_string) * ne) + if elem_ptrs == NULL: + raise MemoryError("Failed to allocate reactant formula element buffer") + coeff_ptrs = malloc(sizeof(cea_real) * ne) + if coeff_ptrs == NULL: + free(elem_ptrs) + raise MemoryError("Failed to allocate reactant formula coefficient buffer") + _element_ptr_buffers.append(elem_ptrs) + _coeff_ptr_buffers.append(coeff_ptrs) + cea_reactants[i].num_elements = ne + cea_reactants[i].elements = elem_ptrs + cea_reactants[i].coefficients = coeff_ptrs + + j = 0 + for fkey, fval in formula_dict.items(): + _encoded = _CString(fkey, "Reactant formula elements") + _reactant_keepalive.append(_encoded) + elem_ptrs[j] = _encoded.ptr + coeff_ptrs[j] = float(fval) + j += 1 + + if reactant.molecular_weight is not None: + cea_reactants[i].has_molecular_weight = 1 + cea_reactants[i].molecular_weight = float(reactant.molecular_weight) + + if reactant.enthalpy is not None: + cea_reactants[i].has_enthalpy = 1 + cea_reactants[i].enthalpy = float(reactant.enthalpy) * float(reactant.molecular_weight) / 1000.0 + cea_reactants[i].enthalpy_units = _enthalpy_units.ptr + + if reactant.temperature is not None: + cea_reactants[i].has_temperature = 1 + cea_reactants[i].temperature = float(reactant.temperature) + cea_reactants[i].temperature_units = _temperature_units.ptr + + if ions: + if products_from_reactants: + ierr = cea_mixture_create_products_from_input_reactants_w_ions(&self.ptr, self.num_species, cea_reactants, nomit, cea_omit) + _check_ierr(ierr, "Mixture.__init__") + ierr = cea_mixture_get_num_species(self.ptr, &num_products) + _check_ierr(ierr, "Mixture.__init__") + self.num_species = num_products + else: + ierr = cea_mixture_create_from_input_reactants_w_ions(&self.ptr, self.num_species, cea_reactants) + _check_ierr(ierr, "Mixture.__init__") else: - # Create the mixture directly from the provided species (products or reactants) - ierr = cea_mixture_create_w_ions(&self.ptr, self.num_species, cea_species) + if products_from_reactants: + ierr = cea_mixture_create_products_from_input_reactants(&self.ptr, self.num_species, cea_reactants, nomit, cea_omit) + _check_ierr(ierr, "Mixture.__init__") + ierr = cea_mixture_get_num_species(self.ptr, &num_products) + _check_ierr(ierr, "Mixture.__init__") + self.num_species = num_products + else: + ierr = cea_mixture_create_from_input_reactants(&self.ptr, self.num_species, cea_reactants) + _check_ierr(ierr, "Mixture.__init__") else: - if products_from_reactants: - # Create the products mixture from the provided reactants - ierr = cea_mixture_create_from_reactants(&self.ptr, self.num_species, cea_species, nomit, cea_omit) - ierr = cea_mixture_get_num_species(self.ptr, &num_products) - self.num_species = num_products + if ions: + if products_from_reactants: + ierr = cea_mixture_create_from_reactants_w_ions(&self.ptr, self.num_species, cea_species, nomit, cea_omit) + _check_ierr(ierr, "Mixture.__init__") + ierr = cea_mixture_get_num_species(self.ptr, &num_products) + _check_ierr(ierr, "Mixture.__init__") + self.num_species = num_products + else: + ierr = cea_mixture_create_w_ions(&self.ptr, self.num_species, cea_species) + _check_ierr(ierr, "Mixture.__init__") else: - # Create the mixture directly from the provided species (products or reactants) - ierr = cea_mixture_create(&self.ptr, self.num_species, cea_species) + if products_from_reactants: + ierr = cea_mixture_create_from_reactants(&self.ptr, self.num_species, cea_species, nomit, cea_omit) + _check_ierr(ierr, "Mixture.__init__") + ierr = cea_mixture_get_num_species(self.ptr, &num_products) + _check_ierr(ierr, "Mixture.__init__") + self.num_species = num_products + else: + ierr = cea_mixture_create(&self.ptr, self.num_species, cea_species) + _check_ierr(ierr, "Mixture.__init__") finally: free(cea_species) free(cea_omit) + if cea_reactants != NULL: + free(cea_reactants) + for i in range(len(_element_ptr_buffers)): + free(_element_ptr_buffers[i]) + for i in range(len(_coeff_ptr_buffers)): + free(_coeff_ptr_buffers[i]) return diff --git a/source/bind/python/cea/samples/rp1311/example5.py b/source/bind/python/cea/samples/rp1311/example5.py index e69de29..b1e533a 100644 --- a/source/bind/python/cea/samples/rp1311/example5.py +++ b/source/bind/python/cea/samples/rp1311/example5.py @@ -0,0 +1,242 @@ +import numpy as np +import cea + +""" +Example 5 from RP-1311 +- HP equilibrium for a solid propellant blend +- Includes one custom reactant (CHOS-Binder) not present in thermo.lib +""" + +pressures = np.array([34.473652, 17.236826, 8.618413, 3.447365, 0.344737]) +weights = np.array([0.7206, 0.1858, 0.09, 0.002, 0.0016], dtype=np.float64) # wt fractions +T_reac = np.array([298.15, 298.15, 298.15, 298.15, 298.15], dtype=np.float64) + +omit_names = [ + "COOH", "C2", "C2H", "CHCO,ketyl", "C2H2,vinylidene", "CH2CO,ketene", "C2H3,vinyl", + "CH3CO,acetyl", "C2H4O,ethylen-o", "CH3CHO,ethanal", "CH3COOH", "(HCOOH)2", + "C2H5", "C2H6", "CH3N2CH3", "CH3OCH3", "C2H5OH", "CCN", "CNC", "C2N2", + "C2O", "C3", "C3H3,propargyl", "C3H4,allene", "C3H4,propyne", "C3H4,cyclo-", + "C3H5,allyl", "C3H6,propylene", "C3H6,cyclo-", "C3H6O", "C3H7,n-propyl", + "C3H7,i-propyl", "C3H8", "C3H8O,1propanol", "C3H8O,2propanol", "C3O2", + "C4", "C4H2", "C4H4,1,3-cyclo-", "C4H6,butadiene", "C4H6,2-butyne", "C4H6,cyclo-", + "C4H8,1-butene", "C4H8,cis2-buten", "C4H8,tr2-butene", "C4H8,isobutene", "C4H8,cyclo-", + "(CH3COOH)2", "C4H9,n-butyl", "C4H9,i-butyl", "C4H9,s-butyl", "C4H9,t-butyl", + "C4H10,isobutane", "C4H10,n-butane", "C4N2", "C5", "C5H6,1,3cyclo-", "C5H8,cyclo-", + "C5H10,1-pentene", "C5H10,cyclo-", "C5H11,pentyl", "C5H11,t-pentyl", "C5H12,n-pentane", + "C5H12,i-pentane", "CH3C(CH3)2CH3", "C6H2", "C6H5,phenyl", "C6H5O,phenoxy", + "C6H6", "C6H5OH,phenol", "C6H10,cyclo-", "C6H12,1-hexene", "C6H12,cyclo-", + "C6H13,n-hexyl", "C7H7,benzyl", "C7H8", "C7H8O,cresol-mx", "C7H14,1-heptene", + "C7H15,n-heptyl", "C7H16,n-heptane", "C8H8,styrene", "C8H10,ethylbenz", + "C8H16,1-octene", "C8H17,n-octyl", "C8H18,isooctane", "C8H18,n-octane", + "C9H19,n-nonyl", "C10H8,naphthale", "C10H21,n-decyl", "C12H9,o-bipheny", "C12H10,biphenyl", + "Jet-A(g)", "HNCO", "HNO", "HNO2", "HNO3", "HCCN", "HCHO,formaldehy", "HCOOH", + "NH", "NH2", "NH2OH", "NCN", "N2H2", "NH2NO2", "N2H4", "H2O2", + "(HCOOH)2", "C6H6(L)", "C7H8(L)", "C8H18(L),n-octa", "Jet-A(L)", "H2O(s)", "H2O(L)", +] + +# Convert h,cal = -2999.082 cal/mol to SI J/kg for Python Reactant input. +chos_binder_mw_kg_per_mol = 14.6652984484e-3 +chos_binder_h_si = cea.units.cal_to_joule(-2999.082) / chos_binder_mw_kg_per_mol + +reactants = [ + "NH4CLO4(I)", + cea.Reactant( + name="CHOS-Binder", + formula={"C": 1.0, "H": 1.86955, "O": 0.031256, "S": 0.008415}, + molecular_weight=14.6652984484, + enthalpy=chos_binder_h_si, + temperature=298.15, + ), + "AL(cr)", + "MgO(cr)", + "H2O(L)", +] + +reac = cea.Mixture(reactants) +prod = cea.Mixture(reactants, products_from_reactants=True, omit=omit_names) +solver = cea.EqSolver(prod, reactants=reac) +solution = cea.EqSolution(solver) + +h0 = reac.calc_property(cea.ENTHALPY, weights, T_reac) +n = len(pressures) +of_ratio_out = np.zeros(n) +P_out = np.zeros(n) +T_out = np.zeros(n) +rho = np.zeros(n) +volume = np.zeros(n) +enthalpy = np.zeros(n) +energy = np.zeros(n) +gibbs = np.zeros(n) +entropy = np.zeros(n) +molecular_weight_M = np.zeros(n) +molecular_weight_MW = np.zeros(n) +gamma_s = np.zeros(n) +cp_eq = np.zeros(n) +cp_fr = np.zeros(n) +cv_eq = np.zeros(n) +cv_fr = np.zeros(n) +mole_fractions = {} +i = 0 + +for p in pressures: + solver.solve(solution, cea.HP, h0 / cea.R, p, weights) + + of_ratio_out[i] = 0.0 + P_out[i] = cea.units.bar_to_atm(p) + if solution.converged: + T_out[i] = solution.T + rho[i] = solution.density * 1.0e-3 + volume[i] = solution.volume * 1.0e3 + enthalpy[i] = cea.units.joule_to_cal(solution.enthalpy) + energy[i] = cea.units.joule_to_cal(solution.energy) + gibbs[i] = cea.units.joule_to_cal(solution.gibbs_energy) + entropy[i] = cea.units.joule_to_cal(solution.entropy) + molecular_weight_M[i] = solution.M + molecular_weight_MW[i] = solution.MW + gamma_s[i] = solution.gamma_s + cp_eq[i] = cea.units.joule_to_cal(solution.cp_eq) + cp_fr[i] = cea.units.joule_to_cal(solution.cp_fr) + cv_eq[i] = cea.units.joule_to_cal(solution.cv_eq) + cv_fr[i] = cea.units.joule_to_cal(solution.cv_fr) + + if i == 0: + for species in solution.mole_fractions: + mole_fractions[species] = np.array([solution.mole_fractions[species]]) + else: + for species in mole_fractions: + mole_fractions[species] = np.append(mole_fractions[species], solution.mole_fractions[species]) + i += 1 + +print("o/f ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(of_ratio_out[i]), end=" ") + else: + print("{0:10.3f}".format(of_ratio_out[i])) + +print("P, atm ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(P_out[i]), end=" ") + else: + print("{0:10.3f}".format(P_out[i])) + +print("T, K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(T_out[i]), end=" ") + else: + print("{0:10.3f}".format(T_out[i])) + +print("Density, g/cc ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3e}".format(rho[i]), end=" ") + else: + print("{0:10.3e}".format(rho[i])) + +print("Volume, cc/g ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3e}".format(volume[i]), end=" ") + else: + print("{0:10.3e}".format(volume[i])) + +print("H, cal/g ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(enthalpy[i]), end=" ") + else: + print("{0:10.3f}".format(enthalpy[i])) + +print("U, cal/g ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(energy[i]), end=" ") + else: + print("{0:10.3f}".format(energy[i])) + +print("G, cal/g ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.1f}".format(gibbs[i]), end=" ") + else: + print("{0:10.1f}".format(gibbs[i])) + +print("S, cal/g-K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(entropy[i]), end=" ") + else: + print("{0:10.3f}".format(entropy[i])) + +print("M, (1/n) ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(molecular_weight_M[i]), end=" ") + else: + print("{0:10.3f}".format(molecular_weight_M[i])) + +print("MW ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.3f}".format(molecular_weight_MW[i]), end=" ") + else: + print("{0:10.3f}".format(molecular_weight_MW[i])) + +print("Gamma_s ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.4f}".format(gamma_s[i]), end=" ") + else: + print("{0:10.4f}".format(gamma_s[i])) + +print("Cp_eq, cal/g-K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.4f}".format(cp_eq[i]), end=" ") + else: + print("{0:10.4f}".format(cp_eq[i])) + +print("Cp_fr, cal/g-K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.4f}".format(cp_fr[i]), end=" ") + else: + print("{0:10.4f}".format(cp_fr[i])) + +print("Cv_eq, cal/g-K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.4f}".format(cv_eq[i]), end=" ") + else: + print("{0:10.4f}".format(cv_eq[i])) + +print("Cv_fr, cal/g-K ", end="") +for i in range(n): + if i < n - 1: + print("{0:10.4f}".format(cv_fr[i]), end=" ") + else: + print("{0:10.4f}".format(cv_fr[i])) + +print() +print("MOLE FRACTIONS") +print("") +trace_species = [] +for species in mole_fractions: + if np.any(mole_fractions[species] > 5e-6): + print("{0:15s}".format(species), end=" ") + for j in range(n): + if j < n - 1: + print("{0:10.5g}".format(mole_fractions[species][j]), end=" ") + else: + print("{0:10.5g}".format(mole_fractions[species][j])) + else: + trace_species.append(species) + +print() +print("TRACE SPECIES:") +max_cols = 10 +nrows = (len(trace_species) + max_cols - 1) // max_cols +for i in range(nrows): + print(" ".join("{0:10s}".format(trace_species[j]) for j in range(i * max_cols, min((i + 1) * max_cols, len(trace_species))))) diff --git a/source/bind/python/cea_def.pxd b/source/bind/python/cea_def.pxd index d131a3e..3d24cc8 100644 --- a/source/bind/python/cea_def.pxd +++ b/source/bind/python/cea_def.pxd @@ -229,6 +229,20 @@ cdef extern from "cea.h": ctypedef struct cea_detonation_solution_t ctypedef cea_detonation_solution_t* cea_detonation_solution + ctypedef struct cea_reactant_input: + cea_string name + cea_int num_elements + const cea_string* elements + const cea_real* coefficients + cea_bool has_molecular_weight + cea_real molecular_weight + cea_bool has_enthalpy + cea_real enthalpy + cea_string enthalpy_units + cea_bool has_temperature + cea_real temperature + cea_string temperature_units + ctypedef struct cea_solver_opts: cea_real trace cea_bool ions @@ -261,6 +275,16 @@ cdef extern from "cea.h": const cea_int nomit, const cea_string omit[]) cpdef cea_err cea_mixture_create_from_reactants_w_ions(cea_mixture *mix, const cea_int nreac, const cea_string reactants[], const cea_int nomit, const cea_string omit[]) + cpdef cea_err cea_mixture_create_from_input_reactants(cea_mixture *mix, const cea_int nreactants, + const cea_reactant_input reactants[]) + cpdef cea_err cea_mixture_create_from_input_reactants_w_ions(cea_mixture *mix, const cea_int nreactants, + const cea_reactant_input reactants[]) + cpdef cea_err cea_mixture_create_products_from_input_reactants(cea_mixture *mix, const cea_int nreactants, + const cea_reactant_input reactants[], + const cea_int nomit, const cea_string omit[]) + cpdef cea_err cea_mixture_create_products_from_input_reactants_w_ions(cea_mixture *mix, const cea_int nreactants, + const cea_reactant_input reactants[], + const cea_int nomit, const cea_string omit[]) cpdef cea_err cea_mixture_destroy(cea_mixture *mix) cpdef cea_err cea_mixture_get_num_species(const cea_mixture mix, cea_int *num_species) cpdef cea_err cea_mixture_get_species_name(const cea_mixture *mix, const cea_int i_species, cea_string *species) diff --git a/source/bind/python/tests/test_rp1311_samples.py b/source/bind/python/tests/test_rp1311_samples.py index 5a0883b..3e4cee6 100644 --- a/source/bind/python/tests/test_rp1311_samples.py +++ b/source/bind/python/tests/test_rp1311_samples.py @@ -34,6 +34,29 @@ "Jet-A(L)", "C6H6(L)", "H2O(s)", "H2O(L)", ] +_OMIT_5 = [ + "COOH", "C2", "C2H", "CHCO,ketyl", "C2H2,vinylidene", "CH2CO,ketene", "C2H3,vinyl", + "CH3CO,acetyl", "C2H4O,ethylen-o", "CH3CHO,ethanal", "CH3COOH", "(HCOOH)2", + "C2H5", "C2H6", "CH3N2CH3", "CH3OCH3", "C2H5OH", "CCN", "CNC", "C2N2", + "C2O", "C3", "C3H3,propargyl", "C3H4,allene", "C3H4,propyne", "C3H4,cyclo-", + "C3H5,allyl", "C3H6,propylene", "C3H6,cyclo-", "C3H6O", "C3H7,n-propyl", + "C3H7,i-propyl", "C3H8", "C3H8O,1propanol", "C3H8O,2propanol", "C3O2", + "C4", "C4H2", "C4H4,1,3-cyclo-", "C4H6,butadiene", "C4H6,2-butyne", "C4H6,cyclo-", + "C4H8,1-butene", "C4H8,cis2-buten", "C4H8,tr2-butene", "C4H8,isobutene", "C4H8,cyclo-", + "(CH3COOH)2", "C4H9,n-butyl", "C4H9,i-butyl", "C4H9,s-butyl", "C4H9,t-butyl", + "C4H10,isobutane", "C4H10,n-butane", "C4N2", "C5", "C5H6,1,3cyclo-", "C5H8,cyclo-", + "C5H10,1-pentene", "C5H10,cyclo-", "C5H11,pentyl", "C5H11,t-pentyl", "C5H12,n-pentane", + "C5H12,i-pentane", "CH3C(CH3)2CH3", "C6H2", "C6H5,phenyl", "C6H5O,phenoxy", + "C6H6", "C6H5OH,phenol", "C6H10,cyclo-", "C6H12,1-hexene", "C6H12,cyclo-", + "C6H13,n-hexyl", "C7H7,benzyl", "C7H8", "C7H8O,cresol-mx", "C7H14,1-heptene", + "C7H15,n-heptyl", "C7H16,n-heptane", "C8H8,styrene", "C8H10,ethylbenz", + "C8H16,1-octene", "C8H17,n-octyl", "C8H18,isooctane", "C8H18,n-octane", + "C9H19,n-nonyl", "C10H8,naphthale", "C10H21,n-decyl", "C12H9,o-bipheny", "C12H10,biphenyl", + "Jet-A(g)", "HNCO", "HNO", "HNO2", "HNO3", "HCCN", "HCHO,formaldehy", "HCOOH", + "NH", "NH2", "NH2OH", "NCN", "N2H2", "NH2NO2", "N2H4", "H2O2", + "(HCOOH)2", "C6H6(L)", "C7H8(L)", "C8H18(L),n-octa", "Jet-A(L)", "H2O(s)", "H2O(L)", +] + @pytest.mark.rp1311 def test_example1_tp_equilibrium(cea_module): @@ -171,6 +194,52 @@ def test_example4_uv_equilibrium_omit(cea_module): assert soln.gamma_s == pytest.approx(1.225707104935455, rel=1e-5) +@pytest.mark.rp1311 +def test_example5_hp_custom_reactant(cea_module): + """RP-1311 Example 5: HP equilibrium with one custom reactant definition.""" + cea = cea_module + pressures = np.array([34.473652, 17.236826, 8.618413, 3.447365, 0.344737], dtype=np.float64) + weights = np.array([0.7206, 0.1858, 0.09, 0.002, 0.0016], dtype=np.float64) + T_reac = np.array([298.15, 298.15, 298.15, 298.15, 298.15], dtype=np.float64) + + chos_binder_mw_kg_per_mol = 14.6652984484e-3 + chos_binder_h_si = cea.units.cal_to_joule(-2999.082) / chos_binder_mw_kg_per_mol + + reactants = [ + "NH4CLO4(I)", + cea.Reactant( + name="CHOS-Binder", + formula={"C": 1.0, "H": 1.86955, "O": 0.031256, "S": 0.008415}, + molecular_weight=14.6652984484, + enthalpy=chos_binder_h_si, + temperature=298.15, + ), + "AL(cr)", + "MgO(cr)", + "H2O(L)", + ] + + reac = cea.Mixture(reactants) + prod = cea.Mixture(reactants, products_from_reactants=True, omit=_OMIT_5) + solver = cea.EqSolver(prod, reactants=reac) + soln = cea.EqSolution(solver) + h0 = reac.calc_property(cea.ENTHALPY, weights, T_reac) + + solver.solve(soln, cea.HP, h0 / cea.R, pressures[0], weights) + assert soln.converged + assert soln.T == pytest.approx(2722.99, rel=5e-5) + assert soln.gamma_s == pytest.approx(1.1928, rel=5e-5) + assert cea.units.joule_to_cal(soln.cp_eq) == pytest.approx(0.5781, rel=5e-4) + assert soln.mole_fractions["CO"] == pytest.approx(0.26456, rel=1e-4) + + solver.solve(soln, cea.HP, h0 / cea.R, pressures[-1], weights) + assert soln.converged + assert soln.T == pytest.approx(2540.82, rel=5e-5) + assert soln.gamma_s == pytest.approx(1.1494, rel=5e-5) + assert cea.units.joule_to_cal(soln.cp_eq) == pytest.approx(0.9585, rel=5e-4) + assert soln.mole_fractions["CO"] == pytest.approx(0.25896, rel=1e-4) + + @pytest.mark.rp1311 def test_example6_detonation_transport(cea_module): """RP-1311 Example 6: CJ detonation + transport, H2/O2 (phi=1.0).""" diff --git a/source/bind/python/tests/test_strings.py b/source/bind/python/tests/test_strings.py index 03d6f2e..d1151b2 100644 --- a/source/bind/python/tests/test_strings.py +++ b/source/bind/python/tests/test_strings.py @@ -7,14 +7,3 @@ def test_mixture_accepts_python_strings(cea_module) -> None: weights = mix.moles_to_weights(np.array([2.0, 1.0, 0.5], dtype=np.float64)) density = mix.calc_property(cea_module.DENSITY, weights, temperature=3000.0, pressure=1.0) assert np.isfinite(density) - - -def test_bytes_inputs_remain_supported(cea_module) -> None: - text_mix = cea_module.Mixture(["H2", "O2"]) - bytes_mix = cea_module.Mixture([b"H2", b"O2"]) - weights = np.array([0.4, 0.6], dtype=np.float64) - - # ENERGY does not require an explicit pressure argument, keeping the comparison simple. - energy_text = text_mix.calc_property(cea_module.ENERGY, weights, temperature=3200.0) - energy_bytes = bytes_mix.calc_property(cea_module.ENERGY, weights, temperature=3200.0) - assert energy_text == pytest.approx(energy_bytes)