Skip to content

Fix mixture entropy/gibbs_energy single-multi scaling mismatches - #132

Open
djkees wants to merge 1 commit into
nasa:mainfrom
djkees:up/mixture-entropy-gibbs-multi-scaling
Open

Fix mixture entropy/gibbs_energy single-multi scaling mismatches#132
djkees wants to merge 1 commit into
nasa:mainfrom
djkees:up/mixture-entropy-gibbs-multi-scaling

Conversation

@djkees

@djkees djkees commented Sep 3, 2026

Copy link
Copy Markdown
Contributor

Summary

Fixes three related numerical bugs in the single-temperature vs multi-temperature overloads of Mixture%calc_entropy and Mixture%calc_gibbs_energy (source/mixture.f90), so that calling either overload with equivalent physical inputs now gives the same, physically correct result.

djkees#66

Changes

  • mixture_calc_entropy_multi divided its result by an extra 1.d3 that mixture_calc_entropy_single did not (the /1.d3 in _single was already commented out). Removed the stray division so both overloads scale identically.
  • mixture_calc_gibbs_energy_multi's nj was never normalized by sum(weights), unlike every other _multi routine (calc_entropy_multi, calc_enthalpy_multi) and unlike calc_gibbs_energy_single (built from calc_enthalpy/calc_entropy, both of which do normalize). Added the missing /total_weight normalization.
  • A third bug found while verifying the above two fixes against the RP-1311 Example 14 reference data already present (but disabled) in mixture_test.pf: both calc_entropy overloads compute log(pressure/...) using the raw pressure argument directly, without dividing by std_pressure (the 1 bar reference pressure) first. calc_gibbs_energy_multi already does this correctly (pr = pressures(j)/std_pressure), but neither calc_entropy overload did. Without this fix, both overloads agree with each other but disagree with the correct physical value. Added /std_pressure to the pressure term in both overloads.
  • Re-enabled the two long-disabled @assertRelativelyEqual assertions in mixture_test.pf's test_calc_properties (calc_gibbs_energy and calc_entropy, both the single- and multi-temperature call forms) — they've been commented out since the initial public release, and now pass with the fix.

Testing

  • Not run via ctest/pFUnit locally: pFUnit's own CMake (vendored via GFE) unconditionally excludes RegexFilter.F90 on Windows (if(NOT WIN32) in pFUnit/src/funit/core/CMakeLists.txt) while FUnit.F90's main driver unconditionally does use pf_RegexFilter, so pFUnit itself fails to build on this Windows/MinGW dev environment regardless of this change. CI's pfunit-tests job builds pFUnit on Ubuntu and should exercise mixture_test.pf and equilibrium_test.pf (including test_sv_derivatives, which uses the now-corrected calc_entropy to build an initial-guess constraint value for an entropy/volume-constrained solve) — please confirm CI is green before merging, especially that test.
  • Verified numerically outside of pFUnit: built cea_core from this branch and linked a standalone Fortran driver against libcea_core.a that replicates test_calc_properties' exact inputs (RP-1311 Example Case 14 mixture: H2O/O2/H2O(L), T=304K, P=5066.25 Pa). Both calc_entropy overloads now produce 9463.0131180265817 and both calc_gibbs_energy overloads now produce -14830015.9..., matching the reference values 9.46301312d+03 and -1.48300159E+07 in mixture_test.pf to 10 significant figures, and matching each other (single vs. multi) exactly to solver tolerance.

Compatibility / Numerical behavior

  • No expected changes to numerical results
  • Expected changes (explain and provide validation)

This intentionally changes the numerical output of Mixture%calc_entropy and Mixture%calc_gibbs_energy (both single- and multi-temperature overloads). Per the audit trail above: (1) neither the _multi overloads (which were never exercised by any real caller — see issue #66) nor calc_entropy_single/calc_gibbs_energy_single (used internally by calc_gibbs_energy_single itself, but not otherwise called in production code outside of this module and one test) had validated numerical behavior to preserve; the only reference data available (mixture_test.pf's RP-1311-sourced expected values) was disabled specifically because these functions didn't produce correct output. (2) The one live caller of calc_entropy outside this module, equilibrium_test.pf's test_sv_derivatives, uses it only to construct an initial-guess entropy value fed into a self-consistent constant-entropy/volume equilibrium solve, not to check against a hardcoded expected value — this should remain valid under the corrected calculation, but I could not confirm convergence locally (see Testing above).


Drafted with Claude's assistance

  • The three bugs and their fixes were derived by direct read of source/mixture.f90:833-967, cross-checked against the normalization pattern used consistently by every other _multi/_single pair in the same file, and against calc_gibbs_energy_multi's existing (correct) std_pressure handling.
  • Independently verified numerically (not taken on faith) by building cea_core from this branch and linking a standalone driver program against it that reproduces mixture_test.pf's test_calc_properties inputs exactly; results matched the RP-1311 reference values to 10 significant figures for both the single- and multi-temperature overloads, both before/after comparison for each fix and cross-checked against a hand-derived formula (S = (H-G)/T using the already-correct calc_gibbs_energy_multi).
  • Was not able to run the repository's own pFUnit test suite locally due to an unrelated, pre-existing pFUnit-on-Windows build issue (see Testing above) — CI should be checked before merging.

Update

The pFUnit-on-Windows build issue referenced in Testing above (tracked as #163) now has a fix: djkees#167. With that fix, cea_core_test was built and run locally on this branch for the first time: ctest -R cea_core_test -V124/124 tests passed, including mixture_test_suite.test_calc_properties (the re-enabled calc_entropy/calc_gibbs_energy assertions this PR's fix targets) and equilibrium_test_suite.test_sv_derivatives (the one live caller of calc_entropy outside this module, flagged above as unconfirmed) — both green.

Also worth noting: this fork's GitHub Actions checks don't appear to run at all (gh pr checks reports nothing on this PR despite .github/workflows/basic_build.yml existing), so the local run above is currently the only automated test evidence for this PR beyond manual review.

entropy_multi had a stray /1.d3, gibbs_energy_multi's nj was never
normalized by sum(weights), and both entropy overloads were missing
/std_pressure in their pressure log term.

Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
(cherry picked from commit 726abee)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant