Steady-State Solver for ZeroD Reactors #1021

paulblum · 2021-04-26T18:08:36Z

Cantera can solve systems of equations representing physically zero-dimensional (i.e., time dependent-only) or one-dimensional (i.e., time + space) systems. The solution of the 1-D problems is done under the assumption of steady state, resulting in a differential-algebraic system of equations to solve. There is a solver implemented in Cantera as part of the 1-D code that solves these problems. 0-D systems can either be transient or steady state, represented by an ODE or algebraic equation respectively; at the moment, the only way to achieve solutions of the steady state problem is to integrate the transient problem until properties stop changing. We would like to modify or implement the existing steady-state 1-D solver for 0-D systems.

Changes proposed in this pull request

Introduce Newton, a simplified version of Cantera's 1-D solver for the solution of 0-D steady-state systems.
Add solveSteady() functionality to ReactorNet in C++ and solve_steady() in Python.
Modify DenseMatrix for optimized Jacobian computation.

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

Fixes #

Checklist

There is a clear use-case for this code change
The commit message has a short title & references relevant issues
Build passes (scons build & scons test) and unit tests address code coverage
The pull request is ready for review

bryanwweber · 2021-04-26T18:18:47Z

@gkogekar FYI, this is the work that @paulblum has been doing on the steady state solver.

codecov · 2021-05-03T05:40:01Z

Codecov Report

Attention: Patch coverage is 0% with 192 lines in your changes missing coverage. Please review.

Project coverage is 65.89%. Comparing base (2d772ab) to head (b1ad331).
Report is 2592 commits behind head on main.

Files with missing lines	Patch %	Lines
src/numerics/Newton.cpp	0.00%	134 Missing ⚠️
src/numerics/DenseMatrix.cpp	0.00%	31 Missing ⚠️
src/zeroD/ReactorNet.cpp	0.00%	21 Missing ⚠️
include/cantera/numerics/Newton.h	0.00%	6 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1021      +/-   ##
==========================================
+ Coverage   65.39%   65.89%   +0.50%     
==========================================
  Files         318      320       +2     
  Lines       46109    47575    +1466     
  Branches    19601    20496     +895     
==========================================
+ Hits        30153    31350    +1197     
- Misses      13452    13686     +234     
- Partials     2504     2539      +35

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

ischoegl · 2021-10-05T17:07:07Z

@paulblum ... as this just came up on the UG, what is the current status of this PR? - and to confirm: this should close Cantera/enhancements/issues/31, correct?

ischoegl · 2021-12-14T17:06:42Z

@paulblum and @bryanwweber … this PR came up again in discussions recently. Are there any plans to push this over the finishing line?

bryanwweber · 2021-12-14T18:09:58Z

@ischoegl Thanks for following up. I don't personally have any plans to work on this, and @paulblum has moved on to a PhD program in a different area, unfortunately, so he probably won't have time either. If you want to pick up the ball from here, please feel free to reuse this code/commits as necessary, or ditch it altogether.

ischoegl · 2021-12-14T20:02:25Z

@bryanwweber ... thanks for the clarification and too bad indeed. It's not something I was intending to work on (and don't have any immediate plans at least at the moment).

speth

I rebased this onto the current main branch so it doesn't get left too far behind. The comments below are partly just to remind me (or whoever finds time to pick this up again) of a couple things that I saw that need some work in the implementation.

Currently, there are no tests, and the simplest thing I could think of to test didn't actually work. This was to use the mix1.py example and replace the call to sim.advance_to_steady_state() with a call to sim.solve_steady(). The solver prints out the message "Failure to converge in 17 steps.", and the solution is clearly not the right steady-state.

speth · 2022-02-04T23:17:23Z

src/numerics/Newton.cpp

+    for(int i = 0; i < MAX; i++) {
+        newtonsolves++;
+        m_config = &m_directsolve_config;
+        if(solve(m_x.data()) > 0) {
+            writelog("\nConverged in {} newton solves, {} timesteps.", newtonsolves, timesteps);
+            return 1;
+        }
+        m_config = &m_timestep_config;
+        copy(m_xsave.begin(), m_xsave.end(), m_x.begin());
+        for(int j = 0; j < MAX; j++) {
+            if (solve(m_x.data()) < 0) {
+                writelog("\nTimestep failure after {} newton solves, {} timesteps.", newtonsolves, timesteps);
+                return -1;
+            }
+            timesteps++;
+        }
+        copy(m_x.begin(), m_x.end(), m_xsave.begin());
+    }
+    writelog("Failure to converge in {} steps.", MAX);
+    return 0;


I think this solver loop only does a much more limited version of the hybrid time-stepping / steady solver structure that's implemented by the MultiNewton class, and is only alternating between solving the steady state and taking time steps of a fixed size. Given the sophistication of the MultiNewton algorithm, it might be worth seeing if we can pose this problem in a way that that solver can use directly.

speth · 2022-02-04T23:19:54Z

src/zeroD/ReactorNet.cpp

+    for (int i = 0; i < m_nv; i++)
+    {
+        m_newt->setBounds(i, -1.0e-3, 1.01);
+    }
+    m_newt->setBounds(0, 1.0e-3, 100);
+    m_newt->setBounds(1, 0, 5);
+    m_newt->setBounds(2, -10000000, 10000000);
+
+    m_newt->setConstants({1});
+    // m_newt->setConstants({0,1,2,11,12});


Settings bounds definitely needs to be moved to the individual reactor classes. I think these are meant to apply only to the base Reactor class, but even then the values for the first (mass) and second (volume) components don't quite make sense.

ischoegl · 2022-06-06T18:24:28Z

I have a relatively broad question here that follows up on #654 (comment), reproduced here:

#1021 and advance_to_steady_state are very different approaches to achieve the same goal. The proposal in #1021 uses the existing steady-state solver for 1-D problems, which is a hybrid transient solver. The idea is, in steady state, the problem reduces to a set of equations that can be solved by Newton iteration... however, Newton iteration can be unstable, in the sense that it is very sensitive to the initial guess, and it can diverge if the initial guess isn't "good" enough. We know that time stepping will always bring the solution closer to the steady state, so #1021 tries a direct Newton solution for the equations, and if it fails, takes a few timesteps, then tries the Newton solution again. This repeats until the Newton solver converges. As it happens, the accuracy of the timesteps doesn't have a strong effect on the Newton solver's final solution, so this solver doesn't use CVODES because it doesn't need all the accuracy that CVODES is capable of providing.

On the other hand, advance_to_steady_state relies on CVODES and repeatedly calls the advance method until some stopping criteria are satisfied. This works because, as I said, timestepping always goes to the steady state (if there is one), but it is potentially much slower than the Newton solution, because time steps take a long time to solve, and you need to take a lot of them to maintain accuracy.

Fundamentally, there's nothing wrong with the 1D approach, but at least from my experience, it is prone to getting stuck for problems that involve marginal (low heating value) mixtures.

But if I am informed correctly, the basic idea for the hybrid solver predates Cantera, and is essentially what was implemented for the original CHEMKIN. My main concern here is that we may be implementing 1980's ideas (which still may be worth keeping, but a lot has happened since) ... why not use CVODES for the transient solver, especially if we will be able to precondition soon? Are there better approaches than a standard newton solver? (e.g. CG etc.). Can we leverage analytical jacobians (which we are close to finally having after #1089 and #1010)? Can we leverage packaged root finders? Can we use some findings here for updated approaches to the 1D solver down the road?

bryanwweber · 2022-06-06T18:53:08Z

why not use CVODES for the transient solver, especially if we will be able to precondition soon? ... Can we leverage analytical jacobians (which we are close to finally having after #1089 and #1010)?

My understanding is that we simply don't need all the overhead of creating CVODES stuff, either the time or memory overhead, since the controls that CVODES provides are simply unnecessary. Likewise for analytical Jacobians.

Can we leverage packaged root finders?

Paul originally attempted KINSOL here, but packaged solvers don't tend to provide control on constraints (e.g., mass fractions) that are required for this application. They also don't provide the hybrid approach, and a Newton algorithm is fairly simple to implement, and it's already done for us and been validated for ~20 years... Again, just my understanding.

Are there better approaches than a standard newton solver? (e.g. CG etc.).

This looks like great enhancement request if you can provide some more details 😄 As this implementation is relatively close to done, I don't think we should close this, unless replacement features are merged first though.

Can we use some findings here for updated approaches to the 1D solver down the road?

Sure, that seems natural...

bryanwweber · 2022-06-06T18:53:49Z

it is prone to getting stuck for problems that involve marginal (low heating value) mixtures.

This sounds like a bug, or something that could be improved more generally. Opening an issue with a sample script that fails would be most welcome!

ischoegl · 2022-06-06T21:59:12Z

@bryanwweber … would Cantera/enhancements#154 help with the KINSOL issues?

bryanwweber · 2022-06-06T22:37:44Z

@ischoegl I don't think so. You still have constraints on the state variables to be positive, so that wouldn't be resolved. You'd also have to run the hybrid loop separately still, KINSOL can't handle that no matter what the state variables are 😞

ischoegl · 2022-06-06T23:16:35Z

@bryanwweber ... I believe KINSOL does support constraints, see one of the examples
https://github.com/LLNL/sundials/blob/e2f29c34f324829302037a1492db480be8828084/examples/kinsol/serial/kinFerTron_dns.c#L47-L48

and you can put constraints on number of iterations, see https://sundials.readthedocs.io/en/latest/kinsol/Usage/index.html#c.KINSetNumMaxIters

One big caveat is that I haven't spent any time with implementations here/

bryanwweber · 2022-06-07T00:43:24Z

@ischoegl I haven't tried it personally either 🤷‍♂️ all I can say is that the people who have tried it (@paulblum and @speth IIRC) have noted that there are enough restrictions on using it that it isn't worth it for our case.

speth · 2022-06-12T16:01:45Z

I think this is the wrong place for this general discussion on the relative merits of different approaches to solving the steady-state problem -- we already have an enhancement issue here: Cantera/enhancements#31.

ischoegl · 2023-04-24T01:03:32Z

Looks like this PR got an ‘honorable’ mention here …

paulblum force-pushed the zeroD_steady branch from 94799ea to 60ab51d Compare May 3, 2021 05:00

This was referenced Jun 9, 2021

Implement a steady-state solver for 0-D systems Cantera/enhancements#31

Open

[WIP] Jet Stirred Reactors pr-omethe-us/PyTeCK#17

Open

paulblum added 9 commits February 4, 2022 16:33

Add initial steady state solver to zeroD

378342d

Add FuncEval-based Newton solver for steadystate

083df60

Remove dependency on external "Jacobian" class

991ecf3

Add constant var capability, cleanup unneeded code

f8515ae

Fix jacobian forced reeval issue to allow jac reuse

0de9ba2

Initial timestep functionality and jac optimizations

3326d31

Solver rewrite

11e2969

Newton solver rewrite part 2

7075f5a

Add configurations, for alt settings btwn solver modes

b1ad331

speth force-pushed the zeroD_steady branch from 21d9b20 to b1ad331 Compare February 4, 2022 23:08

speth reviewed Feb 4, 2022

View reviewed changes

speth added help wanted Reactors labels Feb 4, 2022

ischoegl closed this Feb 27, 2022

ischoegl reopened this Feb 27, 2022

12Chao mentioned this pull request Jun 3, 2022

[Python] Fix ATOL in advance_to_steady_state #1305

Merged

5 tasks

rwest mentioned this pull request Jun 3, 2022

Fix implementation of advance_to_steady_state #654

Open

speth mentioned this pull request May 16, 2024

Improve ReactorNet handling of bulk phase interactions Cantera/enhancements#202

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Steady-State Solver for ZeroD Reactors #1021

Steady-State Solver for ZeroD Reactors #1021

paulblum commented Apr 26, 2021

bryanwweber commented Apr 26, 2021

codecov bot commented May 3, 2021 •

edited

Loading

ischoegl commented Oct 5, 2021 •

edited

Loading

ischoegl commented Dec 14, 2021

bryanwweber commented Dec 14, 2021

ischoegl commented Dec 14, 2021 •

edited

Loading

speth left a comment

speth Feb 4, 2022

speth Feb 4, 2022

ischoegl commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

ischoegl commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

ischoegl commented Jun 6, 2022 •

edited

Loading

bryanwweber commented Jun 7, 2022 •

edited

Loading

speth commented Jun 12, 2022

ischoegl commented Apr 24, 2023

Steady-State Solver for ZeroD Reactors #1021

Are you sure you want to change the base?

Steady-State Solver for ZeroD Reactors #1021

Conversation

paulblum commented Apr 26, 2021

bryanwweber commented Apr 26, 2021

codecov bot commented May 3, 2021 • edited Loading

Codecov Report

ischoegl commented Oct 5, 2021 • edited Loading

ischoegl commented Dec 14, 2021

bryanwweber commented Dec 14, 2021

ischoegl commented Dec 14, 2021 • edited Loading

speth left a comment

Choose a reason for hiding this comment

speth Feb 4, 2022

Choose a reason for hiding this comment

speth Feb 4, 2022

Choose a reason for hiding this comment

ischoegl commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

ischoegl commented Jun 6, 2022

bryanwweber commented Jun 6, 2022

ischoegl commented Jun 6, 2022 • edited Loading

bryanwweber commented Jun 7, 2022 • edited Loading

speth commented Jun 12, 2022

ischoegl commented Apr 24, 2023

codecov bot commented May 3, 2021 •

edited

Loading

ischoegl commented Oct 5, 2021 •

edited

Loading

ischoegl commented Dec 14, 2021 •

edited

Loading

ischoegl commented Jun 6, 2022 •

edited

Loading

bryanwweber commented Jun 7, 2022 •

edited

Loading