Cantera/Python parallelization

[Kyle B]

Thanks to Nick's and Bryan's help, my PSR project from previous posts appears to be working well.  The next issue is running this code in parallel, splitting a large parameter space across multiple cores.  Is there a slick way to do this already included, or should I just use the multiprocessing package?  If the latter, given that these threads will all be running simultaneously, how best to compile/organize the output?

Thanks,

Kyle

[Nick Curtis]

Kyle,

The multiprocessing package is the way to go.  A couple of notes, you want to make sure that each thread has it's own gas/reactors objects.  

The best way to do this is to have a function that starts from, say the mechanism name and initial conditions, and runs from there.

Multiprocessing does support returning data from the called function directly, so for example you could return a tuple of the end residence time / temperature.  

Or you could figure out some way of writing to disk (unique filenames being the key here) and reading in once the runs are completed

Nick

[Bryan W. Weber]

Kyle,

There is no multiprocessing method built-in, other than that if you compile Cantera with MKL it should use multiple cores during a single solution. An example of using the multiprocessing capabilities built in to Python is in the docs for Cantera: http://cantera.github.io/docs/sphinx/html/cython/examples/transport_multiprocessing_viscosity.html

and you can also check out my Monte Carlo uncertainty code: https://github.com/bryanwweber/rcm-temperature-uncertainty/blob/master/monte-carlo-liquid.py for another example (although my example is a little bit harder to parse).

Bryan

On Tuesday, March 29, 2016 at 9:47:54 AM UTC-4, Kyle B wrote:

[Kyle B]

Awesome, thank you both for the help.

On Tuesday, March 29, 2016 at 9:47:54 AM UTC-4, Kyle B wrote:

[Kyle B]

Bryan:

In your code you have all your cases set up in arrays so you can iterate through them, but I don't think I can format it in quite the same way, and I'm at a bit of a loss as to how to pass everything to pool.  For the single-threaded version, I have it set up like this:

# Set parameters to be studied
    # Set pressure, inlet temperature, and equivalence ratio ranges
    p_list = range(1, 2)  # pressure [atm]
    phi_list = [x * 0.1 for x in range(10,11)]
    T_list = [y * 100 for y in range(3,4)]
    res_time_list = [z * 1e-6 for z in range(1,10000)]

    for p in p_list:
        for phi in phi_list:
            reactants = ics.format(phi)
            for Tin in T_list:
                for res_time in res_time_list:
                    outfile.write('{}, {}, {}, {},'.format(p, Tin, phi, res_time))
                    gas_inlet.TPX = Tin, p*ct.one_atm, reactants
                    psr(gas, gas_inlet, reactants, Tin, p, res_time, outfile, integrator)

Obviously this for-loop structure will no longer work, but I'm unsure what the parallel equivalent would be. 

Thanks,

Kyle

[Bryan W. Weber]

Kyle,

There are two methods that will be useful, zip and product.

Zip takes a group of iterables and makes a single iterable from them. In Python, lists (square brackets) and tuples (round brackets) will be the most useful iterables for your case. So using zip, to write a for loop with several variables, it would be:

list_1 = [1,2,3,4]
list_2 = [5,6,7,8]
for i, j in zip(list_1, list_2):
    print(i, j)

1, 5
2, 6
3, 7
4, 8

Product takes the Cartesian product of the iterables, so its equivalent to your nested for-loops:

from itertools import product
for i, j in product(list_1, list_2):
    print(i, j)

1, 5
1, 6
1, 7
1, 8
2, 5
etc.

So what you need to do is create a function that does the setup and runs the simulation, then pass the function and probably a list build from the product

from multiprocessing import Pool
from itertools import product

def run_psr_cases(res_time, reactants, p, Tin, whatever other arguments):
    gas = ct.Solution('gri30.xml')
    gas.TPX = Tin, p, reactants
    # etc.
    # return whatever results you want to collect
    return reac.T, reac.P, netw.time, reac.X

# Set your parameter lists
...
n_processes = 5
with Pool(n_processes) as p:
    results = p.starmap(run_psr_cases, product(all your parameters))

The results will be a list of the results that you returned. You may want to return some of your initial values so you can tell the cases apart. Note that you won't be able to write to output files in the parallel processing unless you're very careful about how you set the filename. Here's a full example, you have to save it to a script and run it like python script_name.py, otherwise it will break in the interactive mode:

from itertools import product
from multiprocessing import Pool


def my_pow(x, y):
    return x, y, "{} to the power of {} equals {}".format(x, y, x**y)

if __name__ == '__main__':
    list_1 = [1, 2, 3, 4]
    list_2 = [5, 6, 7, 8]
    print(list(zip(list_1, list_2)))
    print(list(product(list_1, list_2)))

    with Pool(5) as p:
        a = p.starmap(pow, zip(list_1, list_2))

    print(a)

    with Pool(5) as p:
        b = p.starmap(pow, product(list_1, list_2))

    print(b)

    print(my_pow(3, 2))

    with Pool(5) as p:
        c = p.starmap(my_pow, zip(list_1, list_2))

    print(c)

    with Pool(5) as p:
        d = p.starmap(my_pow, product(list_1, list_2))

    print(d)

Bryan

[Bryan W. Weber]

Kyle,

I should also mention that functions like range, zip, and product return generator objects, so if you try to print them, you won't get anything useful. The benefit of the generator objects is that when you need the next value, it gets computed just at that moment, so it saves memory if your iteration has a lot of values to go through, since it doesn't store all those values in memory. However, once you've iterated through a generator, there's no way to "rewind" them. You can convert a generator to a list by running it through the list function:

print(range(1, 2))
<range object at.... >
print(list(range(1, 2)))
[1, 2]

Lists can iterate through unlimited times.

Bryan

...

[Kyle B]

Bryan:

Thanks, product is definitely what I needed. 

In your monte carlo code you had:

result = np.array(pool.starmap(run_case,send))

So is this array being dynamically filled with the returned compressed temperatures as each process completes?  If so, in my case would I compile the inputs/results into a tuple (so I know the format coming out) and then return that? 

...

[Bryan W. Weber]

Kyle,

I think it is more precise to say that the run_case function returns a temperature, which the starmap function collects into a list that it then returns, and that returned list is passed into the array function where it is converted to an array. So the array is not being dynamically filled, but just filled all at once at the end of the process (actually, Numpy arrays can't be dynamically resized).

The starmap function automatically returns a list whose elements are the return value of the passed in function. So if you return a tuple from your function, you'll get a list of tuples. If you return a list from your function, you'll get a list of lists. By default, returning multiple values from a Python function wraps them in a tuple, but you can manually wrap them in a list and return that instead:

# ... code here
return [val_1, val_2, val_3]

would return a single list whose elements are shown, and that single list would get wrapped into a bigger list by starmap. You can also return whatever you want, a dictionary, a NamedTuple, anything really.

Bryan

[Kyle B]

As an update, this code now runs using a single process in its "parallelized" form.  However, when I start adding processes the code breaks because multiple process are trying to access the same mechanism file simultaneously through the gas = ct.Solution(grimech30.cti) line.  I've tried to create an xml in the parent process to avoid a race condition, but it does not seem like the xml is getting created.  I've attached a stripped copy of the current code. 

Any assistance on resolving this issue would be greatly appreciated.

Thanks,

Kyle

On Tuesday, March 29, 2016 at 9:47:54 AM UTC-4, Kyle B wrote:

[Nick Curtis]

Kyle,

It looks like I was wrong about the race condition (I think the whole cti > xml part stopped being a thing awhile ago, but I haven't touched the Cantera source in some time).

I think the easiest way (which looks like I can't test fully, as the versions of Cantera I have installed are all 'thread-safe') is to use the multiprocess Lock functionality

See the attached, but the basic idea is that we're manually making it so only one Solution object can be created at a time.

Nick

[Kyle B]

The error output from this is below:

C:\Work>python PSR_min.py
Traceback (most recent call last):
  File "PSR_min.py", line 39, in <module>
    gas = ct.Solution('grimech30.cti')
  File "cantera\base.pyx", line 26, in cantera._cantera._SolutionBase.__cinit__
(interfaces\cython\cantera\_cantera.cpp:6785)
  File "cantera\base.pyx", line 47, in cantera._cantera._SolutionBase._init_cti_
xml (interfaces\cython\cantera\_cantera.cpp:7155)
RuntimeError:
***********************************************************************
CanteraError thrown by get_modified_time:
Couldn't open file:./grimech30.cti

[Ray Speth]

Kyle,

It looks like there is a bug in the get_modified_time function on Windows. I have pushed a commit to the development version that I believe fixes the problem. As a workaround, you can use the following (admittedly hackish) "try until you succeed" approach to loading the mechanism file:

import time
while True:
    try:
        gas = ct.Solution('gri30.cti')
        break
    except Exception:
        print('...')
        time.sleep(np.random.random())

My other suggestion is that you should move the creation of Solution objects into the global scope of your script, so that you only have to create the objects once for each worker process, rather than once for each simulation.

Regards,

Ray

[Bryan W. Weber]

Ray,

Kyle and I were working on some offline debugging, where I asked him to move the Solution objects to the local scope (i.e. once per simulation), because I thought there were some memory management issues with using a Solution from the global scope in each simulation. Was I wrong in that memory, or have the problems been fixed?

Thanks!

Bryan

[Ray Speth]

Bryan,

Maybe there's something I'm forgetting, but I don't think creating the Solution objects in the global scope would have ever been problematic unless you were accessing them from multiple threads within the same process. The problems created by the old file-based CTI to XML conversion would have created race conditions no matter where you create the objects, but that should be a complete non-issue now that we just do that in-memory.

Regards,

Ray

[Kyle B]

Ray:

I'll give the workaround a try in the morning, and move the Solution objects back into general scope.  Is the development version the tar.gz/do you want me to give it a try on my end?

Bryan:

Do you still want me to try the version compiled with Boost that we discussed earlier?

On Wednesday, March 30, 2016 at 4:33:04 PM UTC-4, Kyle B wrote:

[Ray Speth]

Kyle,

If you'd like to give the development version a try, you should check out the code using Git and compile the 'master' branch. Instructions are available here if you need them.

Regards,

Ray

[Bryan W. Weber]

Kyle,

Yeah, I'd be interested to know if it works, but more or less just for curiosity. If you want, I'll send you an MSI installer of the master branch this afternoon, or you can follow the instructions Ray linked to if you have Visual Studio on your computer.

Bryan

[Kyle B]

Ray:

When I move the gas objects back into global scope I run into an issue passing those two objects to the psr function via the starmap function, e.g.:

result = np.array(p.starmap(psr, product(gas, gas_inlet, p_list, phi_list, T_list, res_time_list))

What is the appropriate way to pass these objects?

Thanks,

Kyle

[Bryan W. Weber]

Kyle,

If they're in the global scope, you shouldn't have to pass them in.

import cantera as ct

gas = ct.Solution('gri30.xml')

def psr(Tin):
    gas.TP = Tin, None

if __name__ == '__main__':
    print(gas.T)
    psr(600)
    print(gas.T)

doesn't give me any errors and produces the expected output. What is the error message you receive?

Bryan

[Kyle B]

My bad, I had them in __main__, not global.

[Bryan W. Weber]

Kyle,

If you had them under the if statement, that's still in the global scope:

import cantera as ct

def psr(Tin):
    gas.TP = Tin, None

if __name__ == '__main__':
    gas = ct.Solution('gri30.xml')

    print(gas.T)
    psr(600)
    print(gas.T)

also works as expected. The key is that you don't need to pass them in.

Bryan

[Kyle B]

Ray/Bryan:

I am looking into running the development version now, but I am having a recurring issue where Cantera can't find the .cti file.  It happens somewhat randomly, and outputs errors like this:

C:\Work>python PSR_min.py
Traceback (most recent call last):

  File "<string>", line 1, in <module>
  File "C:\Program Files\Python35\lib\multiprocessing\spawn.py", line 106, in sp
awn_main
    exitcode = _main(fd)
  File "C:\Program Files\Python35\lib\multiprocessing\spawn.py", line 115, in _m
ain
    prepare(preparation_data)
  File "C:\Program Files\Python35\lib\multiprocessing\spawn.py", line 226, in pr
epare
    _fixup_main_from_path(data['init_main_from_path'])
  File "C:\Program Files\Python35\lib\multiprocessing\spawn.py", line 278, in _f
ixup_main_from_path
    run_name="__mp_main__")
  File "C:\Program Files\Python35\lib\runpy.py", line 240, in run_path
    pkg_name=pkg_name, script_name=fname)
  File "C:\Program Files\Python35\lib\runpy.py", line 96, in _run_module_code
    mod_name, mod_spec, pkg_name, script_name)
  File "C:\Program Files\Python35\lib\runpy.py", line 85, in _run_code
    exec(code, run_globals)
  File "C:\Work\PSR_min.py", line 8, in <module>

    gas = ct.Solution('grimech30.cti')

  File "interfaces\cython\cantera\base.pyx", line 26, in cantera._cantera._Solut
ionBase.__cinit__ (interfaces\cython\cantera\_cantera.cpp:6934)
  File "interfaces\cython\cantera\base.pyx", line 47, in cantera._cantera._Solut
ionBase._init_cti_xml (interfaces\cython\cantera\_cantera.cpp:7323)
RuntimeError:
***********************************************************************
CanteraError thrown by findInputFile:
Input file grimech30.cti not found in directories
'C:\Program Files\Python35\lib\site-packages\cantera\data',
'.',
'C:\Program Files\Cantera\data',
'C:\Program Files\Cantera\data'
To fix this problem, either:
    a) move the missing files into the local directory;
    b) define environment variable CANTERA_DATA to
         point to the directory containing the file.
***********************************************************************

I get this sporadically on the single-thread version of this code as well, particularly if I make a copy of the .cti file (like attaching it to an email).  Any idea what this might be about?

[Kyle B]

Odd.  It didn't work for me when I had the declarations in that arrangement, but worked fine as soon as I moved them out ...

[Bryan W. Weber]

Kyle,

Oh, you know what, you're right. The difference is how the Pool interacts with the code. So the reason for the if __name__ == '__main__': has to do with the difference between running the code as a script, and importing it. When a Python file is run as a script (i.e., python script_name.py), the __name__ variable has the value '__main__', and the code in the if statement gets run. When the Python file is imported, the __name__ variable has some other value, so the check fails and that code is not run. When the Pool is started, it imports your code, so anything in the if statement will not be run, and this includes the definition of the Solution if you put it there. The reason it works in my sample above is because I only ran it as a script, I never imported it anywhere.

Bryan

[Kyle B]

Ray:

I just tried the development version with your commit from yesterday (thanks for the compile Bryan), and your bug fix seems to have resolved the issue with processes stepping on each others toes.  I can run any number of processes now, no errors.  

Thanks so much for the help.

Kyle

[Kyle B]

Ray/Nick/Bryan:

I've been playing with this PSR code, and I've run across a weird error.  I am setting my residence time with the line:

    res_time_list = [z  for z in np.logspace(-6, -2, 2000)]

When the top end is -2, the code runs fine.  However when I try to increase the end point to -1 and higher, I get a series of these kind of errors:

******************************************************
CanteraError thrown by Phase::setDensity():
density must be positive
******************************************************

Followed by a CVODE integration error.

Any idea what is causing the temperature/density to go negative?  I've attached the full program below.

[Kyle B]

Also, this problem seems to go away when I remove the fuel rich phi = 2.0 condition from the data set.

[Ray Speth]

Kyle,

The "density must be positive" warnings aren't fatal errors, but they are an indicator that the integrator (CVODES) is having trouble finding a solution at a particular time step, so much so that it guessed a negative value for the density.

There are a couple of things that you're doing in your code that make the problem stiffer than it needs to be, which are what I think is leading to the integration errors. First, your method for setting the initial state of the reactor:

gas.TPX = Tinit, p*ct.one_atm, reactants
gas.equilibrate('hp')
gas.TPX = gas.T-300, p*ct.one_atm, reactants

means that you're starting out with an ignitable mixture at a very high temperature, so you're going to get a very rapid ignition of this mixture as the temperature heads toward a value that's much higher than the adiabatic flame temperature for the mixture, even though that's essentially the upper bound for the temperature of the reactor. I think you would be better off just taking the equilibrated mixture based on the original inlet temperature and composition as your initial condition. This has the virtue of being the steady-state solution in the limiting case of long residence times, so should always be reasonably close to the reacting steady-state solution. For this, all you need to do is:

gas.TPX = Tin, p*ct.one_atm, reactants
gas.equilibrate('hp')

Second, you need to pay a bit of attention to the units on the pressure coefficient of the PressureController before picking a value, which are kg/s/Pa. The value you set (100) implies an extra outlet flow rate of 100 kg/s for a 1 Pa increase in the reactor pressure. This makes it very hard for the integrator to solve the mass conservation equation. If you use a much smaller value (say, 0.1) the integrator will have an  easier time, and I think you'll see that it makes no difference to your steady-state results.

Regards,

Ray

[Kyle B]

Hi Ray:

It looks like the pressure controller was the biggest culprit, so I think I'm all set on this front.

The last (I hope) issue is in my validation against Chemkin, which I am using as a "sanity check".  I have included a graph below which shows the Cantera results for 1 atm, phi = 0.5, Tinlet = 300, compared against the Chemkin results for the same conditions (also using the transient solver).  As you can see, the two sets of results vary noticeably (~1.5%) in this case; other conditions show similar kinds of differences.  When I tried to extract the pressure in the reactor using the lines:

        if i == timestep-1:
            return p, phi, Tin, res_time, reactor.T, gas.P

the returned pressure was not the same as the inlet pressure.  This suggested to me that the reactor might not actually be operating at the pressure I want, which would account for the difference in results (I've attached an example output file, pressure is in the last column).  To try to avoid this I tried using the constant pressure reactor, however it turns out that these results lie exactly on top of the previous Cantera-derived results.  At this point I'm at a bit of a loss to explain the differing behavior.  Do you have any thoughts as to what might be the cause?

Thanks,

Kyle

[Ray Speth]

Kyle,

If you run a single case and plot the pressure as a function of time, you'll see that it has not reached a steady state value in those cases where the pressure is not reaching the input value. Given the wide parameter range that you're exploring, what I think you need to do is calculate a value for the pressure controller constant K which gives you a similar time constant for pressure equilibration to the reactor residence time for each simulation.

The other thing that you can also consider doing is instead of integrating to a fixed multiple of the residence time, integrate until certain key reactor parameters reach steady-state values to within some tolerance.

Regards,

Ray

[Kyle B]

Ray:

You nailed it.  Just to run a quick check I set the integration to run to 100 times the residence time, and it comes out right on top of the validation data.  Per your suggestion I switched it into a while-loop format so that it reaches steady state, and it matches up perfectly now.

Thanks,

Kyle

[Kyle B]

As a follow up to this earlier question, I am now trying to modify the upstream of the PSR with output from an equilibrium calculation.  This equilibrium function iterates in order to find the fuel fraction appropriate to a given temperature rise, and then I need to take the output from this function and pass it in to the PSR function.  To do this in parallel I have the following:

    with Pool(nProcs) as p:
        vit_result = np.array(p.starmap(vitiator, product(p_list, Tvit_list, Tvit_target_list, T_tol)), dtype=object)
       
        result = np.array(p.starmap(psr, [vit_result[0,5], product(T_rise_list, Eff_list, T_tol), vit_result[0,3], vit_result[0,6]]))

What I am trying to do is feed the output from the vit_result array into the PSR function as input.  However when I run the above it gives me an error that 'float' is not iterable.  Starmap is apparently trying to iterate over my vit_result input, but really all I want to do is iterate on the stuff inside the product statement, and just feed the other items into the function as regular input.  What is the appropriate way to accomplish this?

Thanks,

Kyle

[Ray Speth]

Kyle,

It's kind of hard to tell what's happening in the provided code, since we don't actually have the definitions for any of your variables. For example, I don't know which are scalars and which are arrays of some dimension. That said, what I think you want to do when calling starmap is pass all of your arguments to product, with arguments that should be the same for each iteration in a list containing just the one value. For example:

>>> f = lambda x,y: (x,y)
>>> list(itertools.starmap(f, itertools.product(range(4), [99])))
[(0, 99), (1, 99), (2, 99), (3, 99)]

Regards,

Ray

[Kyle B]

Thanks for the reply Ray.

I am returning the following values from the "upstream" function

return p, phi, Tvit_in, Tpsr_in, XO2, Ppsr_in, XALL

vit_result(0,3) is Tpsr_in, vit_result(0,5) is Ppsr_in, and vit_result(0,6) is XALL.  The first two are sacalars, however XALL is an array defined by gas[:].X.  My thinking was that product would try to iterate through the values inside XALL, which is not what I want to do, hence putting it outside product.

At the end of the day, what I want to do is get the gas state from the function vitiator to the function psr.  I had tried defining it in global scope thinking that the state would carry over, but that is not the case so far as I can tell.

Thanks,

Kyle

[Kyle B]

Hi Ray:

I've resolved most of the issues with the starmap input, but the one remaining issue is that the XALL list (which contains the species mole fractions) is interpreted by product as a series of input values instead of as a single unit that needs to be given as input to each process.  I've attached the working copy of the code.  Do you have any suggestion as to how to pass this array between functions?

Thanks,

Kyle

On Thursday, May 12, 2016 at 11:20:21 PM UTC-4, Ray Speth wrote:

[Ray Speth]

Kyle,

I think what you want to do is simply pass a list containing XALL as its one element to starmap, e.g.:

result = np.array(p.starmap(psr, product(Ppsr_in, T_rise_list, Eff_list, T_tol, Tpsr_in, [XALL])), dtype=object)

Regards,

Ray