Made pi_p optional in RocketSolver_solve (#36)

markleader · web-flow · commit fda27a8e0960 · 2026-02-19T09:05:36.000-05:00
* made pi_p optional in RocketSolver_solve
diff --git a/.gitignore b/.gitignore
@@ -5,8 +5,10 @@ docs/doctrees/
 docs/doxygen/
 docs/html/
 docs/source/_build/
-docs/*.aux
-docs/*.log
+*.aux
+*.log
+*.toc
+*.out
 extern/gfe
 extern/install
 *.swp
@@ -19,15 +21,20 @@ extern/install
 *.pdf
 *.bat
 *.lib
+*.dylib
+*.so
+*.dll
+*.pyc
+*.pyo
+*.pyd
+*.egg-info/
 !sampleThe/*.inp
 !test/main_interface/*.inp
 !test/main_interface/reference_output/*.out
 __pycache__/
 *.py[cod]
 *$py.class
-source/bind/excel/*.dylib
-source/bind/excel/*.so
-source/bind/excel/*.dll
 source/bind/excel/~$*.xlsm
 source/bind/excel/*.lib
 source/bind/python/samples
+source/bind/python/cea/samples/*.tex
diff --git a/source/bind/c/bindc.F90 b/source/bind/c/bindc.F90
diff --git a/source/bind/c/cea.h b/source/bind/c/cea.h
@@ -432,6 +432,7 @@ extern "C"
       cea_rocket_solution soln,
       const cea_array weights,
       const cea_real pc,
+      // Optional: ignored when n_pi_p == 0 (pi_p may be NULL)
       const cea_array pi_p,
       const cea_int n_pi_p,
       const cea_array subar,
@@ -449,6 +450,7 @@ extern "C"
       cea_rocket_solution soln,
       const cea_array weights,
       const cea_real pc,
+      // Optional: ignored when n_pi_p == 0 (pi_p may be NULL)
       const cea_array pi_p,
       const cea_int n_pi_p,
       const cea_array subar,
diff --git a/source/bind/python/CEA.pyx b/source/bind/python/CEA.pyx
@@ -1830,7 +1830,7 @@ cdef class RocketSolver:
         Raises
         ------
         ValueError
-            If exit conditions not specified, or conflicting parameters provided
+            If conflicting parameters are provided
         """
         cdef cea_err ierr
         cdef int nweights = <int>len(weights)
@@ -1880,12 +1880,15 @@ cdef class RocketSolver:
             else:
                 raise ValueError("RocketSolver.solve: supar must be a list, np.ndarray, or float")
 
-        if ((npi_p + nsubar + nsupar) == 0):
-            raise ValueError("Exit condition not specified.")
-
-        cdef cea_array pi_p_c = <cea_array>malloc(npi_p * sizeof(double))
-        cdef cea_array subar_c = <cea_array>malloc(nsubar * sizeof(double))
-        cdef cea_array supar_c = <cea_array>malloc(nsupar * sizeof(double))
+        cdef cea_array pi_p_c = NULL
+        cdef cea_array subar_c = NULL
+        cdef cea_array supar_c = NULL
+        if npi_p > 0:
+            pi_p_c = <cea_array>malloc(npi_p * sizeof(double))
+        if nsubar > 0:
+            subar_c = <cea_array>malloc(nsubar * sizeof(double))
+        if nsupar > 0:
+            supar_c = <cea_array>malloc(nsupar * sizeof(double))
         if (npi_p > 0 and pi_p_c == NULL) or (nsubar > 0 and subar_c == NULL) or (nsupar > 0 and supar_c == NULL):
             free(wts)
             if pi_p_c != NULL:
diff --git a/source/bind/python/cea_def.pxd b/source/bind/python/cea_def.pxd
@@ -298,12 +298,14 @@ cdef extern from "cea.h":
     cpdef cea_err cea_rocket_solver_create_with_options(cea_rocket_solver *solver, const cea_mixture products,
                                                         const cea_solver_opts opts)
     cpdef cea_err cea_rocket_solver_destroy(cea_rocket_solver *solver)
+    # pi_p is optional when n_pi_p == 0.
     cpdef cea_err cea_rocket_solver_solve_iac(const cea_rocket_solver solver, cea_rocket_solution soln,
                                               const cea_array weights, const cea_real pc, const cea_array pi_p,
                                               const cea_int n_pi_p, const cea_array subar, const cea_int nsubar,
                                               const cea_array supar, const cea_int nsupar, const cea_int n_frz,
                                               const cea_real hc_or_tc, const cea_bool use_hc,
                                               const cea_real tc_est, const cea_bool use_tc_est)
+    # pi_p is optional when n_pi_p == 0.
     cpdef cea_err cea_rocket_solver_solve_fac(const cea_rocket_solver solver, cea_rocket_solution soln,
                                               const cea_array weights, const cea_real pc, const cea_array pi_p,
                                               const cea_int n_pi_p, const cea_array subar, const cea_int nsubar,
diff --git a/source/bind/python/tests/test_rocket_smoke.py b/source/bind/python/tests/test_rocket_smoke.py
@@ -31,3 +31,29 @@ def test_rocket_solver_smoke():
         check=False,
     )
     assert result.returncode == 0, result.stderr + result.stdout
+
+
+@pytest.mark.smoke
+def test_rocket_solver_smoke_without_exit_conditions():
+    code = textwrap.dedent(
+        """
+        import numpy as np
+        import cea
+
+        reactants = ["H2", "O2"]
+        reactants_mix = cea.Mixture(reactants)
+        products_mix = cea.Mixture(reactants, products_from_reactants=True)
+        weights = reactants_mix.moles_to_weights(np.array([2.0, 1.0], dtype=np.float64))
+        solver = cea.RocketSolver(products_mix, reactants=reactants_mix)
+        soln = cea.RocketSolution(solver)
+
+        solver.solve(soln, weights, pc=10.0, tc=3000.0)
+        """
+    )
+    result = subprocess.run(
+        [sys.executable, "-c", code],
+        capture_output=True,
+        text=True,
+        check=False,
+    )
+    assert result.returncode == 0, result.stderr + result.stdout
diff --git a/source/bind/python/tests/test_validation.py b/source/bind/python/tests/test_validation.py
@@ -16,14 +16,6 @@ def test_eqsolver_rejects_non_mixture_reactants(cea_module):
         cea_module.EqSolver(products, reactants=["H2", "O2"])
 
 
-def test_rocket_solver_requires_exit_condition(cea_module):
-    reactants_mix, products_mix, weights = _basic_reactants(cea_module)
-    solver = cea_module.RocketSolver(products_mix, reactants=reactants_mix)
-    soln = cea_module.RocketSolution(solver)
-    with pytest.raises(ValueError):
-        solver.solve(soln, weights, pc=10.0, tc=3000.0)
-
-
 def test_rocket_solver_rejects_conflicting_hc_tc(cea_module):
     reactants_mix, products_mix, weights = _basic_reactants(cea_module)
     solver = cea_module.RocketSolver(products_mix, reactants=reactants_mix)
diff --git a/source/rocket.f90 b/source/rocket.f90
@@ -783,7 +783,7 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar,
         type(RocketSolution), intent(inout) :: soln
         real(dp), intent(in) :: reactant_weights(:)
         real(dp), intent(in) :: pc                   ! Chamber pressure [Pa]
-        real(dp), intent(in) :: pi_p(:)              ! Ratio of chamber pressure to exit pressure [unitless]
+        real(dp), intent(in), optional :: pi_p(:)    ! Ratio of chamber pressure to exit pressure [unitless]
         real(dp), intent(in), optional :: subar(:)   ! Subsonic area ratio [unitless]
         real(dp), intent(in), optional :: supar(:)   ! Supersonic area ratio [unitless]
         integer,  intent(in), optional :: n_frz      ! Station where the composition should be frozen
@@ -806,7 +806,8 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar,
         call log_debug("Starting rocket IAC solve")
 
         ! Set the total number of evaluation points
-        num_pts = 2 + size(pi_p)
+        num_pts = 2  ! infinity + throat
+        if (present(pi_p)) num_pts = num_pts + size(pi_p)
         if (present(subar)) num_pts = num_pts + size(subar)
         if (present(supar)) num_pts = num_pts + size(supar)
 
@@ -893,14 +894,19 @@ subroutine RocketSolver_solve_iac(self, soln, reactant_weights, pc, pi_p, subar,
         ! -----------------------------------------------
         ! Exit conditions: pressure ratio
         ! -----------------------------------------------
-        idx = 3
+        if (present(pi_p)) then
+            idx = 3
 
-        if (frozen .and. idx > n_frz_) then
-            call self%solve_pi_p_frozen(soln, idx, n_frz, pc, pi_p, h_inf, 1)
+            if (frozen .and. idx > n_frz_) then
+                call self%solve_pi_p_frozen(soln, idx, n_frz, pc, pi_p, h_inf, 1)
+            else
+                call self%solve_pi_p(soln, idx, pc, pi_p, h_inf, state1, reactant_weights)
+            end if
+            if (.not. soln%converged) return
         else
-            call self%solve_pi_p(soln, idx, pc, pi_p, h_inf, state1, reactant_weights)
+            ! If pi_p not present, idx stays at 3 for subsequent sections
+            idx = 3
         end if
-        if (.not. soln%converged) return
 
         ! -----------------------------------------------
         ! Exit conditions: subsonic area ratio
@@ -947,7 +953,7 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar,
         type(RocketSolution), intent(inout) :: soln
         real(dp), intent(in) :: reactant_weights(:)
         real(dp), intent(in) :: pc                   ! Injector pressure [Pa]
-        real(dp), intent(in) :: pi_p(:)              ! Ratio of chamber pressure to exit pressure [unitless]
+        real(dp), intent(in), optional :: pi_p(:)    ! Ratio of chamber pressure to exit pressure [unitless]
         real(dp), intent(in), optional :: subar(:)   ! Subsonic area ratio [unitless]
         real(dp), intent(in), optional :: supar(:)   ! Supersonic area ratio [unitless]
         real(dp), intent(in), optional :: ac_at      ! Contraction ratio: ratio of finite chamber area to throat area, Ac/At (FAC only) [unitless]
@@ -994,7 +1000,8 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar,
         call log_debug("Starting rocket FAC solve")
 
         ! Set the total number of evaluation points
-        num_pts = 4 + size(pi_p)
+        num_pts = 4  ! injector + infinity + combustor + throat
+        if (present(pi_p)) num_pts = num_pts + size(pi_p)
         if (present(subar)) num_pts = num_pts + size(subar)
         if (present(supar)) num_pts = num_pts + size(supar)
 
@@ -1206,14 +1213,19 @@ subroutine RocketSolver_solve_fac(self, soln, reactant_weights, pc, pi_p, subar,
         ! -----------------------------------------------
         ! Exit conditions: pressure ratio
         ! -----------------------------------------------
-        idx = 5
+        if (present(pi_p)) then
+            idx = 5
 
-        if (frozen .and. idx > n_frz_) then
-            call self%solve_pi_p_frozen(soln, idx, n_frz_, pc, pi_p, h_inj, 2)
+            if (frozen .and. idx > n_frz_) then
+                call self%solve_pi_p_frozen(soln, idx, n_frz_, pc, pi_p, h_inj, 2)
+            else
+                call self%solve_pi_p(soln, idx, pc, pi_p, h_inj, S_ref, reactant_weights)
+            end if
+            if (.not. soln%converged) return
         else
-            call self%solve_pi_p(soln, idx, pc, pi_p, h_inj, S_ref, reactant_weights)
+            ! If pi_p not present, idx stays at 5 for subsequent sections
+            idx = 5
         end if
-        if (.not. soln%converged) return
 
         ! -----------------------------------------------
         ! Exit conditions: subsonic area ratio
@@ -1264,7 +1276,7 @@ function RocketSolver_solve(self, reactant_weights, pc, pi_p, fac, subar, supar,
         class(RocketSolver), intent(in) :: self
         real(dp), intent(in) :: reactant_weights(:)  !
         real(dp), intent(in) :: pc                   ! Chamber pressure [bar]
-        real(dp), intent(in) :: pi_p(:)              ! Ratio of chamber pressure to exit pressure [unitless]
+        real(dp), intent(in), optional :: pi_p(:)    ! Ratio of chamber pressure to exit pressure [unitless]
         logical,  intent(in), optional :: fac        ! Finite-area combustor flag
         real(dp), intent(in), optional :: subar(:)   ! Subsonic area ratio [unitless]
         real(dp), intent(in), optional :: supar(:)   ! Supersonic area ratio [unitless]
@@ -1289,10 +1301,11 @@ function RocketSolver_solve(self, reactant_weights, pc, pi_p, fac, subar, supar,
 
         ! Call either the FAC or IAC solver
         if (fac_) then
-            call self%solve_fac(soln, reactant_weights, pc, pi_p, subar, supar, ac_at, mdot, n_frz, tc_est, hc, tc)
+            call self%solve_fac(soln, reactant_weights, pc, pi_p=pi_p, subar=subar, supar=supar, &
+                                ac_at=ac_at, mdot=mdot, n_frz=n_frz, tc_est=tc_est, hc=hc, tc=tc)
         else  ! IAC
-            call self%solve_iac(soln, reactant_weights, pc, pi_p, subar=subar, supar=supar, n_frz=n_frz, &
-                tc_est=tc_est, hc=hc, tc=tc)
+            call self%solve_iac(soln, reactant_weights, pc, pi_p=pi_p, subar=subar, supar=supar, &
+                                n_frz=n_frz, tc_est=tc_est, hc=hc, tc=tc)
         end if
 
     end function
diff --git a/source/rocket_test.pf b/source/rocket_test.pf
