diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 1a5ec007b..47fb520bb 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -300,6 +300,11 @@ cdef extern from "cantera/kinetics/Kinetics.h" namespace "Cantera": CxxThermoPhase& thermo(int) + void addPhase(CxxThermoPhase&) except + + void init() except + + void addReaction(shared_ptr[CxxReaction]) except + + void finalize() except + + shared_ptr[CxxReaction] reaction(size_t) except + cbool isReversible(int) except + int reactionType(int) except + diff --git a/interfaces/cython/cantera/base.pyx b/interfaces/cython/cantera/base.pyx index a254c58ad..a58dfe080 100644 --- a/interfaces/cython/cantera/base.pyx +++ b/interfaces/cython/cantera/base.pyx @@ -1,6 +1,7 @@ cdef class _SolutionBase: def __cinit__(self, infile='', phaseid='', phases=(), origin=None, - source=None, thermo=None, species=(), **kwargs): + source=None, thermo=None, species=(), kinetics=None, + reactions=(), **kwargs): # Shallow copy of an existing Solution (for slicing support) cdef _SolutionBase other if origin is not None: @@ -21,7 +22,7 @@ cdef class _SolutionBase: if infile or source: self._init_cti_xml(infile, phaseid, phases, source) elif thermo and species: - self._init_parts(thermo, species) + self._init_parts(thermo, species, kinetics, phases, reactions) else: raise ValueError("Arguments are insufficient to define a phase") @@ -72,7 +73,7 @@ cdef class _SolutionBase: else: self.kinetics = NULL - def _init_parts(self, thermo, species): + def _init_parts(self, thermo, species, kinetics, phases, reactions): """ Instantiate a set of new Cantera C++ objects based on a string defining the model type and a list of Species objects. @@ -84,9 +85,21 @@ cdef class _SolutionBase: self.thermo.addSpecies(S._species) self.thermo.initThermo() + if not kinetics: + kinetics = "none" + + cdef ThermoPhase phase + cdef Reaction reaction if isinstance(self, Kinetics): - # Not yet implemented - self.kinetics = CxxNewKinetics(stringify("none")) + self.kinetics = CxxNewKinetics(stringify(kinetics)) + self.kinetics.addPhase(deref(self.thermo)) + for phase in phases: + self.kinetics.addPhase(deref(phase.thermo)) + self.kinetics.init() + for reaction in reactions: + self.kinetics.addReaction(reaction._reaction) + self.kinetics.finalize() + def __getitem__(self, selection): copy = self.__class__(origin=self) diff --git a/interfaces/cython/cantera/kinetics.pyx b/interfaces/cython/cantera/kinetics.pyx index e756f4513..70f648d5a 100644 --- a/interfaces/cython/cantera/kinetics.pyx +++ b/interfaces/cython/cantera/kinetics.pyx @@ -313,7 +313,7 @@ cdef class InterfaceKinetics(Kinetics): A kinetics manager for heterogeneous reaction mechanisms. The reactions are assumed to occur at an interface between bulk phases. """ - def __init__(self, infile, phaseid='', phases=(), *args, **kwargs): + def __init__(self, infile='', phaseid='', phases=(), *args, **kwargs): super().__init__(infile, phaseid, phases, *args, **kwargs) if self.kinetics.type() not in (kinetics_type_interface, kinetics_type_edge): diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index 3723c6de0..a96e99d09 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -1,6 +1,7 @@ import unittest import numpy as np import re +import itertools import cantera as ct from . import utilities @@ -143,6 +144,86 @@ class TestKinetics(utilities.CanteraTest): self.phase.delta_standard_gibbs) +class KineticsFromReactions(utilities.CanteraTest): + """ + Test for Kinetics objects which are constructed directly from Reaction + objects instead of from CTI/XML files. + """ + def test_idealgas(self): + gas1 = ct.Solution('h2o2.xml') + + S = ct.Species.listFromFile('h2o2.xml') + R = ct.Reaction.listFromFile('h2o2.xml') + gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=S, reactions=R) + + self.assertEqual(gas1.n_reactions, gas2.n_reactions) + gas1.TPY = 800, 2*ct.one_atm, 'H2:0.3, O2:0.7, OH:2e-4, O:1e-3, H:5e-5' + gas2.TPY = gas1.TPY + + self.assertTrue((gas1.reactant_stoich_coeffs() == + gas2.reactant_stoich_coeffs()).all()) + self.assertTrue((gas1.product_stoich_coeffs() == + gas2.product_stoich_coeffs()).all()) + + self.assertArrayNear(gas1.delta_gibbs, + gas2.delta_gibbs) + self.assertArrayNear(gas1.reverse_rate_constants, + gas2.reverse_rate_constants) + self.assertArrayNear(gas1.net_production_rates, + gas2.net_production_rates) + + def test_surface(self): + gas_species = ct.Species.listFromFile('gri30.xml') + surf_species = ct.Species.listFromFile('ptcombust.xml') + reactions = ct.Reaction.listFromFile('ptcombust.xml') + + gas = ct.Solution('ptcombust.xml', 'gas') + surf1 = ct.Interface('ptcombust.xml', 'Pt_surf', [gas]) + + surf2 = ct.Interface(thermo='Surface', kinetics='interface', + species=surf_species, reactions=reactions, + phases=[gas]) + surf1.site_density = surf2.site_density = 5e-9 + gas.TP = surf2.TP = surf1.TP = 900, 2*ct.one_atm + surf2.concentrations = surf1.concentrations + + self.assertEqual(surf1.n_reactions, surf2.n_reactions) + + for k,i in itertools.product(['PT(S)','H2','OH','OH(S)'], + range(surf1.n_species)): + self.assertEqual(surf1.reactant_stoich_coeff(k,i), + surf2.reactant_stoich_coeff(k,i)) + self.assertEqual(surf1.product_stoich_coeff(k,i), + surf2.product_stoich_coeff(k,i)) + + for i in range(surf1.n_reactions): + r1 = surf1.reaction(i) + r2 = surf2.reaction(i) + self.assertEqual(r1.reactants, r2.reactants) + self.assertEqual(r1.products, r2.products) + self.assertEqual(r1.rate.preexponential_factor, + r2.rate.preexponential_factor) + self.assertEqual(r1.rate.temperature_exponent, + r2.rate.temperature_exponent) + self.assertEqual(r1.rate.activation_energy, + r2.rate.activation_energy) + + self.assertArrayNear(surf1.delta_enthalpy, + surf2.delta_enthalpy) + self.assertArrayNear(surf1.forward_rate_constants, + surf2.forward_rate_constants) + self.assertArrayNear(surf1.reverse_rate_constants, + surf2.reverse_rate_constants) + + rop1 = surf1.net_production_rates + rop2 = surf2.net_production_rates + for k in gas.species_names + surf1.species_names: + k1 = surf1.kinetics_species_index(k) + k2 = surf2.kinetics_species_index(k) + self.assertNear(rop1[k1], rop2[k2]) + + class KineticsRepeatability(utilities.CanteraTest): """ Tests to make sure that lazily evaluated of terms in the rate expression