From 65e9ddd6ddfd28841c0e339316569ccf90ab4140 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 21 Dec 2023 15:43:40 -0600 Subject: [PATCH 01/18] Keep cantera in version 2.6.0 for now --- conda-env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conda-env.yml b/conda-env.yml index 9654478d5..46e71b47c 100644 --- a/conda-env.yml +++ b/conda-env.yml @@ -23,7 +23,7 @@ dependencies: - pytest - pylint - pydocstyle -- cantera +- cantera=2.6.0 - h5py * nompi_* # Make sure cantera does not pull in MPI through h5py - gmsh - vtk From 9606842c2f944f1d6ef87d7dbf3d9b0b6055ce4c Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 22 Dec 2023 09:49:59 -0600 Subject: [PATCH 02/18] Update cantera to version 3.0.0 --- conda-env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conda-env.yml b/conda-env.yml index 46e71b47c..9654478d5 100644 --- a/conda-env.yml +++ b/conda-env.yml @@ -23,7 +23,7 @@ dependencies: - pytest - pylint - pydocstyle -- cantera=2.6.0 +- cantera - h5py * nompi_* # Make sure cantera does not pull in MPI through h5py - gmsh - vtk From 0fbf52b63ff9e3b2fad683e72ef86c655a0ca479 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 22 Dec 2023 11:24:26 -0600 Subject: [PATCH 03/18] Update cantera to version 3.0.0 --- conda-env.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/conda-env.yml b/conda-env.yml index 46e71b47c..9654478d5 100644 --- a/conda-env.yml +++ b/conda-env.yml @@ -23,7 +23,7 @@ dependencies: - pytest - pylint - pydocstyle -- cantera=2.6.0 +- cantera - h5py * nompi_* # Make sure cantera does not pull in MPI through h5py - gmsh - vtk From dc8c5aa8605d03b3bed76563deed01dad4f32568 Mon Sep 17 00:00:00 2001 From: Tulio Date: Sat, 6 Jan 2024 18:06:59 -0600 Subject: [PATCH 04/18] Point to cascading branches in pyrometheus --- requirements.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/requirements.txt b/requirements.txt index 9c503a2ff..f006b5dd5 100644 --- a/requirements.txt +++ b/requirements.txt @@ -21,7 +21,7 @@ git+https://github.com/mpi4py/mpi4py#egg=mpi4py --editable git+https://github.com/illinois-ceesd/meshmode.git@production#egg=meshmode --editable git+https://github.com/illinois-ceesd/grudge.git@production#egg=grudge --editable git+https://github.com/illinois-ceesd/pytato.git@production#egg=pytato ---editable git+https://github.com/pyrometheus/pyrometheus.git@tulio/check_reactions#egg=pyrometheus +--editable git+https://github.com/pyrometheus/pyrometheus.git@tulio/entropy#egg=pyrometheus --editable git+https://github.com/illinois-ceesd/logpyle.git#egg=logpyle --editable git+https://github.com/kaushikcfd/feinsum.git#egg=feinsum --editable git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile From 4b52f64d4878e20bd53adf543c1f5ab9c085fe8e Mon Sep 17 00:00:00 2001 From: Tulio Date: Sat, 6 Jan 2024 18:10:50 -0600 Subject: [PATCH 05/18] Cosmetic changes; fix header docs --- mirgecom/transport.py | 2 +- test/test_transport.py | 26 ++++++++++++-------------- 2 files changed, 13 insertions(+), 15 deletions(-) diff --git a/mirgecom/transport.py b/mirgecom/transport.py index 9260488fb..8d7a7137b 100644 --- a/mirgecom/transport.py +++ b/mirgecom/transport.py @@ -339,7 +339,7 @@ class MixtureAveragedTransport(TransportModel): def __init__(self, pyrometheus_mech, alpha=0.6, factor=1.0, lewis=None, epsilon=1e-4, singular_diffusivity=1e-6): - r"""Initialize power law coefficients and parameters. + r"""Initialize mixture averaged transport coefficients and parameters. Parameters ---------- diff --git a/test/test_transport.py b/test/test_transport.py index 8e0414619..da6ede531 100644 --- a/test/test_transport.py +++ b/test/test_transport.py @@ -100,14 +100,15 @@ def inf_norm(x): for pressin in ([0.25, 1.0]): for tempin in ([300.0, 600.0, 900.0, 1200.0, 1500.0, 1800.0, 2100.0]): - print(f"Testing (t,P) = ({tempin}, {pressin})") + cantera_soln.TP = tempin, pressin*cantera.one_atm + print(f"Testing (T, P) = ({cantera_soln.T}, {cantera_soln.P})") # Loop over each individual species by making a single-species mixture for i, name in enumerate(cantera_soln.species_names): - cantera_soln.TP = tempin, pressin*cantera.one_atm cantera_soln.Y = name + ":1" can_t, can_rho, can_y = cantera_soln.TDY + can_p = cantera_soln.P tin = can_t * ones rhoin = can_rho * ones yin = can_y * ones @@ -119,17 +120,18 @@ def inf_norm(x): fluid_state = make_fluid_state(cv, gas_model, tin) + assert inf_norm(fluid_state.temperature - tempin)/tempin < 1e-12 + assert inf_norm(fluid_state.pressure - can_p)/can_p < 1e-12 + # Viscosity mu = fluid_state.tv.viscosity mu_ct = cantera_soln.species_viscosities - err_mu = inf_norm(mu - mu_ct[i]) - assert err_mu < 1.0e-12 + assert inf_norm(mu) - mu_ct[i] < 1.0e-12 # Thermal conductivity kappa = fluid_state.tv.thermal_conductivity kappa_ct = cantera_soln.thermal_conductivity - err_kappa = inf_norm(kappa - kappa_ct) - assert err_kappa < 1.0e-12 + 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 @@ -137,8 +139,7 @@ def inf_norm(x): # individual species diffusivity. However, this tests that # the single-species case is enforced correctly. diff = fluid_state.tv.species_diffusivity - err_diff = inf_norm(diff[i] - sing_diff) - assert err_diff < 1.0e-15 + 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", @@ -170,18 +171,15 @@ def inf_norm(x): # Viscosity mu = fluid_state.tv.viscosity mu_ct = cantera_soln.viscosity - err_mu = inf_norm(mu - mu_ct) - assert err_mu < 1.0e-12 + assert inf_norm(mu - mu_ct) < 1.0e-12 # Thermal conductivity kappa = fluid_state.tv.thermal_conductivity kappa_ct = cantera_soln.thermal_conductivity - err_kappa = inf_norm(kappa - kappa_ct) - assert err_kappa < 1.0e-12 + assert inf_norm(kappa - kappa_ct) < 1.0e-12 # Species diffusivities diff = fluid_state.tv.species_diffusivity diff_ct = cantera_soln.mix_diff_coeffs for i in range(nspecies): - err_diff = inf_norm(diff[i] - diff_ct[i]) - assert err_diff < 1.0e-11 + assert inf_norm(diff[i] - diff_ct[i]) < 1.0e-11 From 3574f99ee70a61585598a5680d85defa6a7c4ad8 Mon Sep 17 00:00:00 2001 From: Tulio Date: Sat, 6 Jan 2024 21:37:54 -0600 Subject: [PATCH 06/18] Force cantera to be 3.0... --- conda-env.yml | 2 +- test/test_transport.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/conda-env.yml b/conda-env.yml index 9654478d5..07e510705 100644 --- a/conda-env.yml +++ b/conda-env.yml @@ -23,7 +23,7 @@ dependencies: - pytest - pylint - pydocstyle -- cantera +- cantera=3.0.0 - h5py * nompi_* # Make sure cantera does not pull in MPI through h5py - gmsh - vtk diff --git a/test/test_transport.py b/test/test_transport.py index da6ede531..e00347305 100644 --- a/test/test_transport.py +++ b/test/test_transport.py @@ -126,7 +126,7 @@ def inf_norm(x): # Viscosity mu = fluid_state.tv.viscosity mu_ct = cantera_soln.species_viscosities - assert inf_norm(mu) - mu_ct[i] < 1.0e-12 + assert inf_norm(mu - mu_ct[i]) < 1.0e-12 # Thermal conductivity kappa = fluid_state.tv.thermal_conductivity From 45f2fa2a25d965904c10c8ed2f0b392eab3929fe Mon Sep 17 00:00:00 2001 From: Tulio Date: Mon, 8 Jan 2024 20:14:07 -0600 Subject: [PATCH 07/18] Comment entropy due to pyro review --- test/test_eos.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_eos.py b/test/test_eos.py index db45d030f..9b695fd63 100644 --- a/test/test_eos.py +++ b/test/test_eos.py @@ -289,7 +289,7 @@ def inf_norm(x): can_p = cantera_soln.P can_e = cantera_soln.int_energy_mass can_h = cantera_soln.enthalpy_mass - can_s = cantera_soln.entropy_mass + # can_s = cantera_soln.entropy_mass can_c = cantera_soln.concentrations ones = dcoll.zeros(actx) + 1.0 @@ -302,7 +302,7 @@ def inf_norm(x): pyro_t = pyro_mechanism.get_temperature(pyro_e, tin, yin) pyro_p = pyro_mechanism.get_pressure(pyro_m, pyro_t, yin) pyro_h = pyro_mechanism.get_mixture_enthalpy_mass(pyro_t, yin) - pyro_s = pyro_mechanism.get_mixture_entropy_mass(pyro_p, pyro_t, yin) + # pyro_s = pyro_mechanism.get_mixture_entropy_mass(pyro_p, pyro_t, yin) pyro_c = pyro_mechanism.get_concentrations(pyro_m, yin) assert inf_norm((pyro_c - can_c)) < 1e-14 @@ -311,7 +311,7 @@ def inf_norm(x): assert inf_norm((pyro_p - can_p) / can_p) < 1e-14 assert inf_norm((pyro_e - can_e) / can_e) < 1e-10 assert inf_norm((pyro_h - can_h) / can_h) < 1e-10 - assert inf_norm((pyro_s - can_s) / can_s) < 1e-10 + # assert inf_norm((pyro_s - can_s) / can_s) < 1e-10 # Test the concentrations zero level y = -ones*y0s From bd38e046a9eceb2a16cf3cf8d8e926776a8608a0 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 12 Jan 2024 09:49:23 -0600 Subject: [PATCH 08/18] Simplify autoignition.py; remove redundant mixture.py --- examples/README.md | 1 - examples/autoignition.py | 86 ++----- examples/mixture.py | 478 --------------------------------------- test/test_eos.py | 6 +- 4 files changed, 17 insertions(+), 554 deletions(-) delete mode 100644 examples/mixture.py diff --git a/examples/README.md b/examples/README.md index e5ab0b4cf..96dfeed76 100644 --- a/examples/README.md +++ b/examples/README.md @@ -7,7 +7,6 @@ The examples and the unique features they exercise are as follows: - `autoignition.py`: Chemistry verification case with Pyrometheus - `heat-source.py`: Diffusion operator - `lump.py`: Lump advection, advection verification case -- `mixture.py`: Mixture EOS with Pyrometheus - `scalar-lump.py`: Scalar component lump advection verification case - `pulse.py`: Acoustic pulse in a box, wall boundary test case - `sod.py`: Sod's shock case: Fluid test case with strong shock diff --git a/examples/autoignition.py b/examples/autoignition.py index d9a49e69c..95105a59b 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -27,40 +27,34 @@ import numpy as np from functools import partial -from meshmode.mesh import BTAG_ALL -from mirgecom.discretization import create_discretization_collection -from grudge.dof_desc import DISCR_TAG_QUAD +from pytools.obj_array import make_obj_array from grudge.shortcuts import make_visualizer - +import grudge.op as op from logpyle import IntervalTimer, set_dt +from mirgecom.discretization import create_discretization_collection from mirgecom.euler import extract_vars_for_logging, units_for_logging -from mirgecom.euler import euler_operator from mirgecom.simutil import ( get_sim_timestep, generate_and_distribute_mesh, write_visfile, - ApplicationOptionsError + ApplicationOptionsError, + check_step, check_naninf_local, check_range_local ) from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point from mirgecom.integrators import rk4_step from mirgecom.steppers import advance_state -from mirgecom.boundary import AdiabaticSlipBoundary from mirgecom.initializers import Uniform from mirgecom.eos import PyrometheusMixture -from mirgecom.gas_model import ( - GasModel, make_fluid_state -) +from mirgecom.gas_model import GasModel, make_fluid_state from mirgecom.utils import force_evaluation from mirgecom.limiter import bound_preserving_limiter from mirgecom.fluid import make_conserved - from mirgecom.logging_quantities import ( initialize_logmgr, logmgr_add_many_discretization_quantities, logmgr_add_cl_device_info, logmgr_add_device_memory_usage, - set_sim_state ) import cantera @@ -164,13 +158,6 @@ def main(actx_class, use_leap=False, use_overintegration=False, nodes = actx.thaw(dcoll.nodes()) ones = dcoll.zeros(actx) + 1.0 - if use_overintegration: - quadrature_tag = DISCR_TAG_QUAD - else: - quadrature_tag = None - - ones = dcoll.zeros(actx) + 1.0 - vis_timer = None if logmgr: @@ -247,15 +234,9 @@ def main(actx_class, use_leap=False, use_overintegration=False, get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) eos = PyrometheusMixture(pyro_mechanism, temperature_guess=temperature_seed) - # {{{ Initialize simple transport model - from mirgecom.transport import MixtureAveragedTransport - transport_model = None - if viscous_terms_on: - transport_model = MixtureAveragedTransport(pyro_mechanism) - # }}} + # {{{ Initialize gas model - gas_model = GasModel(eos=eos, transport=transport_model) - from pytools.obj_array import make_obj_array + gas_model = GasModel(eos=eos) # }}} @@ -268,9 +249,6 @@ def main(actx_class, use_leap=False, use_overintegration=False, initializer = Uniform(dim=dim, pressure=can_p, temperature=can_t, species_mass_fractions=can_y, velocity=velocity) - my_boundary = AdiabaticSlipBoundary() - boundaries = {BTAG_ALL: my_boundary} - from mirgecom.viscous import get_viscous_timestep def get_dt(state): @@ -476,13 +454,11 @@ def my_write_restart(step, t, state, temperature_seed): write_restart_file(actx, rst_data, rst_fname, comm) def my_health_check(cv, dv): - import grudge.op as op health_error = False pressure = dv.pressure temperature = dv.temperature - from mirgecom.simutil import check_naninf_local, check_range_local if check_naninf_local(dcoll, "vol", pressure): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") @@ -518,6 +494,10 @@ def my_health_check(cv, dv): return health_error def my_pre_step(step, t, dt, state): + + if logmgr: + logmgr.tick_before() + cv, tseed = state # update temperature value @@ -526,11 +506,6 @@ def my_pre_step(step, t, dt, state): dv = fluid_state.dv try: - - if logmgr: - logmgr.tick_before() - - from mirgecom.simutil import check_step do_viz = check_step(step=step, interval=nviz) do_restart = check_step(step=step, interval=nrestart) do_health = check_step(step=step, interval=nhealth) @@ -568,34 +543,14 @@ def my_pre_step(step, t, dt, state): # my_write_restart(step=step, t=t, state=fluid_state) raise - return state, dt + return make_obj_array([cv, dv.temperature]), dt def my_post_step(step, t, dt, state): - cv, tseed = state - fluid_state = construct_fluid_state(cv, tseed) - - # Logmgr needs to know about EOS, dt, dim? - # imo this is a design/scope flaw if logmgr: set_dt(logmgr, dt) - set_sim_state(logmgr, dim, cv, gas_model.eos) logmgr.tick_after() - return make_obj_array([cv, fluid_state.temperature]), dt - - from mirgecom.inviscid import ( - inviscid_facial_flux_rusanov, - entropy_stable_inviscid_facial_flux_rusanov - ) - from mirgecom.gas_model import make_operator_fluid_states - from mirgecom.navierstokes import ns_operator - - inv_num_flux_func = entropy_stable_inviscid_facial_flux_rusanov if use_esdg \ - else inviscid_facial_flux_rusanov - - fluid_operator = euler_operator - if viscous_terms_on: - fluid_operator = ns_operator + return state, dt def my_rhs(t, state): cv, tseed = state @@ -604,20 +559,9 @@ def my_rhs(t, state): temperature_seed=tseed, limiter_func=_limit_fluid_cv) - fluid_operator_states = make_operator_fluid_states( - dcoll, fluid_state, gas_model, boundaries=boundaries, - quadrature_tag=quadrature_tag, limiter_func=_limit_fluid_cv) - - fluid_rhs = fluid_operator( - dcoll, state=fluid_state, gas_model=gas_model, time=t, - boundaries=boundaries, operator_states_quad=fluid_operator_states, - quadrature_tag=quadrature_tag, use_esdg=use_esdg, - inviscid_numerical_flux_func=inv_num_flux_func) - chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature) tseed_rhs = fluid_state.temperature - tseed - cv_rhs = fluid_rhs + chem_rhs - return make_obj_array([cv_rhs, tseed_rhs]) + return make_obj_array([chem_rhs, tseed_rhs]) current_dt = get_sim_timestep(dcoll, current_fluid_state, current_t, current_dt, current_cfl, t_final, constant_cfl) diff --git a/examples/mixture.py b/examples/mixture.py deleted file mode 100644 index 4beba862c..000000000 --- a/examples/mixture.py +++ /dev/null @@ -1,478 +0,0 @@ -"""Demonstrate simple gas mixture with Pyrometheus.""" - -__copyright__ = """ -Copyright (C) 2020 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 logging -import numpy as np -from functools import partial - -from meshmode.mesh import BTAG_ALL -from grudge.shortcuts import make_visualizer -from grudge.dof_desc import DISCR_TAG_QUAD - -from mirgecom.discretization import create_discretization_collection -from mirgecom.euler import euler_operator -from mirgecom.simutil import ( - get_sim_timestep, - generate_and_distribute_mesh -) -from mirgecom.io import make_init_message -from mirgecom.mpi import mpi_entry_point - -from mirgecom.integrators import rk4_step -from mirgecom.steppers import advance_state -from mirgecom.boundary import ( - PrescribedFluidBoundary, - AdiabaticSlipBoundary -) -from mirgecom.initializers import Uniform -from mirgecom.eos import PyrometheusMixture - -import cantera - -from logpyle import IntervalTimer, set_dt -from mirgecom.euler import extract_vars_for_logging, units_for_logging -from mirgecom.logging_quantities import ( - initialize_logmgr, - logmgr_add_many_discretization_quantities, - logmgr_add_cl_device_info, - logmgr_add_device_memory_usage, - set_sim_state -) - -logger = logging.getLogger(__name__) - - -class MyRuntimeError(RuntimeError): - """Simple exception to kill the simulation.""" - - pass - - -@mpi_entry_point -def main(actx_class, use_esdg=False, - use_leap=False, casename=None, rst_filename=None, - log_dependent=False, use_overintegration=False): - """Drive example.""" - if casename is None: - casename = "mirgecom" - - from mpi4py import MPI - comm = MPI.COMM_WORLD - rank = comm.Get_rank() - nparts = comm.Get_size() - - from mirgecom.simutil import global_reduce as _global_reduce - global_reduce = partial(_global_reduce, comm=comm) - - logmgr = initialize_logmgr(True, - filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) - - from mirgecom.array_context import initialize_actx, actx_class_is_profiling - actx = initialize_actx(actx_class, comm) - queue = getattr(actx, "queue", None) - use_profiling = actx_class_is_profiling(actx_class) - - # timestepping control - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step - - t_final = 1e-8 - current_cfl = 1.0 - current_dt = 1e-9 - current_t = 0 - current_step = 0 - constant_cfl = False - - # some i/o frequencies - nstatus = 1 - nhealth = 1 - nrestart = 5 - nviz = 100 - - dim = 2 - rst_path = "restart_data/" - rst_pattern = ( - rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" - ) - if rst_filename: # read the grid from restart data - rst_filename = f"{rst_filename}-{rank:04d}.pkl" - from mirgecom.restart import read_restart_data - restart_data = read_restart_data(actx, rst_filename) - local_mesh = restart_data["local_mesh"] - local_nelements = local_mesh.nelements - global_nelements = restart_data["global_nelements"] - assert restart_data["num_parts"] == nparts - else: # generate the grid from scratch - nel_1d = 16 - box_ll = -5.0 - box_ur = 5.0 - from meshmode.mesh.generation import generate_regular_rect_mesh - generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,)*dim, - b=(box_ur,) * dim, nelements_per_axis=(nel_1d,)*dim) - local_mesh, global_nelements = generate_and_distribute_mesh(comm, - generate_mesh) - local_nelements = local_mesh.nelements - - order = 3 - dcoll = create_discretization_collection(actx, local_mesh, order=order) - nodes = actx.thaw(dcoll.nodes()) - - if use_overintegration: - quadrature_tag = DISCR_TAG_QUAD - else: - quadrature_tag = None - - vis_timer = None - - if logmgr: - logmgr_add_cl_device_info(logmgr, queue) - logmgr_add_device_memory_usage(logmgr, queue) - - vis_timer = IntervalTimer("t_vis", "Time spent visualizing") - logmgr.add_quantity(vis_timer) - - logmgr.add_watches([ - ("step.max", "step = {value}, "), - ("t_sim.max", "sim time: {value:1.6e} s\n"), - ("t_step.max", "------- step walltime: {value:6g} s, "), - ("t_log.max", "log walltime: {value:6g} s") - ]) - - if log_dependent: - logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, - extract_vars_for_logging, - units_for_logging) - logmgr.add_watches([ - ("min_pressure", "\n------- P (min, max) (Pa) = ({value:1.9e}, "), - ("max_pressure", "{value:1.9e})\n"), - ("min_temperature", "------- T (min, max) (K) = ({value:7g}, "), - ("max_temperature", "{value:7g})\n")]) - - # Pyrometheus initialization - from mirgecom.mechanisms import get_mechanism_input - mech_input = get_mechanism_input("uiuc_7sp") - sol = cantera.Solution(name="gas", yaml=mech_input) - from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera - pyrometheus_mechanism = \ - get_pyrometheus_wrapper_class_from_cantera(sol)(actx.np) - - nspecies = pyrometheus_mechanism.num_species - eos = PyrometheusMixture(pyrometheus_mechanism, temperature_guess=300) - from mirgecom.gas_model import GasModel, make_fluid_state - gas_model = GasModel(eos=eos) - from pytools.obj_array import make_obj_array - - y0s = np.zeros(shape=(nspecies,)) - for i in range(nspecies-1): - y0s[i] = 1.0 / (10.0 ** (i + 1)) - spec_sum = sum([y0s[i] for i in range(nspecies-1)]) - y0s[nspecies-1] = 1.0 - spec_sum - - # Mixture defaults to STP (p, T) = (1atm, 300K) - velocity = np.zeros(shape=(dim,)) + 1.0 - initializer = Uniform(dim=dim, species_mass_fractions=y0s, velocity=velocity, - pressure=101325.0, temperature=300.0) - - def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): - actx = state_minus.array_context - bnd_discr = dcoll.discr_from_dd(dd_bdry) - nodes = actx.thaw(bnd_discr.nodes()) - return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, - **kwargs), gas_model, - temperature_seed=state_minus.temperature) - - if False: - my_boundary = AdiabaticSlipBoundary() - boundaries = {BTAG_ALL: my_boundary} - else: - boundaries = { - BTAG_ALL: PrescribedFluidBoundary(boundary_state_func=boundary_solution) - } - - if rst_filename: - current_t = restart_data["t"] - current_step = restart_data["step"] - current_cv = restart_data["cv"] - tseed = restart_data["temperature_seed"] - if logmgr: - from mirgecom.logging_quantities import logmgr_set_time - logmgr_set_time(logmgr, current_step, current_t) - else: - # Set the current state from time 0 - current_cv = initializer(x_vec=nodes, eos=eos) - tseed = 300.0 - - def get_fluid_state(cv, tseed): - return make_fluid_state(cv=cv, gas_model=gas_model, - temperature_seed=tseed) - - construct_fluid_state = actx.compile(get_fluid_state) - current_state = construct_fluid_state(current_cv, tseed) - - visualizer = make_visualizer(dcoll) - initname = initializer.__class__.__name__ - eosname = eos.__class__.__name__ - init_message = make_init_message(dim=dim, order=order, - nelements=local_nelements, - global_nelements=global_nelements, - dt=current_dt, t_final=t_final, nstatus=nstatus, - nviz=nviz, cfl=current_cfl, - constant_cfl=constant_cfl, initname=initname, - eosname=eosname, casename=casename) - if rank == 0: - logger.info(init_message) - - def my_write_status(component_errors, dv=None): - status_msg = ( - "------- errors=" - + ", ".join("%.3g" % en for en in component_errors)) - if ((dv is not None) and (not log_dependent)): - temp = dv.temperature - press = dv.pressure - - from grudge.op import nodal_min_loc, nodal_max_loc - tmin = global_reduce(actx.to_numpy(nodal_min_loc(dcoll, "vol", temp)), - op="min") - tmax = global_reduce(actx.to_numpy(nodal_max_loc(dcoll, "vol", temp)), - op="max") - pmin = global_reduce(actx.to_numpy(nodal_min_loc(dcoll, "vol", press)), - op="min") - pmax = global_reduce(actx.to_numpy(nodal_max_loc(dcoll, "vol", press)), - op="max") - dv_status_msg = f"\nP({pmin}, {pmax}), T({tmin}, {tmax})" - status_msg = status_msg + dv_status_msg - - if rank == 0: - logger.info(status_msg) - if rank == 0: - logger.info(status_msg) - - def my_write_viz(step, t, state, dv, exact=None, resid=None): - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, time=t) - if resid is None: - resid = state - exact - viz_fields = [("cv", state), ("dv", dv)] - from mirgecom.simutil import write_visfile - write_visfile(dcoll, viz_fields, visualizer, vizname=casename, - step=step, t=t, overwrite=True, vis_timer=vis_timer, - comm=comm) - - def my_write_restart(step, t, state, tseed): - rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - if rst_fname != rst_filename: - rst_data = { - "local_mesh": local_mesh, - "cv": state, - "temperature_seed": tseed, - "t": t, - "step": step, - "order": order, - "global_nelements": global_nelements, - "num_parts": nparts - } - from mirgecom.restart import write_restart_file - write_restart_file(actx, rst_data, rst_fname, comm) - - def my_health_check(dv, component_errors): - health_error = False - from mirgecom.simutil import check_naninf_local, check_range_local - if check_naninf_local(dcoll, "vol", dv.pressure) \ - or check_range_local(dcoll, "vol", dv.pressure, 1e5, 1.1e5): - health_error = True - logger.info(f"{rank=}: Invalid pressure data found.") - - exittol = .09 - if max(component_errors) > exittol: - health_error = True - if rank == 0: - logger.info("Solution diverged from exact soln.") - - return health_error - - def my_pre_step(step, t, dt, state): - cv, tseed = state - fluid_state = construct_fluid_state(cv, tseed) - dv = fluid_state.dv - - try: - exact = None - component_errors = None - - if logmgr: - logmgr.tick_before() - - from mirgecom.simutil import check_step - do_viz = check_step(step=step, interval=nviz) - do_restart = check_step(step=step, interval=nrestart) - do_health = check_step(step=step, interval=nhealth) - do_status = check_step(step=step, interval=nstatus) - - if do_health: - exact = initializer(x_vec=nodes, eos=eos, time=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(dcoll, cv, exact) - health_errors = global_reduce( - my_health_check(dv, component_errors), op="lor") - if health_errors: - if rank == 0: - logger.info("Fluid solution failed health check.") - raise MyRuntimeError("Failed simulation health check.") - - if do_restart: - my_write_restart(step=step, t=t, state=cv, tseed=tseed) - - if do_viz: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, time=t) - resid = state - exact - my_write_viz(step=step, t=t, state=cv, dv=dv, exact=exact, - resid=resid) - - if do_status: - if component_errors is None: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, time=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(dcoll, cv, exact) - my_write_status(component_errors, dv=dv) - - except MyRuntimeError: - if rank == 0: - logger.info("Errors detected; attempting graceful exit.") - my_write_viz(step=step, t=t, state=cv, dv=dv) - my_write_restart(step=step, t=t, state=cv, tseed=tseed) - raise - - dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, t_final, - constant_cfl) - return state, dt - - def my_post_step(step, t, dt, state): - cv, tseed = state - fluid_state = construct_fluid_state(cv, tseed) - tseed = fluid_state.temperature - # Logmgr needs to know about EOS, dt, dim? - # imo this is a design/scope flaw - if logmgr: - set_dt(logmgr, dt) - set_sim_state(logmgr, dim, cv, eos) - logmgr.tick_after() - return make_obj_array([fluid_state.cv, tseed]), dt - - def my_rhs(t, state): - cv, tseed = state - fluid_state = make_fluid_state(cv, gas_model, temperature_seed=tseed) - return make_obj_array( - [euler_operator(dcoll, state=fluid_state, time=t, - boundaries=boundaries, gas_model=gas_model, - quadrature_tag=quadrature_tag, use_esdg=use_esdg), - 0*tseed]) - - current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, - current_cfl, t_final, constant_cfl) - - current_step, current_t, advanced_state = \ - advance_state(rhs=my_rhs, timestepper=timestepper, - pre_step_callback=my_pre_step, - post_step_callback=my_post_step, dt=current_dt, - state=make_obj_array([current_state.cv, - current_state.temperature]), - t=current_t, t_final=t_final) - - # Dump the final data - if rank == 0: - logger.info("Checkpointing final state ...") - - current_cv, tseed = advanced_state - current_state = make_fluid_state(current_cv, gas_model, temperature_seed=tseed) - final_dv = current_state.dv - final_exact = initializer(x_vec=nodes, eos=eos, time=current_t) - final_resid = current_state.cv - final_exact - my_write_viz(step=current_step, t=current_t, state=current_cv, dv=final_dv, - exact=final_exact, resid=final_resid) - my_write_restart(step=current_step, t=current_t, state=current_state.cv, - tseed=tseed) - - if logmgr: - logmgr.close() - elif use_profiling: - print(actx.tabulate_profiling_data()) - - finish_tol = 1e-16 - assert np.abs(current_t - t_final) < finish_tol - - -if __name__ == "__main__": - import argparse - casename = "uiuc-mixture" - parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}") - parser.add_argument("--lazy", action="store_true", - help="switch to a lazy computation mode") - parser.add_argument("--profiling", action="store_true", - help="turn on detailed performance profiling") - parser.add_argument("--overintegration", action="store_true", - help="use overintegration in the RHS computations") - parser.add_argument("--esdg", action="store_true", - help="use flux-differencing/entropy stable DG for inviscid computations.") - parser.add_argument("--leap", action="store_true", - help="use leap timestepper") - parser.add_argument("--numpy", action="store_true", - help="use numpy-based eager actx.") - parser.add_argument("--restart_file", help="root name of restart file") - parser.add_argument("--casename", help="casename to use for i/o") - args = parser.parse_args() - from warnings import warn - warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") - log_dependent = False - - from mirgecom.simutil import ApplicationOptionsError - if args.esdg: - if not args.lazy and not args.numpy: - raise ApplicationOptionsError("ESDG requires lazy or numpy context.") - if not args.overintegration: - warn("ESDG requires overintegration, enabling --overintegration.") - - from mirgecom.array_context import get_reasonable_array_context_class - actx_class = get_reasonable_array_context_class( - lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) - - logging.basicConfig(format="%(message)s", level=logging.INFO) - if args.casename: - casename = args.casename - rst_filename = None - - if args.restart_file: - rst_filename = args.restart_file - - main(actx_class, use_leap=args.leap, - casename=casename, rst_filename=rst_filename, - use_overintegration=args.overintegration or args.esdg, - use_esdg=args.esdg, log_dependent=log_dependent) - -# vim: foldmethod=marker diff --git a/test/test_eos.py b/test/test_eos.py index 9b695fd63..2d59a1ea5 100644 --- a/test/test_eos.py +++ b/test/test_eos.py @@ -239,10 +239,8 @@ def inf_norm(x): def test_pyrometheus_mechanisms(ctx_factory, mechname): """Test known pyrometheus mechanisms. - This test reproduces a pyrometheus-native test in the MIRGE context. - - Tests that the Pyrometheus mechanism code gets the same thermo properties - as the corresponding mechanism in Cantera. + This test reproduces a pyrometheus-native test in the MIRGE context and + compare thermo properties to the corresponding mechanism in Cantera. """ cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) From a33f66f3e8543c1fd55372d8c1ca56bd77b19411 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 12 Jan 2024 12:16:58 -0600 Subject: [PATCH 09/18] Add a search for mechanisms in local path; Cosmetics in autoignition --- examples/autoignition.py | 65 +++++++++++++++++++++++++-------- mirgecom/mechanisms/__init__.py | 17 ++++++--- 2 files changed, 62 insertions(+), 20 deletions(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index 95105a59b..9141e7a40 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -23,6 +23,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import sys import logging import numpy as np from functools import partial @@ -59,7 +60,31 @@ import cantera + +class SingleLevelFilter(logging.Filter): + def __init__(self, passlevel, reject): + self.passlevel = passlevel + self.reject = reject + + def filter(self, record): + if self.reject: + return (record.levelno != self.passlevel) + else: + return (record.levelno == self.passlevel) + + +h1 = logging.StreamHandler(sys.stdout) +f1 = SingleLevelFilter(logging.INFO, False) +h1.addFilter(f1) +root_logger = logging.getLogger() +root_logger.addHandler(h1) +h2 = logging.StreamHandler(sys.stderr) +f2 = SingleLevelFilter(logging.INFO, True) +h2.addFilter(f2) +root_logger.addHandler(h2) + logger = logging.getLogger(__name__) +logger.setLevel(logging.DEBUG) class MyRuntimeError(RuntimeError): @@ -97,6 +122,8 @@ def main(actx_class, use_leap=False, use_overintegration=False, nel_1d = 8 order = 1 + mech_file = "uiuc_18sp" + # {{{ Time stepping control # This example runs only 3 steps by default (to keep CI ~short) @@ -112,9 +139,9 @@ def main(actx_class, use_leap=False, use_overintegration=False, # Time loop control parameters current_step = 0 - t_final = 1e-8 + t_final = 1e-6 current_cfl = 1.0 - current_dt = 1e-9 + current_dt = 1e-7 current_t = 0 constant_cfl = False @@ -191,18 +218,18 @@ def main(actx_class, use_leap=False, use_overintegration=False, # --- Note: Users may add their own mechanism input file by dropping it into # --- mirgecom/mechanisms alongside the other mech input files. from mirgecom.mechanisms import get_mechanism_input - mech_input = get_mechanism_input("uiuc_7sp") + mech_input = get_mechanism_input(mech_file) cantera_soln = cantera.Solution(name="gas", yaml=mech_input) nspecies = cantera_soln.n_species # Initial temperature, pressure, and mixutre mole fractions are needed to # set up the initial state in Cantera. - temperature_seed = 1500.0 # Initial temperature hot enough to burn + temperature_seed = 1200.0 # Initial temperature hot enough to burn # Parameters for calculating the amounts of fuel, oxidizer, and inert species # which directly sets the species fractions inside cantera - cantera_soln.set_equivalence_ratio(phi=1.0, fuel="C2H4:1,H2:1", + cantera_soln.set_equivalence_ratio(phi=1.0, fuel="C2H4:1", oxidizer={"O2": 1.0, "N2": 3.76}) x = cantera_soln.X @@ -234,6 +261,12 @@ def main(actx_class, use_leap=False, use_overintegration=False, get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) eos = PyrometheusMixture(pyro_mechanism, temperature_guess=temperature_seed) + # print out the mechanism for inspecion + # import pyrometheus as pyro + # with open(f"mechanism.py", "w") as mech_file: + # code = pyro.codegen.python.gen_thermochem_code(cantera_soln) + # print(code, file=mech_file) + # {{{ Initialize gas model gas_model = GasModel(eos=eos) @@ -361,6 +394,7 @@ def get_fluid_state(cv, tseed): temperature_seed = temperature_seed * ones current_cv = force_evaluation(actx, current_cv) + temperature_seed = force_evaluation(actx, temperature_seed) # The temperature_seed going into this function is: # - At time 0: the initial temperature input data (maybe from Cantera) @@ -402,7 +436,7 @@ def get_fluid_state(cv, tseed): f" {eq_pressure=}, {eq_temperature=}," f" {eq_density=}, {eq_mass_fractions=}") - def my_write_status(dt, cfl, dv=None): + def my_write_status(t, dt, cfl, dv=None): status_msg = f"------ {dt=}" if constant_cfl else f"----- {cfl=}" if ((dv is not None) and (not log_dependent)): @@ -418,7 +452,7 @@ def my_write_status(dt, cfl, dv=None): op="min") pmax = global_reduce(actx.to_numpy(nodal_max_loc(dcoll, "vol", press)), op="max") - dv_status_msg = f"\nP({pmin}, {pmax}), T({tmin}, {tmax})" + dv_status_msg = f"\n{t}, P({pmin}, {pmax}), T({tmin}, {tmax})" status_msg = status_msg + dv_status_msg if rank == 0: @@ -463,16 +497,17 @@ def my_health_check(cv, dv): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - if check_range_local(dcoll, "vol", pressure, 1e5, 2.6e5): - health_error = True - logger.info(f"{rank=}: Pressure range violation.") + # if check_range_local(dcoll, "vol", pressure, 1e5, 2.6e5): + # health_error = True + # logger.info(f"{rank=}: Pressure range violation.") if check_naninf_local(dcoll, "vol", temperature): health_error = True logger.info(f"{rank=}: Invalid temperature data found.") - if check_range_local(dcoll, "vol", temperature, 1.498e3, 1.6e3): - health_error = True - logger.info(f"{rank=}: Temperature range violation.") + + # if check_range_local(dcoll, "vol", temperature, 1.198e3, 1.3e3): + # health_error = True + # logger.info(f"{rank=}: Temperature range violation.") # This check is the temperature convergence check # The current *temperature* is what Pyrometheus gets @@ -521,7 +556,7 @@ def my_pre_step(step, t, dt, state): ts_field, cfl, dt = my_get_timestep(t=t, dt=dt, state=fluid_state) if do_status: - my_write_status(dt=dt, cfl=cfl, dv=dv) + my_write_status(t=t, dt=dt, cfl=cfl, dv=dv) if do_restart: my_write_restart(step=step, t=t, state=fluid_state, @@ -587,7 +622,7 @@ def my_rhs(t, state): my_write_viz(step=current_step, t=current_t, dt=dt, state=final_cv, dv=final_dv, production_rates=final_dm, heat_release_rate=final_heat_rls, ts_field=ts_field, cfl=cfl) - my_write_status(dt=dt, cfl=cfl, dv=final_dv) + my_write_status(t=current_dt, dt=dt, cfl=cfl, dv=final_dv) my_write_restart(step=current_step, t=current_t, state=final_fluid_state, temperature_seed=tseed) diff --git a/mirgecom/mechanisms/__init__.py b/mirgecom/mechanisms/__init__.py index 40b986cb1..64249ffaa 100644 --- a/mirgecom/mechanisms/__init__.py +++ b/mirgecom/mechanisms/__init__.py @@ -65,11 +65,18 @@ def import_mechdata(): def get_mechanism_input(mechanism_name: str) -> str: """Get the contents of a mechanism YAML input file.""" - if mechanism_name == "uiuc": - mechanism_name = "uiuc_7sp" - from warnings import warn - warn("The uiuc mechanism was updated in Q3 2023. Use 'uiuc_7sp' instead.", - stacklevel=2) + import os + from pathlib import Path + mech_data = import_mechdata() mech_file = mech_data / get_mechanism_file_name(mechanism_name) + if not os.path.isfile(mech_file): + from warnings import warn + warn("Mechanism not found in default location.", stacklevel=2) + try: + current_path = os.path.abspath(os.getcwd()) + "/" + mech_file = Path(current_path + get_mechanism_file_name(mechanism_name)) + warn("Trying " + current_path, stacklevel=2) + except: + warn("Mechanism not found.", stacklevel=2) return mech_file.read_text() From 2d2d6126728a8c6ab323fb290666b818f1e6e68a Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 12 Jan 2024 12:29:26 -0600 Subject: [PATCH 10/18] Add a search for mechanisms in local path --- mirgecom/mechanisms/__init__.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/mirgecom/mechanisms/__init__.py b/mirgecom/mechanisms/__init__.py index 64249ffaa..5752ea84d 100644 --- a/mirgecom/mechanisms/__init__.py +++ b/mirgecom/mechanisms/__init__.py @@ -73,10 +73,7 @@ def get_mechanism_input(mechanism_name: str) -> str: if not os.path.isfile(mech_file): from warnings import warn warn("Mechanism not found in default location.", stacklevel=2) - try: - current_path = os.path.abspath(os.getcwd()) + "/" - mech_file = Path(current_path + get_mechanism_file_name(mechanism_name)) - warn("Trying " + current_path, stacklevel=2) - except: - warn("Mechanism not found.", stacklevel=2) + current_path = os.path.abspath(os.getcwd()) + "/" + mech_file = Path(current_path + get_mechanism_file_name(mechanism_name)) + warn("Trying " + current_path, stacklevel=2) return mech_file.read_text() From 19d53504d9e27e72ddc01af8fb17fb6ad143d829 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 12 Jan 2024 12:34:29 -0600 Subject: [PATCH 11/18] flake8 --- examples/autoignition.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index 9141e7a40..b73271c38 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -497,17 +497,17 @@ def my_health_check(cv, dv): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - # if check_range_local(dcoll, "vol", pressure, 1e5, 2.6e5): - # health_error = True - # logger.info(f"{rank=}: Pressure range violation.") + if check_range_local(dcoll, "vol", pressure, 1e5, 2.6e5): + health_error = True + logger.info(f"{rank=}: Pressure range violation.") if check_naninf_local(dcoll, "vol", temperature): health_error = True logger.info(f"{rank=}: Invalid temperature data found.") - # if check_range_local(dcoll, "vol", temperature, 1.198e3, 1.3e3): - # health_error = True - # logger.info(f"{rank=}: Temperature range violation.") + if check_range_local(dcoll, "vol", temperature, 1.198e3, 1.3e3): + health_error = True + logger.info(f"{rank=}: Temperature range violation.") # This check is the temperature convergence check # The current *temperature* is what Pyrometheus gets From 8709d24f87c51db59f1065382add42f0de144ddb Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 12 Jan 2024 13:47:56 -0600 Subject: [PATCH 12/18] Put back the 7sp mechanism --- examples/autoignition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index b73271c38..ed4a434f1 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -122,7 +122,7 @@ def main(actx_class, use_leap=False, use_overintegration=False, nel_1d = 8 order = 1 - mech_file = "uiuc_18sp" + mech_file = "uiuc_7sp" # {{{ Time stepping control From 56c7d244538acf41066e5bab296b5f31b894c5fb Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 06:10:05 -0600 Subject: [PATCH 13/18] Do NOT remove mixture.py in this PR --- examples/mixture.py | 478 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 478 insertions(+) create mode 100644 examples/mixture.py diff --git a/examples/mixture.py b/examples/mixture.py new file mode 100644 index 000000000..4beba862c --- /dev/null +++ b/examples/mixture.py @@ -0,0 +1,478 @@ +"""Demonstrate simple gas mixture with Pyrometheus.""" + +__copyright__ = """ +Copyright (C) 2020 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 logging +import numpy as np +from functools import partial + +from meshmode.mesh import BTAG_ALL +from grudge.shortcuts import make_visualizer +from grudge.dof_desc import DISCR_TAG_QUAD + +from mirgecom.discretization import create_discretization_collection +from mirgecom.euler import euler_operator +from mirgecom.simutil import ( + get_sim_timestep, + generate_and_distribute_mesh +) +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point + +from mirgecom.integrators import rk4_step +from mirgecom.steppers import advance_state +from mirgecom.boundary import ( + PrescribedFluidBoundary, + AdiabaticSlipBoundary +) +from mirgecom.initializers import Uniform +from mirgecom.eos import PyrometheusMixture + +import cantera + +from logpyle import IntervalTimer, set_dt +from mirgecom.euler import extract_vars_for_logging, units_for_logging +from mirgecom.logging_quantities import ( + initialize_logmgr, + logmgr_add_many_discretization_quantities, + logmgr_add_cl_device_info, + logmgr_add_device_memory_usage, + set_sim_state +) + +logger = logging.getLogger(__name__) + + +class MyRuntimeError(RuntimeError): + """Simple exception to kill the simulation.""" + + pass + + +@mpi_entry_point +def main(actx_class, use_esdg=False, + use_leap=False, casename=None, rst_filename=None, + log_dependent=False, use_overintegration=False): + """Drive example.""" + if casename is None: + casename = "mirgecom" + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = comm.Get_rank() + nparts = comm.Get_size() + + from mirgecom.simutil import global_reduce as _global_reduce + global_reduce = partial(_global_reduce, comm=comm) + + logmgr = initialize_logmgr(True, + filename=f"{casename}.sqlite", mode="wu", mpi_comm=comm) + + from mirgecom.array_context import initialize_actx, actx_class_is_profiling + actx = initialize_actx(actx_class, comm) + queue = getattr(actx, "queue", None) + use_profiling = actx_class_is_profiling(actx_class) + + # timestepping control + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step + + t_final = 1e-8 + current_cfl = 1.0 + current_dt = 1e-9 + current_t = 0 + current_step = 0 + constant_cfl = False + + # some i/o frequencies + nstatus = 1 + nhealth = 1 + nrestart = 5 + nviz = 100 + + dim = 2 + rst_path = "restart_data/" + rst_pattern = ( + rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" + ) + if rst_filename: # read the grid from restart data + rst_filename = f"{rst_filename}-{rank:04d}.pkl" + from mirgecom.restart import read_restart_data + restart_data = read_restart_data(actx, rst_filename) + local_mesh = restart_data["local_mesh"] + local_nelements = local_mesh.nelements + global_nelements = restart_data["global_nelements"] + assert restart_data["num_parts"] == nparts + else: # generate the grid from scratch + nel_1d = 16 + box_ll = -5.0 + box_ur = 5.0 + from meshmode.mesh.generation import generate_regular_rect_mesh + generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,)*dim, + b=(box_ur,) * dim, nelements_per_axis=(nel_1d,)*dim) + local_mesh, global_nelements = generate_and_distribute_mesh(comm, + generate_mesh) + local_nelements = local_mesh.nelements + + order = 3 + dcoll = create_discretization_collection(actx, local_mesh, order=order) + nodes = actx.thaw(dcoll.nodes()) + + if use_overintegration: + quadrature_tag = DISCR_TAG_QUAD + else: + quadrature_tag = None + + vis_timer = None + + if logmgr: + logmgr_add_cl_device_info(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("t_sim.max", "sim time: {value:1.6e} s\n"), + ("t_step.max", "------- step walltime: {value:6g} s, "), + ("t_log.max", "log walltime: {value:6g} s") + ]) + + if log_dependent: + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, + extract_vars_for_logging, + units_for_logging) + logmgr.add_watches([ + ("min_pressure", "\n------- P (min, max) (Pa) = ({value:1.9e}, "), + ("max_pressure", "{value:1.9e})\n"), + ("min_temperature", "------- T (min, max) (K) = ({value:7g}, "), + ("max_temperature", "{value:7g})\n")]) + + # Pyrometheus initialization + from mirgecom.mechanisms import get_mechanism_input + mech_input = get_mechanism_input("uiuc_7sp") + sol = cantera.Solution(name="gas", yaml=mech_input) + from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera + pyrometheus_mechanism = \ + get_pyrometheus_wrapper_class_from_cantera(sol)(actx.np) + + nspecies = pyrometheus_mechanism.num_species + eos = PyrometheusMixture(pyrometheus_mechanism, temperature_guess=300) + from mirgecom.gas_model import GasModel, make_fluid_state + gas_model = GasModel(eos=eos) + from pytools.obj_array import make_obj_array + + y0s = np.zeros(shape=(nspecies,)) + for i in range(nspecies-1): + y0s[i] = 1.0 / (10.0 ** (i + 1)) + spec_sum = sum([y0s[i] for i in range(nspecies-1)]) + y0s[nspecies-1] = 1.0 - spec_sum + + # Mixture defaults to STP (p, T) = (1atm, 300K) + velocity = np.zeros(shape=(dim,)) + 1.0 + initializer = Uniform(dim=dim, species_mass_fractions=y0s, velocity=velocity, + pressure=101325.0, temperature=300.0) + + def boundary_solution(dcoll, dd_bdry, gas_model, state_minus, **kwargs): + actx = state_minus.array_context + bnd_discr = dcoll.discr_from_dd(dd_bdry) + nodes = actx.thaw(bnd_discr.nodes()) + return make_fluid_state(initializer(x_vec=nodes, eos=gas_model.eos, + **kwargs), gas_model, + temperature_seed=state_minus.temperature) + + if False: + my_boundary = AdiabaticSlipBoundary() + boundaries = {BTAG_ALL: my_boundary} + else: + boundaries = { + BTAG_ALL: PrescribedFluidBoundary(boundary_state_func=boundary_solution) + } + + if rst_filename: + current_t = restart_data["t"] + current_step = restart_data["step"] + current_cv = restart_data["cv"] + tseed = restart_data["temperature_seed"] + if logmgr: + from mirgecom.logging_quantities import logmgr_set_time + logmgr_set_time(logmgr, current_step, current_t) + else: + # Set the current state from time 0 + current_cv = initializer(x_vec=nodes, eos=eos) + tseed = 300.0 + + def get_fluid_state(cv, tseed): + return make_fluid_state(cv=cv, gas_model=gas_model, + temperature_seed=tseed) + + construct_fluid_state = actx.compile(get_fluid_state) + current_state = construct_fluid_state(current_cv, tseed) + + visualizer = make_visualizer(dcoll) + initname = initializer.__class__.__name__ + eosname = eos.__class__.__name__ + init_message = make_init_message(dim=dim, order=order, + nelements=local_nelements, + global_nelements=global_nelements, + dt=current_dt, t_final=t_final, nstatus=nstatus, + nviz=nviz, cfl=current_cfl, + constant_cfl=constant_cfl, initname=initname, + eosname=eosname, casename=casename) + if rank == 0: + logger.info(init_message) + + def my_write_status(component_errors, dv=None): + status_msg = ( + "------- errors=" + + ", ".join("%.3g" % en for en in component_errors)) + if ((dv is not None) and (not log_dependent)): + temp = dv.temperature + press = dv.pressure + + from grudge.op import nodal_min_loc, nodal_max_loc + tmin = global_reduce(actx.to_numpy(nodal_min_loc(dcoll, "vol", temp)), + op="min") + tmax = global_reduce(actx.to_numpy(nodal_max_loc(dcoll, "vol", temp)), + op="max") + pmin = global_reduce(actx.to_numpy(nodal_min_loc(dcoll, "vol", press)), + op="min") + pmax = global_reduce(actx.to_numpy(nodal_max_loc(dcoll, "vol", press)), + op="max") + dv_status_msg = f"\nP({pmin}, {pmax}), T({tmin}, {tmax})" + status_msg = status_msg + dv_status_msg + + if rank == 0: + logger.info(status_msg) + if rank == 0: + logger.info(status_msg) + + def my_write_viz(step, t, state, dv, exact=None, resid=None): + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, time=t) + if resid is None: + resid = state - exact + viz_fields = [("cv", state), ("dv", dv)] + from mirgecom.simutil import write_visfile + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, + step=step, t=t, overwrite=True, vis_timer=vis_timer, + comm=comm) + + def my_write_restart(step, t, state, tseed): + rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "cv": state, + "temperature_seed": tseed, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": nparts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(dv, component_errors): + health_error = False + from mirgecom.simutil import check_naninf_local, check_range_local + if check_naninf_local(dcoll, "vol", dv.pressure) \ + or check_range_local(dcoll, "vol", dv.pressure, 1e5, 1.1e5): + health_error = True + logger.info(f"{rank=}: Invalid pressure data found.") + + exittol = .09 + if max(component_errors) > exittol: + health_error = True + if rank == 0: + logger.info("Solution diverged from exact soln.") + + return health_error + + def my_pre_step(step, t, dt, state): + cv, tseed = state + fluid_state = construct_fluid_state(cv, tseed) + dv = fluid_state.dv + + try: + exact = None + component_errors = None + + if logmgr: + logmgr.tick_before() + + from mirgecom.simutil import check_step + do_viz = check_step(step=step, interval=nviz) + do_restart = check_step(step=step, interval=nrestart) + do_health = check_step(step=step, interval=nhealth) + do_status = check_step(step=step, interval=nstatus) + + if do_health: + exact = initializer(x_vec=nodes, eos=eos, time=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(dcoll, cv, exact) + health_errors = global_reduce( + my_health_check(dv, component_errors), op="lor") + if health_errors: + if rank == 0: + logger.info("Fluid solution failed health check.") + raise MyRuntimeError("Failed simulation health check.") + + if do_restart: + my_write_restart(step=step, t=t, state=cv, tseed=tseed) + + if do_viz: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, time=t) + resid = state - exact + my_write_viz(step=step, t=t, state=cv, dv=dv, exact=exact, + resid=resid) + + if do_status: + if component_errors is None: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, time=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(dcoll, cv, exact) + my_write_status(component_errors, dv=dv) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, state=cv, dv=dv) + my_write_restart(step=step, t=t, state=cv, tseed=tseed) + raise + + dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, t_final, + constant_cfl) + return state, dt + + def my_post_step(step, t, dt, state): + cv, tseed = state + fluid_state = construct_fluid_state(cv, tseed) + tseed = fluid_state.temperature + # Logmgr needs to know about EOS, dt, dim? + # imo this is a design/scope flaw + if logmgr: + set_dt(logmgr, dt) + set_sim_state(logmgr, dim, cv, eos) + logmgr.tick_after() + return make_obj_array([fluid_state.cv, tseed]), dt + + def my_rhs(t, state): + cv, tseed = state + fluid_state = make_fluid_state(cv, gas_model, temperature_seed=tseed) + return make_obj_array( + [euler_operator(dcoll, state=fluid_state, time=t, + boundaries=boundaries, gas_model=gas_model, + quadrature_tag=quadrature_tag, use_esdg=use_esdg), + 0*tseed]) + + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + + current_step, current_t, advanced_state = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, dt=current_dt, + state=make_obj_array([current_state.cv, + current_state.temperature]), + t=current_t, t_final=t_final) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + + current_cv, tseed = advanced_state + current_state = make_fluid_state(current_cv, gas_model, temperature_seed=tseed) + final_dv = current_state.dv + final_exact = initializer(x_vec=nodes, eos=eos, time=current_t) + final_resid = current_state.cv - final_exact + my_write_viz(step=current_step, t=current_t, state=current_cv, dv=final_dv, + exact=final_exact, resid=final_resid) + my_write_restart(step=current_step, t=current_t, state=current_state.cv, + tseed=tseed) + + if logmgr: + logmgr.close() + elif use_profiling: + print(actx.tabulate_profiling_data()) + + finish_tol = 1e-16 + assert np.abs(current_t - t_final) < finish_tol + + +if __name__ == "__main__": + import argparse + casename = "uiuc-mixture" + parser = argparse.ArgumentParser(description=f"MIRGE-Com Example: {casename}") + parser.add_argument("--lazy", action="store_true", + help="switch to a lazy computation mode") + parser.add_argument("--profiling", action="store_true", + help="turn on detailed performance profiling") + parser.add_argument("--overintegration", action="store_true", + help="use overintegration in the RHS computations") + parser.add_argument("--esdg", action="store_true", + help="use flux-differencing/entropy stable DG for inviscid computations.") + parser.add_argument("--leap", action="store_true", + help="use leap timestepper") + parser.add_argument("--numpy", action="store_true", + help="use numpy-based eager actx.") + parser.add_argument("--restart_file", help="root name of restart file") + parser.add_argument("--casename", help="casename to use for i/o") + args = parser.parse_args() + from warnings import warn + warn("Automatically turning off DV logging. MIRGE-Com Issue(578)") + log_dependent = False + + from mirgecom.simutil import ApplicationOptionsError + if args.esdg: + if not args.lazy and not args.numpy: + raise ApplicationOptionsError("ESDG requires lazy or numpy context.") + if not args.overintegration: + warn("ESDG requires overintegration, enabling --overintegration.") + + from mirgecom.array_context import get_reasonable_array_context_class + actx_class = get_reasonable_array_context_class( + lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) + + logging.basicConfig(format="%(message)s", level=logging.INFO) + if args.casename: + casename = args.casename + rst_filename = None + + if args.restart_file: + rst_filename = args.restart_file + + main(actx_class, use_leap=args.leap, + casename=casename, rst_filename=rst_filename, + use_overintegration=args.overintegration or args.esdg, + use_esdg=args.esdg, log_dependent=log_dependent) + +# vim: foldmethod=marker From b2a23bd703f19073e48ec866b26e69a70dfdb158 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 07:41:27 -0600 Subject: [PATCH 14/18] Undo README change --- examples/README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/examples/README.md b/examples/README.md index 96dfeed76..e5ab0b4cf 100644 --- a/examples/README.md +++ b/examples/README.md @@ -7,6 +7,7 @@ The examples and the unique features they exercise are as follows: - `autoignition.py`: Chemistry verification case with Pyrometheus - `heat-source.py`: Diffusion operator - `lump.py`: Lump advection, advection verification case +- `mixture.py`: Mixture EOS with Pyrometheus - `scalar-lump.py`: Scalar component lump advection verification case - `pulse.py`: Acoustic pulse in a box, wall boundary test case - `sod.py`: Sod's shock case: Fluid test case with strong shock From 8fa338cd2477cbcd45faef6371a01d57d933e804 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 08:02:29 -0600 Subject: [PATCH 15/18] typo --- examples/autoignition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index ed4a434f1..d5e6ba926 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -261,7 +261,7 @@ def main(actx_class, use_leap=False, use_overintegration=False, get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) eos = PyrometheusMixture(pyro_mechanism, temperature_guess=temperature_seed) - # print out the mechanism for inspecion + # print out the mechanism for inspection # import pyrometheus as pyro # with open(f"mechanism.py", "w") as mech_file: # code = pyro.codegen.python.gen_thermochem_code(cantera_soln) From 825b7d0362da022b3fc584238018d893c5e2c7e9 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 08:05:10 -0600 Subject: [PATCH 16/18] typo --- examples/autoignition.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index d5e6ba926..dca286789 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -622,7 +622,7 @@ def my_rhs(t, state): my_write_viz(step=current_step, t=current_t, dt=dt, state=final_cv, dv=final_dv, production_rates=final_dm, heat_release_rate=final_heat_rls, ts_field=ts_field, cfl=cfl) - my_write_status(t=current_dt, dt=dt, cfl=cfl, dv=final_dv) + my_write_status(t=current_t, dt=dt, cfl=cfl, dv=final_dv) my_write_restart(step=current_step, t=current_t, state=final_fluid_state, temperature_seed=tseed) From 359215df3cc345dbde6105ee281d14739ab42b86 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 13:11:09 -0600 Subject: [PATCH 17/18] Remove entropy check in --- test/test_eos.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_eos.py b/test/test_eos.py index 2d59a1ea5..bbb9ac51c 100644 --- a/test/test_eos.py +++ b/test/test_eos.py @@ -287,7 +287,6 @@ def inf_norm(x): can_p = cantera_soln.P can_e = cantera_soln.int_energy_mass can_h = cantera_soln.enthalpy_mass - # can_s = cantera_soln.entropy_mass can_c = cantera_soln.concentrations ones = dcoll.zeros(actx) + 1.0 @@ -300,7 +299,6 @@ def inf_norm(x): pyro_t = pyro_mechanism.get_temperature(pyro_e, tin, yin) pyro_p = pyro_mechanism.get_pressure(pyro_m, pyro_t, yin) pyro_h = pyro_mechanism.get_mixture_enthalpy_mass(pyro_t, yin) - # pyro_s = pyro_mechanism.get_mixture_entropy_mass(pyro_p, pyro_t, yin) pyro_c = pyro_mechanism.get_concentrations(pyro_m, yin) assert inf_norm((pyro_c - can_c)) < 1e-14 @@ -309,7 +307,6 @@ def inf_norm(x): assert inf_norm((pyro_p - can_p) / can_p) < 1e-14 assert inf_norm((pyro_e - can_e) / can_e) < 1e-10 assert inf_norm((pyro_h - can_h) / can_h) < 1e-10 - # assert inf_norm((pyro_s - can_s) / can_s) < 1e-10 # Test the concentrations zero level y = -ones*y0s From e84d123c40a4dca273e7df367291aa1c61a69a6c Mon Sep 17 00:00:00 2001 From: Tulio Date: Tue, 6 Feb 2024 10:50:45 -0600 Subject: [PATCH 18/18] Remove tseed_rhs --- examples/autoignition.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/examples/autoignition.py b/examples/autoignition.py index dca286789..925e2a672 100644 --- a/examples/autoignition.py +++ b/examples/autoignition.py @@ -355,8 +355,6 @@ def _limit_fluid_cv(cv, pressure, temperature, dd=None): momentum=mass_lim*cv.velocity, species_mass=mass_lim*spec_lim) - # limit_fluid_cv = actx.compile(_limit_fluid_cv) - def get_temperature_update(cv, temperature): y = cv.species_mass_fractions e = gas_model.eos.internal_energy(cv) / cv.mass @@ -595,8 +593,7 @@ def my_rhs(t, state): limiter_func=_limit_fluid_cv) chem_rhs = eos.get_species_source_terms(cv, fluid_state.temperature) - tseed_rhs = fluid_state.temperature - tseed - return make_obj_array([chem_rhs, tseed_rhs]) + return make_obj_array([chem_rhs, fluid_state.temperature*0.0]) current_dt = get_sim_timestep(dcoll, current_fluid_state, current_t, current_dt, current_cfl, t_final, constant_cfl)