Failure using Cantera's reactor model with modified governing equations

[Gandhali Kogekar]

Hello, I am getting the following error (truncated here) when I run the reactor model from Cantera.

*******************************************************************************
CanteraError thrown by Phase::assignDensity:
density must be positive. density = -1.0509961755470373e-35
*******************************************************************************

Components with largest weighted error estimates:
0: 0
1: 0
2: 0
3: 0
4: 0
5: 0
6: 0
7: 0
8: 0
*******************************************************************************

One thing to note here is that I have modified the sdot term from the evalSurfaces() function to model surface diffusion. A similar standalone code written using scipy.ode package runs successfully. I have attached the plot here. I have compared the initial residual values from both models and they are the same. However, the Cantera reactor model takes larger time steps internally and breaks near t ~ 1e-10 sec. I use the same controls (atol, rtol, maxtimesteps etc) in both models. When I use set_advance_limit() functionality, the time steps are too small. The code takes too long  (~60 mins) to reach just 1e-4 seconds and eventually fails. 

I am surprised to see that the python model successfully runs within a couple of seconds, while Cantera model fails despite using the same governing equations. There are no gas-phase reactions involved, and the python model doesn't solve the gas-phase equations. Are there any other controls/methods from sundials that I am missing here? 

Any help is appreciated. 

Thanks,

Gandhali

[Ray Speth]

Hi Gandhali,

Can you please share your code? There isn't enough information here to speculate on why the solver is struggling or what is causing the density to go negative. The easiest way might be to push your changes to a branch on GitHub and share that link, along with a script for setting up and running the specific simulation in question.

Regards,

Ray

[Gandhali Kogekar]

We are comparing two models - one that uses Cantera's reactor model (with modified equations) and another standalone python code solving the same governing equations. The modified Cantera source code is here. The Python script to run this model is reactor_ct.py. The standalone python notebook is 'surface_diffusion_ode'. Both models use the same yaml file. If we print (dtheta/dt) values in both models, initial values are exactly the same. But Cantera model fails afterward. The immediate jump in the solution at t=0 is mainly because of the rate coefficient value (which is correct). The model works fine if we reduce the rate coefficient value by a factor of 1e6 or so.

[Gandhali Kogekar]

Copying a reply from Ray from another channel: 

You're running into a version of https://github.com/Cantera/enhancements/issues/201. Before you calculate the edge reaction rates, you need to call syncState not only on the edge itself, but also all other surfaces that touch that edge. As long as those surfaces are all defined using independent ThermoPhase objects, it is safe to just do this ahead the other calculations in Reactor::evalSurfaces, for example by just writing:

for(auto S : m_surfaces) {
        S->syncState();
    }

I found this much, much easier to debug by using the lower-level ReactorNet.step method than using ReactorNet.advance . Part of the problem here is that the time scale that you're indicating there to the solver is completely different from the timescale introduced by the reactions. You have dtheta/dt values of around 1e10, giving timescales of ~0.1 ns, but are only asking for the solution at times that are many orders of magnitude larger. Are you sure you've got the units/dimensions of these rates worked out correctly?

UPDATE: 

Using the syncstate before all calculations fixed the issue. With above suggested changes, Cantera model too works perfectly fine using ReactorNet.advance. 

Thanks for your help.