Derivative re-factor and updated Example 7 validation (#110)

* refactored derivatives for cleaner formulation and faster solves

* removed finite differences from smooth derivatives

* updated example7 validation example

* version patch increment

Author: markleader

CHANGELOG.md

@@ -4,6 +4,14 @@ All notable user-visible changes to this project are documented here.
 
 ## [Unreleased]
 
+### Changed
+
+### Fixed
+
+### Added
+
+## [3.3.3] - 2026-08-24
+
 ### Changed
 - Reworked equilibrium total-derivative post-processing as a vector-valued
   response Jacobian and reused one Jacobian factorization across all direct

CMakeLists.txt

@@ -8,7 +8,7 @@ set(CMAKE_MSVC_RUNTIME_LIBRARY "MultiThreaded$<$<CONFIG:Debug>:Debug>")
 set(CMAKE_DISABLE_SOURCE_CHANGES ON)
 set(CMAKE_DISABLE_IN_SOURCE_BUILD ON)
 project(CEA
-    VERSION 3.3.2
+    VERSION 3.3.3
     LANGUAGES Fortran
 )
 

Doxyfile

@@ -48,7 +48,7 @@ PROJECT_NAME           = "CEA"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 3.3.2
+PROJECT_NUMBER         = 3.3.3
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewers a

docs/source/conf.py

@@ -22,7 +22,7 @@
 copyright = ''
 author = 'Mark Leader'
 version = '3.3'
-release = '3.3.2'
+release = '3.3.3'
 
 # -- General configuration ---------------------------------------------------
 # https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration

project.md

@@ -1,6 +1,6 @@
 ---
 project: CEA
-version: 3.3.2
+version: 3.3.3
 src_dir: ./source
 output_dir: ./docs
 project_github: https://github.com/nasa/cea
@@ -18,4 +18,4 @@ search: true
 print_creation_date: false
 extra_filetypes:
 exclude: *.pf
----
+---

source/bind/python/cea/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "3.3.2"
+__version__ = "3.3.3"
 
 # initialize libcea, loading in the default data files
 from cea.lib.libcea import init as libcea_init

source/main.f90

@@ -556,7 +556,7 @@ subroutine run_shock_problem(prob, thermo, solver, solutions, branch_codes)
         type(ShockSolution) :: solution
         real(dp) :: T0, P0, u1, mach1
         real(dp), allocatable :: weights(:)
-        logical :: incident, input_reflected, frozen, equilibrium, incident_frozen, reflected_frozen, reflected, use_mach
+        logical :: input_reflected, frozen, equilibrium, incident_frozen, reflected_frozen, reflected, use_mach
         logical :: stop_after_branch, branch_failed
         integer :: i, k, idx_P, idx_T, npts, num_u1, num_P, num_T
         integer, parameter :: BR_INCD_EQ = 1
@@ -585,21 +585,12 @@ subroutine run_shock_problem(prob, thermo, solver, solutions, branch_codes)
         products = Mixture(thermo, product_names, ions=prob%problem%include_ions)
 
         ! Get the problem flags
-        if (prob%problem%shk_incident) then
-            incident = .true.
-        else
-            incident = .false.
-        end if
-
         if (prob%problem%shk_reflected) then
             input_reflected = .true.
             reflected = .true.
         else
             input_reflected = .false.
             reflected = .false.
-            if (.not. incident) then
-                incident = .true.  ! Incident by default
-            end if
         end if
 
         if (prob%problem%equilibrium) then
@@ -617,29 +608,16 @@ subroutine run_shock_problem(prob, thermo, solver, solutions, branch_codes)
             end if
         end if
 
-        ! Match the SHCK branch sequencing. When the input explicitly requests
-        ! incident-shock output (`inc`), CEA2 reports only that incident branch,
-        ! with the reflected state attached if requested. The extra mixed
-        ! equilibrium/frozen permutations are only explored for reflected-only
-        ! request patterns.
-        if (incident) then
-            npts = 1
-        else
-            npts = 0
-            if (equilibrium) then
-                npts = npts + 1
-                if (input_reflected .and. frozen) then
-                    npts = npts + 1
-                else if (input_reflected) then
-                    ! Equilibrium-only reflected shock prints incident and reflected states together.
-                end if
-            end if
-            if (frozen) then
-                npts = npts + 1
-                if (input_reflected .and. equilibrium) npts = npts + 1
-            end if
-            if (npts == 0) npts = 1
-        end if
+        ! Match the SHCK branch sequencing. Each requested composition model
+        ! produces an incident branch. A reflected request applies each
+        ! requested model to every incident branch, so requesting both
+        ! equilibrium and frozen models produces four incident/reflected
+        ! combinations.
+        npts = 0
+        if (equilibrium) npts = npts + 1
+        if (frozen) npts = npts + 1
+        if (npts == 0) npts = 1
+        if (input_reflected) npts = npts*npts
 
         ! Get the loop sizes
         num_u1 = 1
@@ -740,7 +718,7 @@ subroutine run_shock_problem(prob, thermo, solver, solutions, branch_codes)
                 if (.not. solution%converged .and. input_reflected) branch_failed = .true.
             end do
 
-            if ((.not. incident) .and. input_reflected .and. frozen) then
+            if (input_reflected .and. frozen) then
                 k = k + 1
                 branch_codes(k) = BR_REFL_EQ_FROM_EQ
                 incident_frozen = .false.
@@ -784,7 +762,7 @@ subroutine run_shock_problem(prob, thermo, solver, solutions, branch_codes)
             if (input_reflected .and. branch_failed) stop_after_branch = .true.
         end if
 
-        if (frozen .and. .not. stop_after_branch .and. (.not. incident .or. .not. equilibrium)) then
+        if (frozen .and. .not. stop_after_branch) then
             k = k + 1
             if (input_reflected) then
                 branch_codes(k) = BR_INCD_FRZ_REFL_FRZ
@@ -1101,7 +1079,9 @@ subroutine shock_output(ioout, prob, solver, solutions, branch_codes)
             have_reflected_state = .false.
             do i = 1, m
                 if (solutions(i, 1, k)%eq_soln(2)%T > 0.0d0) have_incident_state = .true.
-                if (solutions(i, 1, k)%eq_soln(3)%T > 0.0d0) have_reflected_state = .true.
+                if (size(solutions(i, 1, k)%eq_soln) >= 3) then
+                    if (solutions(i, 1, k)%eq_soln(3)%T > 0.0d0) have_reflected_state = .true.
+                end if
             end do
             if (.not. have_incident_state) then
                 write_incd_frz = .false.

test/main_interface/example7.inp

@@ -6,22 +6,21 @@
 !          reactants at the temperature given (300 K) using the thermo.lib
 !          coefficients.
 !     (c)  Reactants are given in moles (moles = ...).
-!     (d)  Initial gas pressures are 10 and 20 mm Hg (p,mmhg=10,20,)
-!     (e)  Seven initial gas velocities are assigned (u1=1000,1100,1200,
-!          1250,1300,1350,1400,).  Note units of u1 are always m/s.
+!     (d)  The initial gas pressure is 10 mm Hg (p,mmhg=10,).
+!     (e)  Six initial gas velocities are assigned (u1=1100,1200,1250,
+!          1300,1350,1400,).  Note units of u1 are always m/s.
 !     (f)  Equilibrium calculations are to be performed for incident shock
 !          conditions (incd  eql).
 !     (g)  Frozen calculations are to be performed for incident shock
 !          conditions (incd  froz).
-!     (h)  No 'outp' dataset is given since the default values of the
-!          the parameters have the desired  values (e.g. SI units).
+!     (h)  Long-format output is requested for validation.
 
 reac   name= H2  moles= 0.050  t(k) 300.00
        name= O2  moles= 0.050  t(k) 300.00
        name= Ar  moles= 0.900  t(k) 300.00
 
 problem  case=7 p,mmhg=10, shock u1=1100,1200,1250,1300,1350,1400,
-  incd refl eql
+  incd eql froz
 output long
 
-end
+end

test/main_interface/parse_output.py

@@ -75,6 +75,8 @@ def parse_output(fname):
     trace_next = False  # Trace species get parsed on the next line
     frozen = False  # Transport property type
     count_blank = 0  # Number of blank lines in a row
+    shock_mode = "unknown"
+    shock_state = "unknown"
 
     # Parse the file
     f = open(fname, "r")
@@ -94,6 +96,17 @@ def parse_output(fname):
         if vals[0] == "!":
             continue
 
+        if vals[0].startswith("WARNING") or vals[0] == "ANSWERS":
+            continue
+
+        upper_line = line.upper()
+        if "EQUILIBRIUM COMPOSITION FOR INCIDENT SHOCKED" in upper_line:
+            shock_mode = "equilibrium"
+            continue
+        if "FROZEN COMPOSITION FOR INCIDENT SHOCKED" in upper_line:
+            shock_mode = "frozen"
+            continue
+
         # Check for data headings
         if vals[0] == "THERMODYNAMIC":
             if len(vals) > 1:
@@ -105,8 +118,18 @@ def parse_output(fname):
             parse_flag = Parse.THERMO
             continue
 
-        if (vals[0] == "INITIAL") or (vals[0] == "SHOCKED"):
+        if vals[0] == "INITIAL":
             parse_flag = Parse.SHOCK
+            shock_state = "initial"
+            continue
+
+        if vals[0] == "SHOCKED":
+            parse_flag = Parse.SHOCK
+            shock_state = "reflected" if "REFLECTED" in upper_line else "incident"
+            if "EQUILIBRIUM" in upper_line:
+                shock_mode = "equilibrium"
+            elif "FROZEN" in upper_line:
+                shock_mode = "frozen"
             continue
 
         if (vals[0] == "UNBURNED") or (vals[0] == "BURNED"):
@@ -200,11 +223,13 @@ def parse_output(fname):
             for i in range(len(vals)):
                 vals[i] = parse_numeric(vals[i])
 
+            property_key = f"{shock_mode}:{shock_state}:{property_name}"
+
             # Check if these values are already in the results dict; we are just adding more points
-            if (property_name in shock):  # Adding new entries to existing values
-                shock[property_name]["vals"] = np.append(shock[property_name]["vals"], np.array(vals))
+            if (property_key in shock):  # Adding new entries to existing values
+                shock[property_key]["vals"] = np.append(shock[property_key]["vals"], np.array(vals))
             else:  # Adding a new property type
-                shock[property_name] = {"units":unit_name, "vals":np.array(vals)}
+                shock[property_key] = {"units":unit_name, "vals":np.array(vals)}
 
         if parse_flag == Parse.DETON:
             # Parse the name and the units
@@ -250,6 +275,12 @@ def parse_output(fname):
                 elif (name_vals[0] == "SHOCKED"):
                     trace_next = False
                     parse_flag = Parse.SHOCK
+                    upper_line = line.upper()
+                    shock_state = "reflected" if "REFLECTED" in upper_line else "incident"
+                    if "EQUILIBRIUM" in upper_line:
+                        shock_mode = "equilibrium"
+                    elif "FROZEN" in upper_line:
+                        shock_mode = "frozen"
                     continue
                 else:
                     for species_name in name_vals: