fix: seed FAC supersonic-area-ratio Newton iteration from throat gamma_s - #28
Merged
Merged
Conversation
Uses soln%throat_idx instead of hardcoded station 2, matching RP-1311 Eq. 6.22 and the pattern already used elsewhere in rocket.f90. For IAC problems throat_idx=2 so behavior is unchanged; for FAC problems throat_idx=4, so this now seeds from the throat instead of infinity. Fixes #10
djkees
marked this pull request as ready for review
September 1, 2026 17:25
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
RocketSolver_solve_suparandRocketSolver_solve_supar_frozen(source/rocket.f90:896,990) hardcodesoln%eq_partials(2)%gamma_sas the throat-station seed for the FAC/IAC supersonic-area-ratio Newton iteration. For IAC problemsthroat_idx=2so this happened to be correct, but for FAC problemsthroat_idx=4(stations: 1=injector, 2=infinity, 3=combustor end, 4=throat), so the infinity station'sgamma_swas used instead of the throat's.Changes
Replace the hardcoded
soln%eq_partials(2)%gamma_swithsoln%eq_partials(soln%throat_idx)%gamma_sat both sites, matching the pattern already used elsewhere in the file (rocket.f90:661,724) and matching RP-1311 Eq. 6.22, which explicitly states "the value of γₛ is that determined for throat conditions."Testing
ctest --test-dir build-dev --output-on-failure— all 15 tests pass, including the FAC-moderocket_test_suite.test_rocket_fac_only_area_ratios.python test/main_interface/test_main.py— 14/14 RP-1311 book-example comparisons pass against published reference values (rtol=1e-2), includingexample9/example10(FAC,supar=25,50,75), which exercise exactly this code path.example9/example10are unchanged before/after the fix (2, 2, 3 iterations) — for these particular cases the infinity- and throat-stationgamma_shappen to be close enough that the seed choice doesn't change convergence speed.gamma_sgap — H2/O2 at O/F 3–20, RP-1/O2, chamber pressures 10–68 atm, andsuparfrom 2.01 to 150 (largest gap found: ~1.03% relative, at O/F=3, pc=68 atm, supar=100,gamma2=1.23075vsgammaT=1.24338). Pre-fix and post-fix iteration counts and converged results were identical in every case tried. So while this is a genuine correctness fix — it replaces a physically wrong index with the one RP-1311 Eq. 6.22 specifies — no input found so far (across 8 cases total) demonstrates it changing convergence speed or outcome. The "could hit the 10-iteration cap unconverged" motivation is a theoretical concern based on the code being wrong, not something reproduced numerically; flagging that distinction here rather than overstating the empirical case for this being a robustness fix.Compatibility / Numerical behavior
Only the initial-guess seed for a Newton iteration changes; the converged fixed point is unaffected for a converged case (see benchmark results above). Only FAC problems using
supar >= 2.0are affected — IAC problems (throat_idx=2) see no change at all.Drafted with Claude's assistance
rocket.f90, cross-checked against RP-1311 (NASA-CR-19950013764) Eq. 6.21/6.22, which was fetched and read directly (see screenshot in issue history).ctestsuite and the RP-1311 benchmark comparison script (test/main_interface/test_main.py) against published reference values, not just re-stated from the issue.#10