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