Negative forward and reverse rates of progress

[markus.p...@gmail.com]

Dear User Group members,

I recently found that Cantera calculates (slightly) negative forward and reverse rates of progress (O(1e-30)). I used Cantera 2.2.1 with Python 3.5 interface and simulated constant-volume auto-ignition with GRI-Mech 3.0 mechanism. Can anybody tell me how this can be?

Thanks,
Markus

[Ray Speth]

Markus,

I think you will find that the negative forward/reverse rates of progress correspond to cases where you also have negative mass fractions for some species. These negative mass fractions should generally be on the order of the absolute integrator tolerance, and can be ignored.

Regards,

Ray

[markus.p...@gmail.com]

Dear Ray,

does this mean that in Cantera there is no lower limit (namely zero) imposed on the calculated species fractions? Is it somehow possible to impose such a limit?

Thanks,
Markus

[Ray Speth]

Markus,

Correct, there is no lower limit imposed, and there is no mechanism in Cantera for doing so. In general, attempting to impose such a constraint makes the integration more likely to fail. Cantera follows the CVODES manual's "advice on controlling unphysical negative values" (see page 34), which states among other things: "Remember that a small negative value in y returned by cvodes, with magnitude comparable to abstol or less, is equivalent to zero as far as the computation is concerned.".

Regards,

Ray

[Weiqi Ji]

Hi Ray,

When there is a negative mass fraction, will the mass fraction be clipped to a non-negative value during evaluating the rate of progress? Such as q = k[A][B] and [A] = -1E-12. Whether Cantera actually use [A] = -1E-12 or [A] = 0?

For reaction orders being integer, both might work well. But for fractional reaction order, like [A]^0.5, a negative concentration will not be allowed.

Thanks,

Weiqi

[Steven DeCaluwe]

Hi Weiqi,

Within the Cantera objects, the quantity of a species is indeed clipped to be >= 0. So the integrator does not enforce it, but the properties used by the object/class methods do in fact impose a minima of zero.

For example, in the python interface:

>>> import cantera as ct

>>> gas = ct.Solution('air.yaml’)

>>> gas.X = [0, 1, -1, 0, 1, 0, 0,3]

>>> gas.X

array([0. , 0.2, 0. , 0. , 0.2, 0. , 0. , 0.6])

In the Matlab toolbox, one can specify the `nonorm’ option, which prevents normalization such that the mole or mass fractions do not necessarily sum to one, but I would be very surprised if this circumvented the minimum value of zero.  But since I do not have Matlab and cannot currently check, I thought I would at least mention it.

Best,

Steven 

--
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/2dd64519-780f-423f-878c-003718e805f0n%40googlegroups.com.

[Weiqi Ji]

Hi Steven

I thought the concentration is not clipped to be >= 0 as the clipping may induce numerically issue, according to Ray's discussion above. 

Below is the suggestion that Ray refers to

> 

1. The user's RHS function should never change a negative value in the solution vector y to a non-negative value, as a "solution" to this problem. This can cause instability. If the RHS routine cannot tolerate a zero or negative value (e.g. because there is a square root or log of it), then the offending value should be changed to zero or a tiny positive number in a temporary variable (not in the input y vector) for the purposes of computing f(t,y). [http://sundials.wikidot.com/integrating-over-discontinuities]

Based on the above suggestion, if the negative value is not allowed in calculating rates, I guess the concentrations are clipped?

But if a negative value is allowed, I am not sure if it is clipped?

Thanks,

Weiqi

[Ray Speth]

Hi all,

Using the "unnormalized" methods for setting composition (which are used in both Cantera's 0D and 1D solvers) does indeed permit the concentration to be set to negative values. There are a few cases that have to be handled when it comes to reaction rates, in order to avoid problematic behavior.

The simplest case, of an elementary reaction where one species has a negative value for a concentration, and provides part of the motivation for this behavior. Here a small negative concentration will give a small negative forward rate of progress, and drive the concentration of this species back towards zero.

If a reaction involves two species with negative concentrations, using the negative concentrations would result in a positive forward rate constant, driving those species concentrations to even more negative values. This is of course not desirable, so we clip the minimum value for the rate to zero. A similar logic is applied for reactions involving 3 species.

For reactions with non-integer orders, or more than 3 reactants, the rate is set to zero if any reactant has a negative concentration.

You can see the implementations of all these cases in the StoichManager.h header file.

Regards,

Ray

[Steven DeCaluwe]

Thanks, Ray.  One other situation I wonder about: what if a reaction has a reactant with negative concentration, but the product(s) have zero concentration? 

I know it’s not possible to chase down / anticipate / prevent every single scenario and/or side-branch, but just something that crossed my mind.

Best,

Steven

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/f099fd48-f8a7-4b9d-986d-beeaec2d2433n%40googlegroups.com.

[Ray Speth]

Steven,

Cantera doesn't apply any special handling based on the product concentrations, so it would indeed drive those concentrations to negative values. I think this isn't a problem though, because the behavior is self-limiting, with the reactant concentration is being driven towards zero. At the very least, I've never seen a case where this situation led to integration failures.

Regards,

Ray

[Steven DeCaluwe]

Ray,

Yeah, I didn’t really think it would be a problem that would cause major issues - to reach that state in the first place (negative reactant concentration and zero product concentration) would most likely require some other reaction to drive it into that state.  If not, the two sets of species (reactants and products) in that case would likely oscillate between zero and slightly negative values.  Like I said, was more just curious.  

Thanks for the insights,

Steven

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/254b7bd5-ec33-4697-b8c7-f06c8cd38a9en%40googlegroups.com.

[Ingmar Schoegl]

All,

Curiously, one of my students stumbled across a related issue last week. As mentioned above, small negative values in species concentrations can result from integration and aren't an issue per se. These solutions can obviously saved without problem. Interestingly, when one attempts to set a the gas state to previously attained values (from Python), the negative values are truncated, so it's not possible to recreate the same state. I have not decided yet on whether this is an issue or not, as the negative values are unphysical. On the other hand, the importers for SolutionArray (and FlameBase) are affected by this, which means that information is lost.

-ingmar-

[Ingmar Schoegl]

PS: I stand corrected, as it is possible to set negative values via set_unnormalized_mole_fractions; the issue arising in the context of SolutionArray (and FlameBase) remains.

[Bryan Weber]

Hi Ingmar,

That's an interesting problem. Can you use set_unnormalized_mole_fractions on the SolutionArray? If I recall correctly, that forwards most functions to the underlying Solution instance. Which still leaves FlameBase.

Bryan

[Shepherd, Joseph E.]

Folks,

You can obtain negative concentration values during time integration even with a sophisticated ODE solver, both in python and MATLAB. I have also observed this in the past when the using CHEMKIN and backward difference solvers with adaptive step size control. It is a generic problem with time integration of stiff problems that have internal boundary layers, regions of large rates of decrease of one variable embedded in a slowly varying base flow.

In chemical kinetics, the problem often takes place when computing flow or batch reactions without diffusion. One or species rapidly decreases at the end of the reaction exothermic reaction region and the solver time-control algorithm has been taking increasing large steps up until this point. The algorithm is unable to sufficiently decrease the step size to avoid overpredicting the species concentration, which then becomes negative. This is a particular problem with certain fuel-oxidizer computations. For example, computing the ZND model of methane-oxygen reaction using almost any mechanism (GRI for example) will result in this issue unless the maximum step size and tolerances are tweaked to prevent this. These is an extremely rapid drop in mass fraction from O(E-1) to < O(E-13) in in CH3O, CH3, and CH2 at the end of the energy release reaction zone. In flame or flow reactor computations with diffusion, this same change would happen over a sufficiently large time or distance that the usual integrator can adjust over multiple time steps to avoid creating negative concentrations, but this doesn't happen in simple flow reactor models. The consequence can be that the simulation will throw an error which originates from Cantera attempting to evaluate log(concentration) or similar action that is a numerical crime with negative concentrations. Superficially the problem appears to be with Cantera but it is a solver issue. The same sort of symptom can happen in some cases with the equilibrate function when using certain options that involve entropy and can be fixed by careful attention to the tolerances and choice of equilibrium solver. VCS tends to be more robust than Gibbs minimization in this regard.

The solution is to monitor the species concentrations and adjust the maximum step size and tolerances to enable the solver to resolve these internal boundary layers of species.

Best,

Joe

You can see the implementations of all these cases in the StoichManager.h <https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgithub.com%2FCantera%2Fcantera%2Fblob%2Fmain%2Finclude%2Fcantera%2Fkinetics%2FStoichManager.h&data=04%7C01%7Cdecaluwe%40mines.edu%7C880e1f8cd9ff4aeb1ee708d902a79d27%7C997209e009b346239a4d76afa44a675c%7C0%7C0%7C637543750513947715%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=E%2FJoO764ZpoI6b36UEMfS%2Fx5HIjyeyd1hNsKuH1FskM%3D&reserved=0> header file.

Regards,
Ray


On Sunday, April 18, 2021 at 2:28:51 PM UTC-4 Weiqi Ji wrote:


Hi Steven

I thought the concentration is not clipped to be >= 0 as the clipping may induce numerically issue, according to Ray's discussion above.

Below is the suggestion that Ray refers to
>

1. The user's RHS function should never change a negative value in the solution vector y to a non-negative value, as a "solution" to this problem. This can cause instability. If the RHS routine cannot tolerate a zero or negative value (e.g. because there is a square root or log of it), then the offending value should be changed to zero or a tiny positive number in a temporary variable (not in the input y vector) for the purposes of computing f(t,y). [http://sundials.wikidot.com/integrating-over-discontinuities <https://nam04.safelinks.protection.outlook.com/?url=http%3A%2F%2Fsundials.wikidot.com%2Fintegrating-over-discontinuities&data=04%7C01%7Cdecaluwe%40mines.edu%7C880e1f8cd9ff4aeb1ee708d902a79d27%7C997209e009b346239a4d76afa44a675c%7C0%7C0%7C637543750513957717%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=fp5swGDrw2b%2BmqrAQ6oL87zJYHhWNZR0CsxyxkKSaUQ%3D&reserved=0> ]

To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/2dd64519-780f-423f-878c-003718e805f0n%40googlegroups.com <https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fcantera-users%2F2dd64519-780f-423f-878c-003718e805f0n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7Cdecaluwe%40mines.edu%7C880e1f8cd9ff4aeb1ee708d902a79d27%7C997209e009b346239a4d76afa44a675c%7C0%7C0%7C637543750513967710%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=r0yaqb07UAS61urfii366%2FTx4DdR0H8XHbC%2FoibxLS8%3D&reserved=0> .

--
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/f099fd48-f8a7-4b9d-986d-beeaec2d2433n%40googlegroups.com <https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fcantera-users%2Ff099fd48-f8a7-4b9d-986d-beeaec2d2433n%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7Cdecaluwe%40mines.edu%7C880e1f8cd9ff4aeb1ee708d902a79d27%7C997209e009b346239a4d76afa44a675c%7C0%7C0%7C637543750513977705%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=QDpm5gl%2Fqe%2FsohBdI8Hxrwux3yJnZUwtz3M%2BY1NvXuw%3D&reserved=0> .

--
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/254b7bd5-ec33-4697-b8c7-f06c8cd38a9en%40googlegroups.com <https://nam04.safelinks.protection.outlook.com/?url=https%3A%2F%2Fgroups.google.com%2Fd%2Fmsgid%2Fcantera-users%2F254b7bd5-ec33-4697-b8c7-f06c8cd38a9en%2540googlegroups.com%3Futm_medium%3Demail%26utm_source%3Dfooter&data=04%7C01%7Cdecaluwe%40mines.edu%7C880e1f8cd9ff4aeb1ee708d902a79d27%7C997209e009b346239a4d76afa44a675c%7C0%7C0%7C637543750513987699%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=bTc0c66NmXj92eMV9Z%2Bmsf9c64k2OLb6Wy5NHDwselk%3D&reserved=0> .

--
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 <mailto:cantera-user...@googlegroups.com> .
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/b5123624-8d2c-4f9e-8032-f56430145973n%40googlegroups.com <https://groups.google.com/d/msgid/cantera-users/b5123624-8d2c-4f9e-8032-f56430145973n%40googlegroups.com?utm_medium=email&utm_source=footer> .

[Weiqi Ji]

Hi All,

Add one piece of experience to it. Such negative concentration is linked to stiffness. The rate constants are usually within a certain range, for example, the collision limits. Such physics constraints allow us to estimate the stiffness of the system can suggest common tolerance for the solver. It is OK if we only need to compute a couple of simulations by closely monitoring the simulation.

However, when doing optimization by searching a range of kinetic parameters, it is a headache problem as systems with different rate constants can behavior pretty different. A solver abortion could destroy the whole optimization pipeline. This is some issue people seldomly faced in the non-stiff system. One solution is to impose constraints on A, b, Ea. But for a complex reaction network, it is not so straightforward to determine the range, as in simple problems with a handful of parameters.

Does anyone have suggestions on dealing with such situations?

[Ingmar Schoegl]

All,

I believe the possibility of getting spurious negative values in chemistry simulations (and other stiff systems) is relatively well known, so it is certainly not specific to Cantera. I believe limiting step sizes close to zero would be a great idea (and may be relatively simple to implement - at least if it is not strictly enforced; a somewhat related approach already exists with `set_advance_limits` that is currently used to refine portions with large temperature gradients (see https://cantera.org/examples/python/reactors/reactor1.py.html)

As an interesting aside, it turns out that Cantera's reactor network examples currently truncate those negative values (due to a bug in SolutionArray.append that one of my students just pointed out to me earlier today, see https://github.com/Cantera/cantera/issues/1016). This would mean that the issue is likely a lot more common than it appears.

-ingmar-

You received this message because you are subscribed to a topic in the Google Groups "Cantera Users' Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/cantera-users/EkInq8sKyA4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to cantera-user...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cantera-users/c5ca36c9-73e7-4f91-8682-bde5ffaa46f3n%40googlegroups.com.