1from collections import OrderedDict
2import numpy as np
3import gc
4
5import cantera as ct
6from . import utilities
7
8
9class TestThermoPhase(utilities.CanteraTest):
10    def setUp(self):
11        self.phase = ct.Solution('h2o2.yaml', transport_model=None)
12
13    def test_source(self):
14        self.assertEqual(self.phase.source, 'h2o2.yaml')
15
16    def test_missing_phases_key(self):
17        yaml = '''
18        species:
19        - name: H
20          composition: {H: 1}
21          thermo:
22            model: NASA7
23            temperature-ranges: [200.0, 1000.0, 6000.0]
24            data:
25            - [2.5, 0.0, 0.0, 0.0, 0.0, 2.547366e+04, -0.44668285]
26            - [2.5, 0.0, 0.0, 0.0, 0.0, 2.547366e+04, -0.44668285]
27            note: L6/94
28        '''
29        with self.assertRaisesRegex(ct.CanteraError, "Key 'phases' not found"):
30            _ = ct.Solution(yaml=yaml)
31
32    def test_base_attributes(self):
33        self.assertIsInstance(self.phase.name, str)
34        self.assertIsInstance(self.phase.phase_of_matter, str)
35        self.assertIsInstance(self.phase.thermo_model, str)
36        self.assertIsInstance(self.phase.kinetics_model, str)
37        self.assertIsInstance(self.phase.transport_model, str)
38        self.assertIsInstance(self.phase.composite, tuple)
39        self.assertEqual(len(self.phase.composite), 3)
40        self.assertEqual(self.phase.composite,
41                         (self.phase.thermo_model,
42                          self.phase.kinetics_model,
43                          self.phase.transport_model))
44        self.phase.name = 'spam'
45        self.assertEqual(self.phase.name, 'spam')
46        with self.assertRaises(AttributeError):
47            self.phase.type = 'eggs'
48
49    def test_phases(self):
50        self.assertEqual(self.phase.n_phases, 1)
51        self.assertEqual(self.phase.phase_of_matter, "gas")
52
53    def test_states(self):
54        self.assertEqual(self.phase._native_state, ('T', 'D', 'Y'))
55        self.assertIn('TPY', self.phase._full_states.values())
56        self.assertIn('TD', self.phase._partial_states.values())
57
58    def test_species(self):
59        self.assertEqual(self.phase.n_species, 10)
60        for i,name in enumerate(self.phase.species_names):
61            self.assertEqual(name, self.phase.species_name(i))
62            self.assertEqual(i, self.phase.species_index(name))
63            self.assertEqual(i, self.phase.species_index(i))
64        with self.assertRaisesRegex(ct.CanteraError, 'IndexError thrown by Phase::checkSpeciesIndex'):
65            self.phase.species(self.phase.n_species)
66
67    def test_elements(self):
68        self.assertEqual(self.phase.n_elements, 4)
69        for i,symbol in enumerate(self.phase.element_names):
70            self.assertEqual(symbol, self.phase.element_name(i))
71            self.assertEqual(i, self.phase.element_index(symbol))
72            self.assertEqual(i, self.phase.element_index(i))
73
74    def test_n_atoms(self):
75        data = [(1, 'O', 'O'), (2, 'O', 'O2'), (1, b'H', b'OH'),
76                (2, 'H', 'H2O'), (2, 'O', 'H2O2'), (1, 'Ar', 'AR'),
77                (0, 'O', 'H'), (0, 'H', 'AR'), (0, 'Ar', 'HO2')]
78        for (n, elem, species) in data:
79            self.assertEqual(self.phase.n_atoms(species, elem), n)
80            mElem = self.phase.element_index(elem)
81            kSpec = self.phase.species_index(species)
82            self.assertEqual(self.phase.n_atoms(kSpec, mElem), n)
83
84        with self.assertRaisesRegex(ValueError, 'No such species'):
85            self.phase.n_atoms('C', 'H2')
86        with self.assertRaisesRegex(ValueError, 'No such element'):
87            self.phase.n_atoms('H', 'CH4')
88
89    def test_elemental_mass_fraction(self):
90        self.phase.Y = 'H2O:0.5, O2:0.5'
91        Zo = self.phase.elemental_mass_fraction('O')
92        Zh = self.phase.elemental_mass_fraction('H')
93        Zar = self.phase.elemental_mass_fraction('Ar')
94
95        mO = self.phase.element_index('O')
96        self.assertEqual(Zo, self.phase.elemental_mass_fraction(mO))
97        self.assertNear(Zo, 0.5 + 0.5 * (15.999 / 18.015))
98        self.assertNear(Zh, 0.5 * (2.016 / 18.015))
99        self.assertEqual(Zar, 0.0)
100
101        with self.assertRaisesRegex(ValueError, 'No such element'):
102            self.phase.elemental_mass_fraction('C')
103        with self.assertRaisesRegex(ValueError, 'No such element'):
104            self.phase.elemental_mass_fraction(5)
105
106    def test_elemental_mole_fraction(self):
107        self.phase.X = 'H2O:0.5, O2:0.5'
108        Zo = self.phase.elemental_mole_fraction('O')
109        Zh = self.phase.elemental_mole_fraction('H')
110        Zar = self.phase.elemental_mole_fraction('Ar')
111
112        mO = self.phase.element_index('O')
113        self.assertEqual(Zo, self.phase.elemental_mole_fraction(mO))
114        self.assertNear(Zo, (0.5 + 1) / (0.5*3 + 0.5*2))
115        self.assertNear(Zh, (2*0.5) / (0.5*3 + 0.5*2))
116        self.assertEqual(Zar, 0.0)
117
118        with self.assertRaisesRegex(ValueError, 'No such element'):
119            self.phase.elemental_mole_fraction('C')
120        with self.assertRaisesRegex(ValueError, 'No such element'):
121            self.phase.elemental_mole_fraction(5)
122
123    def test_elemental_mass_mole_fraction(self):
124        # expected relationship between elemental mass and mole fractions
125        comps = ['H2O:0.5, O2:0.5', 'H2:0.1, O2:0.4, H2O2:0.3, AR:0.2',
126                 'O2:0.1, H2:0.9']
127        for comp in comps:
128            self.phase.X = comp
129
130            denom = sum(self.phase.elemental_mole_fraction(i)
131                        * self.phase.atomic_weight(i)
132                        for i in range(self.phase.n_elements))
133
134            for i in range(self.phase.n_elements):
135                self.assertNear(self.phase.elemental_mass_fraction(i),
136                                self.phase.elemental_mole_fraction(i)
137                                * self.phase.atomic_weight(i) / denom)
138
139    def test_weights(self):
140        atomic_weights = self.phase.atomic_weights
141        molecular_weights = self.phase.molecular_weights
142        self.assertEqual(self.phase.n_elements, len(atomic_weights))
143        self.assertEqual(self.phase.n_species, len(molecular_weights))
144
145        for i,mw in enumerate(molecular_weights):
146            test_weight = 0.0
147            for j,aw in enumerate(atomic_weights):
148                test_weight += aw * self.phase.n_atoms(i,j)
149            self.assertNear(test_weight, mw)
150
151    def test_charges(self):
152        gas = ct.Solution('ch4_ion.yaml')
153        charges = gas.charges
154        test = {'E': -1., 'N2': 0., 'H3O+': 1.}
155        for species, charge in test.items():
156            self.assertIn(species, gas.species_names)
157            index = gas.species_index(species)
158            self.assertEqual(charges[index], charge)
159
160    def test_setComposition(self):
161        X = np.zeros(self.phase.n_species)
162        X[2] = 1.0
163        self.phase.X = X
164        Y = self.phase.Y
165
166        self.assertEqual(list(X), list(Y))
167
168    def test_setComposition_singleton(self):
169        X = np.zeros((1,self.phase.n_species,1))
170        X[0,2,0] = 1.0
171        self.phase.X = X
172        Y = self.phase.Y
173
174        self.assertEqual(list(X[0,:,0]), list(Y))
175
176    def test_setCompositionString(self):
177        self.phase.X = 'h2:1.0, o2:1.0'
178        X = self.phase.X
179        self.assertNear(X[0], 0.5)
180        self.assertNear(X[3], 0.5)
181
182        with self.assertRaisesRegex(ct.CanteraError, 'Unknown species'):
183            self.phase.X = 'H2:1.0, CO2:1.5'
184
185    def test_setCompositionStringBad(self):
186        X0 = self.phase.X
187        with self.assertRaisesRegex(ct.CanteraError, 'Trouble processing'):
188            self.phase.X = 'H2:1.0, O2:asdf'
189        self.assertArrayNear(X0, self.phase.X)
190
191        with self.assertRaisesRegex(ct.CanteraError, 'Trouble processing'):
192            self.phase.X = 'H2:1e-x4'
193        self.assertArrayNear(X0, self.phase.X)
194
195        with self.assertRaisesRegex(ct.CanteraError, 'decimal point in exponent'):
196            self.phase.X = 'H2:1e-1.4'
197        self.assertArrayNear(X0, self.phase.X)
198
199        with self.assertRaisesRegex(ct.CanteraError, 'Duplicate key'):
200            self.phase.X = 'H2:0.5, O2:1.0, H2:0.1'
201        self.assertArrayNear(X0, self.phase.X)
202
203    def test_setCompositionDict(self):
204        self.phase.X = {b'H2':1.0, b'O2':3.0}
205        X = self.phase.X
206        self.assertNear(X[0], 0.25)
207        self.assertNear(X[3], 0.75)
208
209        self.phase.Y = {'H2':1.0, 'O2':3.0}
210        Y = self.phase.Y
211        self.assertNear(Y[0], 0.25)
212        self.assertNear(Y[3], 0.75)
213
214    def test_getCompositionDict(self):
215        self.phase.X = 'oh:1e-9, O2:0.4, AR:0.6'
216        self.assertEqual(len(self.phase.mole_fraction_dict(1e-7)), 2)
217        self.assertEqual(len(self.phase.mole_fraction_dict()), 3)
218
219        self.phase.Y = 'O2:0.4, AR:0.6'
220        Y1 = self.phase.mass_fraction_dict()
221        self.assertNear(Y1['O2'], 0.4)
222        self.assertNear(Y1['AR'], 0.6)
223
224    def test_setCompositionNoNorm(self):
225        X = np.zeros(self.phase.n_species)
226        X[2] = 1.0
227        X[0] = 0.01
228        self.phase.set_unnormalized_mole_fractions(X)
229        self.assertArrayNear(self.phase.X, X)
230        self.assertNear(sum(X), 1.01)
231
232        Y = np.zeros(self.phase.n_species)
233        Y[2] = 1.0
234        Y[0] = 0.01
235        self.phase.set_unnormalized_mass_fractions(Y)
236        self.assertArrayNear(self.phase.Y, Y)
237        self.assertNear(sum(Y), 1.01)
238
239    def test_setCompositionNoNormBad(self):
240        X = np.zeros(self.phase.n_species - 1)
241        with self.assertRaisesRegex(ValueError, 'incorrect length'):
242            self.phase.set_unnormalized_mole_fractions(X)
243
244        with self.assertRaisesRegex(ValueError, 'incorrect length'):
245            self.phase.set_unnormalized_mass_fractions([1,2,3])
246
247    def test_setCompositionDict_bad1(self):
248        with self.assertRaisesRegex(ct.CanteraError, 'Unknown species'):
249            self.phase.X = {'H2':1.0, 'HCl':3.0}
250
251    def test_setCompositionDict_bad2(self):
252        with self.assertRaises(TypeError):
253            self.phase.Y = {'H2':1.0, 'O2':'xx'}
254
255    def test_setCompositionSlice(self):
256        self.phase['h2', 'o2'].X = 0.1, 0.9
257        X = self.phase.X
258        self.assertNear(X[0], 0.1)
259        self.assertNear(X[3], 0.9)
260
261    def test_setCompositionSingleSpecies(self):
262        gas = ct.Solution('argon.xml')
263        gas.X = [1]
264        gas.Y = np.array([[1.001]])
265        self.assertEqual(gas.Y[0], 1.0)
266
267    def test_setCompositionSlice_bad(self):
268        X0 = self.phase.X
269        with self.assertRaisesRegex(ValueError, 'incorrect length'):
270            self.phase['H2','O2'].Y = [0.1, 0.2, 0.3]
271        self.assertArrayNear(self.phase.X, X0)
272
273    def test_setCompositionEmpty_bad(self):
274        X0 = self.phase.X
275        with self.assertRaisesRegex(ValueError, 'incorrect length'):
276            self.phase.Y = np.array([])
277        self.assertArrayNear(self.phase.X, X0)
278
279    @utilities.slow_test
280    def test_set_equivalence_ratio_stoichiometric(self):
281        gas = ct.Solution('gri30.yaml', transport_model=None)
282        for fuel in ('C2H6', 'H2:0.7, CO:0.3', 'NH3:0.4, CH3OH:0.6'):
283            for oxidizer in ('O2:1.0, N2:3.76', 'H2O2:1.0'):
284                gas.set_equivalence_ratio(1.0, fuel, oxidizer)
285                gas.equilibrate('TP')
286                # Almost everything should end up as CO2, H2O and N2
287                self.assertGreater(sum(gas['H2O','CO2','N2'].X), 0.999999)
288
289    def test_set_equivalence_ratio_lean(self):
290        gas = ct.Solution('gri30.yaml', transport_model=None)
291        excess = 0
292        for phi in np.linspace(0.9, 0, 5):
293            gas.set_equivalence_ratio(phi, 'CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76')
294            gas.equilibrate('TP')
295            self.assertGreater(gas['O2'].X[0], excess)
296            excess = gas['O2'].X[0]
297        self.assertNear(sum(gas['O2','N2'].X), 1.0)
298
299    def test_set_equivalence_ratio_sulfur(self):
300        sulfur_species = [k for k in ct.Species.listFromFile('nasa_gas.xml') if k.name in ("SO", "SO2")]
301        gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
302                          species=ct.Species.listFromFile('gri30.xml') + sulfur_species,
303                          reactions=ct.Reaction.listFromFile('gri30.xml'))
304        fuel = 'CH3:0.5, SO:0.25, OH:0.125, N2:0.125'
305        ox = 'O2:0.5, SO2:0.25, CO2:0.125, CH:0.125'
306
307        def test_sulfur_results(gas, fuel, ox, basis):
308            gas.set_equivalence_ratio(2.0, fuel, ox, basis)
309            Z = gas.mixture_fraction(fuel, ox, basis)
310            self.assertNear(gas.stoich_air_fuel_ratio(fuel, ox, basis)/((1.0-Z)/Z),  2.0)
311            gas.set_mixture_fraction(Z, fuel, ox, basis)
312            self.assertNear(gas['SO2'].X[0], 31.0/212.0)
313            self.assertNear(gas['O2'].X[0],  31.0/106.0)
314            self.assertNear(gas['SO'].X[0],  11.0/106.0)
315            self.assertNear(gas['CO2'].X[0], 31.0/424.0)
316            self.assertNear(gas['CH3'].X[0], 11.0/53.0)
317            self.assertNear(gas['N2'].X[0],  11.0/212.0)
318            self.assertNear(gas['CH'].X[0],  31.0/424.0)
319            self.assertNear(gas['OH'].X[0],  11.0/212.0)
320            self.assertNear(gas.equivalence_ratio(fuel, ox, basis),  2.0)
321
322        test_sulfur_results(gas, fuel, ox, 'mole')
323
324        gas.TPX = None, None, fuel
325        fuel = gas.Y
326        gas.TPX = None, None, ox
327        ox = gas.Y
328        test_sulfur_results(gas, fuel, ox, 'mass')
329
330    def test_equivalence_ratio(self):
331        gas = ct.Solution('gri30.yaml', transport_model=None)
332        for phi in np.linspace(0.5, 2.0, 5):
333            gas.set_equivalence_ratio(phi, 'CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76')
334            self.assertNear(phi, gas.equivalence_ratio('CH4:0.8, CH3OH:0.2', 'O2:1.0, N2:3.76'))
335        # Check sulfur species
336        sulfur_species = [k for k in ct.Species.listFromFile('nasa_gas.xml') if k.name in ("SO", "SO2")]
337        gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics',
338                          species=ct.Species.listFromFile('gri30.xml') + sulfur_species)
339        for phi in np.linspace(0.5, 2.0, 5):
340            gas.set_equivalence_ratio(phi, 'CH3:0.5, SO:0.25, OH:0.125, N2:0.125', 'O2:0.5, SO2:0.25, CO2:0.125')
341            self.assertNear(phi, gas.equivalence_ratio('CH3:0.5, SO:0.25, OH:0.125, N2:0.125', 'O2:0.5, SO2:0.25, CO2:0.125'))
342        gas.X = 'CH4:1' # pure fuel
343        self.assertEqual(gas.equivalence_ratio(), np.inf)
344
345    def test_get_set_equivalence_ratio_functions(self):
346        fuel = "CH4:0.2,O2:0.02,N2:0.1,CO:0.05,CO2:0.02"
347        ox = "O2:0.21,N2:0.79,CO:0.04,CH4:0.01,CO2:0.03"
348
349        gas = ct.Solution('gri30.yaml', transport_model=None)
350        gas.TPX = 300, 1e5, fuel
351        Y_Cf = gas.elemental_mass_fraction("C")
352        Y_Of = gas.elemental_mass_fraction("O")
353        gas.TPX = 300, 1e5, ox
354        Y_Co = gas.elemental_mass_fraction("C")
355        Y_Oo = gas.elemental_mass_fraction("O")
356
357        def test_equil_results(gas, fuel, ox, Y_Cf, Y_Of, Y_Co, Y_Oo, basis):
358            gas.TP = 300, 1e5
359            gas.set_equivalence_ratio(1.3, fuel, ox, basis)
360            T = gas.T
361
362            # set mixture to burnt state to make sure that equivalence ratio and
363            # mixture fraction are independent of reaction progress
364            gas.equilibrate("HP");
365
366            phi = gas.equivalence_ratio(fuel, ox, basis)
367            phi_loc = gas.equivalence_ratio()
368            mf = gas.mixture_fraction(fuel, ox, basis)
369            mf_C = gas.mixture_fraction(fuel, ox, basis, element="C")
370            mf_O = gas.mixture_fraction(fuel, ox, basis, element="O")
371            l = gas.stoich_air_fuel_ratio(fuel, ox, basis)
372
373            gas.set_mixture_fraction(mf, fuel,ox, basis)
374            phi2 = gas.equivalence_ratio(fuel, ox, basis)
375
376            self.assertNear(phi, 1.3)
377            self.assertNear(phi2, 1.3)
378            self.assertNear(phi_loc, 1.1726068608195617)
379            self.assertNear(mf, 0.13415725911057605)
380            self.assertNear(mf_C, (gas.elemental_mass_fraction("C")-Y_Co)/(Y_Cf-Y_Co))
381            self.assertNear(mf_O, (gas.elemental_mass_fraction("O")-Y_Oo)/(Y_Of-Y_Oo))
382            self.assertNear(l, 8.3901204498353561)
383            self.assertNear(gas.P, 1e5)
384            self.assertNear(T, 300.0)
385
386        test_equil_results(gas, fuel, ox, Y_Cf, Y_Of, Y_Co, Y_Oo, 'mole')
387
388        # do the same for mass-based functions
389
390        gas.TPX = None, None, fuel
391        fuel = gas.Y
392        Y_Cf = gas.elemental_mass_fraction("C")
393        Y_Of = gas.elemental_mass_fraction("O")
394        gas.TPX = None, None, ox
395        ox = gas.Y
396        Y_Co = gas.elemental_mass_fraction("C")
397        Y_Oo = gas.elemental_mass_fraction("O")
398
399        test_equil_results(gas, fuel, ox, Y_Cf, Y_Of, Y_Co, Y_Oo, 'mass')
400
401    def test_full_report(self):
402        report = self.phase.report(threshold=0.0)
403        self.assertIn(self.phase.name, report)
404        self.assertIn('temperature', report)
405        self.assertNotIn('minor', report)
406        for name in self.phase.species_names:
407            self.assertIn(name, report)
408
409    def test_default_report(self):
410        self.phase.X = 'H2:0.1, O2:0.9, HO2:1e-10, H2O2:1e-20'
411        report = self.phase.report()
412        self.assertIn('minor', report)
413        for name in (' H2 ', ' O2 ', ' HO2 '):
414            self.assertIn(name, report)
415        for name in (' H2O2 ', ' OH ', ' AR '):
416            self.assertNotIn(name, report)
417
418    def test_name(self):
419        self.phase.name = 'something'
420        self.assertEqual(self.phase.name, 'something')
421        self.assertIn('something', self.phase.report())
422
423    def test_badLength(self):
424        X = np.zeros(5)
425        with self.assertRaisesRegex(ValueError, 'incorrect length'):
426            self.phase.X = X
427        with self.assertRaisesRegex(ValueError, 'incorrect length'):
428            self.phase.Y = X
429
430    def test_mass_basis(self):
431        self.assertEqual(self.phase.basis, 'mass')
432        self.assertNear(self.phase.density_mass, self.phase.density)
433        self.assertNear(self.phase.enthalpy_mass, self.phase.h)
434        self.assertNear(self.phase.entropy_mass, self.phase.s)
435        self.assertNear(self.phase.int_energy_mass, self.phase.u)
436        self.assertNear(self.phase.volume_mass, self.phase.v)
437        self.assertNear(self.phase.cv_mass, self.phase.cv)
438        self.assertNear(self.phase.cp_mass, self.phase.cp)
439
440    def test_molar_basis(self):
441        self.phase.basis = 'molar'
442        self.assertEqual(self.phase.basis, 'molar')
443        self.assertNear(self.phase.density_mole, self.phase.density)
444        self.assertNear(self.phase.enthalpy_mole, self.phase.h)
445        self.assertNear(self.phase.entropy_mole, self.phase.s)
446        self.assertNear(self.phase.int_energy_mole, self.phase.u)
447        self.assertNear(self.phase.volume_mole, self.phase.v)
448        self.assertNear(self.phase.cv_mole, self.phase.cv)
449        self.assertNear(self.phase.cp_mole, self.phase.cp)
450
451    def check_setters(self, T1, rho1, Y1):
452        T0, rho0, Y0 = self.phase.TDY
453        self.phase.TDY = T1, rho1, Y1
454        X1 = self.phase.X
455        P1 = self.phase.P
456        h1 = self.phase.h
457        s1 = self.phase.s
458        u1 = self.phase.u
459        v1 = self.phase.v
460
461        def check_state(T, rho, Y):
462            self.assertNear(self.phase.T, T)
463            self.assertNear(self.phase.density, rho)
464            self.assertArrayNear(self.phase.Y, Y)
465
466        self.phase.TDY = T0, rho0, Y0
467        self.phase.TPY = T1, P1, Y1
468        check_state(T1, rho1, Y1)
469
470        self.phase.TDY = T0, rho0, Y0
471        self.phase.UVY = u1, v1, Y1
472        check_state(T1, rho1, Y1)
473
474        self.phase.TDY = T0, rho0, Y0
475        self.phase.HPY = h1, P1, Y1
476        check_state(T1, rho1, Y1)
477
478        self.phase.TDY = T0, rho0, Y0
479        self.phase.SPY = s1, P1, Y1
480        check_state(T1, rho1, Y1)
481
482        self.phase.TDY = T0, rho0, Y0
483        self.phase.TPX = T1, P1, X1
484        check_state(T1, rho1, Y1)
485
486        self.phase.TDY = T0, rho0, Y0
487        self.phase.UVX = u1, v1, X1
488        check_state(T1, rho1, Y1)
489
490        self.phase.TDY = T0, rho0, Y0
491        self.phase.HPX = h1, P1, X1
492        check_state(T1, rho1, Y1)
493
494        self.phase.TDY = T0, rho0, Y0
495        self.phase.SPX = s1, P1, X1
496        check_state(T1, rho1, Y1)
497
498        self.phase.TDY = T0, rho0, Y0
499        self.phase.SVX = s1, v1, X1
500        check_state(T1, rho1, Y1)
501
502        self.phase.TDY = T0, rho0, Y0
503        self.phase.SVY = s1, v1, Y1
504        check_state(T1, rho1, Y1)
505
506        self.phase.TDY = T0, rho0, Y0
507        self.phase.DPX = rho1, P1, X1
508        check_state(T1, rho1, Y1)
509
510        self.phase.TDY = T0, rho0, Y0
511        self.phase.DPY = rho1, P1, Y1
512        check_state(T1, rho1, Y1)
513
514    def test_setState_mass(self):
515        self.check_setters(T1 = 500.0, rho1 = 1.5,
516                           Y1 = [0.1, 0.0, 0.0, 0.1, 0.4, 0.2, 0.0, 0.0, 0.2, 0.0])
517
518    def test_setState_mole(self):
519        self.phase.basis = 'molar'
520        self.check_setters(T1 = 750.0, rho1 = 0.02,
521                           Y1 = [0.2, 0.1, 0.0, 0.3, 0.1, 0.0, 0.0, 0.2, 0.1, 0.0])
522
523    def test_setters_hold_constant(self):
524        props = ('T','P','s','h','u','v','X','Y')
525        pairs = [('TP', 'T', 'P'), ('SP', 's', 'P'),
526                 ('UV', 'u', 'v')]
527
528        self.phase.TDX = 1000, 1.5, 'H2O:0.1, O2:0.95, AR:3.0'
529        values = {}
530        for p in props:
531            values[p] = getattr(self.phase, p)
532
533        for pair, first, second in pairs:
534            self.phase.TDX = 500, 2.5, 'H2:0.1, O2:1.0, AR:3.0'
535            first_val = getattr(self.phase, first)
536            second_val = getattr(self.phase, second)
537
538            setattr(self.phase, pair, (values[first], None))
539            self.assertNear(getattr(self.phase, first), values[first])
540            self.assertNear(getattr(self.phase, second), second_val)
541
542            self.phase.TDX = 500, 2.5, 'H2:0.1, O2:1.0, AR:3.0'
543            setattr(self.phase, pair, (None, values[second]))
544            self.assertNear(getattr(self.phase, first), first_val)
545            self.assertNear(getattr(self.phase, second), values[second])
546
547            self.phase.TDX = 500, 2.5, 'H2:0.1, O2:1.0, AR:3.0'
548            setattr(self.phase, pair + 'X', (None, None, values['X']))
549            self.assertNear(getattr(self.phase, first), first_val)
550            self.assertNear(getattr(self.phase, second), second_val)
551
552            self.phase.TDX = 500, 2.5, 'H2:0.1, O2:1.0, AR:3.0'
553            setattr(self.phase, pair + 'Y', (None, None, values['Y']))
554            self.assertNear(getattr(self.phase, first), first_val)
555            self.assertNear(getattr(self.phase, second), second_val)
556
557    def test_setter_errors(self):
558        with self.assertRaises(TypeError):
559            self.phase.TD = 400
560
561        with self.assertRaisesRegex(AssertionError, 'incorrect number'):
562            self.phase.TP = 300, 101325, 'CH4:1.0'
563
564        with self.assertRaisesRegex(AssertionError, 'incorrect number'):
565            self.phase.HPY = 1.2e6, 101325
566
567        with self.assertRaisesRegex(AssertionError, 'incorrect number'):
568            self.phase.UVX = -4e5, 4.4, 'H2:1.0', -1
569
570    def test_invalid_property(self):
571        x = self.phase
572        with self.assertRaises(AttributeError):
573            x.foobar = 300
574        with self.assertRaises(AttributeError):
575            x.foobar
576
577    def check_getters(self):
578        T,D,X = self.phase.TDX
579        self.assertNear(T, self.phase.T)
580        self.assertNear(D, self.phase.density)
581        self.assertArrayNear(X, self.phase.X)
582
583        T,D,Y = self.phase.TDY
584        self.assertNear(T, self.phase.T)
585        self.assertNear(D, self.phase.density)
586        self.assertArrayNear(Y, self.phase.Y)
587
588        T,D = self.phase.TD
589        self.assertNear(T, self.phase.T)
590        self.assertNear(D, self.phase.density)
591
592        T,P,X = self.phase.TPX
593        self.assertNear(T, self.phase.T)
594        self.assertNear(P, self.phase.P)
595        self.assertArrayNear(X, self.phase.X)
596
597        T,P,Y = self.phase.TPY
598        self.assertNear(T, self.phase.T)
599        self.assertNear(P, self.phase.P)
600        self.assertArrayNear(Y, self.phase.Y)
601
602        T,P = self.phase.TP
603        self.assertNear(T, self.phase.T)
604        self.assertNear(P, self.phase.P)
605
606        H,P,X = self.phase.HPX
607        self.assertNear(H, self.phase.h)
608        self.assertNear(P, self.phase.P)
609        self.assertArrayNear(X, self.phase.X)
610
611        H,P,Y = self.phase.HPY
612        self.assertNear(H, self.phase.h)
613        self.assertNear(P, self.phase.P)
614        self.assertArrayNear(Y, self.phase.Y)
615
616        H,P = self.phase.HP
617        self.assertNear(H, self.phase.h)
618        self.assertNear(P, self.phase.P)
619
620        U,V,X = self.phase.UVX
621        self.assertNear(U, self.phase.u)
622        self.assertNear(V, self.phase.v)
623        self.assertArrayNear(X, self.phase.X)
624
625        U,V,Y = self.phase.UVY
626        self.assertNear(U, self.phase.u)
627        self.assertNear(V, self.phase.v)
628        self.assertArrayNear(Y, self.phase.Y)
629
630        U,V = self.phase.UV
631        self.assertNear(U, self.phase.u)
632        self.assertNear(V, self.phase.v)
633
634        S,P,X = self.phase.SPX
635        self.assertNear(S, self.phase.s)
636        self.assertNear(P, self.phase.P)
637        self.assertArrayNear(X, self.phase.X)
638
639        S,P,Y = self.phase.SPY
640        self.assertNear(S, self.phase.s)
641        self.assertNear(P, self.phase.P)
642        self.assertArrayNear(Y, self.phase.Y)
643
644        S,P = self.phase.SP
645        self.assertNear(S, self.phase.s)
646        self.assertNear(P, self.phase.P)
647
648        S,V,X = self.phase.SVX
649        self.assertNear(S, self.phase.s)
650        self.assertNear(V, self.phase.v)
651        self.assertArrayNear(X, self.phase.X)
652
653        S,V,Y = self.phase.SVY
654        self.assertNear(S, self.phase.s)
655        self.assertNear(V, self.phase.v)
656        self.assertArrayNear(Y, self.phase.Y)
657
658        S,V = self.phase.SV
659        self.assertNear(S, self.phase.s)
660        self.assertNear(V, self.phase.v)
661
662        D,P,X = self.phase.DPX
663        self.assertNear(D, self.phase.density)
664        self.assertNear(P, self.phase.P)
665        self.assertArrayNear(X, self.phase.X)
666
667        D,P,Y = self.phase.DPY
668        self.assertNear(D, self.phase.density)
669        self.assertNear(P, self.phase.P)
670        self.assertArrayNear(Y, self.phase.Y)
671
672        D,P = self.phase.DP
673        self.assertNear(D, self.phase.density)
674        self.assertNear(P, self.phase.P)
675
676    def test_getState_mass(self):
677        self.phase.TDY = 350.0, 0.7, 'H2:0.1, H2O2:0.1, AR:0.8'
678        self.check_getters()
679
680    def test_getState_mole(self):
681        self.phase.basis = 'molar'
682        self.phase.TDX = 350.0, 0.01, 'H2:0.1, O2:0.3, AR:0.6'
683        self.check_getters()
684
685    def test_getState(self):
686        self.assertNear(self.phase.P, ct.one_atm)
687        self.assertNear(self.phase.T, 300)
688
689    def test_partial_molar(self):
690        self.phase.TDY = 350.0, 0.6, 'H2:0.1, H2O2:0.1, AR:0.8'
691        self.assertNear(sum(self.phase.partial_molar_enthalpies * self.phase.X),
692                        self.phase.enthalpy_mole)
693
694        self.assertNear(sum(self.phase.partial_molar_entropies * self.phase.X),
695                        self.phase.entropy_mole)
696
697        self.assertNear(sum(self.phase.partial_molar_int_energies * self.phase.X),
698                        self.phase.int_energy_mole)
699
700        self.assertNear(sum(self.phase.chemical_potentials * self.phase.X),
701                        self.phase.gibbs_mole)
702
703        self.assertNear(sum(self.phase.partial_molar_cp * self.phase.X),
704                        self.phase.cp_mole)
705
706    def test_nondimensional(self):
707        self.phase.TDY = 850.0, 0.2, 'H2:0.1, H2O:0.6, AR:0.3'
708        H = (sum(self.phase.standard_enthalpies_RT * self.phase.X) *
709             ct.gas_constant * self.phase.T)
710        self.assertNear(H, self.phase.enthalpy_mole)
711
712        U = (sum(self.phase.standard_int_energies_RT * self.phase.X) *
713             ct.gas_constant * self.phase.T)
714        self.assertNear(U, self.phase.int_energy_mole)
715
716        cp = sum(self.phase.standard_cp_R * self.phase.X) * ct.gas_constant
717        self.assertNear(cp, self.phase.cp_mole)
718
719    def test_activities(self):
720        self.phase.TDY = 850.0, 0.2, 'H2:0.1, H2O:0.6, AR:0.3'
721        self.assertArrayNear(self.phase.X, self.phase.activities)
722
723        self.assertArrayNear(self.phase.activity_coefficients,
724                             np.ones(self.phase.n_species))
725
726    def test_isothermal_compressibility(self):
727        self.assertNear(self.phase.isothermal_compressibility, 1.0/self.phase.P)
728
729    def test_thermal_expansion_coeff(self):
730        self.assertNear(self.phase.thermal_expansion_coeff, 1.0/self.phase.T)
731
732    def test_ref_info(self):
733        self.assertNear(self.phase.reference_pressure, ct.one_atm)
734        self.assertNear(self.phase.min_temp, 300.0)
735        self.assertNear(self.phase.max_temp, 3500.0)
736
737    def test_unpicklable(self):
738        import pickle
739        with self.assertRaises(NotImplementedError):
740            pickle.dumps(self.phase)
741
742    def test_uncopyable(self):
743        import copy
744        with self.assertRaises(NotImplementedError):
745            copy.copy(self.phase)
746
747    def test_add_species(self):
748        ref = ct.Solution('gri30.yaml', transport_model=None)
749        n_orig = self.phase.n_species
750        self.phase.add_species(ref.species('CO2'))
751        self.phase.add_species(ref.species('CO'))
752
753        self.assertEqual(self.phase.n_species, n_orig + 2)
754        self.assertIn('CO2', self.phase.species_names)
755        self.assertIn('CO', self.phase.species_names)
756
757        state = 400, 2e5, 'H2:0.7, CO2:0.2, CO:0.1'
758        ref.TPY = state
759        self.phase.TPY = state
760        self.assertNear(self.phase.enthalpy_mass, ref.enthalpy_mass)
761        self.assertNear(self.phase.entropy_mole, ref.entropy_mole)
762        self.assertArrayNear(ref[self.phase.species_names].partial_molar_cp,
763                             self.phase.partial_molar_cp)
764
765    def test_add_species_disabled(self):
766        ref = ct.Solution('gri30.yaml', transport_model=None)
767
768        reactor = ct.IdealGasReactor(self.phase)
769        with self.assertRaisesRegex(ct.CanteraError, 'Cannot add species'):
770            self.phase.add_species(ref.species('CH4'))
771        del reactor
772        gc.collect()
773        self.phase.add_species(ref.species('CH4'))
774
775        flame = ct.FreeFlame(self.phase, width=0.1)
776        with self.assertRaisesRegex(ct.CanteraError, 'Cannot add species'):
777            self.phase.add_species(ref.species('CO'))
778        del flame
779        gc.collect()
780        self.phase.add_species(ref.species('CO'))
781
782        mix = ct.Mixture([(self.phase, 2.0)])
783        with self.assertRaisesRegex(ct.CanteraError, 'Cannot add species'):
784            self.phase.add_species(ref.species('CH2O'))
785        del mix
786        gc.collect()
787        self.phase.add_species(ref.species('CH2O'))
788
789    def test_add_species_duplicate(self):
790        species = self.phase.species('H2O2')
791        with self.assertRaisesRegex(ct.CanteraError, 'already contains'):
792            self.phase.add_species(species)
793
794
795class TestThermo(utilities.CanteraTest):
796    def setUp(self):
797        self.gas = ct.ThermoPhase('h2o2.xml')
798        self.gas.TPX = 450, 2e5, 'H2:1.0, O2:0.4, AR:3, H2O:0.1'
799
800    def test_setSV_lowT(self):
801        """
802        Set state in terms of (s,v) when the end temperature is below the
803        phase's nominal temperature limit.
804        """
805
806        self.gas.TPX = 450, 1e5, 'H2:1.0, O2:0.4, AR:3'
807        s1, v1 = self.gas.SV
808        self.gas.SV = s1, 3 * v1
809
810        self.assertNear(self.gas.s, s1)
811        self.assertNear(self.gas.v, 3 * v1)
812        self.assertTrue(self.gas.T < self.gas.min_temp)
813
814    def test_setSV_low_invalid(self):
815        self.gas.TPX = 450, 1e5, 'H2:1.0, O2:0.4, AR:3'
816        self.gas.SV = 4600, None
817        with self.assertRaises(ct.CanteraError):
818            self.gas.SV = -1000, None
819
820    def test_setSV_highT(self):
821        """
822        Set state in terms of (s,v) when the end temperature is above the
823        phase's nominal temperature limit.
824        """
825
826        self.gas.TPX = 2900, 1e5, 'H2:1.0, O2:0.4, AR:3'
827        s1, v1 = self.gas.SV
828        self.gas.SV = s1, 0.3 * v1
829
830        self.assertNear(self.gas.s, s1)
831        self.assertNear(self.gas.v, 0.3 * v1)
832        self.assertTrue(self.gas.T > self.gas.max_temp)
833
834    def test_setHP_lowT(self):
835        """
836        Set state in terms of (s,v) when the end temperature is below the
837        phase's nominal temperature limit.
838        """
839
840        self.gas.TPX = 450, 1e5, 'H2:1.0, O2:0.4, AR:3'
841        deltaH = 1.25e5
842        h1, p1 = self.gas.HP
843        self.gas.HP = h1 - deltaH, None
844
845        self.assertNear(self.gas.h, h1 - deltaH)
846        self.assertNear(self.gas.P, p1)
847        self.assertTrue(self.gas.T < self.gas.min_temp)
848
849    def test_setHP_low_invalid(self):
850        """
851        Set state in terms of (h,p) when the enthalpy would imply a negative
852        temperature
853        """
854
855        self.gas.TPX = 300, 101325, 'H2:1.0'
856        with self.assertRaises(ct.CanteraError):
857            self.gas.HP = -4e6, 101325
858
859    def test_setHP_highT(self):
860        """
861        Set state in terms of (s,v) when the end temperature is above the
862        phase's nominal temperature limit.
863        """
864
865        self.gas.TPX = 2800, 1e5, 'H2:1.0, O2:0.4, AR:3'
866        deltaH = 8.25e5
867        h1, p1 = self.gas.HP
868        self.gas.HP = h1 + deltaH, None
869
870        self.assertNear(self.gas.h, h1 + deltaH)
871        self.assertNear(self.gas.P, p1)
872        self.assertTrue(self.gas.T > self.gas.max_temp)
873
874    def test_volume(self):
875        """
876        This phase should follow the ideal gas law
877        """
878        g = self.gas
879        self.assertNear(g.P, g.density_mole * ct.gas_constant * g.T)
880
881        self.assertNear(
882            g.P / g.density,
883            ct.gas_constant / g.mean_molecular_weight * g.T)
884
885        self.assertNear(g.density, 1.0 / g.volume_mass)
886
887    def test_energy(self):
888        g = self.gas
889        mmw = g.mean_molecular_weight
890        self.assertNear(g.enthalpy_mass, g.enthalpy_mole / mmw)
891        self.assertNear(g.int_energy_mass, g.int_energy_mole / mmw)
892        self.assertNear(g.gibbs_mass, g.gibbs_mole / mmw)
893        self.assertNear(g.entropy_mass, g.entropy_mole / mmw)
894
895        self.assertNear(g.cv_mass, g.cv_mole / mmw)
896        self.assertNear(g.cp_mass, g.cp_mole / mmw)
897        self.assertNear(g.cv_mole + ct.gas_constant, g.cp_mole)
898
899    def test_nondimensional(self):
900        g = self.gas
901        R = ct.gas_constant
902
903        self.assertNear(np.dot(g.standard_cp_R, g.X), g.cp_mole / R)
904        self.assertNear(np.dot(g.standard_enthalpies_RT, g.X),
905                        g.enthalpy_mole / (R*g.T))
906
907        Smix_R = - np.dot(g.X, np.log(g.X+1e-20))
908        self.assertNear(np.dot(g.standard_entropies_R, g.X) + Smix_R,
909                        g.entropy_mole / R)
910        self.assertNear(np.dot(g.standard_gibbs_RT, g.X) - Smix_R,
911                        g.gibbs_mole / (R*g.T))
912
913
914class TestInterfacePhase(utilities.CanteraTest):
915    def setUp(self):
916        self.gas = ct.Solution('diamond.xml', 'gas')
917        self.solid = ct.Solution('diamond.xml', 'diamond')
918        self.interface = ct.Interface('diamond.xml', 'diamond_100',
919                                      (self.gas, self.solid))
920
921    def test_properties(self):
922        self.interface.site_density = 100
923        self.assertNear(self.interface.site_density, 100)
924
925    def test_coverages_array(self):
926        C = np.zeros(self.interface.n_species)
927        C[1] = 0.25
928        C[3] = 0.125
929        C[4] = 0.125
930        self.interface.coverages = C
931        C = self.interface.coverages
932        # should now be normalized
933        self.assertNear(C[1], 0.5)
934        self.assertNear(C[3], 0.25)
935        self.assertNear(C[4], 0.25)
936        self.assertNear(sum(C), 1.0)
937
938    def test_coverages_string(self):
939        self.interface.coverages = 'c6HM:0.2, c6H*:0.8'
940        C = self.interface.coverages
941        self.assertNear(C[self.interface.species_index('c6HM')], 0.2)
942        self.assertNear(C[self.interface.species_index('c6H*')], 0.8)
943
944    def test_coverages_dict(self):
945        self.interface.coverages = {'c6**':1.0, 'c6*M':3.0}
946        C = self.interface.coverages
947        self.assertNear(C[self.interface.species_index('c6**')], 0.25)
948        self.assertNear(C[self.interface.species_index('c6*M')], 0.75)
949
950
951class ImportTest(utilities.CanteraTest):
952    """
953    Test the various ways of creating a Solution object
954    """
955    def check(self, gas, phase, T, P, nSpec, nElem):
956        self.assertEqual(gas.name, phase)
957        self.assertNear(gas.T, T)
958        self.assertNear(gas.P, P)
959        self.assertEqual(gas.n_species, nSpec)
960        self.assertEqual(gas.n_elements, nElem)
961
962    def test_import_phase_cti(self):
963        gas1 = ct.Solution('air-no-reactions.cti', 'air')
964        self.check(gas1, 'air', 300, 101325, 8, 3)
965
966        gas2 = ct.Solution('air-no-reactions.cti', 'notair')
967        self.check(gas2, 'notair', 900, 5*101325, 7, 2)
968
969    def test_import_phase_cti2(self):
970        # This should import the first phase, i.e. 'air'
971        gas = ct.Solution('air-no-reactions.cti')
972        self.check(gas, 'air', 300, 101325, 8, 3)
973
974    def test_import_phase_xml(self):
975        gas1 = ct.Solution('air-no-reactions.xml', 'air')
976        self.check(gas1, 'air', 300, 101325, 8, 3)
977
978        gas2 = ct.Solution('air-no-reactions.xml', 'notair')
979        self.check(gas2, 'notair', 900, 5*101325, 7, 2)
980
981    def test_import_phase_cti_text(self):
982        cti_def = """
983ideal_gas(name='spam', elements='O H',
984          species='gri30: all',
985          options='skip_undeclared_elements',
986          initial_state=state(temperature=350, pressure=2e6))
987"""
988        gas = ct.Solution(source=cti_def)
989        self.check(gas, 'spam', 350, 2e6, 8, 2)
990
991    def test_import_phase_xml_text(self):
992        xml_def = """
993<?xml version="1.0"?>
994<ctml>
995  <validate reactions="yes" species="yes"/>
996  <phase dim="3" id="spam">
997    <elementArray datasrc="elements.xml">O</elementArray>
998    <speciesArray datasrc="gri30.xml#species_data">all
999      <skip element="undeclared"/>
1000    </speciesArray>
1001    <state>
1002      <temperature units="K">350.0</temperature>
1003      <pressure units="Pa">2000000.0</pressure>
1004    </state>
1005    <thermo model="IdealGas"/>
1006    <kinetics model="GasKinetics"/>
1007    <transport model="None"/>
1008  </phase>
1009</ctml>"""
1010        gas = ct.Solution(source=xml_def)
1011        self.check(gas, 'spam', 350, 2e6, 2, 1)
1012
1013    def test_import_from_species(self):
1014        gas1 = ct.Solution('h2o2.yaml', transport_model=None)
1015        gas1.TPX = 350, 101325, 'H2:0.3, O2:0.7'
1016        gas1.equilibrate('HP')
1017
1018        species = ct.Species.listFromFile('h2o2.yaml')
1019        gas2 = ct.ThermoPhase(thermo='IdealGas', species=species)
1020        gas2.TPX = 350, 101325, 'H2:0.3, O2:0.7'
1021        gas2.equilibrate('HP')
1022        self.assertEqual(gas1.n_elements, gas2.n_elements)
1023        self.assertEqual(gas1.species_names, gas2.species_names)
1024        self.assertNear(gas1.T, gas2.T)
1025        self.assertArrayNear(gas1.X, gas2.X)
1026
1027    def test_checkReactionBalance(self):
1028        with self.assertRaisesRegex(ct.CanteraError, 'reaction is unbalanced'):
1029            ct.Solution('h2o2_unbalancedReaction.xml')
1030
1031    def test_yaml_ideal_gas_simple(self):
1032        gas = ct.ThermoPhase('ideal-gas.yaml', 'simple')
1033        self.check(gas, 'simple', 500, 10 * ct.one_atm, 3, 2)
1034
1035    def test_yaml_ideal_gas_remote_species(self):
1036        gas = ct.ThermoPhase('ideal-gas.yaml', 'species-remote')
1037        self.check(gas, 'species-remote', 300, ct.one_atm, 4, 2)
1038
1039    def test_yaml_duplicate(self):
1040        with self.assertRaisesRegex(ct.CanteraError, 'duplicate'):
1041            gas = ct.ThermoPhase('ideal-gas.yaml', 'duplicate-species')
1042
1043
1044class TestSpecies(utilities.CanteraTest):
1045    def setUp(self):
1046        self.gas = ct.Solution('h2o2.yaml', transport_model=None)
1047
1048    def test_standalone(self):
1049        s = ct.Species('CH4', {'C':1, 'H':4})
1050
1051        self.assertEqual(s.name, 'CH4')
1052        c = s.composition
1053        self.assertEqual(len(c), 2)
1054        self.assertEqual(c['C'], 1)
1055        self.assertEqual(c['H'], 4)
1056
1057    def test_defaults(self):
1058        s = ct.Species('H2')
1059        self.assertEqual(s.size, 1.0)
1060        self.assertEqual(s.charge, 0.0)
1061
1062        self.assertIsNone(s.thermo)
1063        self.assertIsNone(s.transport)
1064
1065    def test_index_accessor(self):
1066        for k in range(self.gas.n_species):
1067            s = self.gas.species(k)
1068            self.assertEqual(s.name, self.gas.species_name(k))
1069
1070            for m,n in s.composition.items():
1071                self.assertEqual(n, self.gas.n_atoms(k,m))
1072
1073    def test_species_noargs(self):
1074        for k,s in enumerate(self.gas.species()):
1075            self.assertEqual(s.name, self.gas.species_name(k))
1076
1077    def test_name_accessor(self):
1078        for name in self.gas.species_names:
1079            s = self.gas.species(name)
1080            self.assertEqual(s.name, name)
1081
1082    def test_fromCti(self):
1083        h2_cti = """
1084            species(
1085                name="H2",
1086                atoms="H:2",
1087                thermo=(
1088                    NASA([200.00, 1000.00],
1089                         [2.344331120E+00, 7.980520750E-03, -1.947815100E-05,
1090                          2.015720940E-08, -7.376117610E-12, -9.179351730E+02,
1091                          6.830102380E-01]),
1092                    NASA([1000.00, 3500.00],
1093                         [3.337279200E+00, -4.940247310E-05, 4.994567780E-07,
1094                         -1.795663940E-10, 2.002553760E-14, -9.501589220E+02,
1095                         -3.205023310E+00])
1096                ),
1097                transport=gas_transport(geom="linear",
1098                                        diam=2.92,
1099                                        well_depth=38.00,
1100                                        polar=0.79,
1101                                        rot_relax=280.00),
1102                note = "TPIS78"
1103            )"""
1104        s1 = self.gas.species('H2')
1105        s2 = ct.Species.fromCti(h2_cti)
1106        self.assertEqual(s2.name, 'H2')
1107        self.assertEqual(s1.composition, s2.composition)
1108        self.assertEqual(s1.thermo.cp(350), s2.thermo.cp(350))
1109
1110    def test_fromXml(self):
1111        import xml.etree.ElementTree as ET
1112        root = ET.parse(self.cantera_data_path / "h2o2.xml").getroot()
1113        h2_node = root.find('.//species[@name="H2"]')
1114        h2_string = ET.tostring(h2_node)
1115
1116        s1 = self.gas.species('H2')
1117        s2 = ct.Species.fromXml(h2_string)
1118
1119        self.assertEqual(s2.name, 'H2')
1120        self.assertEqual(s1.composition, s2.composition)
1121        self.assertEqual(s1.thermo.cp(350), s2.thermo.cp(350))
1122
1123    def test_listFromFile_cti(self):
1124        S = ct.Species.listFromFile('h2o2.cti')
1125        self.assertEqual(S[3].name, self.gas.species_name(3))
1126
1127    def test_listFromFile_xml(self):
1128        S = ct.Species.listFromFile('h2o2.xml')
1129        self.assertEqual(S[3].name, self.gas.species_name(3))
1130
1131    def test_listfromFile_yaml(self):
1132        S = ct.Species.listFromFile('h2o2.yaml')
1133        self.assertEqual({sp.name for sp in S}, set(self.gas.species_names))
1134
1135    def test_listFromCti(self):
1136        S = ct.Species.listFromCti((self.cantera_data_path / "h2o2.cti").read_text())
1137        self.assertEqual(S[3].name, self.gas.species_name(3))
1138
1139    def test_listFomYaml(self):
1140        yaml = '''
1141        - name: H2O
1142          composition: {H: 2, O: 1}
1143          thermo: {model: constant-cp, h0: 100}
1144        - name: HO2
1145          composition: {H: 1, O: 2}
1146          thermo: {model: constant-cp, h0: 200}
1147        '''
1148        species = ct.Species.listFromYaml(yaml)
1149        self.assertEqual(species[0].name, 'H2O')
1150        self.assertEqual(species[1].composition, {'H': 1, 'O': 2})
1151        self.assertNear(species[0].thermo.h(300), 100)
1152
1153    def test_listFromYaml_section(self):
1154        species = ct.Species.listFromYaml(
1155            (self.test_data_path / "ideal-gas.yaml").read_text(),
1156            'species')
1157
1158        self.assertEqual(species[0].name, 'O2')
1159        self.assertEqual(species[1].composition, {'N': 1, 'O': 1})
1160
1161    def test_listFromXml(self):
1162        S = ct.Species.listFromXml((self.cantera_data_path / "h2o2.xml").read_text())
1163        self.assertEqual(S[4].name, self.gas.species_name(4))
1164
1165    def test_modify_thermo(self):
1166        S = {sp.name: sp for sp in ct.Species.listFromFile('h2o2.xml')}
1167        self.gas.TPX = 400, 2*ct.one_atm, 'H2:1.0'
1168        g0 = self.gas.gibbs_mole
1169
1170        self.gas.TPX = None, None, 'O2:1.0'
1171        self.assertNotAlmostEqual(g0, self.gas.gibbs_mole)
1172        # Replace O2 thermo with the data from H2
1173        S['O2'].thermo = S['H2'].thermo
1174        self.gas.modify_species(self.gas.species_index('O2'), S['O2'])
1175        self.assertNear(g0, self.gas.gibbs_mole)
1176
1177    def test_modify_thermo_invalid(self):
1178        S = {sp.name: sp for sp in ct.Species.listFromFile('h2o2.xml')}
1179
1180        orig = S['H2']
1181        thermo = orig.thermo
1182        copy = ct.Species('foobar', orig.composition)
1183        copy.thermo = thermo
1184        with self.assertRaisesRegex(ct.CanteraError, 'modifySpecies'):
1185            self.gas.modify_species(self.gas.species_index('H2'), copy)
1186
1187        copy = ct.Species('H2', {'H': 3})
1188        copy.thermo = thermo
1189        with self.assertRaisesRegex(ct.CanteraError, 'modifySpecies'):
1190            self.gas.modify_species(self.gas.species_index('H2'), copy)
1191
1192        copy = ct.Species('H2', orig.composition)
1193        copy.thermo = ct.ConstantCp(thermo.min_temp, thermo.max_temp,
1194            thermo.reference_pressure, [300, 123, 456, 789])
1195        with self.assertRaisesRegex(ct.CanteraError, 'modifySpecies'):
1196            self.gas.modify_species(self.gas.species_index('H2'), copy)
1197
1198        copy = ct.Species('H2', orig.composition)
1199        copy.thermo = ct.NasaPoly2(thermo.min_temp+200, thermo.max_temp,
1200            thermo.reference_pressure, thermo.coeffs)
1201        with self.assertRaisesRegex(ct.CanteraError, 'modifySpecies'):
1202            self.gas.modify_species(self.gas.species_index('H2'), copy)
1203
1204    def test_alias(self):
1205        self.gas.add_species_alias('H2', 'hydrogen')
1206        self.assertEqual(self.gas.species_index('hydrogen'), 0)
1207        self.gas.X = 'hydrogen:.5, O2:.5'
1208        self.assertNear(self.gas.X[0], 0.5)
1209        with self.assertRaisesRegex(ct.CanteraError, 'Invalid alias'):
1210            self.gas.add_species_alias('H2', 'O2')
1211        with self.assertRaisesRegex(ct.CanteraError, 'Unable to add alias'):
1212            self.gas.add_species_alias('spam', 'eggs')
1213
1214    def test_isomers(self):
1215        gas = ct.Solution('nDodecane_Reitz.yaml')
1216        iso = gas.find_isomers({'C':4, 'H':9, 'O':2})
1217        self.assertEqual(len(iso), 2)
1218        iso = gas.find_isomers('C:4, H:9, O:2')
1219        self.assertEqual(len(iso), 2)
1220        iso = gas.find_isomers({'C':7, 'H':15})
1221        self.assertEqual(len(iso), 1)
1222        iso = gas.find_isomers({'C':7, 'H':16})
1223        self.assertEqual(len(iso), 0)
1224
1225
1226class TestSpeciesThermo(utilities.CanteraTest):
1227    h2o_coeffs = [
1228        1000.0, 3.03399249E+00, 2.17691804E-03, -1.64072518E-07,
1229        -9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00,
1230        4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09,
1231        1.77197817E-12, -3.02937267E+04, -8.49032208E-01
1232    ]
1233    def setUp(self):
1234        self.gas = ct.Solution('h2o2.yaml', transport_model=None)
1235        self.gas.X = 'H2O:1.0'
1236
1237    def test_create(self):
1238        st = ct.NasaPoly2(300, 3500, 101325, self.h2o_coeffs)
1239
1240        for T in [300, 500, 900, 1200, 2000]:
1241            self.gas.TP = T, 101325
1242            self.assertNear(st.cp(T), self.gas.cp_mole)
1243            self.assertNear(st.h(T), self.gas.enthalpy_mole)
1244            self.assertNear(st.s(T), self.gas.entropy_mole)
1245
1246    def test_invalid(self):
1247        with self.assertRaisesRegex(ValueError, 'incorrect length'):
1248            # not enough coefficients
1249            st = ct.NasaPoly2(300, 3500, 101325,
1250                              [1000.0, 3.03399249E+00, 2.17691804E-03])
1251
1252    def test_wrap(self):
1253        st = self.gas.species('H2O').thermo
1254
1255        self.assertIsInstance(st, ct.NasaPoly2)
1256
1257        for T in [300, 500, 900, 1200, 2000]:
1258            self.gas.TP = T, 101325
1259            self.assertNear(st.cp(T), self.gas.cp_mole)
1260            self.assertNear(st.h(T), self.gas.enthalpy_mole)
1261            self.assertNear(st.s(T), self.gas.entropy_mole)
1262
1263    def test_coeffs(self):
1264        st = ct.NasaPoly2(300, 3500, 101325, self.h2o_coeffs)
1265        self.assertEqual(st.min_temp, 300)
1266        self.assertEqual(st.max_temp, 3500)
1267        self.assertEqual(st.reference_pressure, 101325)
1268        self.assertArrayNear(self.h2o_coeffs, st.coeffs)
1269        self.assertEqual(st.n_coeffs, len(st.coeffs))
1270        self.assertTrue(st._check_n_coeffs(st.n_coeffs))
1271
1272    def test_nasa9_load(self):
1273        gas = ct.Solution('airNASA9.cti')
1274        st = gas.species(3).thermo
1275        self.assertIsInstance(st, ct.Nasa9PolyMultiTempRegion)
1276        self.assertEqual(st.n_coeffs, len(st.coeffs))
1277        self.assertTrue(st._check_n_coeffs(st.n_coeffs))
1278
1279    def test_nasa9_create(self):
1280        gas = ct.Solution('airNASA9.cti')
1281        st = gas.species(3).thermo
1282        t_min = st.min_temp
1283        t_max = st.max_temp
1284        p_ref = st.reference_pressure
1285        coeffs = st.coeffs
1286        st2 = ct.Nasa9PolyMultiTempRegion(t_min, t_max, p_ref, coeffs)
1287        self.assertIsInstance(st2, ct.Nasa9PolyMultiTempRegion)
1288        self.assertEqual(st.min_temp, t_min)
1289        self.assertEqual(st.max_temp, t_max)
1290        self.assertEqual(st.reference_pressure, p_ref)
1291        for T in range(300, 20000, 1000):
1292            self.assertNear(st.cp(T), st2.cp(T))
1293            self.assertNear(st.h(T), st2.h(T))
1294            self.assertNear(st.s(T), st2.s(T))
1295
1296    def test_shomate_load(self):
1297        sol = ct.Solution('thermo-models.yaml', 'molten-salt-Margules')
1298        st = sol.species(0).thermo
1299        self.assertIsInstance(st, ct.ShomatePoly2)
1300        self.assertEqual(st.n_coeffs, len(st.coeffs))
1301        self.assertTrue(st._check_n_coeffs(st.n_coeffs))
1302
1303    def test_shomate_create(self):
1304        sol = ct.Solution('thermo-models.yaml', 'molten-salt-Margules')
1305        st = sol.species(0).thermo
1306        t_min = st.min_temp
1307        t_max = st.max_temp
1308        p_ref = st.reference_pressure
1309        coeffs = st.coeffs
1310        st2 = ct.ShomatePoly2(t_min, t_max, p_ref, coeffs)
1311        self.assertIsInstance(st2, ct.ShomatePoly2)
1312        self.assertEqual(st.min_temp, t_min)
1313        self.assertEqual(st.max_temp, t_max)
1314        self.assertEqual(st.reference_pressure, p_ref)
1315        for T in [300, 500, 700, 900]:
1316            self.assertNear(st.cp(T), st2.cp(T))
1317            self.assertNear(st.h(T), st2.h(T))
1318            self.assertNear(st.s(T), st2.s(T))
1319
1320    def test_piecewise_gibbs_load(self):
1321        sol = ct.Solution('thermo-models.yaml', 'HMW-NaCl-electrolyte')
1322        st = sol.species(1).thermo
1323        self.assertIsInstance(st, ct.Mu0Poly)
1324        self.assertEqual(st.n_coeffs, len(st.coeffs))
1325        self.assertTrue(st._check_n_coeffs(st.n_coeffs))
1326
1327    def test_piecewise_gibbs_create1(self):
1328        # use OH- ion data from test/thermo/phaseConstructors.cpp
1329        h298 = -230.015e6
1330        T1 = 298.15
1331        mu1 = -91.50963
1332        T2 = 333.15
1333        mu2 = -85
1334        pref = 101325
1335        coeffs = [2, h298, T1, mu1*ct.gas_constant*T1, T2, mu2*ct.gas_constant*T2]
1336        st2 = ct.Mu0Poly(200, 3500, pref, coeffs)
1337        self.assertIsInstance(st2, ct.Mu0Poly)
1338        self.assertEqual(st2.n_coeffs, len(coeffs))
1339        self.assertEqual(st2.n_coeffs, len(st2.coeffs))
1340
1341    def test_piecewise_gibbs_create2(self):
1342        sol = ct.Solution('thermo-models.yaml', 'HMW-NaCl-electrolyte')
1343        st = sol.species(1).thermo
1344        t_min = st.min_temp
1345        t_max = st.max_temp
1346        p_ref = st.reference_pressure
1347        coeffs = st.coeffs
1348        st2 = ct.Mu0Poly(t_min, t_max, p_ref, coeffs)
1349        self.assertIsInstance(st2, ct.Mu0Poly)
1350        self.assertEqual(st.min_temp, t_min)
1351        self.assertEqual(st.max_temp, t_max)
1352        self.assertEqual(st.reference_pressure, p_ref)
1353        for T in [300, 500, 700, 900]:
1354            self.assertNear(st.cp(T), st2.cp(T))
1355            self.assertNear(st.h(T), st2.h(T))
1356            self.assertNear(st.s(T), st2.s(T))
1357
1358
1359class TestQuantity(utilities.CanteraTest):
1360    @classmethod
1361    def setUpClass(cls):
1362        utilities.CanteraTest.setUpClass()
1363        cls.gas = ct.Solution('gri30.yaml', transport_model=None)
1364
1365    def setUp(self):
1366        self.gas.TPX = 300, 101325, 'O2:1.0, N2:3.76'
1367
1368    def test_mass_moles(self):
1369        q1 = ct.Quantity(self.gas, mass=5)
1370        self.assertNear(q1.mass, 5)
1371        self.assertNear(q1.moles, 5 / q1.mean_molecular_weight)
1372
1373        q1.mass = 7
1374        self.assertNear(q1.moles, 7 / q1.mean_molecular_weight)
1375
1376        q1.moles = 9
1377        self.assertNear(q1.moles, 9)
1378        self.assertNear(q1.mass, 9 * q1.mean_molecular_weight)
1379
1380    def test_extensive(self):
1381        q1 = ct.Quantity(self.gas, mass=5)
1382        self.assertNear(q1.mass, 5)
1383
1384        self.assertNear(q1.volume * q1.density, q1.mass)
1385        self.assertNear(q1.V * q1.density, q1.mass)
1386        self.assertNear(q1.int_energy, q1.moles * q1.int_energy_mole)
1387        self.assertNear(q1.enthalpy, q1.moles * q1.enthalpy_mole)
1388        self.assertNear(q1.entropy, q1.moles * q1.entropy_mole)
1389        self.assertNear(q1.gibbs, q1.moles * q1.gibbs_mole)
1390        self.assertNear(q1.int_energy, q1.U)
1391        self.assertNear(q1.enthalpy, q1.H)
1392        self.assertNear(q1.entropy, q1.S)
1393        self.assertNear(q1.gibbs, q1.G)
1394
1395    def test_multiply(self):
1396        q1 = ct.Quantity(self.gas, mass=5)
1397        q2 = q1 * 2.5
1398        self.assertNear(q1.mass * 2.5, q2.mass)
1399        self.assertNear(q1.moles * 2.5, q2.moles)
1400        self.assertNear(q1.entropy * 2.5, q2.entropy)
1401        self.assertArrayNear(q1.X, q2.X)
1402
1403    def test_multiply_HP(self):
1404        self.gas.TPX = 500, 101325, 'CH4:1.0, O2:0.4'
1405        q1 = ct.Quantity(self.gas, mass=2, constant='HP')
1406        q2 = ct.Quantity(self.gas, mass=1, constant='HP')
1407        q2.equilibrate('HP')
1408        q3 = 0.2 * q1 + q2 * 0.4
1409        self.assertNear(q1.P, q3.P)
1410        self.assertNear(q1.enthalpy_mass, q3.enthalpy_mass)
1411        self.assertNear(q2.enthalpy_mass, q3.enthalpy_mass)
1412
1413    def test_iadd(self):
1414        q0 = ct.Quantity(self.gas, mass=5)
1415        q1 = ct.Quantity(self.gas, mass=5)
1416        q2 = ct.Quantity(self.gas, mass=5)
1417        q2.TPX = 500, 101325, 'CH4:1.0'
1418
1419        q1 += q2
1420        self.assertNear(q0.mass + q2.mass, q1.mass)
1421        # addition is at constant UV
1422        self.assertNear(q0.U + q2.U, q1.U)
1423        self.assertNear(q0.V + q2.V, q1.V)
1424        self.assertArrayNear(q0.X*q0.moles + q2.X*q2.moles, q1.X*q1.moles)
1425
1426    def test_add(self):
1427        q1 = ct.Quantity(self.gas, mass=5)
1428        q2 = ct.Quantity(self.gas, mass=5)
1429        q2.TPX = 500, 101325, 'CH4:1.0'
1430
1431        q3 = q1 + q2
1432        self.assertNear(q1.mass + q2.mass, q3.mass)
1433        # addition is at constant UV
1434        self.assertNear(q1.U + q2.U, q3.U)
1435        self.assertNear(q1.V + q2.V, q3.V)
1436        self.assertArrayNear(q1.X*q1.moles + q2.X*q2.moles, q3.X*q3.moles)
1437
1438    def test_equilibrate(self):
1439        self.gas.TPX = 300, 101325, 'CH4:1.0, O2:0.2, N2:1.0'
1440        q1 = ct.Quantity(self.gas)
1441        self.gas.equilibrate('HP')
1442        T2 = self.gas.T
1443
1444        self.assertNear(q1.T, 300)
1445        q1.equilibrate('HP')
1446        self.assertNear(q1.T, T2)
1447
1448    def test_incompatible(self):
1449        gas2 = ct.Solution('h2o2.yaml', transport_model=None)
1450        q1 = ct.Quantity(self.gas)
1451        q2 = ct.Quantity(gas2)
1452
1453        with self.assertRaisesRegex(ValueError, 'different phase definitions'):
1454            q1+q2
1455
1456
1457class TestMisc(utilities.CanteraTest):
1458    def test_stringify_bad(self):
1459        with self.assertRaises(AttributeError):
1460            ct.Solution(3)
1461
1462    def test_case_sensitive_names(self):
1463        gas = ct.Solution('h2o2.yaml', transport_model=None)
1464        self.assertFalse(gas.case_sensitive_species_names)
1465        self.assertEqual(gas.species_index('h2'), 0)
1466        gas.X = 'h2:.5, o2:.5'
1467        self.assertNear(gas.X[0], 0.5)
1468        gas.Y = 'h2:.5, o2:.5'
1469        self.assertNear(gas.Y[0], 0.5)
1470
1471        gas.case_sensitive_species_names = True
1472        self.assertTrue(gas.case_sensitive_species_names)
1473        with self.assertRaises(ValueError):
1474            gas.species_index('h2')
1475        with self.assertRaisesRegex(ct.CanteraError, 'Unknown species'):
1476            gas.X = 'h2:1.0, o2:1.0'
1477        with self.assertRaisesRegex(ct.CanteraError, 'Unknown species'):
1478            gas.Y = 'h2:1.0, o2:1.0'
1479
1480        gas_cti = """ideal_gas(
1481            name="gas",
1482            elements=" S C Cs ",
1483            species=" nasa: all ",
1484            options=["skip_undeclared_elements"],
1485            initial_state=state(temperature=300, pressure=(1, "bar"))
1486        )"""
1487        ct.suppress_thermo_warnings(True)
1488        gas = ct.Solution(source=gas_cti)
1489        with self.assertRaisesRegex(ct.CanteraError, 'is not unique'):
1490            gas.species_index('cs')
1491        gas.case_sensitive_species_names = True
1492        with self.assertRaises(ValueError):
1493            gas.species_index('cs')
1494
1495
1496class TestElement(utilities.CanteraTest):
1497    @classmethod
1498    def setUpClass(cls):
1499        utilities.CanteraTest.setUpClass()
1500        cls.ar_sym = ct.Element('Ar')
1501        cls.ar_name = ct.Element('argon')
1502        cls.ar_num = ct.Element(18)
1503
1504    def test_element_multiple_possibilities(self):
1505        carbon = ct.Element('Carbon')
1506        self.assertEqual(carbon.name, 'carbon')
1507        self.assertEqual(carbon.symbol, 'C')
1508
1509    def test_element_weight(self):
1510        self.assertNear(self.ar_sym.weight, 39.95)
1511        self.assertNear(self.ar_name.weight, 39.95)
1512        self.assertNear(self.ar_num.weight, 39.95)
1513
1514    def test_element_symbol(self):
1515        self.assertEqual(self.ar_sym.symbol, 'Ar')
1516        self.assertEqual(self.ar_name.symbol, 'Ar')
1517        self.assertEqual(self.ar_num.symbol, 'Ar')
1518
1519    def test_element_name(self):
1520        self.assertEqual(self.ar_sym.name, 'argon')
1521        self.assertEqual(self.ar_name.name, 'argon')
1522        self.assertEqual(self.ar_num.name, 'argon')
1523
1524    def test_element_atomic_number(self):
1525        self.assertEqual(self.ar_sym.atomic_number, 18)
1526        self.assertEqual(self.ar_name.atomic_number, 18)
1527        self.assertEqual(self.ar_num.atomic_number, 18)
1528
1529    def test_element_name_not_present(self):
1530        with self.assertRaisesRegex(ct.CanteraError, 'element not found'):
1531            ct.Element('I am not an element')
1532
1533    def test_element_atomic_number_small(self):
1534        with self.assertRaisesRegex(ct.CanteraError, 'IndexError'):
1535            ct.Element(0)
1536
1537    def test_element_atomic_number_big(self):
1538        num_elements = ct.Element.num_elements_defined
1539        with self.assertRaisesRegex(ct.CanteraError, 'IndexError'):
1540            ct.Element(num_elements + 1)
1541
1542    def test_element_no_weight(self):
1543        with self.assertRaisesRegex(ct.CanteraError, 'no stable isotopes'):
1544            ct.Element('Tc')
1545
1546    def test_element_bad_input(self):
1547        with self.assertRaisesRegex(TypeError, 'input argument to Element'):
1548            ct.Element(1.2345)
1549
1550    def test_get_isotope(self):
1551        d_sym = ct.Element('D')
1552        self.assertEqual(d_sym.atomic_number, 1)
1553        self.assertNear(d_sym.weight, 2.0141017781)
1554        self.assertEqual(d_sym.name, 'deuterium')
1555        self.assertEqual(d_sym.symbol, 'D')
1556
1557        d_name = ct.Element('deuterium')
1558        self.assertEqual(d_name.atomic_number, 1)
1559        self.assertNear(d_name.weight, 2.0141017781)
1560        self.assertEqual(d_name.name, 'deuterium')
1561        self.assertEqual(d_name.symbol, 'D')
1562
1563    def test_elements_lists(self):
1564        syms = ct.Element.element_symbols
1565        names = ct.Element.element_names
1566        num_elements = ct.Element.num_elements_defined
1567        self.assertEqual(len(syms), num_elements)
1568        self.assertEqual(len(names), num_elements)
1569
1570
1571class TestSolutionArray(utilities.CanteraTest):
1572    @classmethod
1573    def setUpClass(cls):
1574        utilities.CanteraTest.setUpClass()
1575        cls.gas = ct.Solution('h2o2.yaml')
1576
1577    def test_passthrough(self):
1578        states = ct.SolutionArray(self.gas, 3)
1579        self.assertEqual(states.n_species, self.gas.n_species)
1580        self.assertEqual(states.reaction_equation(10),
1581                         self.gas.reaction_equation(10))
1582
1583    def test_meta(self):
1584        meta = {'foo': 'bar', 'spam': 'eggs'}
1585        states = ct.SolutionArray(self.gas, 3, meta=meta)
1586        self.assertEqual(states.meta['foo'], 'bar')
1587        self.assertEqual(states.meta['spam'], 'eggs')
1588
1589    def test_get_state(self):
1590        states = ct.SolutionArray(self.gas, 4)
1591        H, P = states.HP
1592        self.assertEqual(H.shape, (4,))
1593        self.assertEqual(P.shape, (4,))
1594
1595        S, P, Y = states.SPY
1596        self.assertEqual(S.shape, (4,))
1597        self.assertEqual(P.shape, (4,))
1598        self.assertEqual(Y.shape, (4, self.gas.n_species))
1599
1600    def test_idealgas_getters(self):
1601        N = 11
1602        states = ct.SolutionArray(self.gas, N)
1603        getters = 'TDPUVHSXY'  # omit getters that contain Q
1604
1605        # obtain setters from thermo objects
1606        all_getters = [k for k in dir(self.gas)
1607                       if not set(k) - set(getters) and len(k)>1]
1608
1609        # ensure that getters do not raise attribute errors
1610        for g in all_getters:
1611            out = getattr(states, g)
1612            self.assertEqual(len(out), len(g))
1613
1614    def test_properties_onedim(self):
1615        N = 11
1616        states = ct.SolutionArray(self.gas, N)
1617        T = np.linspace(300, 2200, N)
1618        P = np.logspace(3, 8, N)
1619        X = 'H2:0.5, O2:0.4, AR:0.1, H2O2:0.01, OH:0.001'
1620        states.TPX = T, P, X
1621
1622        self.assertArrayNear(states.T, T)
1623        self.assertArrayNear(states.P, P)
1624
1625        h = states.enthalpy_mass
1626        ropr = states.reverse_rates_of_progress
1627        Dkm = states.mix_diff_coeffs
1628        for i in range(N):
1629            self.gas.TPX = T[i], P[i], X
1630            self.assertNear(self.gas.enthalpy_mass, h[i])
1631            self.assertArrayNear(self.gas.reverse_rates_of_progress, ropr[i])
1632            self.assertArrayNear(self.gas.mix_diff_coeffs, Dkm[i])
1633
1634    def test_properties_ndim(self):
1635        states = ct.SolutionArray(self.gas, (2,3,5))
1636
1637        T = np.linspace(300, 2200, 5)
1638        P = np.logspace(3, 8, 2)[:,np.newaxis, np.newaxis]
1639        X = np.random.random((3,1,self.gas.n_species))
1640        states.TPX = T, P, X
1641
1642        TT, PP = states.TP
1643
1644        h = states.enthalpy_mass
1645        ropr = states.reverse_rates_of_progress
1646        Dkm = states.mix_diff_coeffs
1647
1648        self.assertEqual(h.shape, (2,3,5))
1649        self.assertEqual(ropr.shape, (2,3,5,self.gas.n_reactions))
1650        self.assertEqual(Dkm.shape, (2,3,5,self.gas.n_species))
1651
1652        for i,j,k in np.ndindex(TT.shape):
1653            self.gas.TPX = T[k], P[i], X[j]
1654            self.assertNear(self.gas.enthalpy_mass, h[i,j,k])
1655            self.assertArrayNear(self.gas.reverse_rates_of_progress, ropr[i,j,k])
1656            self.assertArrayNear(self.gas.mix_diff_coeffs, Dkm[i,j,k])
1657
1658    def test_slicing_onedim(self):
1659        states = ct.SolutionArray(self.gas, 5)
1660        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1661        T0 = states.T
1662        H0 = states.enthalpy_mass
1663
1664        # Verify that original object is updated when slices change
1665        state = states[1]
1666        state.TD = 300, 0.5
1667        self.assertNear(states.T[0], 500)
1668        self.assertNear(states.T[1], 300)
1669        self.assertNear(states.P[2], 2e5)
1670        self.assertNear(states.density[1], 0.5)
1671
1672        # Verify that the slices are updated when the original object changes
1673        states.TD = 900, None
1674        self.assertNear(state.T, 900)
1675        self.assertNear(states.density[1], 0.5)
1676
1677    def test_slicing_ndim(self):
1678        states = ct.SolutionArray(self.gas, (2,5))
1679        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1680        T0 = states.T
1681        H0 = states.enthalpy_mass
1682
1683        # Verify that original object is updated when slices change
1684        row2 = states[1]
1685        row2.TD = 300, 0.5
1686        T = states.T
1687        D = states.density
1688        self.assertArrayNear(T[0], T0[0])
1689        self.assertArrayNear(T[1], 300*np.ones(5))
1690        self.assertArrayNear(D[1], 0.5*np.ones(5))
1691
1692        col3 = states[:,2]
1693        col3.TD = 400, 2.5
1694        T = states.T
1695        D = states.density
1696        self.assertArrayNear(T[:,2], 400*np.ones(2))
1697        self.assertArrayNear(D[:,2], 2.5*np.ones(2))
1698
1699        # Verify that the slices are updated when the original object changes
1700        states.TP = 900, None
1701        self.assertArrayNear(col3.T, 900*np.ones(2))
1702        self.assertArrayNear(row2.T, 900*np.ones(5))
1703
1704    def test_extra_create_by_dict(self):
1705        extra = OrderedDict([('grid', np.arange(10)),
1706                             ('velocity', np.random.rand(10))])
1707        states = ct.SolutionArray(self.gas, 10, extra=extra)
1708        keys = list(states._extra.keys())
1709        self.assertEqual(keys[0], 'grid')
1710        self.assertArrayNear(states.grid, np.arange(10))
1711
1712    def test_extra_no_shape(self):
1713        # The shape of the value for "prop" here is (), which is falsey
1714        # and causes the use of np.full()
1715        states = ct.SolutionArray(self.gas, 3, extra={"prop": 1})
1716        self.assertEqual(states.prop.shape, (3,))
1717        self.assertArrayNear(states.prop, np.array((1, 1, 1)))
1718
1719        # Check a multidimensional SolutionArray
1720        states = ct.SolutionArray(self.gas, (2, 2), extra={"prop": 3})
1721        self.assertEqual(states.prop.shape, (2, 2))
1722        self.assertArrayNear(states.prop, np.array(((3, 3), (3, 3))))
1723
1724    def test_extra_not_empty(self):
1725        """Test that a non-empty SolutionArray raises a ValueError if
1726           initial values for properties are not supplied.
1727        """
1728        with self.assertRaisesRegex(ValueError, "Initial values for extra properties"):
1729            ct.SolutionArray(self.gas, 3, extra=["prop"])
1730        with self.assertRaisesRegex(ValueError, "Initial values for extra properties"):
1731            ct.SolutionArray(self.gas, 3, extra=np.array(["prop", "prop2"]))
1732
1733    def test_extra_create_multidim(self):
1734        # requires matching first dimensions
1735        extra_val = [[1, 2, 3] for i in range(5)]
1736        states = ct.SolutionArray(self.gas, 5, extra={"prop": extra_val})
1737        self.assertEqual(states.prop.shape, (5, 3,))
1738        states = ct.SolutionArray(self.gas, 5, extra={"prop": np.array(extra_val)})
1739        self.assertEqual(states.prop.shape, (5, 3,))
1740        states = ct.SolutionArray(self.gas, (3, 4,), extra={"prop": np.ones((3, 4, 5,))})
1741        self.assertEqual(states.prop.shape, (3, 4, 5,))
1742        states = ct.SolutionArray(self.gas, 1, extra={"prop": extra_val[0]})
1743        self.assertEqual(states.prop.shape, (1, 3,))
1744        states = ct.SolutionArray(self.gas, 1, extra={"prop": [2]})
1745        self.assertEqual(states.prop.shape, (1,))
1746        with self.assertRaisesRegex(ValueError, "Unable to map"):
1747            ct.SolutionArray(self.gas, (3, 3), extra={"prop": np.arange(3)})
1748
1749    def test_extra_create_by_iterable(self):
1750        states = ct.SolutionArray(self.gas, extra=("prop1"))
1751        self.assertEqual(states.prop1.shape, (0,))
1752
1753        # An integer is not an iterable, and only bare strings are
1754        # turned into iterables
1755        with self.assertRaisesRegex(ValueError, "Extra properties"):
1756            ct.SolutionArray(self.gas, extra=2)
1757
1758    def test_extra_not_string(self):
1759        with self.assertRaisesRegex(TypeError, "is not a string"):
1760            ct.SolutionArray(self.gas, extra=[1])
1761
1762    def test_extra_no_objects(self):
1763        with self.assertRaisesRegex(ValueError, "not supported"):
1764            prop = np.array([0, [1, 2], (3, 4)], dtype=object)
1765            states = ct.SolutionArray(self.gas, 3, extra={"prop": prop})
1766
1767    def test_extra_reserved_names(self):
1768        with self.assertRaisesRegex(ValueError, "name is already used"):
1769            ct.SolutionArray(self.gas, extra=["creation_rates"])
1770
1771        with self.assertRaisesRegex(ValueError, "name is already used"):
1772            ct.SolutionArray(self.gas, extra={"creation_rates": 0})
1773
1774    def test_extra_create_by_string(self):
1775        states = ct.SolutionArray(self.gas, extra="prop")
1776        self.assertEqual(states.prop.shape, (0,))
1777
1778    def test_extra_setattr(self):
1779        states = ct.SolutionArray(self.gas, 7, extra={'prop': range(7)})
1780        states.prop = 0
1781        self.assertArrayNear(states.prop, np.zeros((7,)))
1782        mod_array = np.linspace(0, 10, 7)
1783        states.prop = mod_array
1784        self.assertArrayNear(states.prop, mod_array)
1785        with self.assertRaisesRegex(ValueError, "Incompatible shapes"):
1786            states.prop = [1, 2]
1787
1788    def test_assign_to_slice(self):
1789        states = ct.SolutionArray(self.gas, 7, extra={'prop': range(7)})
1790        array = np.arange(7)
1791        self.assertArrayNear(states.prop, array)
1792        states.prop[1] = -5
1793        states.prop[3:5] = [0, 1]
1794        array_mod = np.array([0, -5, 2, 0, 1, 5, 6])
1795        self.assertArrayNear(states.prop, array_mod)
1796        # assign to multi-dimensional extra
1797        extra_val = [[1, 2, 3] for i in range(5)]
1798        states = ct.SolutionArray(self.gas, 5, extra={"prop": extra_val})
1799        states.prop[:, 1] = -1
1800        self.assertArrayNear(states.prop[:, 1], -1 * np.ones((5,)))
1801        states.prop[2, :] = -2
1802        self.assertArrayNear(states.prop[2, :], -2 * np.ones((3,)))
1803
1804    def test_extra_create_by_ndarray(self):
1805        properties_array = np.array(["prop1", "prop2", "prop3"])
1806        states = ct.SolutionArray(self.gas, shape=(0,), extra=properties_array)
1807        self.assertEqual(states.prop1.shape, (0,))
1808        self.assertEqual(states.prop2.shape, (0,))
1809        self.assertEqual(states.prop3.shape, (0,))
1810        # Ensure that a 2-dimensional array is flattened
1811        properties_array = np.array((["prop1"], ["prop2"]))
1812        states = ct.SolutionArray(self.gas, extra=properties_array)
1813        self.assertEqual(states.prop1.shape, (0,))
1814        self.assertEqual(states.prop2.shape, (0,))
1815
1816    def test_append(self):
1817        states = ct.SolutionArray(self.gas, 5)
1818        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1819        self.assertEqual(states.cp_mass.shape, (5,))
1820
1821        states.append(T=1100, P=3e5, X='AR:1.0')
1822        self.assertEqual(states.cp_mass.shape, (6,))
1823        self.assertNear(states.P[-1], 3e5)
1824        self.assertNear(states.T[-1], 1100)
1825
1826        self.gas.TPX = 1200, 5e5, 'O2:0.3, AR:0.7'
1827        states.append(self.gas.state)
1828        self.assertEqual(states.cp_mass.shape, (7,))
1829        self.assertNear(states.P[-1], 5e5)
1830        self.assertNear(states.X[-1, self.gas.species_index('AR')], 0.7)
1831
1832        self.gas.TPX = 300, 1e4, 'O2:0.5, AR:0.5'
1833        HPY = self.gas.HPY
1834        self.gas.TPX = 1200, 5e5, 'O2:0.3, AR:0.7'  # to make sure it gets changed
1835        states.append(HPY=HPY)
1836        self.assertEqual(states.cp_mass.shape, (8,))
1837        self.assertNear(states.P[-1], 1e4)
1838        self.assertNear(states.T[-1], 300)
1839
1840    def test_append_with_extra(self):
1841        states = ct.SolutionArray(self.gas, 5, extra={"prop": "value"})
1842        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1843        self.assertEqual(states._shape, (5,))
1844        states.append(T=1100, P=3e5, X="AR:1.0", prop="value2")
1845        self.assertEqual(states.prop[-1], "value2")
1846        self.assertEqual(states.prop.shape, (6,))
1847        states.append(T=1100, P=3e5, X="AR:1.0", prop=100)
1848        # NumPy converts to the existing type of the array
1849        self.assertEqual(states.prop[-1], "100")
1850        self.assertEqual(states.prop.shape, (7,))
1851        # two-dimensional input array
1852        states = ct.SolutionArray(self.gas, 1, extra={"prop": [1, 2, 3]})
1853        states.append(T=1100, P=3e5, X="AR:1.0", prop=['a', 'b', 'c'])
1854        self.assertEqual(states._shape, (2,))
1855
1856    def test_append_failures(self):
1857        states = ct.SolutionArray(self.gas, 5, extra={"prop": "value"})
1858        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1859        self.assertEqual(states._shape, (5,))
1860
1861        with self.assertRaisesRegex(TypeError, "Missing keyword arguments for extra"):
1862            states.append(T=1100, P=3e5, X="AR:1.0")
1863        # Failing to append a state shouldn't change the size
1864        self.assertEqual(states._shape, (5,))
1865
1866        with self.assertRaisesRegex(KeyError, "does not specify"):
1867            # I is not a valid property
1868            states.append(TPI=(1100, 3e5, "AR:1.0"), prop="value2")
1869        # Failing to append a state shouldn't change the size
1870        self.assertEqual(states._shape, (5,))
1871
1872        with self.assertRaisesRegex(KeyError, "is not a valid"):
1873            # I is not a valid property
1874            states.append(T=1100, P=3e5, I="AR:1.0", prop="value2")
1875        # Failing to append a state shouldn't change the size
1876        self.assertEqual(states._shape, (5,))
1877
1878        with self.assertRaisesRegex(ValueError, "incompatible value"):
1879            # prop has incompatible dimensions
1880            states.append(T=1100, P=3e5, X="AR:1.0", prop=[1, 2, 3])
1881        # Failing to append a state shouldn't change the size
1882        self.assertEqual(states._shape, (5,))
1883
1884        states = ct.SolutionArray(self.gas, 1, extra={"prop": [1, 2, 3]})
1885        with self.assertRaisesRegex(ValueError, "does not match"):
1886            # prop has incorrect dimensions
1887            states.append(T=1100, P=3e5, X="AR:1.0", prop=['a', 'b'])
1888        # Failing to append a state shouldn't change the size
1889        self.assertEqual(states._shape, (1,))
1890
1891    def test_purefluid(self):
1892        water = ct.Water()
1893        states = ct.SolutionArray(water, 5)
1894        states.TQ = 400, np.linspace(0, 1, 5)
1895
1896        P = states.P
1897        for i in range(1, 5):
1898            self.assertNear(P[0], P[i])
1899
1900        states.TP = np.linspace(400, 500, 5), 101325
1901        self.assertArrayNear(states.Q.squeeze(), np.ones(5))
1902
1903    def test_phase_of_matter(self):
1904        water = ct.Water()
1905        states = ct.SolutionArray(water, 5)
1906        T = [300, 500, water.critical_temperature*2, 300]
1907        P = [101325, 101325, 101325, water.critical_pressure*2]
1908        states[:4].TP = T, P
1909        states[4].TQ = 300, .4
1910        pom = ['liquid', 'gas', 'supercritical', 'supercritical', 'liquid-gas-mix']
1911        self.assertEqual(list(states.phase_of_matter), pom)
1912
1913    def test_purefluid_getters(self):
1914        N = 11
1915        water = ct.Water()
1916        states = ct.SolutionArray(water, N)
1917        getters = 'TDPUVHSQ'  # omit getters that contain X or Y
1918
1919        # obtain setters from thermo objects
1920        all_getters = [k for k in dir(water)
1921                       if not set(k) - set(getters) and len(k)>1]
1922
1923        # ensure that getters do not raise attribute errors
1924        for g in all_getters:
1925            out = getattr(states, g)
1926            self.assertEqual(len(out), len(g))
1927
1928    def test_purefluid_yaml_state(self):
1929        yaml = """
1930        phases:
1931        - name: water
1932          thermo: pure-fluid
1933          species: [{{liquidvapor.yaml/species: [H2O]}}]
1934          state: {{T: {T}, P: 101325, Q: {Q}}}
1935          pure-fluid-name: water
1936        """
1937        w = ct.PureFluid(yaml=yaml.format(T=373.177233, Q=0.5))
1938        self.assertNear(w.Q, 0.5)
1939
1940        with self.assertRaisesRegex(ct.CanteraError, "setState"):
1941            ct.PureFluid(yaml=yaml.format(T=373, Q=0.5))
1942
1943        w = ct.PureFluid(yaml=yaml.format(T=370, Q=0.0))
1944        self.assertNear(w.P, 101325)
1945
1946        with self.assertRaisesRegex(ct.CanteraError, "setState"):
1947            ct.PureFluid(yaml=yaml.format(T=370, Q=1.0))
1948
1949    def test_sort(self):
1950        np.random.seed(0)
1951        t = np.random.random(101)
1952        T = np.linspace(300., 1000., 101)
1953        P = ct.one_atm * (1. + 10.*np.random.random(101))
1954
1955        states = ct.SolutionArray(self.gas, 101, extra={'t': t})
1956        states.TP = T, P
1957
1958        states.sort('t')
1959        self.assertTrue((states.t[1:] - states.t[:-1] > 0).all())
1960        self.assertFalse((states.T[1:] - states.T[:-1] > 0).all())
1961        self.assertFalse(np.allclose(states.P, P))
1962
1963        states.sort('T')
1964        self.assertFalse((states.t[1:] - states.t[:-1] > 0).all())
1965        self.assertTrue((states.T[1:] - states.T[:-1] > 0).all())
1966        self.assertArrayNear(states.P, P)
1967
1968        states.sort('T', reverse=True)
1969        self.assertTrue((states.T[1:] - states.T[:-1] < 0).all())
1970
1971    def test_set_equivalence_ratio(self):
1972        states = ct.SolutionArray(self.gas, 8)
1973        phi = np.linspace(.5, 2., 8)
1974        args = 'H2:1.0', 'O2:1.0'
1975        states.set_equivalence_ratio(phi, *args)
1976        states.set_equivalence_ratio(phi[0], *args)
1977        states.set_equivalence_ratio(list(phi), *args)
1978
1979        with self.assertRaises(ValueError):
1980            states.set_equivalence_ratio(phi[:-1], *args)
1981
1982        states = ct.SolutionArray(self.gas, (2,4))
1983        states.set_equivalence_ratio(phi.reshape((2,4)), *args)
1984
1985        with self.assertRaises(ValueError):
1986            states.set_equivalence_ratio(phi, *args)
1987        with self.assertRaises(ValueError):
1988            states.set_equivalence_ratio(phi.reshape((4,2)), *args)
1989
1990    def test_species_slicing(self):
1991        states = ct.SolutionArray(self.gas, (2,5))
1992        states.TPX = np.linspace(500, 1000, 5), 2e5, 'H2:0.5, O2:0.4'
1993        states.equilibrate('HP')
1994        self.assertArrayNear(states('H2').X.squeeze(),
1995                             states.X[...,self.gas.species_index('H2')])
1996
1997        kk = (self.gas.species_index('OH'), self.gas.species_index('O'))
1998        self.assertArrayNear(states('OH','O').partial_molar_cp,
1999                             states.partial_molar_cp[...,kk])
2000
2001    def test_slice_SolutionArray(self):
2002        soln = ct.SolutionArray(self.gas, 10)
2003        arr = soln[2:9:3]
2004        self.assertEqual(len(arr.T), 3)
2005
2006    def test_zero_length_slice_SolutionArray(self):
2007        states = ct.SolutionArray(self.gas, 4)
2008        arr1 = states[3:3]
2009        self.assertEqual(len(arr1.T), 0)
2010        self.assertEqual(arr1.X.shape, (0,10))
2011        self.assertEqual(arr1.n_reactions, 29)
2012
2013        states.TP = [100,300,900,323.23], ct.one_atm
2014        arr2 = states[slice(0)]
2015        self.assertEqual(len(arr2.T), 0)
2016