Handling negative mass fractions

[H V]

Hi all 

I'm having an issue with my reactor network where I see negative concentrations for  a significant number of the species. The concentrations are a few orders of magnitude higher than the solver tolerances, and this is causing the integration to fail. 

I looked at these two helpful threads

https://groups.google.com/g/cantera-users/c/EkInq8sKyA4/m/TEMuwCU8AQAJ

https://github.com/Cantera/cantera/issues/1117

As I understood, it's a common problem but I'm wondering what the best way to get around it is. I thought of manually setting the concentrations to 0 after each step but I see this not advisable. Any thoughts / comments on how to tackle this would be great!

I'm interested in seeing the species evolution at a fixed rate of temperature rise. I see that the integration only fails when the rate is quite slow ~1K/s whereas fast rates ~1K/ns work. 

[H V]

Alos should mention that I checked that the rates do not violate the collisional limit

[Steven DeCaluwe]

Hi H V,

Two solutions that I have used, in the past:

1. Cheap / easy, but maybe not fundamentally sound: 

a. When you feed the concentrations to Cantera, set a min of zero.  You are not manually setting the solution vector to zero, just monkeying with how Cantera sets the state.  Since these are presumably very small concentrations, I would hope the error introduced is negligible.

b. Similarly, send Cantera the absolute value of the mass/mole fractions.

2. Is your state variable currently molar / mass concentration, or mole / mass fraction?  If the latter, reconfiguring to solve for concentration might help matters.  Since the fractions are not actually conserved quantities, there are sometimes odd artifacts that lead to problems like this.

This is all a little speculative, but hopefully there’s an idea in there, somewhere.

Cheers,

Steven DeCaluwe

--
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 visit https://groups.google.com/d/msgid/cantera-users/aef8e014-cd6f-4fc9-9a9f-5865f3b750f5n%40googlegroups.com.

[H V]

Hi Steven

Thanks for the reply!

For 1. How do I set the concentrations to a minimum of zero? Is there somewhere I can specify this? 

For 2.I current set the state using gas.TPX. I assume this means cantera sets the state variable to mole fractions. Would I be correct in using set_unnormalized_mass_fractions() instead ?

[Steven DeCaluwe]

Hi H V,

Apologies - I didn’t look at your code until just now.  I had assumed, since you were asking about overwriting mole fractions, that you were integrating your own set of equations, rather than using one of Cantera’s pre-defined reactors.

Since you’re running a constant pressure reactor with a fixed dTdt, at this point you could just write your own integration function and use something like scipy’s solve_ivp to integrate it.

Something like this would work, as your residual function.  There is likely a way to use the extensible reactor to do the same thing, if you don’t want to bother with the overhead of setting up an integrator.

def residual(t, y, args):

#initialize resudial

ydot = np.zeros_like(y)

# Assume Temperature index is 0 (and user-defined dTdt is in the 0 position in ‘args’):

T = y[0]

ydot[0] = args[0]

# Assume user-defined pressure is args[1]

p = args[1]

# Species mole fractions are y[1:]

x_k = np.maximum(y[1:, np.zeros_like(y[1:]))

#alternately, use:

#x_k = np.abs(y[1:])

gas.TPX = T, P, x_k

omega = gas.net_producton_rates

gas.basis = ‘molar’. #there is likely a way to set this once, when you create the phase

c = gas.density

dot[:1] = omega / c

On Jul 13, 2025, at 7:00 AM, H V <bjs...@gmail.com> wrote:

Hi Steve

Thanks for the reply!

For 1. How do I set the concentrations to a minimum of zero? Is there somewhere I can specify this? 

For 2.I current set the state using gas.TPX. I assume this means cantera sets the state variable to mole fractions. Would I be correct in using set_unnormalized_mass_fractions() instead ?

On Friday, 11 July 2025 at 17:58:48 UTC+2 S. DeCaluwe wrote:

To view this discussion visit https://groups.google.com/d/msgid/cantera-users/3d5e61e9-e875-45b5-95a8-2d17dc0119f8n%40googlegroups.com.

[Ray Speth]

Hi,

Your mechanism most certainly does have rate constants that exceed the collisional rate limit. At 300 K, the highest reverse rate constant is 7.14e162, and there are dozens of others in excess of 1.0e30.

One other problem, which might seem a bit paradoxical, is that the integration limits you are setting are too stringent. It is not really possible to solve using double precision arithmetic to a relative tolerance of 1e-16. If you stick to the default tolerances (rtol of 1e-9), the solver fares quite a bit better, though the divergent behavior does eventually become a problem.

To help set things back on course before they go too far astray, you can reset the negative species, as Steven suggested. Specifically, you can rewrite the integration loop as:

n_reset = 0 while r.thermo.T < T_max: net.step() states.append(TPX=r.thermo.TPX, t=net.time, normalize=False) state = net.get_state() Y = state[2:] print(f't = {net.time}, Ypos = {sum(Y >= 0)}, Yneg = {sum(Y < 0)}, Ymin = {min(Y)}') if min(Y) < -1e-10: print('Resetting negative mass fractions') gas.TPY = None, None, Y r.syncState() net.reinitialize() n_reset += 1

Along with leaving the tolerances alone, I find that this gives a reasonable looking solution to the problem and only needs to reset the integrator once.

Regards,
Ray

​

[H V]

Dear Ray, Dear DeSteven, 

Thanks for the useful tips! 

For the collisional rate limit I did indeed not check the collisional limit for the reverse reactions. 

I was using the rate limit proposed here -> 10.1016/j.combustflame.2017.08.005 to check the collisonal limit. 

For a single species is the appropriate value the self collisional limit? 

With respect to resetting the negative species, is there an advantage to resetting the negative values in the extensible reactor rather than in the loop? 

class ImposedTemperatureReactor(ct.ExtensibleIdealGasConstPressureReactor):
def __init__(self, *args, dTdt, **kwargs):
super().__init__(*args, **kwargs)
self.dTdt = dTdt
self.kEnergy = self.component_index('temperature')

def after_eval(self, t, LHS, RHS):
LHS[self.kEnergy] = 1 #dt
RHS[self.kEnergy] = self.dTdt #dTdt/dt

# Fix 1: Handle negative concentrations and derivatives
n_species = len(self.thermo.species_names)
species_start_idx = 2 # Skip mass and energy equations

# Check if any species has negative concentration
concentrations = self.thermo.concentrations.copy()
has_negative = np.any(concentrations < 0)

if has_negative:
# Set negative concentrations to zero
concentrations_fixed = np.maximum(concentrations, 0.0)
# Update the reactor state using concentrations
self.thermo.concentrations = concentrations_fixed

# Also clamp derivatives for very small concentrations
for i in range(n_species):
species_idx = species_start_idx + i
# If concentration is very small and derivative is negative, clamp derivative to zero
if self.thermo.Y[i] < 0.0 and RHS[species_idx] < 0:
RHS[species_idx] = 0.0
self.thermo.Y[i] = 0.0

# re-normalize the mass fractions
self.thermo.Y = self.thermo.Y/ sum(self.thermo.Y)

This makes the integration quite slow but seems more stable. 

Thanks!!