Using the residual in python

[Remi Roncen]

Dear all,

I hope this message finds you well. 

After a 7+ years break in chemistry and combustion, I am starting again in this field, and I have greatly appreciated Cantera and the available tutorials, they really kick-started my learning again !

I have mostly dabbled with the python interface, and I am now trying to dive deeper to make sure I have everything right. The goal I set myself with is to reproduce in python the case of the infinite adiabatic flame, using Cantera only  “for the chemistry” as a black box. This means re-writing the residual calculation, and the Newton solver (re-calculate the Jacobian). I think there would be value for some Cantera users in having such a tutorial. I know I would have liked it 🙂

This approach of re-writing one’s own residual is similar to some of the available tutorials (i.e., 1D_packbed.ipynb), except that a DAE integrator was used to advance the solution in space. 

I was able to access the residual calculated by Cantera, in order to compare it with my own. There are major differences in their amplitudes, which tells me that something is very wrong with my calculations. 

I have attached three python files, in case someone can help:

• my_adiabatic_flame.py, in which I use Cantera’s solver to obtain a solution for the flame, and where I display the different terms composing the residual that I re-calculated, based on that solution. I also compare the total residual with that of Cantera (re-scaled, and multiplied by 5e5 to be of comparable order with my residuals, for the plot). 

• my_tools.py, where I coded the spatial derivatives. I use a central difference adapted for non-uniform grids for interior grid points. At the inlet, an upward scheme is used, and at the outlet, a backward scheme. I (quickly) tested the implementation of the first and second order of these derivatives in test_derivatives.py.

Before continuing with the Jacobian implementation and then the Newton solver, I was wondering if someone experienced in the flux calculations could have a look in my_adiabatic_flame.py, near lines  90-107, and in the writing of the residual terms, lines 110-147, to see if it seems ok to you. Is there any other reason why the residual would be so different from  the one calculated by Cantera (maybe the way I import it) ? I suspect a problem of units somewhere.

Any help is appreciated,

Regards, and a happy new year in advance to you :)

Rémi

[Ray Speth]

Hi Rémi,

There are a lot of choices that have to be made in the finite difference discretization of the nonlinear equations for the laminar flame. Some of them are important for obtaining a solution with the desired stability and accuracy, while many others have minor effects on these properties. However, some of these minor choices will have significant effects on the local solution that lead to large values in the residual from a state vector that satisfied a different solver implementation.

I think what this means is that you're probably only going to be able to check the correctness of your residual function after implementing the hybrid steady/time-stepping damped Newton solver, which has its own share of challenging details.

Regards,

Ray

[Ronny Mönnig]

Very nice, informative correspondence 

Btw, I would be interested in a m-file of this stated problem....does it exist already?

Thank you in advance 

Best,

Ronny 

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/c3c155c3-f59f-4c41-a332-d7d4e95c60ean%40googlegroups.com.

[Remi Roncen]

Dear Ray and Ronny,

thanks a lot for your answers.

Ray, I currently use Cantera's solution, and only recalculate the residual from this converged state. I believe that this approach does not require coding the Newton step and Jacobian at this stage, even though it is still on my todo list. I have also coded the divHeatFlux calculation in Cantera's way, to reproduce the discretization method exactly, and found very similar results to my own (I did make a mistake in the derivative calculation of the initial script in this post). Overall, the mismatch points to an error not on the discretization (which did not affect the residual) or the Newton solver (since I directly used Cantera's solution), but something else. There was a problem in the normalization of the residual that I fixed, but which did not really improve the match between the shape of my residual and that of Cantera. The search continue, and I will try to update the initial script with my progress.

Ronny, I would also be interested if this already exists in Matlab, simply to verify my own integration since both languages are so similar :) Please let me know if you come across anything.

Best,

Rémi

[Remi Roncen]

Dear all,

I have recently been able to make some progress regarding my initial goal, which was to match the residual of Cantera using a python implementation. This was made possible thanks to the help of T. Zirwes, who pointed out to me that some transport properties were not evaluated at the grid nodes, but between green nodes (setGasAtMidpoint function baffled me at first, see StFlow.cpp line 231).

I intend to provide the community with a short code were residuals are matched. In the meantime, I have started to tackle the problem of actually solving the equations, using the Newton procedure. I am starting to get a good grasp of what's happening behind the scene when calling the "solve" function, but there is still a thing or two troubling me. I seek your help on the following:

Once a steady-state Newton solve has failed, a time-stepping method is used (for 10 time steps). This advances the solution in the hope to get closer to the radius of convergence of the steady-state Newton solver.

To efficiently time-step, the Jacobian is updated by subtracting "A*rdt" to it. A is a diagonal matrix, used to "mask" all the algebraic equations composing the residuals (which can  not be advanced in time directly). rdt is the inverse of the time step. 

Looking at OneDim.cpp near line 226, I see that the solve function seems to hard-code the call to the residual with rdt = 0 (line 230 : m_jac->eval(x, xnew, 0.0) ). This means that in the residual function, defined in StFlow.cpp, all the terms that include "rdt" and the "_prev" terms are null. My question is: is this intended ? And if yes, what is the purpose of terms like "- rdt*(Y(x,k,j) - Y_prev(k,j))" ? Because there are a lot of "solve" in the code, I am not sure whether I missed one with a value of rdt that is not 0. Are the terms with "rdt" in the residual redundant with the Jacobian transient update ?

Best,

Rémi

[Thorsten Zirwes]

Hi, thanks for the shout out :)

Regarding rdt: If I remember correctly, the eval function is called both for the Newton step and the pseudo time step. If it is called in “Newton” (steady-state) mode, rdt is set to zero, so the time derivative rdt*(Y(x,k,j) - Y_prev(k,j)) is ignored. For the peuso time step (transient) mode, Y and Y_prev as well as rdt are set and the time derivative is included in the residuals. I guess this is to avoid having to write the same code twice. Here, the residual function is called with the actual time step: https://github.com/Cantera/cantera/blob/main/src/oneD/OneDim.cpp#L274 as a result of this https://github.com/Cantera/cantera/blob/main/src/oneD/OneDim.cpp#L233, where the eval function of the specific 1D domain is called. I haven't looked at the code too closely but this is my takeaway. Hope this helps!

[Remi Roncen]

Hi Thorsten,

thanks for the help, I'll keep the rdt terms.

I have made some more tests recently, showing I had a problem with the calculation of a Newton step when using a sparse matrix in python. 

I still need to figure out the situation with the "diag" matrix, as I don't see its values being set at the left boundary condition in StFlow.cpp (near lines 418). I guess this question relates to a more general one about the setting of boundary condition "objects" on each side of a domain. Any idea where I should look for the Free Flame BC for instance ?

 
Best,

Rémi

[Thorsten Zirwes]

The boundaries are set at two different places in the code:

First, a "general" boundary is set in the residual evaluation of the domain itself (here for the left boundary): https://github.com/Cantera/cantera/blob/main/src/oneD/StFlow.cpp#L470-L503 (for the right boundary, a separate function exists: https://github.com/Cantera/cantera/blob/main/src/oneD/StFlow.cpp#L964-L993).

Second, inside the specific boundary implementation which overwrites the general boundary. For the free flame, the left boundary should be Inlet1D: https://github.com/Cantera/cantera/blob/main/src/oneD/Boundary1D.cpp#L165-L224

[Remi Roncen]

Dear all,

progress has been made since the time of my initial message, and I would now like to share it with the community. Those who are interested can have a look a pyromancy, a code that attempts to only call Cantera for the chemistry part, and that re-introduces its own residual calculation, Jacobian calculation and Newton solver. Doing so was quite educational, even though I am still stuck.

In its current state, the code I developed is not really useful, as it does not allow one to solve for a free flame. At least not yet ! The reason is still a bit unclear to me, but I currently suspect that solving the linear system is the key step I am still not doing properly. In Cantera, this step relies on LAPACK solvers that can handle the band matrix structure of the Jacobian, with an LU step for the factorization. While I was able to call similar LAPACK routines from python’s scipy package, I found that they would not “tolerate'' the Jacobian I calculated, as it is too ill-conditioned. It is thus possible that errors remain in my Residual or Jacobian calculation. Or that I am missing a key conditioning step.

For those of you who want to dig in the code, let me know :) Hopefully I’m not the only one with an interest in this python partial “rewrite”. 

In terms of structure, pyromancy somewhat follows Cantera’s, with some changes. I am still not fully set on them, and things might change in the future. Currently, only a single domain is supported, and only the equations of a free flame have been coded.

• I have created an “examples” folder with a script in which I compare the residual calculation done by Cantera and pyromancy. There is also a “start_here.py” script to get you started, should you wish to 🙂

• I have also created a “tests” folder where I stored the current tests I am developing to see what works, and especially what doesn’t (it’s not a unit-test folder).

I am currently working on the script “tests/test_FreeFlame_convergence.py”, where I create a flame (which can be initialized by Cantera’s solution), a Jacobian object (I compare my Jacobian implementation with the default implementation that does not make use of the banded structure), and I calculate a Newton step in a few different ways. I then attempt to solve the problem, to no avail at the moment. I created a SolidConduction class to code a (simpler) residual for the heat equation (with more well-conditioned Jacobian). The resolution seems to work, even though I obtain small differences between the Newton solver and the time-stepping solution.

I am short of ideas in terms of what to try next, especially since I do not know how to return the Jacobian matrix directly from Cantera in order to compare it with my own. So I’ll consider any input 🙂

Best,

Remi

[Ronny Mönnig]

Regarding solvers I can recommend EIGEN library of C++

Not sure, if’s applicable here, but maybe give it a try within your code 

[Remi Roncen]

Thanks for the solver idea Ronny, will check it out !

[Remi Roncen]

Dear all,

I finally was able to solve my problem. For future reference, I will explain what was happening.

To anchor the flame in the mesh, the FreeFlame needs to have its temperature "anchored" at one location. This means that at this chosen grid point, the residual is set as the difference between the current temperature, and the one fixed in advance.

Now, the trick here is that this algebraic equation (which is thus not included in the time-stepping procedure, i.e., "diag=0") is not applied on the energy equation with "offset_T", but on the continuity equation with "offset_U", in line 973 of StFlow.cpp.

I had been mixing things up a bit in my code, which ended in the condition being coded twice, once in the energy equation and once in the continuity equation. This resulted in a singular Jacobian. 

Part of my confusion came from line 474 of StFlow.cpp where, if the energy equation is disabled, the residual for offset_T is also the difference between the temperature and a fixed temperature (on all grid points this time). Also, if the energy equation is disabled, the continuity equation at the "anchoring" location is changed from a temperature difference to a momentum difference, in order to avoid a singular Jacobian I presume.

I think Cantera's way of doing things is the correct one, i.e., to anchor the flame with a temperature difference on the continuity equation residual. I suspect that the problem can be quite more sensitive to the energy equation residual than it is to the continuity residual. It is best to avoid "wasting" the information contained in the residual of the energy equation.

I will now proceed to clean-up the code, add the grid refinement (already coded somewhere else), improve the comments, the tests, and write a documentation that details the different steps & struggles I went through in re-coding this small portion of Cantera. Can't wait :D

See y'all,

Rémi

ps : the link to access pyromancy --> https://bitbucket.org/rroncen/pyromancy/