Steady-state reactor network solver

[speth]

Changes proposed in this pull request

• Generalize the 1D flame solver into a steady-state solver for generic ODE/DAE systems. This is the new SteadyStateSolver class, which is also used as a base class for OneDim.
• Use this solver to implement a new ReactorNet.solveSteady method as an alternative to the time-marching approach used in the Python-only advance_to_steady_state method.
• Make steady-state solver more robust by allowing it to continue in time-stepping mode after encountering certain errors while trying to solve the steady-state problem. These include both Jacobian factorization errors and errors evaluating the residual function with a "bad" state vector.

This solver isn't currently a full replacement for advance_to_steady_state, as there are some complicated interactions with equations that don't translate directly into the steady-state problem, implicit constraints that are automatically satisfied when integrating the normal ODEs but are not currently imposed in the steady solver, and possibly other mathematical issues that I haven't worked out yet.

A simple example of the first case is the equation for the volume of a control volume reactor, where the ODE reduces to
$dV/dt = 0$. This is fine for an ODE solver, but for the steady solver, it is a problem because there's no information that can be used to solve for $V$. Instead, this equation needs to be replaced with the algebraic equation $V(t) - V_0 = 0$. This is handled in the implementation here by a new steadyConstraints method of each reactor class that indicates which equations need to be transformed into such constraints.

Unfortunately, this only works in the case where the constraint corresponds to an existing equation. Other cases, such as the constraint that the mass of a constant pressure reactor should be constant work fine in the case of some formulations (ConstPressureReactor, where mass is one of the state variables) but not others (ConstPressureMoleReactor, where the mass is a derived quantity based on the moles of each species).

The cases where the solver does not work have been configured to raise exceptions when initializing the steady-state solver, and are indicated in docstring for the solveSteady method. I think we will eventually be able to resolve at least some of these limitations, though this will require a somewhat more flexible approach to modifying the set of equations to be solved.

If applicable, fill in the issue number this pull request is fixing

• Supersedes Steady-State Solver for ZeroD Reactors #1021.
• Completes Implement a steady-state solver for 0-D systems enhancements#31

If applicable, provide an example illustrating new features this pull request is introducing

This solver seems to work particularly well for cases like combustor.py, where the interest is in solving for a number of nearby steady-state conditions:

residence_time = 0.1  # starting residence time
while combustor.T > 500:
    sim.initial_time = 0.0  # reset the integrator
    sim.solve_steady()
    print('tres = {:.2e}; T = {:.1f}'.format(residence_time, combustor.T))
    states.append(combustor.thermo.state, tres=residence_time)
    residence_time *= 0.9  # decrease the residence time for the next iteration

At least in this case, it is a drop-in replacement for sim.advance_to_steady_state.

Checklist

• The pull request includes a clear description of this code change
• Commit messages have short titles and reference relevant issues
• Build passes (scons build & scons test) and unit tests address code coverage
• Style & formatting of contributed code follows contributing guidelines
• The pull request is ready for review

[bryanwweber]

This is so cool Ray, thanks for picking this work up and pushing it through 🎉🎉🎉🎉

[codecov]

Codecov Report

Attention: Patch coverage is 76.00000% with 138 lines in your changes missing coverage. Please review.

Project coverage is 74.22%. Comparing base (305d5a1) to head (d848a9d).

Files with missing lines	Patch %	Lines
src/numerics/SteadyStateSystem.cpp	83.33%	22 Missing and 5 partials ⚠️
src/zeroD/ReactorNet.cpp	85.12%	7 Missing and 11 partials ⚠️
src/oneD/MultiNewton.cpp	74.57%	12 Missing and 3 partials ⚠️
include/cantera/numerics/SteadyStateSystem.h	65.71%	12 Missing ⚠️
src/oneD/OneDim.cpp	73.33%	12 Missing ⚠️
src/zeroD/ConstPressureMoleReactor.cpp	0.00%	11 Missing ⚠️
src/zeroD/ConstPressureReactor.cpp	60.86%	6 Missing and 3 partials ⚠️
src/zeroD/IdealGasConstPressureMoleReactor.cpp	0.00%	8 Missing ⚠️
src/zeroD/Reactor.cpp	72.41%	4 Missing and 4 partials ⚠️
src/oneD/Sim1D.cpp	71.42%	2 Missing and 2 partials ⚠️
... and 6 more

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1907      +/-   ##
==========================================
+ Coverage   74.11%   74.22%   +0.10%     
==========================================
  Files         445      448       +3     
  Lines       55454    55746     +292     
  Branches     9121     9190      +69     
==========================================
+ Hits        41101    41377     +276     
- Misses      11262    11263       +1     
- Partials     3091     3106      +15     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.