[bug] Equilibrium solver doesn't converge for detonation; CEA2 converges

[dfnca]

Summary
Some detonation problems that run fine in CEA2 fail with an unconverged equilibrium solver.

Reproduction

• Interface: confirmed with legacy, python interfaces
• Minimal input / steps (attach .inp if relevant):
  reac
  oxid O2 wt=1
  fuel C2H4 wt=1
  prob det phi=0.4 t,k=298.15 p,bar=4
  outp short
  end

Expected behavior
Compute a solution consistent with CEA2, within solver tolerances.

Actual behavior
ERROR: DETON: Equilibrium solver did not converge
WARNING: Singular update matrix encountered at iteration 2
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
WARNING: Singular update matrix encountered at iteration 1
CRITICAL: EqSolver_solve: Too many singular matrices encountered.
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_OVERFLOW_FLAG IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
STOP 1

Environment

• OS: Ubuntu + gfortran 13.3, Win 11 + Intel oneAPI 2025.3
• CEA version/commit: 3.0.3 / 358065a

Additional context
There is some weirdness:

1. Building with gfortran using the core preset, the example problem gives multiple "ERROR: DETON: Equilibrium solver did not converge" messages but ultimately converges but the solution doesn't quite match CEA2 in some respects.
2. Building with gfortran using the dev preset, the example problem produces the result listed above in "What actually happened". Intel builds on Windows also produce this result, with the exception of the fpe output.
3. The attached testDetIssueHPPT.py replicates the initial guess HP problem and the TP problem that fails to converge on the first pass through the iteration loop in detonation.f90, using values that correspond to the sample problem listed in the Reproduction section. But the python HPPT example converges, unlike detonation.f90. Calling a detonation problem from the python interface produces the same unconverged equilibrium solution as does the legacy interface.
4. Zeroing soln%dln_n and soln%dln_T in the EqSolver_solve initialization section after line 1950 of equilibrium.f90 appears to fix the issue. The sample problem converges with results that match CEA2, as do other problems that produced similar behavior. I'm not necessarily proposing this as a fix because I don't understand why it works.

testDetIssueHPPT.py

[dfnca]

There are other instances where the EQ solver does or does not converge depending on history. The attached python code illustrates an example where condensed phases are present. The example computes TP solutions for a matrix of T, P, and o/f combinations. The basic example does not converge for some T. P, o/f combinations (testEQtable.py), but creating a new instance of EqSolution prior to each solve call converges for all cases (testEQtableResetSoln.py). I presume that this could be addressed by re-initializing some critical EqSolution variables within the solve subroutine, but I haven't tested it.

The detonation, equilibrium, and maybe shock and rocket solver robustness is affected by some variables that are not reiniailized on each new solver call, when reusing a solution object.

testEQtable.py
testEQtableResetSoln.py

[markleader]

Both of these are fixed now by PR #47 - I'll merge it shortly.

Your second comment here was addressed by updating the warm-start logic in equilibrium.f90. There are still some failures in testEQtable.py while all cases in testEQtableResetSoln.py still pass. But the behavior is now consistent with CEA2. I see some non-converging cases when I run this subset of the problems in CEA2:

   problem case=test-table tp  p(bar)=.1,1,10,t(k)=3000,4000,
             o/f=5,20,35
   reac
     fuel= AL  wtfrac=.5
     fuel= C2H4  wtfrac=.5
     oxid= Air wtfrac=1
   output si
   end

But when I run each of those 18 total cases one-by-one in CEA2, they all pass. PR #47 is consistent here.

In CEA2, the solutions use the previous solution as the initial condition for the next solve. This almost always improves the solver robustness and reduces the number of iterations required. There is something about this problem where the previous solution causes the next one to fail, while starting from scratch works fine. The specific fix in the PR saves a seed state for the last converged iteration, so if a EqSolution is passed with converged=.false., the initial conditions revert to the saved seed state - this is consistent with CEA2.

A potential improvement could be: if a problem doesn't converged, and it started from initial conditions different than the default, then reset to default initial conditions and re-try the solve one more time. I'm not ready to commit to that right now, but it could make sense as a future enhancement.