diff --git a/include/cantera/cython/wrappers.h b/include/cantera/cython/wrappers.h index 4025d8c07..0b50207d3 100644 --- a/include/cantera/cython/wrappers.h +++ b/include/cantera/cython/wrappers.h @@ -37,6 +37,13 @@ public: } }; +// Function for assigning elements of Array2D, since Cython has trouble +// with assigning to the reference returned by operator() +void CxxArray2D_set(Cantera::Array2D& array, size_t i, size_t j, double value) +{ + array(i,j) = value; +} + // Function which populates a 1D array #define ARRAY_FUNC(PREFIX, CLASS_NAME, FUNC_NAME) \ void PREFIX ## _ ## FUNC_NAME(Cantera::CLASS_NAME* object, double* data) \ diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index 4b52c8206..68422604c 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -61,6 +61,13 @@ cdef extern from "cantera/base/smart_ptr.h": T* get() void reset(T*) +cdef extern from "cantera/base/Array.h" namespace "Cantera": + cdef cppclass CxxArray2D "Cantera::Array2D": + CxxArray2D() + CxxArray2D(size_t, size_t) + void resize(size_t, size_t) + double operator()(size_t, size_t) + cdef extern from "cantera/thermo/SpeciesThermoInterpType.h": cdef cppclass CxxSpeciesThermo "Cantera::SpeciesThermoInterpType": CxxSpeciesThermo() @@ -283,6 +290,7 @@ cdef extern from "cantera/kinetics/Reaction.h" namespace "Cantera": CxxPlog rate cdef cppclass CxxChebyshevRate "Cantera::ChebyshevRate": + CxxChebyshevRate(double, double, double, double, CxxArray2D) double Tmin() double Tmax() double Pmin() @@ -686,6 +694,9 @@ cdef extern from "cantera/cython/wrappers.h": cdef void CxxSetLogger "setLogger" (CxxPythonLogger*) + # workaround for Cython assignment limitations + cdef void CxxArray2D_set(CxxArray2D, size_t, size_t, double) + # ThermoPhase composition cdef void thermo_getMassFractions(CxxThermoPhase*, double*) except + cdef void thermo_setMassFractions(CxxThermoPhase*, double*) except + diff --git a/interfaces/cython/cantera/reaction.pyx b/interfaces/cython/cantera/reaction.pyx index 93674f3ea..39b67761b 100644 --- a/interfaces/cython/cantera/reaction.pyx +++ b/interfaces/cython/cantera/reaction.pyx @@ -384,6 +384,20 @@ cdef class ChebyshevReaction(Reaction): c = np.fromiter(r.rate.coeffs(), np.double) return c.reshape((r.rate.nTemperature(), r.rate.nPressure())) + def set_parameters(self, Tmin, Tmax, Pmin, Pmax, coeffs): + cdef CxxChebyshevReaction* r = self.reaction + + cdef CxxArray2D data + data.resize(len(coeffs), len(coeffs[0])) + cdef double value + cdef int i + cdef int j + for i,row in enumerate(coeffs): + for j,value in enumerate(row): + CxxArray2D_set(data, i, j, value) + + r.rate = CxxChebyshevRate(Tmin, Tmax, Pmin, Pmax, data) + cdef class CoverageDepenency: cdef public double a diff --git a/interfaces/cython/cantera/test/test_kinetics.py b/interfaces/cython/cantera/test/test_kinetics.py index f904e43c7..f1ed19c2a 100644 --- a/interfaces/cython/cantera/test/test_kinetics.py +++ b/interfaces/cython/cantera/test/test_kinetics.py @@ -718,3 +718,28 @@ class TestReaction(utilities.CanteraTest): gas1.forward_rate_constants[0]) self.assertNear(gas2.net_rates_of_progress[0], gas1.net_rates_of_progress[0]) + + def test_chebyshev(self): + gas1 = ct.Solution('pdep-test.cti') + species = ct.Species.listFromFile('pdep-test.cti') + + r = ct.ChebyshevReaction() + r.reactants = {'R5': 1, 'H': 1} + r.products = {'P5A': 1, 'P5B': 1} + r.set_parameters(Tmin=300.0, Tmax=2000.0, Pmin=1000, Pmax=10000000, + coeffs=[[ 5.28830e+00, -1.13970e+00, -1.20590e-01, 1.60340e-02], + [ 1.97640e+00, 1.00370e+00, 7.28650e-03, -3.04320e-02], + [ 3.17700e-01, 2.68890e-01, 9.48060e-02, -7.63850e-03], + [-3.12850e-02, -3.94120e-02, 4.43750e-02, 1.44580e-02]]) + + gas2 = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', + species=species, reactions=[r]) + + gas2.X = gas1.X = 'R5:0.3, P5A:0.6, H:0.1' + + for T,P in itertools.product([300, 500, 1500], [1e4, 4e5, 3e6]): + gas1.TP = gas2.TP = T, P + self.assertNear(gas2.forward_rate_constants[0], + gas1.forward_rate_constants[4]) + self.assertNear(gas2.net_rates_of_progress[0], + gas1.net_rates_of_progress[4])