Speeding up 0D reactor network

[Unni Kurumbail]

Hi,

Hope you're all doing well. I'm working on a PFR simulation in Python using Cantera. I adapted the Lagrangian version of a PFR simulation from the helpful Cantera example file located here: https://cantera.org/examples/python/reactors/pfr.py.html

However, in my situation, I am imposing an experimental temperature distribution on the reactor profile and therefore need to (I think) reinitialize the reactor network during each iteration to update the thermodynamics state of my Solution object:

(relevant code snippets below)

        r = ct.IdealGasConstPressureReactor(gas,energy='off')
        sim = ct.ReactorNet([r])

for n in range(num_steps):

            thermo = r.thermo
            thermo.TPX = new_T, thermo.TPX[1],thermo.TPX[2]
            r.insert(thermo)
            sim.reinitialize()
            sim.advance(t_i)

            ...

I was surprised to observe that just adding the sim.reinitialize() code significantly increases the simulation time. For a small example with 1000 iterations the execution time increased from 2s to 60s. Would anyone know a faster way to apply this temperature distribution without reinitializing each time? Many thanks!

Sincerely,

Unni

[Ray Speth]

Hi Unni,

Using the sim.reinitialize() function frequently prevents the Sundials integrator from operating efficiently, which it can only do when it is running uninterrupted. I would recommend implementing your imposed temperature distribution using the ExtensibleReactor class. This way, you can set the energy_enabled option on the Reactor object to False, eliminating the energy equation from the system, and then implementing a modification to update_state that will set the temperature according to your input. This way, the temperature will be updated as part of the internal calculations initiated by Sundials, and the solver will be to run with reasonable time steps.

Regards,
Ray

​

[Unni Kurumbail]

Ray,

Thanks for the idea, I gave it a go but seem to be running into some difficulties. I implemented a simple ExtensibleReactor object that just changes the temperature based on an externally determined variable 'curr_T':

        class ExtendedFlowReactor(ct.ExtensibleIdealGasConstPressureReactor):
            def before_update_state(self,y):
                #y is the state vector. Element 1 is temperature. Replace with correct temperature
                y[1] = curr_T

Surprisingly, running the extended version of the reactor takes even more time than reinitializing. Whereas one loop of reinitializing/advancing took about 50ms before, it takes about 100ms using the extended reactor. Furthermore, I observed that after a few iterations the integrator throws a CVODE error (-3) which it did not throw before. As such I'm wondering if there's more that needs to be added to properly update the state? Looking at the C++ implementation of this code it does seem that simply swapping in the new temperature should be a one-line fix. But maybe there is another function that is better to modify (e.g. set_state)?

Many thanks,

Unni

[Ray Speth]

Hi Unni,

Sorry, I maybe should have clarified, the update_state method is not a place to modify the state vector y that is provided by the ODE solver. There’s actually no point at which you should modify the function inputs that come from Sundials, as I believe it risks messing up the internal state of the solver.

I looked at this a bit more, and I unfortunately don’t see a consistent way to achieve this right now, without having to fully reimplement the eval method as well. The difficulty is that the update_state is responsible for several things, including updating the ThermoPhase object, saving that state, and setting consistent values for the mass and volume of the reactor based on that, but methods for all of those pieces are not currently exposed in Python. In particular, following the C++ implementation, I think what you would want to do (in theory) is use an after_update_state method that would look something like:

def after_update_state(self, y):
    # composition and pressure will be set by base method
    self.gas.TP = some_function_of_t(self.net.time), None
    self.volume = self.mass / self.gas.density
    self.update_connected(false)  # saves the current thermo state

However, there’s currently no way to call the base update_connected method from the Python delegate class, and I’m not sure of a good workaround.

Regards,
Ray

​

[Unni Kurumbail]

Ray,

 

Thanks for the very helpful information, I will consider other options then.

 

Sincerely,

Unni

--
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/0IeoYgEr-Ng/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/4a7678c6-c25d-4f95-8342-e8f573cab4e6n%40googlegroups.com.