diff --git a/mirgecom/eos.py b/mirgecom/eos.py index 9ad289767..20ad1cd81 100644 --- a/mirgecom/eos.py +++ b/mirgecom/eos.py @@ -722,11 +722,14 @@ def get_enthalpy(self, temperature: DOFArray, def get_species_molecular_weights(self): """Get the species molecular weights.""" return self._pyrometheus_mech.wts + # return self._pyrometheus_mech.molecular_weights def species_enthalpies(self, cv: ConservedVars, temperature: DOFArray) -> DOFArray: """Get the species specific enthalpies.""" spec_r = self._pyrometheus_mech.gas_constant/self._pyrometheus_mech.wts + # spec_r = (self._pyrometheus_mech.gas_constant + # / self._pyrometheus_mech.molecular_weights) return (spec_r * temperature * self._pyrometheus_mech.get_species_enthalpies_rt(temperature)) diff --git a/test/test_chemistry.py b/test/test_chemistry.py new file mode 100644 index 000000000..5c8862a7f --- /dev/null +++ b/test/test_chemistry.py @@ -0,0 +1,355 @@ +"""Test the chemistry source terms.""" + +__copyright__ = """ +Copyright (C) 2024 University of Illinois Board of Trustees +""" + + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" +import numpy as np +import pyopencl as cl +import pytest +import cantera +from pytools.obj_array import make_obj_array + +from grudge import op + +from meshmode.array_context import PyOpenCLArrayContext +from meshmode.mesh.generation import generate_regular_rect_mesh +from meshmode.array_context import ( # noqa + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests) + +from mirgecom.fluid import make_conserved +from mirgecom.eos import PyrometheusMixture +from mirgecom.discretization import create_discretization_collection +from mirgecom.mechanisms import get_mechanism_input +from mirgecom.thermochemistry import ( + get_pyrometheus_wrapper_class_from_cantera, + get_pyrometheus_wrapper_class +) + + +@pytest.mark.parametrize(("mechname", "fuel", "rate_tol"), + [("uiuc_7sp", "C2H4", 1e-11), + ("sandiego", "H2", 1e-9)]) +@pytest.mark.parametrize("reactor_type", + ["IdealGasReactor", "IdealGasConstPressureReactor"]) +@pytest.mark.parametrize(("pressure", "nsteps"), + [(25000.0, 100), (101325.0, 50)]) +def test_pyrometheus_kinetics(ctx_factory, mechname, fuel, rate_tol, reactor_type, + pressure, nsteps, output_mechanism=True): + """Test known pyrometheus reaction mechanisms. + + This test reproduces a pyrometheus-native test in the MIRGE context. + + Tests that the Pyrometheus mechanism code gets the same reaction rates as + the corresponding mechanism in Cantera. The reactions are integrated in + time and verified against homogeneous reactors in Cantera. + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + dim = 1 + nel_1d = 4 + order = 4 + + mesh = generate_regular_rect_mesh(a=(-0.5,) * dim, b=(0.5,) * dim, + nelements_per_axis=(nel_1d,) * dim) + + dcoll = create_discretization_collection(actx, mesh, order=order) + ones = dcoll.zeros(actx) + 1.0 + + # Pyrometheus initialization + mech_input = get_mechanism_input(mechname) + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + + if output_mechanism: + import pyrometheus + # write then load the mechanism file to yield better readable pytest coverage + with open(f"./{mechname}.py", "w") as mech_file: + code = pyrometheus.codegen.python.gen_thermochem_code(cantera_soln) + print(code, file=mech_file) + + import importlib + pyromechlib = importlib.import_module(f"{mechname}") + pyro_obj = get_pyrometheus_wrapper_class( + pyro_class=pyromechlib.Thermochemistry)(actx.np) + else: + pyro_obj = get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) + + nspecies = pyro_obj.num_species + print(f"PyrometheusMixture::NumSpecies = {nspecies}") + + tempin = 1200.0 + pressin = pressure + print(f"Testing (t,P) = ({tempin}, {pressin})") + + # Homogeneous reactor to get test data + cantera_soln.set_equivalence_ratio(phi=1.0, fuel=fuel+":1", + oxidizer="O2:1.0,N2:3.76") + cantera_soln.TP = tempin, pressin + + # constant density, variable pressure + if reactor_type == "IdealGasReactor": + reactor = cantera.IdealGasReactor(cantera_soln, # pylint: disable=no-member + name="Batch Reactor") + + # constant pressure, variable density + if reactor_type == "IdealGasConstPressureReactor": + reactor = cantera.IdealGasConstPressureReactor( # pylint: disable=no-member + cantera_soln, name="Batch Reactor") + + sim = cantera.ReactorNet([reactor]) # pylint: disable=no-member + + def inf_norm(x): + return actx.to_numpy(op.norm(dcoll, x, np.inf)) + + def get_mixture_entropy_mass(pressure, temperature, mass_fractions): + mmw = pyro_obj.get_mix_molecular_weight(mass_fractions) + + return 1.0/mmw * get_mixture_entropy_mole(pressure, temperature, + mass_fractions) + + def get_mole_average_property(mass_fractions, spec_property): + mmw = pyro_obj.get_mix_molecular_weight(mass_fractions) + mole_fracs = pyro_obj.get_mole_fractions(mmw, mass_fractions) + return sum(mole_fracs[i] * spec_property[i] for i in range(nspecies)) + + # def get_mixture_enthalpy_mole(temperature, mass_fractions): + # h0_rt = pyro_obj.get_species_enthalpies_rt(temperature) + # hmix = get_mole_average_property(mass_fractions, h0_rt) + # return pyro_obj.gas_constant * temperature * hmix + + def get_mixture_entropy_mole(pressure, temperature, mass_fractions): + mmw = pyro_obj.get_mix_molecular_weight(mass_fractions) + # necessary to avoid nans in the log function below + x = actx.np.where( + actx.np.less(pyro_obj.get_mole_fractions(mmw, mass_fractions), 1e-16), + 1e-16, pyro_obj.get_mole_fractions(mmw, mass_fractions)) + s0_r = pyro_obj.get_species_entropies_r(pressure, temperature) + s_t_mix = get_mole_average_property(mass_fractions, s0_r) + s_x_mix = get_mole_average_property(mass_fractions, actx.np.log(x)) + return pyro_obj.gas_constant * (s_t_mix - s_x_mix) + + time = 0.0 + dt = 2.5e-6 + for _ in range(nsteps): + time += dt + sim.advance(time) + + # Get state from Cantera + can_t = reactor.T + can_p = cantera_soln.P + can_rho = reactor.density + can_y = reactor.Y + # print(f"can_y = {can_y}") + + tin = can_t * ones + pin = can_p * ones + rhoin = can_rho * ones + yin = can_y * ones + + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + # assert properties used internaly in the chemistry + + for i in range(nspecies): + # species entropy + can_s = cantera_soln.standard_entropies_R[i] + spec_entropy = pyro_obj.get_species_entropies_r(pin, tin)[i] + abs_error_s = inf_norm(spec_entropy - can_s) + assert abs_error_s < 1.0e-13 + + # species Gibbs energy + can_g = cantera_soln.standard_gibbs_RT[i] + spec_gibbs = pyro_obj.get_species_gibbs_rt(pin, tin)[i] + abs_error_g = inf_norm(spec_gibbs - can_g) + assert abs_error_g < 1.0e-13 + + # mixture entropy mole + can_s_mix_mole = cantera_soln.entropy_mole + s_mix_mole = get_mixture_entropy_mole(pin, tin, yin) + abs_error_s_mix_mole = inf_norm(s_mix_mole - can_s_mix_mole) + assert abs_error_s_mix_mole/can_s_mix_mole < 2.0e-11 + assert abs_error_s_mix_mole < 5.0e-6 + + # mixture entropy mass + can_s_mix_mass = cantera_soln.entropy_mass + s_mix_mass = get_mixture_entropy_mass(pin, tin, yin) + abs_error_s_mix_mass = inf_norm(s_mix_mass - can_s_mix_mass) + assert abs_error_s_mix_mass/can_s_mix_mass < 2.0e-11 + assert abs_error_s_mix_mass < 5.0e-6 + + # delta enthalpy + can_delta_h = cantera_soln.delta_enthalpy/(pyro_obj.gas_constant*tin) + nu = cantera_soln.product_stoich_coeffs - cantera_soln.reactant_stoich_coeffs + delta_h = nu.T@pyro_obj.get_species_enthalpies_rt(tin) + abs_error_delta_h = inf_norm(can_delta_h - delta_h) + assert abs_error_delta_h < 1e-13 + + # delta entropy + # zero or negative mole fractions values are troublesome due to the log + # see CHEMKIN manual for more details + mmw = pyro_obj.get_mix_molecular_weight(yin) + _x = pyro_obj.get_mole_fractions(mmw, yin) + mole_fracs = actx.np.where(actx.np.less(_x, 1e-15), 1e-15, _x) # noqa + delta_s = nu.T@(pyro_obj.get_species_entropies_r(pin, tin) + - actx.np.log(mole_fracs)) + # exclude meaningless check on entropy for irreversible reaction + for i, reaction in enumerate(cantera_soln.reactions()): + if reaction.reversible: + can_delta_s = cantera_soln.delta_entropy[i]/pyro_obj.gas_constant + assert inf_norm(can_delta_s - delta_s[i]) < 1e-13 + + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + pyro_c = pyro_obj.get_concentrations(rhoin, yin) + # print(f"pyro_conc = {pyro_c}") + + # forward rates + kfw_pm = pyro_obj.get_fwd_rate_coefficients(tin, pyro_c) + kfw_ct = cantera_soln.forward_rate_constants + for i, _ in enumerate(cantera_soln.reactions()): + assert inf_norm((kfw_pm[i] - kfw_ct[i]) / kfw_ct[i]) < 1.0e-13 + + # equilibrium rates + keq_pm = actx.np.exp(-1.*pyro_obj.get_equilibrium_constants(pin, tin)) + keq_ct = cantera_soln.equilibrium_constants + for i, reaction in enumerate(cantera_soln.reactions()): + if reaction.reversible: # skip irreversible reactions + assert inf_norm((keq_pm[i] - keq_ct[i]) / keq_ct[i]) < 1.0e-13 + + # reverse rates + krv_pm = pyro_obj.get_rev_rate_coefficients(pin, tin, pyro_c) + krv_ct = cantera_soln.reverse_rate_constants + for i, reaction in enumerate(cantera_soln.reactions()): + if reaction.reversible: # skip irreversible reactions + assert inf_norm((krv_pm[i] - krv_ct[i]) / krv_ct[i]) < 1.0e-13 + + # reaction progress + rates_pm = pyro_obj.get_net_rates_of_progress(pin, tin, pyro_c) + rates_ct = cantera_soln.net_rates_of_progress + for i, rates in enumerate(rates_ct): + assert inf_norm(rates_pm[i] - rates) < rate_tol + + # species production/destruction + omega_pm = pyro_obj.get_net_production_rates(rhoin, tin, yin) + omega_ct = cantera_soln.net_production_rates + for i, omega in enumerate(omega_ct): + assert inf_norm(omega_pm[i] - omega) < rate_tol + + # check that the reactions progress far enough (and remains stable) + assert can_t > 1800.0 + assert can_t < 3200.0 + + +@pytest.mark.parametrize(("mechname", "fuel", "rate_tol"), + [("uiuc_7sp", "C2H4", 1e-11), + ("sandiego", "H2", 1e-9)]) +@pytest.mark.parametrize("reactor_type", + ["IdealGasReactor", "IdealGasConstPressureReactor"]) +@pytest.mark.parametrize(("pressure", "nsteps"), + [(25000.0, 100), (101325.0, 50)]) +def test_mirgecom_kinetics(ctx_factory, mechname, fuel, rate_tol, reactor_type, + pressure, nsteps): + """Test of known pyrometheus reaction mechanisms in the MIRGE context. + + Tests that the Pyrometheus mechanism code gets the same reaction rates as + the corresponding mechanism in Cantera. The reactions are integrated in + time and verified against a homogeneous reactor in Cantera. + """ + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + dim = 1 + nel_1d = 1 + order = 4 + + mesh = generate_regular_rect_mesh(a=(-0.5,) * dim, b=(0.5,) * dim, + nelements_per_axis=(nel_1d,) * dim) + + dcoll = create_discretization_collection(actx, mesh, order=order) + zeros = dcoll.zeros(actx) + ones = dcoll.zeros(actx) + 1.0 + + def inf_norm(x): + return actx.to_numpy(op.norm(dcoll, x, np.inf)) + + # Pyrometheus initialization + mech_input = get_mechanism_input(mechname) + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + + pyro_obj = get_pyrometheus_wrapper_class_from_cantera( + cantera_soln, temperature_niter=5)(actx.np) + + nspecies = pyro_obj.num_species + print(f"PyrometheusMixture::NumSpecies = {nspecies}") + + tempin = 1200.0 + print(f"Testing (t,P) = ({tempin}, {pressure})") + + # Homogeneous reactor to get test data + cantera_soln.set_equivalence_ratio(phi=1.0, fuel=fuel+":1", + oxidizer="O2:1.0,N2:3.76") + cantera_soln.TP = tempin, pressure + + eos = PyrometheusMixture(pyro_obj, temperature_guess=tempin) + + # constant density, variable pressure + if reactor_type == "IdealGasReactor": + reactor = cantera.IdealGasReactor( # pylint: disable=no-member + cantera_soln, name="Batch Reactor") + + # constant pressure, variable density + if reactor_type == "IdealGasConstPressureReactor": + reactor = cantera.IdealGasConstPressureReactor( # pylint: disable=no-member + cantera_soln, name="Batch Reactor") + + net = cantera.ReactorNet([reactor]) # pylint: disable=no-member + + time = 0.0 + dt = 2.5e-6 + for _ in range(nsteps): + time += dt + net.advance(time) + + tin = reactor.T * ones + rhoin = reactor.density * ones + yin = reactor.Y * ones + ein = rhoin * eos.get_internal_energy(temperature=tin, + species_mass_fractions=yin) + + cv = make_conserved(dim=dim, mass=rhoin, energy=ein, + momentum=make_obj_array([zeros]), species_mass=rhoin*yin) + + temp = eos.temperature(cv=cv, temperature_seed=tin) + + # do NOT test anything else. If this match, then everything is ok + omega_mc = eos.get_production_rates(cv, temp) + omega_ct = cantera_soln.net_production_rates + for i in range(cantera_soln.n_species): + assert inf_norm((omega_mc[i] - omega_ct[i])) < rate_tol + + # check that the reactions progress far enough and are stable + assert reactor.T > 1800.0 + assert reactor.T < 3200.0 diff --git a/test/test_eos.py b/test/test_eos.py index bbb9ac51c..8b627e662 100644 --- a/test/test_eos.py +++ b/test/test_eos.py @@ -42,9 +42,6 @@ pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests) -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests, -) from mirgecom.fluid import make_conserved from mirgecom.eos import IdealSingleGas, PyrometheusMixture from mirgecom.gas_model import GasModel, make_fluid_state @@ -59,7 +56,8 @@ @pytest.mark.parametrize("mechname", ["air_3sp", "uiuc_7sp", "sandiego", "uiuc_8sp_phenol", "uiuc_4sp_oxidation"]) @pytest.mark.parametrize("dim", [1, 2, 3]) -def test_mixture_dependent_properties(ctx_factory, mechname, dim): +@pytest.mark.parametrize("pressure", [25000.0, 101325.0]) +def test_mixture_dependent_properties(ctx_factory, mechname, dim, pressure): """Test MixtureEOS functionality.""" cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -91,26 +89,26 @@ def test_mixture_dependent_properties(ctx_factory, mechname, dim): def inf_norm(x): return actx.to_numpy(op.norm(dcoll, x, np.inf)) - # first check each species individually for a fixed temperature + # ~~~ First check each species individually for a fixed temperature tempin = 600.0 - pressin = 101325.0 eos = PyrometheusMixture(pyro_obj, temperature_guess=tempin) gas_model = GasModel(eos=eos) for i in range(nspecies): x = np.zeros(nspecies,) x[i] = 1.0 - cantera_soln.TPX = tempin, pressin, x + cantera_soln.TPX = tempin, pressure, x can_t, can_rho, can_y = cantera_soln.TDY - can_p = pressin + can_p = cantera_soln.P + pin = can_p * ones tin = can_t * ones yin = can_y * ones rhoin = can_rho * ones # First, check density - mass = eos.get_density(pressin*ones, tempin*ones, yin) + mass = eos.get_density(pin, tin, yin) abs_err_m = inf_norm(mass - can_rho) assert abs_err_m/can_rho < 1.0e-14 assert abs_err_m < 1.0e-10 @@ -119,7 +117,7 @@ def inf_norm(x): cv = make_conserved(dim=2, mass=rhoin, momentum=make_obj_array([zeros, zeros]), energy=rhoin*gas_model.eos.get_internal_energy(tin, yin), - species_mass=can_rho * can_y * ones) + species_mass=rhoin*yin) fluid_state = make_fluid_state(cv, gas_model, tin) @@ -158,8 +156,9 @@ def inf_norm(x): assert abs_err_e/np.abs(can_e) < 1.0e-12 assert abs_err_e < 1.0e-6 + # enthalpy can_h = cantera_soln.enthalpy_mass - enthalpy = fluid_state.dv.species_enthalpies[i] + enthalpy = fluid_state.dv.species_enthalpies[i] # pylint: disable=no-member abs_err_h = inf_norm(enthalpy - can_h) assert abs_err_h/np.abs(can_h) < 1.0e-12 assert abs_err_h < 1.0e-6 @@ -170,20 +169,18 @@ def inf_norm(x): for tempin in ([300.0, 600.0, 900.0, 1200.0]): - print(f"Testing (t,P) = ({tempin}, {pressin})") - - eos = PyrometheusMixture(pyro_obj, temperature_guess=tempin) - gas_model = GasModel(eos=eos) + print(f"Testing (t,P) = ({tempin}, {pressure})") - cantera_soln.TPX = tempin, pressin, x + cantera_soln.TPX = tempin, pressure, x can_t, can_rho, can_y = cantera_soln.TDY - can_p = pressin + can_p = cantera_soln.P + pin = can_p * ones tin = can_t * ones - rhoin = can_rho * ones yin = can_y * ones + rhoin = can_rho * ones - mass = eos.get_density(pressin*ones, tempin*ones, yin) + mass = eos.get_density(pin, tin, yin) abs_err_m = inf_norm(mass - can_rho) assert abs_err_m/can_rho < 1.0e-14 assert abs_err_m < 1.0e-10 @@ -191,7 +188,7 @@ def inf_norm(x): cv = make_conserved(dim=2, mass=rhoin, momentum=make_obj_array([zeros, zeros]), energy=rhoin*gas_model.eos.get_internal_energy(tin, yin), - species_mass=can_rho * can_y * ones) + species_mass=rhoin*yin) fluid_state = make_fluid_state(cv, gas_model, tin) @@ -236,7 +233,7 @@ def inf_norm(x): @pytest.mark.parametrize("mechname", ["air_3sp", "uiuc_7sp", "sandiego", "uiuc_8sp_phenol", "uiuc_4sp_oxidation"]) -def test_pyrometheus_mechanisms(ctx_factory, mechname): +def test_pyrometheus_mechanisms(ctx_factory, mechname, output_mechanism=True): """Test known pyrometheus mechanisms. This test reproduces a pyrometheus-native test in the MIRGE context and @@ -264,8 +261,23 @@ def inf_norm(x): # Pyrometheus initialization mech_input = get_mechanism_input(mechname) - sol = cantera.Solution(name="gas", yaml=mech_input) - pyro_mechanism = get_pyrometheus_wrapper_class_from_cantera(sol)(actx.np) + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + + if output_mechanism: + import pyrometheus + # write then load the mechanism file to yield better readable pytest coverage + with open(f"./{mechname}.py", "w") as mech_file: + code = pyrometheus.codegen.python.gen_thermochem_code(cantera_soln) + print(code, file=mech_file) + + import importlib + from mirgecom.thermochemistry import get_pyrometheus_wrapper_class + pyromechlib = importlib.import_module(f"{mechname}") + pyro_mechanism = get_pyrometheus_wrapper_class( + pyro_class=pyromechlib.Thermochemistry)(actx.np) + else: + pyro_mechanism = \ + get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) nspecies = pyro_mechanism.num_species print(f"PyrometheusMixture::NumSpecies = {nspecies}") @@ -274,12 +286,11 @@ def inf_norm(x): temp0 = 300.0 y0s = np.ones(shape=(nspecies,))/float(nspecies) - for fac in range(1, 11): + for fac in [0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]: pressin = fac * press0 tempin = fac * temp0 print(f"Testing (t,P) = ({tempin}, {pressin})") - cantera_soln = cantera.Solution(name="gas", yaml=mech_input) cantera_soln.TPY = tempin, pressin, y0s can_t, can_m, can_y = cantera_soln.TDY @@ -328,7 +339,6 @@ def inf_norm(x): assert inf_norm(conc[spec]) < 1e-14 -# TODO remove this test.. It is already covered in the other ones @pytest.mark.parametrize("mechname", ["uiuc_7sp", "sandiego"]) @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("y0", [0, 1]) @@ -356,14 +366,23 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): dcoll = create_discretization_collection(actx, mesh, order=order) nodes = actx.thaw(dcoll.nodes()) + ones = dcoll.zeros(actx) + 1.0 + + def inf_norm(x): + return actx.to_numpy(op.norm(dcoll, x, np.inf)) + # Pyrometheus initialization mech_input = get_mechanism_input(mechname) - sol = cantera.Solution(name="gas", yaml=mech_input) - prometheus_mechanism = get_pyrometheus_wrapper_class_from_cantera(sol)(actx.np) + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + pyro_mechanism = get_pyrometheus_wrapper_class_from_cantera( + cantera_soln, temperature_niter=5)(actx.np) + + nspecies = pyro_mechanism.num_species + print(f"PyrometheusMixture::Mechanism = {mechname}") + print(f"PyrometheusMixture::NumSpecies = {nspecies}") - nspecies = prometheus_mechanism.num_species - print(f"PrometheusMixture::Mechanism = {mechname}") - print(f"PrometheusMixture::NumSpecies = {nspecies}") + eos = PyrometheusMixture(pyro_mechanism) + gas_model = GasModel(eos=eos) press0 = 101500.0 temp0 = 300.0 @@ -373,28 +392,28 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): y0s[0] = 1.0 - np.sum(y0s[1:]) velocity = vel * np.ones(shape=(dim,)) - for fac in range(1, 7): + for fac in [0.125, 0.25, 0.5, 1.0, 2.0, 4.0, 8.0]: tempin = fac * temp0 pressin = fac * press0 print(f"Testing {mechname}(t,P) = ({tempin}, {pressin})") - ones = dcoll.zeros(actx) + 1.0 tin = tempin * ones pin = pressin * ones yin = y0s * ones - tguess = 300.0 - pyro_rho = prometheus_mechanism.get_density(pin, tin, yin) - pyro_e = prometheus_mechanism.get_mixture_internal_energy_mass(tin, yin) - pyro_t = prometheus_mechanism.get_temperature(pyro_e, tguess, yin) - pyro_p = prometheus_mechanism.get_pressure(pyro_rho, pyro_t, yin) + # avoid starting too far from the actual temperature + tguess = tempin + ones*100.0*np.random.random() + assert inf_norm(tguess) > 0.0 + + pyro_rho = pyro_mechanism.get_density(pin, tin, yin) + pyro_e = pyro_mechanism.get_mixture_internal_energy_mass(tin, yin) + pyro_t = pyro_mechanism.get_temperature(pyro_e, tguess, yin) + pyro_p = pyro_mechanism.get_pressure(pyro_rho, pyro_t, yin) - print(f"prom(rho, y, p, t, e) = ({pyro_rho}, {y0s}, " - f"{pyro_p}, {pyro_t}, {pyro_e})") + # print(f"prom(rho, y, p, t, e) = ({pyro_rho}, {y0s}, " + # f"{pyro_p}, {pyro_t}, {pyro_e})") - eos = PyrometheusMixture(prometheus_mechanism) - gas_model = GasModel(eos=eos) initializer = Uniform( dim=dim, pressure=pyro_p, temperature=pyro_t, species_mass_fractions=y0s, velocity=velocity) @@ -408,15 +427,12 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): y = cv.species_mass_fractions rho = cv.mass - print(f"pyro_y = {y}") - print(f"pyro_eos.p = {p}") - print(f"pyro_eos.temp = {temperature}") - print(f"pyro_eos.e = {internal_energy}") - - def inf_norm(x): - return actx.to_numpy(op.norm(dcoll, x, np.inf)) + # print(f"pyro_y = {y}") + # print(f"pyro_eos.p = {p}") + # print(f"pyro_eos.temp = {temperature}") + # print(f"pyro_eos.e = {internal_energy}") - tol = 1e-14 + tol = 5e-14 assert inf_norm((cv.mass - pyro_rho) / pyro_rho) < tol assert inf_norm((temperature - pyro_t) / pyro_t) < tol assert inf_norm((internal_energy - pyro_e) / pyro_e) < tol @@ -424,245 +440,28 @@ def inf_norm(x): # Test the concentrations zero level y = -1.0*y - print(f"{y=}") - conc = prometheus_mechanism.get_concentrations(rho, y) - print(f"{conc=}") + # print(f"{y=}") + conc = pyro_mechanism.get_concentrations(rho, y) + # print(f"{conc=}") for spec in range(nspecies): assert max(conc[spec]).all() >= 0 zlev = 1e-3 test_mech = \ - get_pyrometheus_wrapper_class_from_cantera(sol, + get_pyrometheus_wrapper_class_from_cantera(cantera_soln, zero_level=zlev)(actx.np) y = 0*y + zlev - print(f"{y=}") + # print(f"{y=}") conc = test_mech.get_concentrations(rho, y) - print(f"{conc=}") + # print(f"{conc=}") for spec in range(nspecies): assert max(conc[spec]).all() == 0 -@pytest.mark.parametrize(("mechname", "fuel", "rate_tol"), - [("uiuc_7sp", "C2H4", 1e-11), - ("sandiego", "H2", 1e-9)]) -@pytest.mark.parametrize("reactor_type", - ["IdealGasReactor", "IdealGasConstPressureReactor"]) -def test_pyrometheus_kinetics(ctx_factory, mechname, fuel, rate_tol, reactor_type): - """Test known pyrometheus reaction mechanisms. - - This test reproduces a pyrometheus-native test in the MIRGE context. - - Tests that the Pyrometheus mechanism code gets the same reaction rates as - the corresponding mechanism in Cantera. The reactions are integrated in - time and verified against homogeneous reactors in Cantera. - """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) - - dim = 1 - nel_1d = 4 - - mesh = generate_regular_rect_mesh( - a=(-0.5,) * dim, b=(0.5,) * dim, nelements_per_axis=(nel_1d,) * dim - ) - - order = 4 - - logger.info(f"Number of elements {mesh.nelements}") - - dcoll = create_discretization_collection(actx, mesh, order=order) - ones = dcoll.zeros(actx) + 1.0 - - # Pyrometheus initialization - mech_input = get_mechanism_input(mechname) - cantera_soln = cantera.Solution(name="gas", yaml=mech_input) - - pyro_obj = get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) - - nspecies = pyro_obj.num_species - print(f"PyrometheusMixture::NumSpecies = {nspecies}") - - tempin = 1200.0 - pressin = cantera.one_atm - print(f"Testing (t,P) = ({tempin}, {pressin})") - - # Homogeneous reactor to get test data - cantera_soln.set_equivalence_ratio(phi=1.0, fuel=fuel+":1", - oxidizer="O2:1.0,N2:3.76") - cantera_soln.TP = tempin, pressin - - # constant density, variable pressure - if reactor_type == "IdealGasReactor": - reactor = cantera.IdealGasReactor(cantera_soln, name="Batch Reactor") - - # constant pressure, variable density - if reactor_type == "IdealGasConstPressureReactor": - reactor = cantera.IdealGasConstPressureReactor(cantera_soln, - name="Batch Reactor") - - sim = cantera.ReactorNet([reactor]) - - time = 0.0 - dt = 2e-6 - for _ in range(50): - time += dt - sim.advance(time) - - # Get state from Cantera - can_t = reactor.T - can_p = cantera_soln.P - can_rho = reactor.density - can_y = reactor.Y - print(f"can_y = {can_y}") - - tin = can_t * ones - pin = can_p * ones - rhoin = can_rho * ones - yin = can_y * ones - - pyro_c = pyro_obj.get_concentrations(rhoin, yin) - print(f"pyro_conc = {pyro_c}") - - # Print - def inf_norm(x): - return actx.to_numpy(op.norm(dcoll, x, np.inf)) - - # forward rates - kfw_pm = pyro_obj.get_fwd_rate_coefficients(tin, pyro_c) - kfw_ct = cantera_soln.forward_rate_constants - for i, _ in enumerate(cantera_soln.reactions()): - assert inf_norm((kfw_pm[i] - kfw_ct[i]) / kfw_ct[i]) < 1.0e-13 - - # equilibrium rates - keq_pm = actx.np.exp(-1.*pyro_obj.get_equilibrium_constants(pin, tin)) - keq_ct = cantera_soln.equilibrium_constants - for i, reaction in enumerate(cantera_soln.reactions()): - if reaction.reversible: # skip irreversible reactions - assert inf_norm((keq_pm[i] - keq_ct[i]) / keq_ct[i]) < 1.0e-13 - - # reverse rates - krv_pm = pyro_obj.get_rev_rate_coefficients(pin, tin, pyro_c) - krv_ct = cantera_soln.reverse_rate_constants - for i, reaction in enumerate(cantera_soln.reactions()): - if reaction.reversible: # skip irreversible reactions - assert inf_norm((krv_pm[i] - krv_ct[i]) / krv_ct[i]) < 1.0e-13 - - # reaction progress - rates_pm = pyro_obj.get_net_rates_of_progress(pin, tin, pyro_c) - rates_ct = cantera_soln.net_rates_of_progress - for i, rates in enumerate(rates_ct): - assert inf_norm(rates_pm[i] - rates) < rate_tol - - # species production/destruction - omega_pm = pyro_obj.get_net_production_rates(rhoin, tin, yin) - omega_ct = cantera_soln.net_production_rates - for i, omega in enumerate(omega_ct): - assert inf_norm(omega_pm[i] - omega) < rate_tol - - # check that the reactions progress far enough - assert can_t > 2000.0 - assert can_t < 4000.0 - - -@pytest.mark.parametrize(("mechname", "fuel", "rate_tol"), - [("uiuc_7sp", "C2H4", 1e-11), - ("sandiego", "H2", 1e-9)]) -@pytest.mark.parametrize("reactor_type", - ["IdealGasReactor", "IdealGasConstPressureReactor"]) -def test_mirgecom_kinetics(ctx_factory, mechname, fuel, rate_tol, reactor_type): - """Test of known pyrometheus reaction mechanisms in the MIRGE context. - - Tests that the Pyrometheus mechanism code gets the same reaction rates as - the corresponding mechanism in Cantera. The reactions are integrated in - time and verified against a homogeneous reactor in Cantera. - """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) - - dim = 1 - nel_1d = 1 - - mesh = generate_regular_rect_mesh( - a=(-0.5,) * dim, b=(0.5,) * dim, nelements_per_axis=(nel_1d,) * dim - ) - - order = 4 - - logger.info(f"Number of elements {mesh.nelements}") - - dcoll = create_discretization_collection(actx, mesh, order=order) - zeros = dcoll.zeros(actx) - ones = dcoll.zeros(actx) + 1.0 - - def inf_norm(x): - return actx.to_numpy(op.norm(dcoll, x, np.inf)) - - # Pyrometheus initialization - mech_input = get_mechanism_input(mechname) - cantera_soln = cantera.Solution(name="gas", yaml=mech_input) - - pyro_obj = get_pyrometheus_wrapper_class_from_cantera( - cantera_soln, temperature_niter=5)(actx.np) - - nspecies = pyro_obj.num_species - print(f"PyrometheusMixture::NumSpecies = {nspecies}") - - tempin = 1200.0 - pressin = cantera.one_atm - print(f"Testing (t,P) = ({tempin}, {pressin})") - - # Homogeneous reactor to get test data - cantera_soln.set_equivalence_ratio(phi=1.0, fuel=fuel+":1", - oxidizer="O2:1.0,N2:3.76") - cantera_soln.TP = tempin, pressin - - eos = PyrometheusMixture(pyro_obj, temperature_guess=tempin) - - # constant density, variable pressure - if reactor_type == "IdealGasReactor": - reactor = cantera.IdealGasReactor(cantera_soln, name="Batch Reactor") - - # constant pressure, variable density - if reactor_type == "IdealGasConstPressureReactor": - reactor = cantera.IdealGasConstPressureReactor(cantera_soln, - name="Batch Reactor") - - net = cantera.ReactorNet([reactor]) - - time = 0.0 - dt = 2e-6 - for _ in range(50): - time += dt - net.advance(time) - - can_t = reactor.T - tin = can_t * ones - rhoin = reactor.density * ones - yin = reactor.Y * ones - ein = rhoin * eos.get_internal_energy(temperature=tin, - species_mass_fractions=yin) - - cv = make_conserved(dim=dim, mass=rhoin, energy=ein, - momentum=make_obj_array([zeros]), species_mass=rhoin*yin) - - temp = eos.temperature(cv=cv, temperature_seed=tin) - - omega_mc = eos.get_production_rates(cv, temp) - omega_ct = cantera_soln.net_production_rates - for i in range(cantera_soln.n_species): - assert inf_norm((omega_mc[i] - omega_ct[i])) < rate_tol - - # check that the reactions progress far enough - assert can_t > 2000.0 - assert can_t < 4000.0 - - @pytest.mark.parametrize("mechname", ["uiuc_7sp_const_gamma"]) def test_temperature_constant_cv(ctx_factory, mechname): - """.""" + """Test mechanism with constant heat capacity.""" cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) actx = PyOpenCLArrayContext(queue) @@ -694,8 +493,8 @@ def inf_norm(x): pyro_obj = get_pyrometheus_wrapper_class_from_cantera( cantera_soln, temperature_niter=0)(actx.np) - pressin = cantera.one_atm - eos = PyrometheusMixture(pyro_obj, temperature_guess=0.) # XXX dummy + pressin = cantera.one_atm # pylint: disable=no-member + eos = PyrometheusMixture(pyro_obj, temperature_guess=666.) for tin in ([300.0, 600.0, 900.0, 1200.0, 1500.0, 1800.0, 2100.0]): cantera_soln.TP = tin, pressin diff --git a/test/test_transport.py b/test/test_transport.py index e00347305..186d204b3 100644 --- a/test/test_transport.py +++ b/test/test_transport.py @@ -51,15 +51,21 @@ from mirgecom.gas_model import GasModel, make_fluid_state from mirgecom.discretization import create_discretization_collection from mirgecom.mechanisms import get_mechanism_input -from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera +from mirgecom.thermochemistry import ( + get_pyrometheus_wrapper_class_from_cantera, + get_pyrometheus_wrapper_class +) +import pyrometheus logger = logging.getLogger(__name__) -@pytest.mark.parametrize("mechname", ["sandiego"]) +@pytest.mark.parametrize("mechname", ["uiuc_7sp"]) @pytest.mark.parametrize("dim", [2]) @pytest.mark.parametrize("order", [1, 3, 5]) -def test_pyrometheus_transport(ctx_factory, mechname, dim, order): +@pytest.mark.parametrize("use_lewis", [True, False]) +def test_pyrometheus_transport(ctx_factory, mechname, dim, order, use_lewis, + output_mechanism=True): """Test mixture-averaged transport properties.""" cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) @@ -82,25 +88,41 @@ def inf_norm(x): # Pyrometheus initialization mech_input = get_mechanism_input(mechname) - cantera_soln = cantera.Solution(name="gas", yaml=mech_input) - pyro_obj = get_pyrometheus_wrapper_class_from_cantera( - cantera_soln, temperature_niter=3)(actx.np) + ct_transport_model = "unity-Lewis-number" if use_lewis else "mixture-averaged" + cantera_soln = cantera.Solution(name="gas", yaml=mech_input, + transport_model=ct_transport_model) + + if output_mechanism: + # write then load the mechanism file to yield better readable pytest coverage + with open(f"./{mechname}.py", "w") as mech_file: + code = pyrometheus.codegen.python.gen_thermochem_code(cantera_soln) + print(code, file=mech_file) + + import importlib + pyromechlib = importlib.import_module(f"{mechname}") + pyro_obj = get_pyrometheus_wrapper_class( + pyro_class=pyromechlib.Thermochemistry)(actx.np) + else: + pyro_obj = get_pyrometheus_wrapper_class_from_cantera( + cantera_soln, temperature_niter=3)(actx.np) nspecies = pyro_obj.num_species print(f"PyrometheusMixture::NumSpecies = {nspecies}") # Transport data initialization sing_diff = 1e-6 + lewis = np.ones(nspecies,) if use_lewis else None transport_model = MixtureAveragedTransport(pyro_obj, epsilon=1e-4, - singular_diffusivity=sing_diff) + singular_diffusivity=sing_diff, + lewis=lewis) eos = PyrometheusMixture(pyro_obj, temperature_guess=666.) gas_model = GasModel(eos=eos, transport=transport_model) for pressin in ([0.25, 1.0]): for tempin in ([300.0, 600.0, 900.0, 1200.0, 1500.0, 1800.0, 2100.0]): - cantera_soln.TP = tempin, pressin*cantera.one_atm + cantera_soln.TP = tempin, pressin*101325.0 print(f"Testing (T, P) = ({cantera_soln.T}, {cantera_soln.P})") # Loop over each individual species by making a single-species mixture @@ -113,7 +135,7 @@ def inf_norm(x): rhoin = can_rho * ones yin = can_y * ones - cv = make_conserved(dim=2, mass=rhoin, + cv = make_conserved(dim=dim, mass=rhoin, momentum=make_obj_array([zeros, zeros]), energy=rhoin*gas_model.eos.get_internal_energy(tin, yin), species_mass=rhoin*yin) @@ -133,13 +155,14 @@ def inf_norm(x): kappa_ct = cantera_soln.thermal_conductivity assert inf_norm(kappa - kappa_ct) < 1.0e-12 - # NOTE: Individual species are exercised in Pyrometheus. - # Since the transport model enforce a singular-species case - # to avoid numerical issues when Yi -> 1, we can not test the - # individual species diffusivity. However, this tests that - # the single-species case is enforced correctly. - diff = fluid_state.tv.species_diffusivity - assert inf_norm(diff[i] - sing_diff) < 1.0e-15 + if not use_lewis: + # NOTE: Individual species are exercised in Pyrometheus. + # Since the transport model enforce a singular-species case + # to avoid numerical issues when Yi -> 1, we can not test the + # individual species diffusivity. However, this tests that + # the single-species case is enforced correctly. + diff = fluid_state.tv.species_diffusivity + assert inf_norm(diff[i] - sing_diff) < 1.0e-15 # prescribe an actual mixture cantera_soln.set_equivalence_ratio(phi=1.0, fuel="H2:1", @@ -154,7 +177,7 @@ def inf_norm(x): rhoin = can_rho * ones yin = can_y * ones - cv = make_conserved(dim=2, mass=rhoin, + cv = make_conserved(dim=dim, mass=rhoin, momentum=make_obj_array([zeros for _ in range(dim)]), energy=rhoin*gas_model.eos.get_internal_energy(tin, yin), species_mass=rhoin*yin) @@ -182,4 +205,4 @@ def inf_norm(x): diff = fluid_state.tv.species_diffusivity diff_ct = cantera_soln.mix_diff_coeffs for i in range(nspecies): - assert inf_norm(diff[i] - diff_ct[i]) < 1.0e-11 + assert inf_norm(diff[i] - diff_ct[i]) < 2.0e-11