Problem Creating Custom Extensible Reactor

[Ian Walters]

I was trying to create a ZND solver for laminar detonation calculations using the extensible reactor framework in Python, as I've found that the built-in Cantera reactors that use CVODE for integrating the ODEs are faster than custom implementations that use the integrators included with scipy - eg. LSODA. However, I was having trouble getting the reactor to work correctly, and I suspect it's because I'm missing some simple aspect of the extensible reactor class.

The specific issue is that the state does not appear to be updating. I've used replace_get_state, replace_update_state, and replace_eval to change all of the variables and equations being solved. Since the ZND equations require 3 + n_species state variables, I've used the standard ExtensibleReactor as the base class. I assume there wouldn't be any benefit to using the ExtensibleIdealGasReactor instead, since again I'm replacing all of the equations?

I've attached a MWE of the reactor and hard-coded the initial condition to remove any dependence to the Shock & Detonation Toolbox.

Please let me know if you need any additional information to understand the issue - thanks in advance for any help!

Best regards,

Ian

[Ray Speth]

Hi Ian,

It looks like you’ve run into a snag with the behavior of the Reactor.thermo property, which automatically calls Reactor.syncState. Since you haven’t overridden, syncState, it retains its default implementation and is restoring an initial state that never gets updated. I think the most efficient fix for this would be to store a different reference to the Solution object in your reactor, e.g. as self.gas, and access it through that property rather than as self.thermo.

One other change you’ll probably want to make is to make sure that you unnormalized mass fractions in replace_update_state so the finite difference Jacobian will be computed correctly, e.g.:

def replace_update_state(self, y): self.gas.set_unnormalized_mass_fractions(y[3:]) self.gas.TD = y[1], y[0] self.x = y[2]

Regards,
Ray

​

[Ian Walters]

Ray,

Thank you so much for the response - sure enough, that resolved my problem. As I originally expected, I was missing something simple.

I also appreciate the recommendation to use the unnormalized mass fraction setting routine - I've encountered issues with not using that before, but had forgotten here.

During my preliminary profiling the ExtensibleReactor version appears about 5x faster than using scipy's built in integrators, which is promising.

Thank you once again!

Ian