Jacobian, Timescales

[Siddhant Deshmukh]

Hello,

I was wondering if the Jacobian of the supplied network to KROME is stored anywhere / is accessible at runtime. I would like to compute its eigenvalues to see the chemical timescales of the system at various densities / temperatures. As a follow-up, does KROME have some way of determining a characteristic timescale for the network to proceed to steady-state, or the slowest rate in the network (at a given temperature)?

Cheers,

Siddhant

[tommaso grassi]

Hi Siddhant,

the Jacobian is automatically created by DLSODES by using the MF=222 options, more details are in the comments of krome/solver/opkdmain.f before the DLSODES call.

For constant temperature models the Jacobian is algebraically created by KROME (see krome_ode.f90 after using the preprocessor).

Determining the characteristic time-scales is a missing features in KROME, but clearly it is of great interest, hence if you plan to do something along this line we are interested in the future developments.

cheers,

-tg

[Siddhant Deshmukh]

Hi Tommasso,

Great, thanks for the info. Is the algebraically created Jacobian for constant temperature models available as a matrix in PyKROME? I am essentially doing chemical post-processing on model outputs where I have the gas density and temperature for each grid cell. For this, I would like to evaluate the eigenvalues of the Jacobian as the (inverse) timescales of the reaction network. Hence it would be quite useful if, for a given set of number densities and a temperature, I could get the Jacobian matrix in PyKROME.

I am following the statement in https://tcg.mae.cornell.edu/pubs/Pope_CTM_97.pdf on page 2 (paragraph above equation 5) which relates the Jacobian of the chemical reaction network to the timescales. However, I do not yet know what is sufficient to yield a "characteristic timescale" and whether this is simply the longest one. I will keep you posted if this leads anywhere.

Cheers,

Siddhant

[tommaso grassi]

This is not available in PyKROME now, I am working on it. I'll let you know when ready.

Please let me know if the timescale thing work somehow, I am really interested in this :)

[Siddhant Deshmukh]

Ok great to hear. For now, I just copied the analytical Jacobian over into a Python file. The timescale seems to work but I would hesitate to blindly recommend this method. Perhaps it is also because of my re action network + temperatures/densities since I am looking at stellar atmospheres, but some of the eigenvalues are really large values, and some others are positive. That being said, if I 'prune' these to remove the ones that seem unlikely to be accurate, I get timescales that roughly correspond with intuition. I tried to use this method for very simple reaction networks and it works quite well, leading me to believe that the network I'm using is the problem, not the Jacobian analysis itself.

I plan to look at this in more detail (and work on adding KROME to our code) next month, so I will keep you posted on other methods I try.

[tommaso grassi]

Note that the Jacobian as implemented in KROME is not very correct if you have variable temperature.

This is because in principle one should also take into account that in d(dx_i/dt)/dT also the rate coefficients are variables, and the finite differences should be applied also there and not only to the d(dT/dt)/dx_j as it is now.

Or (but this is more complicated) one should include the analytical derivative of the rate coefficients w.r.t. the temperature.

A simple example is a network like
 A -> B

where let's say you have just the rate coefficient k(T) = T^2

The time derivative of B is then dB/dt = k(T)*A = T^2 * A

Hence the Jacobian terms for dB/dt are

 d(dB/dt)/dA = T^2

 d(dB/dt)/dB = 0

However, if you also have variable temperature there is another set of terms (that are not included in KROME if you use the "algebraic" Jacobian)

 d(dB/dt)/dT = 2*T*A

The "automatic" Jacobian that use finite differences for all the terms already takes into account this by construction.

I hope this is not too confusing :)

[Siddhant Deshmukh]

Yes, understood, thanks. What I did was take the algebraic Jacobian and run KROME on a grid where each cell has fixed temperature. So it is not physically accurate of course, and I should be including the variable temperature. But I was just doing some testing to see if the first-order approximation of computing a timescale from a given gas density + temperature was feasible.  But I wonder if adding the temperature dependence here will fix the issues.

[tommaso grassi]

In principle if the temperature doesn't change, the temperature is not a problem, because it's not a variable and therefore it doesn't go into the Jacobian.

But you might have some chemical rate that depends on other parameters that are also variables (or that are linked to variables), like for example some density-dependent rate, where the ratio between the critical density and the actual density comes into play (i.e. the rate is not constant during a call to KROME, or the solver).

Additionally, it might happen that even if you have dT/dt = 0, but T is anyway a variable in the solver, this might cause some problem when the solver "explores" the solution, but I am not sure if this is an actual issue or not.

btw, a couple of days ago this paper was on arxiv https://arxiv.org/abs/2105.12070 and it might be in part helpful.