Updating Cantera species thermo data without writing to xml file (working from Python)

[Nwike]

I am running a parameter fitting simulation (in python) where I need to update the value of the standard enthalpy and entropy for Cantera species (in an ideal solution phase). The xml species description looks something like this:

<species name="DummySpecies">

      <atomArray> C:1 H:1</atomArray>

      <thermo>

        <const_cp Tmax="300.0" Tmin="298.0">

          <t0 units="K">298.15</t0>

          <h0 units="J/mol">0.0</h0>

          <s0 units="J/mol/K"> 0.0 </s0>

          <cp0 units="J/mol/K"> 0.0</cp0>

        </const_cp>

      </thermo>

      <standardState model="constant_incompressible">

        <molarVolume> 0.0 </molarVolume>

      </standardState>

    </species>

 

And I am interested in updating the values of ‘h0’ and ‘s0’. To run this simulation, I load the phase data from xml and use it to create a mixture object like so:

phase      = cantera.import_phases(XmlFile,['Phase Names'])

mixture   = cantera.Mixture(phase)

I then pass the mixture object to the fitting code that performs an iterative calculation. The issue is that at each iteration, I need to update the enthalpy and/or entropy information for the species. Currently, I write the updated values back onto the xml file, and read it again using the ‘import_phases()’ function, but this is inefficient. I tried unsuccessfully to directly update the species data, but this is not possible once the phase object has been created:

phase[phase ID].species()[Species ID].thermo.coeffs[Coefficient index for h0] = New_h0_value *********** does not work

Is there a way to update these values without first writing to, and rereading from the input xml file?

[Bryan W. Weber]

Hi,

What version of Cantera are you using, from what interface, on what OS? Can you please include the complete script that you run and the complete XML or CTI file?

Best,

Bryan

[Nwike]

Hello Bryan

Here's the information you requested:

Cantera Version: 2.3.0

Interface: Python 3.6.1

OS: Windows 10

Python scripts:

fit.py: main script for running the parameter fit and plotting results

meth.py: support function definitions (the update step takes place in function 'obj' in the 'meth.py' file

twophase.xml: main xml file

elementz.xml: updated xml elements file to include inert diluent

data.csv: experimental data table

Thanks

Nwike

[Nick Curtis]

Hi Nwike,

Here's one way to utilize an updated species thermo coeff via python (w/o rewriting the xml):

import cantera as ct
gas = ct.Solution('gri30.cti')
h2 = gas.species(0)
h2_new = ct.Species(h2.name, h2.composition, h2.charge, h2.size)
new_coeffs = h2.thermo.coeffs.copy()
new_coeffs[0] = 1
h2_new.thermo = ct.NasaPoly2(h2.thermo.min_temp, h2.thermo.max_temp, h2.thermo.reference_pressure, new_coeffs)
gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=[h2_new] + gas.species()[1:], reactions=gas.reactions())
assert gas.species(0).thermo.coeffs[0] == 1

Constructing the gas directly from the pre-constructed species / reaction objects seems to be faster on my system (at least, you don't have to re-parse the XML for large files).

Nick

[Nwike]

Hello Nick

 

Thanks a lot for your comments. I implemented the steps you suggested for my system and everything looked good until I got to the last line (gas = ct.Solution...), which threw new errors.

 

Case 1: Error: CanteraError thrown by getElementWeight: element not found: dummy

‘dummy’ refers to a fictitious element I created to model the inert species, dodecane, in elementz.xml. Below’s the script that throws this error:

 

import cantera as ct

p = ct.Solution('yesdummy.xml','PC88A_liquid')

h2 = p.species(0)

h2_new = ct.Species(h2.name, h2.composition, h2.charge, h2.size)

new_coeff = h2.thermo.coeffs.copy()

new_coeff[1] = -2350000

h2_new.thermo = ct.ConstantCp(h2.thermo.min_temp, h2.thermo.max_temp, h2.thermo.reference_pressure, new_coeff)

p = ct.Solution(thermo='IdealSolidSolution', species=[h2_new]+p.species()[1:], datasrc='elementz.xml')

 

 

Case 2: Error: CanteraError thrown by IdealSolidSolnPhase::setDensity: Density is not an independent variable

For this case, I  removed the ‘dodecane’ species from the phase to temporarily get rid of the first problem, but then run into this new exception. I do not know if anyone has seen this before (I did a quick search of topics on this forum). Here’s the corresponding script:

 

import cantera as ct

p = ct.Solution('nodummy.xml','PC88A_liquid')

h2 = p.species(0)

h2_new = ct.Species(h2.name, h2.composition, h2.charge, h2.size)

new_coeff = p.species(0).thermo.coeffs.copy()

new_coeff[1] = -2350000

h2_new.thermo = ct.ConstantCp(h2.thermo.min_temp, h2.thermo.max_temp, h2.thermo.reference_pressure, new_coeff)

p = ct.Solution(thermo='IdealSolidSolution', species=[h2_new]+p.species()[1:])

 

I’ve also attached the xml files. Thanks a mil.

Nwike

[Nick Curtis]

Hi Nwike,

Sorry, I didn't look closely enough at your scripts!  I didn't realize you weren't working with ideal gasses (which is 99% of what I use Cantera for).

You will likely have to solve both of these problems to get your script to work, starting with the second:

This error appears to be thrown from here, so the good news is that the Solution class is at least trying to initialize the IdealSolidSolution (which, neat! I didn't know it would do that!), and the docs say:

>In this equation of state implementation, the density is a function only of the mole fractions.

>Therefore, it can't be an independent variable. Instead, the pressure is used as the independent variable. Functions which try to set the thermodynamic state by calling setDensity() may cause an exception to be thrown.

It appears to me that the phase is being created using this call, which assumes a formGC=0.

I'm not quite able to figure out exactly where the incorrect set density call is coming from at the moment.

The best way to do so is likely:

1. Compile a debug version of cantera

2. Save the second script as 'test.py'

3. run 'gdb python'

4. inside gdb, type 'run test.py'

and post the backtrace that GDB spits out.

I'll try to take a look at it later if possible as well (or Ray or Bryan might have a better idea of exactly how the ThermoPhase initialization from python works, I admit I'm only casually familiar with it)

Best,

Nick

[Nwike]

 Thanks, Nick. I'm trying to setup the 'debug run' but I've not quite figured out how to install gdb on my windows machine. I'll try to I get that figured out first.

Nwike

[Nick Curtis]

Nwike,

I made at least some progress here...

First, on the dev-branch of cantera (i.e., compiled from be011b4) the second script fails with a different (and more obvious) error, namely it complains that the molar volume of the copied species is not specified.

One path I see to make this possible is to define the AnyMap class in _cantera.pyd, such that the field extra can be accessed in cython / python.

Then, you could utilize the kwargs of the python Species object to initialize the extra map.

At this point it may be worth opening an issue on the tracker such that species with non-default arguments (typically read in from the extra map) can be initialized via python.

Additionally, it would be a good idea to get Ray Speth's feedback here, to make sure that the solution I see is one that's worth pursuing.

Nick

[Nwike]

Hi Nick

This sounds promising. I'm not sure what it takes to make this update, but if it is straightforward, that'll be awesome. We can   Ray's thoughts on this.

Nwike

[Ray Speth]

Nwike,

The reason your initial attempt to modify the Species object in place did not work is because changes to Species objects do not take effect until the modify_species function is called. (as noted in a recent update to the documentation). The calls to make should look something like this:

k = p.species_index(name)
S = p.species(k)
S.thermo = ct.ConstantCp(...)
p.modify_species(k, S)

I don’t really understand the purpose of introducing the “dummy” element in your other example. In any case, it is not currently possible to add elements explicitly using the Python interface — they must come from an input file, or are added automatically from the set of actual elements. Assuming the species with the dummy element is not the one you would like to change properties of, you could create the XML definition of a phase containing only that species, and then add the other species using the add_species method of ThermoPhase.

All of these capabilities for dynamically adding and modifying species are mostly targeted at ideal gases at the moment. If they work for anything else, it’s a happy accident. There have been some changes under the hood, such as the introduction of the ‘extra’ field that Nick noted, that will eventually provide a mechanism for supplying such properties from Python, but that is still some way off.

Regards,
Ray

[Nwike]

Hello Ray

 It works now once I add the modify_species() call.

I actually added the dummy element to a copy of the default elements.xml file (renamed elementz.xml), and load it with the simulation - so I wasn't trying to add a dummy element from the python interface. I brought this up because of the 'element not found' cantera error when I tried the ct.Solution(..) line using Nick's sample code. 

Anyway, thanks a lot, Nick & Ray for your help with this.

Nwike

[William Michalak]

Hi Nwike,

This is a useful thread. Now that you have worked out the kinks, would you be so kind to post your final working version?

Best,

William

[Nwike]

Sure William

So a quick recap: the problem was that I had a parameter fitting model which updates standard enthalpy values for species at every iteration. The original implementation approached it this way:

At each iteration:

1. update standard enthalpy value for species k of phase, p
2. write this back to the species k data in an xml file 
3. read back the species k data from the xml file
4. proceed to fitting
   

The current implementation (based on recommendations from Nick and Ray) replaces steps 2 and 3 with the direct update, so there is no need to write to/read from physical memory at each iteration. I have included the updated implementation in the attached files. 

fit.py: main script for running the parameter fit and plotting results

meth.py: support function definitions (the update step takes place in function 'obj' in the 'meth.py' file

twophase.xml: main xml file

elementz.xml: updated xml elements file to include inert diluent

data.csv: experimental data table

Regards

Nwike

[William Michalak]

Hi Nwike,

Thanks for posting the final solution.

I am interesting in extending this approach to optimizing a reaction network, where I want to tune the A and E values to minimize the least squares on product composition. Nick pointed out how to access thermo properties and coefficients, but I don't seem to be able to find the reaction coefficients. Does anyone know how I can access these in the same fashion?

Thanks,

William

[Nwike]

Hi William

Sorry for the delay in responding to this. I have not done something similar before, but I think you can handle the update step using the Cantera Arrhenius object:

import cantera as ct

g = ct.Solution('gri30.xml') # source file

n = 2 # index of reaction of interest

A,b,E = 1,2,3 # Arhenius kinetic parameters
 g.reaction(n).rate =  ct.Arrhenius(A,b,E) # update step

Let me know if this is not what you're looking for.

Nwike

[Ray Speth]

Hi,

As in the case of modifying the thermo, this will not work as written. You need to take the Reaction object, modify the rate, and then call the modify_reaction function, as noted in the documentation.

Regards,
Ray