From 8a2ce4c580a46d9386dc6848383813253d4caf69 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 11 Jan 2024 10:42:24 -0600 Subject: [PATCH 01/16] Add files --- examples/flame1d/flame1d.py | 900 +++++++++++++++++++++++++++++++++++ examples/flame1d/x_025um.dat | 715 ++++++++++++++++++++++++++++ examples/flame1d/y_025um.dat | 4 + 3 files changed, 1619 insertions(+) create mode 100644 examples/flame1d/flame1d.py create mode 100644 examples/flame1d/x_025um.dat create mode 100644 examples/flame1d/y_025um.dat diff --git a/examples/flame1d/flame1d.py b/examples/flame1d/flame1d.py new file mode 100644 index 000000000..8f0b28c2f --- /dev/null +++ b/examples/flame1d/flame1d.py @@ -0,0 +1,900 @@ +"""mirgecom driver for the 1D flame demonstration.""" + +__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 yaml +import logging +import sys +import numpy as np +import pyopencl as cl +import pyopencl.tools as cl_tools +import cantera +from functools import partial +from pytools.obj_array import make_obj_array + +from arraycontext import thaw, freeze + +from grudge.dof_desc import DOFDesc, BoundaryDomainTag +from grudge.dof_desc import DD_VOLUME_ALL +from grudge.shortcuts import compiled_lsrk45_step +from grudge.shortcuts import make_visualizer +from grudge import op + +from logpyle import IntervalTimer, set_dt +from mirgecom.profiling import PyOpenCLProfilingArrayContext +from mirgecom.navierstokes import ns_operator, grad_cv_operator, grad_t_operator +from mirgecom.simutil import ( + check_step, + get_sim_timestep, + generate_and_distribute_mesh, + write_visfile, + check_naninf_local, + check_range_local, + global_reduce +) +from mirgecom.utils import force_evaluation +from mirgecom.restart import write_restart_file +from mirgecom.io import make_init_message +from mirgecom.mpi import mpi_entry_point +from mirgecom.steppers import advance_state +from mirgecom.fluid import make_conserved +from mirgecom.eos import PyrometheusMixture +from mirgecom.gas_model import GasModel, make_fluid_state, make_operator_fluid_states +from mirgecom.logging_quantities import ( + initialize_logmgr, logmgr_add_cl_device_info, logmgr_set_time, + logmgr_add_device_memory_usage) + +####################################################################################### + +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): + pass + + +class _FluidOpStatesTag: + pass + + +class _FluidGradCVTag: + pass + + +class _FluidGradTempTag: + pass + + +def sponge_func(cv, cv_ref, sigma): + return sigma*(cv_ref - cv) + + +class InitSponge: + r""" + .. automethod:: __init__ + .. automethod:: __call__ + """ + def __init__(self, *, x_min=None, x_max=None, y_min=None, y_max=None, x_thickness=None, y_thickness=None, amplitude): + r"""Initialize the sponge parameters. + Parameters + ---------- + x0: float + sponge starting x location + thickness: float + sponge extent + amplitude: float + sponge strength modifier + """ + + self._x_min = x_min + self._x_max = x_max + self._y_min = y_min + self._y_max = y_max + self._x_thickness = x_thickness + self._y_thickness = y_thickness + self._amplitude = amplitude + + def __call__(self, x_vec, *, time=0.0): + """Create the sponge intensity at locations *x_vec*. + Parameters + ---------- + x_vec: numpy.ndarray + Coordinates at which solution is desired + time: float + Time at which solution is desired. The strength is (optionally) + dependent on time + """ + xpos = x_vec[0] + ypos = x_vec[1] + actx = xpos.array_context + zeros = 0*xpos + + sponge = xpos*0.0 + + if (self._x_max is not None): + x0 = (self._x_max - self._x_thickness) + dx = +((xpos - x0)/self._x_thickness) + sponge = sponge + self._amplitude * actx.np.where( + actx.np.greater(xpos, x0), + actx.np.where(actx.np.greater(xpos, self._x_max), + 1.0, 3.0*dx**2 - 2.0*dx**3), + 0.0 + ) + + if (self._x_min is not None): + x0 = (self._x_min + self._x_thickness) + dx = -((xpos - x0)/self._x_thickness) + sponge = sponge + self._amplitude * actx.np.where( + actx.np.less(xpos, x0), + actx.np.where(actx.np.less(xpos, self._x_min), + 1.0, 3.0*dx**2 - 2.0*dx**3), + 0.0 + ) + + return sponge + + +@mpi_entry_point +def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, + use_leap=False, use_profiling=False, casename=None, lazy=False, + restart_file=None, use_overintegration=False): + + from mpi4py import MPI + comm = MPI.COMM_WORLD + rank = 0 + rank = comm.Get_rank() + nparts = comm.Get_size() + + logmgr = initialize_logmgr(use_logmgr, filename=(f"{casename}.sqlite"), + mode="wo", 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) + + # ~~~~~~~~~~~~~~~~~~ + + import time + t_start = time.time() + t_shutdown = 720*60 + + # ~~~~~~~~~~~~~~~~~~ + + rst_path = "restart_data/" + viz_path = "viz_data/" + vizname = viz_path+casename + rst_pattern = rst_path+"{cname}-{step:06d}-{rank:04d}.pkl" + + # default i/o frequencies + ngarbage = 10 + nviz = 25000 + nrestart = 25000 + nhealth = 1 + nstatus = 100 + + order = 4 + my_mechanism = "uiuc_7sp" + cantera_file = "adiabatic_flame_uiuc_7sp_phi1.00_p1.00_E:H0.0_Y.csv" + original_rst_file = "flame_1D_025um_C2H4_p=4_uiuc_7sp-000000-0000.pkl" +# transport = "power-law" +# transport = "mix-lewis" + transport = "mix" + + # default timestepping control + integrator = "compiled_lsrk45" + current_dt = 2.0e-9 + t_final = 1.0 + + niter = 2000000 + +###################################################### + + local_dt = False + constant_cfl = True + current_cfl = 0.4 + + dim = 2 + +############################################################################## + + def _compiled_stepper_wrapper(state, t, dt, rhs): + return compiled_lsrk45_step(actx, state, t, dt, rhs) + + force_eval_stepper = True + if integrator == "compiled_lsrk45": + timestepper = _compiled_stepper_wrapper + force_eval_stepper = False + + if rank == 0: + print("\n#### Simulation control data: ####") + print(f"\tnviz = {nviz}") + print(f"\tnrestart = {nrestart}") + print(f"\tnhealth = {nhealth}") + print(f"\tnstatus = {nstatus}") + if (constant_cfl == False): + print(f"\tcurrent_dt = {current_dt}") + print(f"\tt_final = {t_final}") + else: + print(f"\tconstant_cfl = {constant_cfl}") + print(f"\tcurrent_cfl = {current_cfl}") + print(f"\tniter = {niter}") + print(f"\torder = {order}") + print(f"\tTime integration = {integrator}") + +############################################################################## + + restart_step = None + if restart_file is None: + + xx = np.loadtxt('x_025um.dat') + yy = np.loadtxt('y_025um.dat') + + coords = tuple((xx,yy)) + + from meshmode.mesh.generation import generate_box_mesh + generate_mesh = partial(generate_box_mesh, + axis_coords=coords, + periodic=(False, True), + boundary_tag_to_face={ + "inlet": ["-x"], + "outlet": ["+x"]}) + + local_mesh, global_nelements = ( + generate_and_distribute_mesh(comm, generate_mesh)) + local_nelements = local_mesh.nelements + + else: # Restart + from mirgecom.restart import read_restart_data + restart_data = read_restart_data(actx, restart_file) + restart_step = restart_data["step"] + local_mesh = restart_data["local_mesh"] + local_nelements = local_mesh.nelements + global_nelements = restart_data["global_nelements"] + restart_order = int(restart_data["order"]) + + assert comm.Get_size() == restart_data["num_parts"] + + from mirgecom.discretization import create_discretization_collection + dcoll = create_discretization_collection(actx, local_mesh, order, + mpi_communicator=comm) + + nodes = actx.thaw(dcoll.nodes()) + + dd_vol = DD_VOLUME_ALL + + zeros = nodes[0]*0.0 + ones = nodes[0]*0.0 + 1.0 + + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + quadrature_tag = DISCR_TAG_QUAD if use_overintegration else DISCR_TAG_BASE + +########################################################################## + + # {{{ Set up initial state using Cantera + + # Use Cantera for initialization + print("Using Cantera", cantera.__version__) + + # import os + # current_path = os.path.abspath(os.getcwd()) + "/" + # mechanism_file = current_path + my_mechanism + mechanism_file = "uiuc_7sp" + + from mirgecom.mechanisms import get_mechanism_input + mech_input = get_mechanism_input(mechanism_file) + + cantera_soln = cantera.Solution(name="gas", yaml=mech_input) + nspecies = cantera_soln.n_species + + cantera_soln.set_equivalence_ratio(phi=1.0, fuel="C2H4:1,H2:1", + oxidizer={"O2": 1.0, "N2": 3.76}) + + # Initial temperature, pressure, and mixture mole fractions are needed to + # set up the initial state in Cantera. + if False: + cantera_data = np.genfromtxt(cantera_file, skip_header=1, delimiter=',') + + temp_unburned = cantera_data[0,2] + temp_burned = cantera_data[-1,2] + + y_unburned = cantera_data[0,4:] + y_burned = cantera_data[-1,4:] + + vel_unburned = cantera_data[0,1] + vel_burned = cantera_data[-1,1] + + mass_unburned = cantera_data[0,3] + mass_burned = cantera_data[-1,3] + + if True: + temp_unburned = 300.0 + + cantera_soln.TP = temp_unburned, 101325.0 + y_unburned = cantera_soln.Y + mass_unburned = cantera_soln.density + + cantera_soln.equilibrate("HP") + temp_burned = cantera_soln.T + y_burned = cantera_soln.Y + mass_burned = cantera_soln.density + + vel_unburned = 0.5 + vel_burned = vel_unburned*mass_unburned/mass_burned + + # Uncomment next line to make pylint fail when it can't find cantera.one_atm + one_atm = cantera.one_atm + pres_unburned = one_atm + pres_burned = one_atm + + # }}} + + # {{{ Create Pyrometheus thermochemistry object & EOS + + # Import Pyrometheus EOS + from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera + pyrometheus_mechanism = \ + get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) + eos = PyrometheusMixture(pyrometheus_mechanism, temperature_guess=1350.0) + + species_names = pyrometheus_mechanism.species_names + + # }}} + + if transport == "power-law": + from mirgecom.transport import PowerLawTransport + transport_model = PowerLawTransport(lewis=np.ones((nspecies,)), beta=4.093e-7) + + if transport == "mix-lewis": + from mirgecom.transport import MixtureAveragedTransport + transport_model = MixtureAveragedTransport(pyrometheus_mechanism, + lewis=np.ones(nspecies,)) + if transport == "mix": + from mirgecom.transport import MixtureAveragedTransport + transport_model = MixtureAveragedTransport(pyrometheus_mechanism) + + gas_model = GasModel(eos=eos, transport=transport_model) + + print(f"Pyrometheus mechanism species names {species_names}") + print(f"Unburned:") + print(f"T = {temp_unburned}") + print(f"D = {mass_unburned}") + print(f"Y = {y_unburned}") + print(f"U = {vel_unburned}") + print(f" ") + print(f"Burned:") + print(f"T = {temp_burned}") + print(f"D = {mass_burned}") + print(f"Y = {y_burned}") + print(f"U = {vel_burned}") + print(f" ") + +############################################################################### + + from mirgecom.limiter import bound_preserving_limiter + + def _limit_fluid_cv(cv, pressure, temperature, dd=None): + + # limit species + spec_lim = make_obj_array([ + bound_preserving_limiter(dcoll, cv.species_mass_fractions[i], + mmin=0.0, mmax=1.0, modify_average=True, dd=dd) + for i in range(nspecies)]) + + # normalize to ensure sum_Yi = 1.0 + aux = cv.mass*0.0 + for i in range(0, nspecies): + aux = aux + spec_lim[i] + spec_lim = spec_lim/aux + + # recompute density + mass_lim = eos.get_density(pressure=pressure, + temperature=temperature, species_mass_fractions=spec_lim) + + # recompute energy + energy_lim = mass_lim*( + gas_model.eos.get_internal_energy(temperature, + species_mass_fractions=spec_lim) + + 0.5*np.dot(cv.velocity, cv.velocity)) + + # make a new CV with the limited variables + return make_conserved(dim=dim, mass=mass_lim, energy=energy_lim, + momentum=mass_lim*cv.velocity, species_mass=mass_lim*spec_lim) + + def _get_fluid_state(cv, temp_seed): + return make_fluid_state(cv=cv, gas_model=gas_model, + temperature_seed=temp_seed, + limiter_func=_limit_fluid_cv + ) + + get_fluid_state = actx.compile(_get_fluid_state) + +############################################################################## + + logmgr = initialize_logmgr(use_logmgr, filename=(f"{casename}.sqlite"), + mode="wo", mpi_comm=comm) + + vis_timer = None + if logmgr: + logmgr_add_cl_device_info(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("dt.max", "dt: {value:1.6e} s, "), + ("t_sim.max", "sim time: {value:1.6e} s, "), + ("t_step.max", "------- step walltime: {value:6g} s\n") + ]) + + try: + logmgr.add_watches(["memory_usage_python.max", + "memory_usage_gpu.max"]) + except KeyError: + pass + + if use_profiling: + logmgr.add_watches(["pyopencl_array_time.max"]) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + gc_timer = IntervalTimer("t_gc", "Time spent garbage collecting") + logmgr.add_quantity(gc_timer) + +############################################################################## + + if restart_file is None: + + current_t = 0.0 + current_step = 0 + first_step = 0 + + flame_start_loc = 0.0 + + tseed = 1234.56789 + + from mirgecom.initializers import PlanarDiscontinuity + # use the burned conditions with a lower temperature + bulk_init = PlanarDiscontinuity(dim=dim, disc_location=flame_start_loc, + sigma=0.0005, nspecies=nspecies, + temperature_right=temp_burned, temperature_left=temp_unburned, + pressure_right=pres_burned, pressure_left=pres_unburned, + velocity_right=make_obj_array([vel_burned, 0.0]), + velocity_left=make_obj_array([vel_unburned, 0.0]), + species_mass_right=y_burned, species_mass_left=y_unburned) + + if rank == 0: + logging.info("Initializing soln.") + # for Discontinuity initial conditions + current_cv = bulk_init(x_vec=nodes, eos=eos, time=0.) + + else: + + if local_dt: + current_t = restart_data["step"] + else: + current_t = restart_data["t"] + current_step = restart_step + first_step = restart_step + 0 + + if restart_order != order: + print('Wrong order...') + sys.exit() + else: + current_cv = restart_data["cv"] + tseed = restart_data["temperature_seed"] + + if logmgr: + logmgr_set_time(logmgr, current_step, current_t) + + current_state = get_fluid_state(current_cv, tseed) + current_state = force_evaluation(actx, current_state) + +############################################################################## + + # initialize the sponge field + sponge_x_thickness = 0.065 + sponge_amp = 10000.0 + + xMaxLoc = +0.100 + xMinLoc = -0.100 + + sponge_init = InitSponge(x_max=xMaxLoc, + x_min=xMinLoc, + x_thickness=sponge_x_thickness, + amplitude=sponge_amp) + + sponge_sigma = sponge_init(x_vec=nodes) + + if False: + cantera_data = read_restart_data(actx, original_rst_file) + ref_cv = force_evaluation(actx, cantera_data["cv"]) + if True: + ref_cv = force_evaluation(actx, bulk_init(x_vec=nodes, eos=eos, time=0.)) + +################################################################# + + from mirgecom.boundary import ( + PrescribedFluidBoundary,PressureOutflowBoundary, + LinearizedOutflow2DBoundary, LinearizedInflow2DBoundary) + + inflow_cv_cond = op.project(dcoll, dd_vol, dd_vol.trace("inlet"), ref_cv) + def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): + return make_fluid_state(cv=inflow_cv_cond, gas_model=gas_model, + temperature_seed=300.0) + + inflow_bnd = PrescribedFluidBoundary( + boundary_state_func=inlet_bnd_state_func) + +# inflow_bnd = LinearizedInflow2DBoundary( +# free_stream_density=mass_unburned, free_stream_pressure=101325.0, +# free_stream_velocity=make_obj_array([vel_unburned, 0.0]), +# free_stream_species_mass_fractions=y_unburned) + + outflow_bnd = LinearizedOutflow2DBoundary( + free_stream_density=mass_burned, free_stream_pressure=101325.0, + free_stream_velocity=make_obj_array([vel_burned, 0.0]), + free_stream_species_mass_fractions=y_burned) +# outflow_bnd = PressureOutflowBoundary(boundary_pressure=101325.0) + + boundaries = {BoundaryDomainTag("inlet"): inflow_bnd, + BoundaryDomainTag("outlet"): outflow_bnd} + +#################################################################################### + + visualizer = make_visualizer(dcoll) + + initname = "flame1D" + eosname = eos.__class__.__name__ + init_message = make_init_message(dim=dim, order=order, + nelements=local_nelements, global_nelements=global_nelements, + t_initial=current_t, 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 get_production_rates(cv, temperature): +# return make_obj_array([eos.get_production_rates(cv, temperature)]) +# compute_production_rates = actx.compile(get_production_rates) + + def my_write_viz(step, t, dt, state): + + y = state.cv.species_mass_fractions + gas_const = gas_model.eos.gas_const(cv=state.cv) + + # reaction_rates, = compute_production_rates(state.cv, state.temperature) + viz_fields = [("CV_rho", state.cv.mass), + ("CV_rhoU", state.cv.momentum), + ("CV_rhoE", state.cv.energy), + ("DV_P", state.pressure), + ("DV_T", state.temperature), + ("DV_U", state.velocity[0]), + ("DV_V", state.velocity[1]), + # ("reaction_rates", reaction_rates), + ("sponge", sponge_sigma), + ("R", gas_const), + ("gamma", eos.gamma(state.cv, state.temperature)), + ("dt", dt), + ("mu", state.tv.viscosity), + ("kappa", state.tv.thermal_conductivity), + ] + + # species mass fractions + viz_fields.extend(("Y_"+species_names[i], y[i]) for i in range(nspecies)) + + # species diffusivity + viz_fields.extend( + ("diff_"+species_names[i], state.tv.species_diffusivity[i]) + for i in range(nspecies)) + + print('Writing solution file...') + write_visfile(dcoll, viz_fields, visualizer, vizname=vizname, + step=step, t=t, overwrite=True, comm=comm) + + def my_write_restart(step, t, cv, tseed): + restart_fname = rst_pattern.format(cname=casename, step=step, rank=rank) + if restart_fname != restart_file: + rst_data = { + "local_mesh": local_mesh, + "cv": cv, + "temperature_seed": tseed, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": nparts + } + + write_restart_file(actx, rst_data, restart_fname, comm) + +################################################################## + + def my_health_check(cv, dv): + health_error = False + pressure = force_evaluation(actx, dv.pressure) + temperature = force_evaluation(actx, dv.temperature) + + if check_naninf_local(dcoll, "vol", pressure): + health_error = True + logger.info(f"{rank=}: NANs/Infs in pressure data.") + + if check_naninf_local(dcoll, "vol", temperature): + health_error = True + logger.info(f"{rank=}: NANs/Infs in temperature data.") + + return health_error + +############################################################################## + + def my_pre_step(step, t, dt, state): + + if logmgr: + logmgr.tick_before() + + cv, tseed = state + cv = force_evaluation(actx, cv) + tseed = force_evaluation(actx, tseed) + + fluid_state = get_fluid_state(cv, tseed) + + if local_dt: + t = force_evaluation(actx, t) + dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, + gas_model, constant_cfl=constant_cfl, local_dt=local_dt) + dt = force_evaluation(actx, dt) + #dt = force_evaluation(actx, actx.np.minimum(dt, current_dt)) + else: + if constant_cfl: + dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, + t_final, constant_cfl, local_dt) + + try: + 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_garbage = check_step(step=step, interval=ngarbage) + + t_elapsed = time.time() - t_start + if t_shutdown - t_elapsed < 300.0: + my_write_restart(step=step, t=t, cv=fluid_state.cv, + tseed=tseed) + sys.exit() + + if do_garbage: + with gc_timer.start_sub_timer(): + from warnings import warn + warn("Running gc.collect() to work around memory growth issue ") + import gc + gc.collect() + + if do_health: + health_errors = global_reduce( + my_health_check(fluid_state.cv, fluid_state.dv), 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, cv=fluid_state.cv, + tseed=tseed) + + if do_viz: + my_write_viz(step=step, t=t, dt=dt, state=fluid_state) + + except MyRuntimeError: + if rank == 0: + logger.info("Errors detected; attempting graceful exit.") + my_write_viz(step=step, t=t, dt=dt, state=fluid_state) + raise + + return make_obj_array([fluid_state.cv, fluid_state.temperature]), dt + + def my_rhs(t, state): + cv, tseed = state + + fluid_state = make_fluid_state(cv=cv, gas_model=gas_model, + temperature_seed=tseed, limiter_func=_limit_fluid_cv) + + operator_states_quad = make_operator_fluid_states( + dcoll, fluid_state, gas_model, boundaries, quadrature_tag, + comm_tag=_FluidOpStatesTag, limiter_func=_limit_fluid_cv) + + grad_cv = grad_cv_operator( + dcoll, gas_model, boundaries, fluid_state, time=t, + quadrature_tag=quadrature_tag, comm_tag=_FluidGradCVTag, + limiter_func=_limit_fluid_cv, + operator_states_quad=operator_states_quad) + + grad_t = grad_t_operator( + dcoll, gas_model, boundaries, fluid_state, time=t, + quadrature_tag=quadrature_tag, comm_tag=_FluidGradTempTag, + limiter_func=_limit_fluid_cv, + operator_states_quad=operator_states_quad) + + ns_rhs = ns_operator(dcoll, gas_model, fluid_state, boundaries, time=t, + quadrature_tag=quadrature_tag, grad_cv=grad_cv, grad_t=grad_t, + operator_states_quad=operator_states_quad) + + chem_rhs = eos.get_species_source_terms(fluid_state.cv, + fluid_state.temperature) + + #sponge_rhs = sponge_func(cv=fluid_state.cv, cv_ref=ref_cv, + # sigma=sponge_sigma) + rhs = ns_rhs + chem_rhs #+ sponge_rhs + + return make_obj_array([rhs, zeros]) + + def my_post_step(step, t, dt, state): + if step == first_step + 1: + with gc_timer.start_sub_timer(): + import gc + gc.collect() + # Freeze the objects that are still alive so they will not + # be considered in future gc collections. + logger.info("Freezing GC objects to reduce overhead of " + "future GC collections") + gc.freeze() + + min_dt = np.min(actx.to_numpy(dt)) if local_dt else dt + if logmgr: + set_dt(logmgr, min_dt) + logmgr.tick_after() + + return state, dt + +############################################################################## + + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) + current_dt = force_evaluation(actx, current_dt) + + if rank == 0: + logging.info("Stepping.") + + (current_step, current_t, stepper_state) = \ + advance_state(rhs=my_rhs, timestepper=timestepper, + pre_step_callback=my_pre_step, + post_step_callback=my_post_step, + istep=current_step, dt=current_dt, t=current_t, + t_final=t_final, max_steps=niter, local_dt=local_dt, + force_eval=force_eval_stepper, + state=make_obj_array([current_state.cv, tseed]), + # compile_rhs=False + ) + current_cv, tseed = stepper_state + current_state = make_fluid_state(current_cv, gas_model, tseed) + + # Dump the final data + if rank == 0: + logger.info("Checkpointing final state ...") + + my_write_viz(step=current_step, t=current_t, dt=current_dt, + state=current_state) + my_write_restart(step=current_step, t=current_t, cv=current_state.cv, + temperature_seed=tseed) + + if logmgr: + logmgr.close() + elif use_profiling: + print(actx.tabulate_profiling_data()) + + exit() + + +if __name__ == "__main__": + logging.basicConfig(format="%(message)s", level=logging.INFO) + + import argparse + parser = argparse.ArgumentParser(description="MIRGE-Com 1D Flame Driver") + parser.add_argument("-r", "--restart_file", type=ascii, + dest="restart_file", nargs="?", action="store", + help="simulation restart file") + parser.add_argument("-i", "--input_file", type=ascii, + dest="input_file", nargs="?", action="store", + help="simulation config file") + parser.add_argument("-c", "--casename", type=ascii, + dest="casename", nargs="?", action="store", + help="simulation case name") + parser.add_argument("--profiling", action="store_true", default=False, + help="enable kernel profiling [OFF]") + parser.add_argument("--log", action="store_true", default=True, + help="enable logging profiling [ON]") + parser.add_argument("--overintegration", action="store_true", default=False, + help="enable overintegration [OFF]") + parser.add_argument("--esdg", action="store_true", + help="use entropy stable DG for inviscid computations.") + parser.add_argument("--lazy", action="store_true", default=False, + help="enable lazy evaluation [OFF]") + parser.add_argument("--numpy", action="store_true", + help="use numpy-based eager actx.") + + args = parser.parse_args() + + # for writing output + casename = "flame1D" + if args.casename: + print(f"Custom casename {args.casename}") + casename = (args.casename).replace("'", "") + else: + print(f"Default casename {casename}") + + restart_file = None + if args.restart_file: + restart_file = (args.restart_file).replace("'", "") + print(f"Restarting from file: {restart_file}") + + input_file = None + if args.input_file: + input_file = (args.input_file).replace("'", "") + print(f"Reading user input from {args.input_file}") + else: + print("No user input file, using default values") + + print(f"Running {sys.argv[0]}\n") + + from warnings import warn + 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) + + main(actx_class, use_logmgr=args.log, casename=casename, + restart_file=restart_file, use_overintegration=args.overintegration) + +# vim: foldmethod=marker diff --git a/examples/flame1d/x_025um.dat b/examples/flame1d/x_025um.dat new file mode 100644 index 000000000..cc7ce824b --- /dev/null +++ b/examples/flame1d/x_025um.dat @@ -0,0 +1,715 @@ +-2.041730741244119907e-01 +-1.983270622123499416e-01 +-1.926513224919013390e-01 +-1.871408955788444461e-01 +-1.817909665370416317e-01 +-1.765968606712136624e-01 +-1.715540394422544745e-01 +-1.666580965015173887e-01 +-1.619047538406075959e-01 +-1.572898580533165225e-01 +-1.528093767064319930e-01 +-1.484593948162528376e-01 +-1.442361114277293799e-01 +-1.401358362932405843e-01 +-1.361549866481058246e-01 +-1.322900840800138311e-01 +-1.285377514896332563e-01 +-1.248947101397491904e-01 +-1.213577767903471910e-01 +-1.179238609171413682e-01 +-1.145899620111162920e-01 +-1.113531669567230098e-01 +-1.082106474864382717e-01 +-1.051596577094627977e-01 +-1.021975317123992311e-01 +-9.932168122981324676e-02 +-9.652959338264238687e-02 +-9.381882848247651008e-02 +-9.118701789979119510e-02 +-8.863186199427147693e-02 +-8.615112810541737665e-02 +-8.374264860167553171e-02 +-8.140431898639219344e-02 +-7.913409605893263754e-02 +-7.692999612936025911e-02 +-7.479009328511522503e-02 +-7.271251770817829807e-02 +-7.069545404124924493e-02 +-6.873713980151230363e-02 +-6.683586384060265229e-02 +-6.498996484942823337e-02 +-6.319782990654045085e-02 +-6.145789306878531727e-02 +-5.976863400300363310e-02 +-5.812857665758452280e-02 +-5.653628797271160328e-02 +-5.499037662817479050e-02 +-5.348949182765361243e-02 +-5.203232211840975313e-02 +-5.061759424535746538e-02 +-4.924407203851058346e-02 +-4.791055533283399698e-02 +-4.661587891955575741e-02 +-4.535891152802348780e-02 +-4.413855483721546014e-02 +-4.295374251604261451e-02 +-4.180343929160296179e-02 +-4.068664004457417382e-02 +-3.960236893095399335e-02 +-3.854967852938100126e-02 +-3.752764901329071895e-02 +-3.653538734718364545e-02 +-3.557202650630299223e-02 +-3.463672471904022387e-02 +-3.372866473140646859e-02 +-3.284705309292709763e-02 +-3.199111946333547257e-02 +-3.116011593946010558e-02 +-3.035331640171702941e-02 +-2.957001587963637226e-02 +-2.880952993586874400e-02 +-2.807119406813318382e-02 +-2.735436312858409783e-02 +-2.665841076008983873e-02 +-2.598272884893036358e-02 +-2.532672699343572684e-02 +-2.468983198810112759e-02 +-2.407148732272772901e-02 +-2.347115269615161412e-02 +-2.288830354413596765e-02 +-2.232243058101398073e-02 +-2.177303935468195492e-02 +-2.123964981455377446e-02 +-2.072179589209922901e-02 +-2.021902509359967126e-02 +-1.973089810476514813e-02 +-1.925698840686755273e-02 +-1.879688190405435361e-02 +-1.835017656151726601e-02 +-1.791648205419970469e-02 +-1.749541942573605222e-02 +-1.708662075732473792e-02 +-1.668972884624579309e-02 +-1.630439689374196463e-02 +-1.593028820199067497e-02 +-1.556707587990204499e-02 +-1.521444255748590109e-02 +-1.487208010853818818e-02 +-1.453968938140448586e-02 +-1.421697993758535769e-02 +-1.390366979795513644e-02 +-1.359948519637239721e-02 +-1.330416034046682551e-02 +-1.301743717939345471e-02 +-1.273906517835134665e-02 +-1.246880109966968901e-02 +-1.220640879027002149e-02 +-1.195165897531888857e-02 +-1.170432905789060493e-02 +-1.146420292446508754e-02 +-1.123107075609079859e-02 +-1.100472884504779918e-02 +-1.078497941685071203e-02 +-1.057163045743606498e-02 +-1.036449554538300920e-02 +-1.016339368902081930e-02 +-9.968149168280829309e-03 +-9.778591381154624895e-03 +-9.594554694624330207e-03 +-9.415878299934723680e-03 +-9.242406072080737098e-03 +-9.073986433387545500e-03 +-8.910472221064058562e-03 +-8.751720558614070958e-03 +-8.597592730992724072e-03 +-8.447954063399184046e-03 +-8.302673803599631280e-03 +-8.161625007677735660e-03 +-8.024684429112787648e-03 +-7.891732411088566226e-03 +-7.762652781938836691e-03 +-7.637332753638127747e-03 +-7.515662823249089769e-03 +-7.397536677240314738e-03 +-7.282851098591018937e-03 +-7.171505876601411145e-03 +-7.063403719329947386e-03 +-6.958450168580953589e-03 +-6.856553517368337881e-03 +-6.757624729783274012e-03 +-6.661577363195833774e-03 +-6.568327492722590699e-03 +-6.477793637894199187e-03 +-6.389896691458867538e-03 +-6.304559850259516818e-03 +-6.221708548124225101e-03 +-6.141270390711320824e-03 +-6.063175092252190770e-03 +-5.987354414136530911e-03 +-5.913742105286375582e-03 +-5.842273844266806897e-03 +-5.772887183082759874e-03 +-5.705521492612811510e-03 +-5.640117909632279269e-03 +-5.576619285379335288e-03 +-5.514970135619195676e-03 +-5.455116592162749696e-03 +-5.397006355797268583e-03 +-5.340588650588063736e-03 +-5.285814179511166139e-03 +-5.232635081378255934e-03 +-5.181004889016207343e-03 +-5.130878488664703495e-03 +-5.082212080556446944e-03 +-5.034963140645518608e-03 +-4.989090383450442250e-03 +-4.944553725979494371e-03 +-4.901314252706729171e-03 +-4.859334181568122392e-03 +-4.818576830948115633e-03 +-4.779006587627720494e-03 +-4.740588875666171959e-03 +-4.703290126188940101e-03 +-4.667077748055705431e-03 +-4.631920099382662130e-03 +-4.597786459894270643e-03 +-4.564646158254326638e-03 +-4.532464991306797870e-03 +-4.501202516133172814e-03 +-4.470810848816979786e-03 +-4.441234589365540379e-03 +-4.412411543636324089e-03 +-4.384273970275990792e-03 +-4.356750139496174540e-03 +-4.329766046085505721e-03 +-4.303247167161224594e-03 +-4.277120192851587704e-03 +-4.251314683807731754e-03 +-4.225764622822316199e-03 +-4.200409829529845422e-03 +-4.175197198664148783e-03 +-4.150081705905538793e-03 +-4.125027103969023355e-03 +-4.100006208953808906e-03 +-4.075000657206401776e-03 +-4.049999999999989372e-03 +-4.024999999999989524e-03 +-3.999999999999989675e-03 +-3.974999999999989826e-03 +-3.949999999999989977e-03 +-3.924999999999990129e-03 +-3.899999999999990280e-03 +-3.874999999999990431e-03 +-3.849999999999990583e-03 +-3.824999999999990734e-03 +-3.799999999999990885e-03 +-3.774999999999991036e-03 +-3.749999999999991188e-03 +-3.724999999999991339e-03 +-3.699999999999991490e-03 +-3.674999999999991641e-03 +-3.649999999999991793e-03 +-3.624999999999991944e-03 +-3.599999999999992095e-03 +-3.574999999999992246e-03 +-3.549999999999992398e-03 +-3.524999999999992549e-03 +-3.499999999999992700e-03 +-3.474999999999992852e-03 +-3.449999999999993003e-03 +-3.424999999999993154e-03 +-3.399999999999993305e-03 +-3.374999999999993457e-03 +-3.349999999999993608e-03 +-3.324999999999993759e-03 +-3.299999999999993910e-03 +-3.274999999999994062e-03 +-3.249999999999994213e-03 +-3.224999999999994364e-03 +-3.199999999999994515e-03 +-3.174999999999994667e-03 +-3.149999999999994818e-03 +-3.124999999999994969e-03 +-3.099999999999995121e-03 +-3.074999999999995272e-03 +-3.049999999999995423e-03 +-3.024999999999995574e-03 +-2.999999999999995726e-03 +-2.974999999999995877e-03 +-2.949999999999996028e-03 +-2.924999999999996179e-03 +-2.899999999999996331e-03 +-2.874999999999996482e-03 +-2.849999999999996633e-03 +-2.824999999999996785e-03 +-2.799999999999996936e-03 +-2.774999999999997087e-03 +-2.749999999999997238e-03 +-2.724999999999997390e-03 +-2.699999999999997541e-03 +-2.674999999999997692e-03 +-2.649999999999997843e-03 +-2.624999999999997995e-03 +-2.599999999999998146e-03 +-2.574999999999998297e-03 +-2.549999999999998448e-03 +-2.524999999999998600e-03 +-2.499999999999998751e-03 +-2.474999999999998902e-03 +-2.449999999999999054e-03 +-2.424999999999999205e-03 +-2.399999999999999356e-03 +-2.374999999999999507e-03 +-2.349999999999999659e-03 +-2.324999999999999810e-03 +-2.299999999999999961e-03 +-2.275000000000000112e-03 +-2.250000000000000264e-03 +-2.225000000000000415e-03 +-2.200000000000000566e-03 +-2.175000000000000717e-03 +-2.150000000000000869e-03 +-2.125000000000001020e-03 +-2.100000000000001171e-03 +-2.075000000000001323e-03 +-2.050000000000001474e-03 +-2.025000000000001625e-03 +-2.000000000000001776e-03 +-1.975000000000001928e-03 +-1.950000000000001862e-03 +-1.925000000000001796e-03 +-1.900000000000001731e-03 +-1.875000000000001665e-03 +-1.850000000000001600e-03 +-1.825000000000001534e-03 +-1.800000000000001469e-03 +-1.775000000000001403e-03 +-1.750000000000001337e-03 +-1.725000000000001272e-03 +-1.700000000000001206e-03 +-1.675000000000001141e-03 +-1.650000000000001075e-03 +-1.625000000000001010e-03 +-1.600000000000000944e-03 +-1.575000000000000878e-03 +-1.550000000000000813e-03 +-1.525000000000000747e-03 +-1.500000000000000682e-03 +-1.475000000000000616e-03 +-1.450000000000000551e-03 +-1.425000000000000485e-03 +-1.400000000000000419e-03 +-1.375000000000000354e-03 +-1.350000000000000288e-03 +-1.325000000000000223e-03 +-1.300000000000000157e-03 +-1.275000000000000092e-03 +-1.250000000000000026e-03 +-1.224999999999999960e-03 +-1.199999999999999895e-03 +-1.174999999999999829e-03 +-1.149999999999999764e-03 +-1.124999999999999698e-03 +-1.099999999999999633e-03 +-1.074999999999999567e-03 +-1.049999999999999501e-03 +-1.024999999999999436e-03 +-9.999999999999993703e-04 +-9.749999999999993047e-04 +-9.499999999999993476e-04 +-9.249999999999993904e-04 +-8.999999999999994333e-04 +-8.749999999999994761e-04 +-8.499999999999995190e-04 +-8.249999999999995618e-04 +-7.999999999999996047e-04 +-7.749999999999996475e-04 +-7.499999999999996904e-04 +-7.249999999999997332e-04 +-6.999999999999997760e-04 +-6.749999999999998189e-04 +-6.499999999999998617e-04 +-6.249999999999999046e-04 +-5.999999999999999474e-04 +-5.749999999999999903e-04 +-5.500000000000000331e-04 +-5.250000000000000760e-04 +-5.000000000000001188e-04 +-4.750000000000001617e-04 +-4.500000000000001503e-04 +-4.250000000000001390e-04 +-4.000000000000001276e-04 +-3.750000000000001162e-04 +-3.500000000000001049e-04 +-3.250000000000000935e-04 +-3.000000000000000821e-04 +-2.750000000000000708e-04 +-2.500000000000000594e-04 +-2.250000000000000481e-04 +-2.000000000000000367e-04 +-1.750000000000000253e-04 +-1.500000000000000140e-04 +-1.250000000000000026e-04 +-1.000000000000000048e-04 +-7.500000000000000698e-05 +-5.000000000000000240e-05 +-2.500000000000000120e-05 +-0.000000000000000000e+00 +2.500000000000000120e-05 +5.000000000000000240e-05 +7.500000000000000698e-05 +1.000000000000000048e-04 +1.250000000000000026e-04 +1.500000000000000140e-04 +1.750000000000000253e-04 +2.000000000000000367e-04 +2.250000000000000481e-04 +2.500000000000000594e-04 +2.750000000000000708e-04 +3.000000000000000821e-04 +3.250000000000000935e-04 +3.500000000000001049e-04 +3.750000000000001162e-04 +4.000000000000001276e-04 +4.250000000000001390e-04 +4.500000000000001503e-04 +4.750000000000001617e-04 +5.000000000000001188e-04 +5.250000000000000760e-04 +5.500000000000000331e-04 +5.749999999999999903e-04 +5.999999999999999474e-04 +6.249999999999999046e-04 +6.499999999999998617e-04 +6.749999999999998189e-04 +6.999999999999997760e-04 +7.249999999999997332e-04 +7.499999999999996904e-04 +7.749999999999996475e-04 +7.999999999999996047e-04 +8.249999999999995618e-04 +8.499999999999995190e-04 +8.749999999999994761e-04 +8.999999999999994333e-04 +9.249999999999993904e-04 +9.499999999999993476e-04 +9.749999999999993047e-04 +9.999999999999993703e-04 +1.024999999999999436e-03 +1.049999999999999501e-03 +1.074999999999999567e-03 +1.099999999999999633e-03 +1.124999999999999698e-03 +1.149999999999999764e-03 +1.174999999999999829e-03 +1.199999999999999895e-03 +1.224999999999999960e-03 +1.250000000000000026e-03 +1.275000000000000092e-03 +1.300000000000000157e-03 +1.325000000000000223e-03 +1.350000000000000288e-03 +1.375000000000000354e-03 +1.400000000000000419e-03 +1.425000000000000485e-03 +1.450000000000000551e-03 +1.475000000000000616e-03 +1.500000000000000682e-03 +1.525000000000000747e-03 +1.550000000000000813e-03 +1.575000000000000878e-03 +1.600000000000000944e-03 +1.625000000000001010e-03 +1.650000000000001075e-03 +1.675000000000001141e-03 +1.700000000000001206e-03 +1.725000000000001272e-03 +1.750000000000001337e-03 +1.775000000000001403e-03 +1.800000000000001469e-03 +1.825000000000001534e-03 +1.850000000000001600e-03 +1.875000000000001665e-03 +1.900000000000001731e-03 +1.925000000000001796e-03 +1.950000000000001862e-03 +1.975000000000001928e-03 +2.000000000000001776e-03 +2.025000000000001625e-03 +2.050000000000001474e-03 +2.075000000000001323e-03 +2.100000000000001171e-03 +2.125000000000001020e-03 +2.150000000000000869e-03 +2.175000000000000717e-03 +2.200000000000000566e-03 +2.225000000000000415e-03 +2.250000000000000264e-03 +2.275000000000000112e-03 +2.299999999999999961e-03 +2.324999999999999810e-03 +2.349999999999999659e-03 +2.374999999999999507e-03 +2.399999999999999356e-03 +2.424999999999999205e-03 +2.449999999999999054e-03 +2.474999999999998902e-03 +2.499999999999998751e-03 +2.524999999999998600e-03 +2.549999999999998448e-03 +2.574999999999998297e-03 +2.599999999999998146e-03 +2.624999999999997995e-03 +2.649999999999997843e-03 +2.674999999999997692e-03 +2.699999999999997541e-03 +2.724999999999997390e-03 +2.749999999999997238e-03 +2.774999999999997087e-03 +2.799999999999996936e-03 +2.824999999999996785e-03 +2.849999999999996633e-03 +2.874999999999996482e-03 +2.899999999999996331e-03 +2.924999999999996179e-03 +2.949999999999996028e-03 +2.974999999999995877e-03 +2.999999999999995726e-03 +3.024999999999995574e-03 +3.049999999999995423e-03 +3.074999999999995272e-03 +3.099999999999995121e-03 +3.124999999999994969e-03 +3.149999999999994818e-03 +3.174999999999994667e-03 +3.199999999999994515e-03 +3.224999999999994364e-03 +3.249999999999994213e-03 +3.274999999999994062e-03 +3.299999999999993910e-03 +3.324999999999993759e-03 +3.349999999999993608e-03 +3.374999999999993457e-03 +3.399999999999993305e-03 +3.424999999999993154e-03 +3.449999999999993003e-03 +3.474999999999992852e-03 +3.499999999999992700e-03 +3.524999999999992549e-03 +3.549999999999992398e-03 +3.574999999999992246e-03 +3.599999999999992095e-03 +3.624999999999991944e-03 +3.649999999999991793e-03 +3.674999999999991641e-03 +3.699999999999991490e-03 +3.724999999999991339e-03 +3.749999999999991188e-03 +3.774999999999991036e-03 +3.799999999999990885e-03 +3.824999999999990734e-03 +3.849999999999990583e-03 +3.874999999999990431e-03 +3.899999999999990280e-03 +3.924999999999990129e-03 +3.949999999999989977e-03 +3.974999999999989826e-03 +3.999999999999989675e-03 +4.024999999999989524e-03 +4.049999999999989372e-03 +4.075000657206401776e-03 +4.100006208953808906e-03 +4.125027103969023355e-03 +4.150081705905538793e-03 +4.175197198664148783e-03 +4.200409829529845422e-03 +4.225764622822316199e-03 +4.251314683807731754e-03 +4.277120192851587704e-03 +4.303247167161224594e-03 +4.329766046085505721e-03 +4.356750139496174540e-03 +4.384273970275990792e-03 +4.412411543636324089e-03 +4.441234589365540379e-03 +4.470810848816979786e-03 +4.501202516133172814e-03 +4.532464991306797870e-03 +4.564646158254326638e-03 +4.597786459894270643e-03 +4.631920099382662130e-03 +4.667077748055705431e-03 +4.703290126188940101e-03 +4.740588875666171959e-03 +4.779006587627720494e-03 +4.818576830948115633e-03 +4.859334181568122392e-03 +4.901314252706729171e-03 +4.944553725979494371e-03 +4.989090383450442250e-03 +5.034963140645518608e-03 +5.082212080556446944e-03 +5.130878488664703495e-03 +5.181004889016207343e-03 +5.232635081378255934e-03 +5.285814179511166139e-03 +5.340588650588063736e-03 +5.397006355797268583e-03 +5.455116592162749696e-03 +5.514970135619195676e-03 +5.576619285379335288e-03 +5.640117909632279269e-03 +5.705521492612811510e-03 +5.772887183082759874e-03 +5.842273844266806897e-03 +5.913742105286375582e-03 +5.987354414136530911e-03 +6.063175092252190770e-03 +6.141270390711320824e-03 +6.221708548124225101e-03 +6.304559850259516818e-03 +6.389896691458867538e-03 +6.477793637894199187e-03 +6.568327492722590699e-03 +6.661577363195833774e-03 +6.757624729783274012e-03 +6.856553517368337881e-03 +6.958450168580953589e-03 +7.063403719329947386e-03 +7.171505876601411145e-03 +7.282851098591018937e-03 +7.397536677240314738e-03 +7.515662823249089769e-03 +7.637332753638127747e-03 +7.762652781938836691e-03 +7.891732411088566226e-03 +8.024684429112787648e-03 +8.161625007677735660e-03 +8.302673803599631280e-03 +8.447954063399184046e-03 +8.597592730992724072e-03 +8.751720558614070958e-03 +8.910472221064058562e-03 +9.073986433387545500e-03 +9.242406072080737098e-03 +9.415878299934723680e-03 +9.594554694624330207e-03 +9.778591381154624895e-03 +9.968149168280829309e-03 +1.016339368902081930e-02 +1.036449554538300920e-02 +1.057163045743606498e-02 +1.078497941685071203e-02 +1.100472884504779918e-02 +1.123107075609079859e-02 +1.146420292446508754e-02 +1.170432905789060493e-02 +1.195165897531888857e-02 +1.220640879027002149e-02 +1.246880109966968901e-02 +1.273906517835134665e-02 +1.301743717939345471e-02 +1.330416034046682551e-02 +1.359948519637239721e-02 +1.390366979795513644e-02 +1.421697993758535769e-02 +1.453968938140448586e-02 +1.487208010853818818e-02 +1.521444255748590109e-02 +1.556707587990204499e-02 +1.593028820199067497e-02 +1.630439689374196463e-02 +1.668972884624579309e-02 +1.708662075732473792e-02 +1.749541942573605222e-02 +1.791648205419970469e-02 +1.835017656151726601e-02 +1.879688190405435361e-02 +1.925698840686755273e-02 +1.973089810476514813e-02 +2.021902509359967126e-02 +2.072179589209922901e-02 +2.123964981455377446e-02 +2.177303935468195492e-02 +2.232243058101398073e-02 +2.288830354413596765e-02 +2.347115269615161412e-02 +2.407148732272772901e-02 +2.468983198810112759e-02 +2.532672699343572684e-02 +2.598272884893036358e-02 +2.665841076008983873e-02 +2.735436312858409783e-02 +2.807119406813318382e-02 +2.880952993586874400e-02 +2.957001587963637226e-02 +3.035331640171702941e-02 +3.116011593946010558e-02 +3.199111946333547257e-02 +3.284705309292709763e-02 +3.372866473140646859e-02 +3.463672471904022387e-02 +3.557202650630299223e-02 +3.653538734718364545e-02 +3.752764901329071895e-02 +3.854967852938100126e-02 +3.960236893095399335e-02 +4.068664004457417382e-02 +4.180343929160296179e-02 +4.295374251604261451e-02 +4.413855483721546014e-02 +4.535891152802348780e-02 +4.661587891955575741e-02 +4.791055533283399698e-02 +4.924407203851058346e-02 +5.061759424535746538e-02 +5.203232211840975313e-02 +5.348949182765361243e-02 +5.499037662817479050e-02 +5.653628797271160328e-02 +5.812857665758452280e-02 +5.976863400300363310e-02 +6.145789306878531727e-02 +6.319782990654045085e-02 +6.498996484942823337e-02 +6.683586384060265229e-02 +6.873713980151230363e-02 +7.069545404124924493e-02 +7.271251770817829807e-02 +7.479009328511522503e-02 +7.692999612936025911e-02 +7.913409605893263754e-02 +8.140431898639219344e-02 +8.374264860167553171e-02 +8.615112810541737665e-02 +8.863186199427147693e-02 +9.118701789979119510e-02 +9.381882848247651008e-02 +9.652959338264238687e-02 +9.932168122981324676e-02 +1.021975317123992311e-01 +1.051596577094627977e-01 +1.082106474864382717e-01 +1.113531669567230098e-01 +1.145899620111162920e-01 +1.179238609171413682e-01 +1.213577767903471910e-01 +1.248947101397491904e-01 +1.285377514896332563e-01 +1.322900840800138311e-01 +1.361549866481058246e-01 +1.401358362932405843e-01 +1.442361114277293799e-01 +1.484593948162528376e-01 +1.528093767064319930e-01 +1.572898580533165225e-01 +1.619047538406075959e-01 +1.666580965015173887e-01 +1.715540394422544745e-01 +1.765968606712136624e-01 +1.817909665370416317e-01 +1.871408955788444461e-01 +1.926513224919013390e-01 +1.983270622123499416e-01 +2.041730741244119907e-01 diff --git a/examples/flame1d/y_025um.dat b/examples/flame1d/y_025um.dat new file mode 100644 index 000000000..28c76dc59 --- /dev/null +++ b/examples/flame1d/y_025um.dat @@ -0,0 +1,4 @@ +0.000000000000000000e+00 +2.500000000000000120e-05 +5.000000000000000240e-05 +7.500000000000000698e-05 From b096617eb51ecabcc1695a597252233b7ec6f027 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 07:03:23 -0600 Subject: [PATCH 02/16] flake8 and pylint --- examples/flame1d/flame1d.py | 374 +++++++++++++--------------- examples/mixture.py | 478 ------------------------------------ 2 files changed, 169 insertions(+), 683 deletions(-) delete mode 100644 examples/mixture.py diff --git a/examples/flame1d/flame1d.py b/examples/flame1d/flame1d.py index 8f0b28c2f..07dda71fc 100644 --- a/examples/flame1d/flame1d.py +++ b/examples/flame1d/flame1d.py @@ -23,26 +23,21 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ -import yaml import logging -import sys +import sys +import gc +from warnings import warn +from functools import partial import numpy as np -import pyopencl as cl -import pyopencl.tools as cl_tools import cantera -from functools import partial from pytools.obj_array import make_obj_array -from arraycontext import thaw, freeze - -from grudge.dof_desc import DOFDesc, BoundaryDomainTag -from grudge.dof_desc import DD_VOLUME_ALL +from grudge.dof_desc import BoundaryDomainTag, DD_VOLUME_ALL from grudge.shortcuts import compiled_lsrk45_step from grudge.shortcuts import make_visualizer from grudge import op from logpyle import IntervalTimer, set_dt -from mirgecom.profiling import PyOpenCLProfilingArrayContext from mirgecom.navierstokes import ns_operator, grad_cv_operator, grad_t_operator from mirgecom.simutil import ( check_step, @@ -54,18 +49,18 @@ global_reduce ) from mirgecom.utils import force_evaluation -from mirgecom.restart import write_restart_file +from mirgecom.restart import write_restart_file, read_restart_data from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point from mirgecom.steppers import advance_state from mirgecom.fluid import make_conserved from mirgecom.eos import PyrometheusMixture -from mirgecom.gas_model import GasModel, make_fluid_state, make_operator_fluid_states +from mirgecom.gas_model import ( + GasModel, make_fluid_state, make_operator_fluid_states) from mirgecom.logging_quantities import ( initialize_logmgr, logmgr_add_cl_device_info, logmgr_set_time, logmgr_add_device_memory_usage) -####################################################################################### class SingleLevelFilter(logging.Filter): def __init__(self, passlevel, reject): @@ -74,20 +69,19 @@ def __init__(self, passlevel, reject): def filter(self, record): if self.reject: - return (record.levelno != self.passlevel) - else: - return (record.levelno == self.passlevel) + return record.levelno != self.passlevel + 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) +# 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) @@ -110,6 +104,7 @@ class _FluidGradTempTag: def sponge_func(cv, cv_ref, sigma): + """Apply the sponge.""" return sigma*(cv_ref - cv) @@ -118,8 +113,9 @@ class InitSponge: .. automethod:: __init__ .. automethod:: __call__ """ - def __init__(self, *, x_min=None, x_max=None, y_min=None, y_max=None, x_thickness=None, y_thickness=None, amplitude): + def __init__(self, x_min, x_max, x_thickness, amplitude): r"""Initialize the sponge parameters. + Parameters ---------- x0: float @@ -132,56 +128,44 @@ def __init__(self, *, x_min=None, x_max=None, y_min=None, y_max=None, x_thicknes self._x_min = x_min self._x_max = x_max - self._y_min = y_min - self._y_max = y_max self._x_thickness = x_thickness - self._y_thickness = y_thickness self._amplitude = amplitude - def __call__(self, x_vec, *, time=0.0): + def __call__(self, x_vec): """Create the sponge intensity at locations *x_vec*. + Parameters ---------- x_vec: numpy.ndarray Coordinates at which solution is desired - time: float - Time at which solution is desired. The strength is (optionally) - dependent on time """ xpos = x_vec[0] - ypos = x_vec[1] actx = xpos.array_context - zeros = 0*xpos sponge = xpos*0.0 - if (self._x_max is not None): - x0 = (self._x_max - self._x_thickness) - dx = +((xpos - x0)/self._x_thickness) - sponge = sponge + self._amplitude * actx.np.where( - actx.np.greater(xpos, x0), - actx.np.where(actx.np.greater(xpos, self._x_max), - 1.0, 3.0*dx**2 - 2.0*dx**3), - 0.0 - ) - - if (self._x_min is not None): - x0 = (self._x_min + self._x_thickness) - dx = -((xpos - x0)/self._x_thickness) - sponge = sponge + self._amplitude * actx.np.where( - actx.np.less(xpos, x0), - actx.np.where(actx.np.less(xpos, self._x_min), - 1.0, 3.0*dx**2 - 2.0*dx**3), - 0.0 - ) + x0 = self._x_max - self._x_thickness + dx = +(xpos - x0)/self._x_thickness + sponge = sponge + self._amplitude * actx.np.where( + actx.np.greater(xpos, x0), + actx.np.where( + actx.np.greater(xpos, self._x_max), 1., 3.*dx**2 - 2.*dx**3), + 0.) + + x0 = self._x_min + self._x_thickness + dx = -(xpos - x0)/self._x_thickness + sponge = sponge + self._amplitude * actx.np.where( + actx.np.less(xpos, x0), + actx.np.where( + actx.np.less(xpos, self._x_min), 1., 3.*dx**2 - 2.*dx**3), + 0.) return sponge @mpi_entry_point -def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, - use_leap=False, use_profiling=False, casename=None, lazy=False, - restart_file=None, use_overintegration=False): +def main(actx_class, use_esdg=False, use_overintegration=False, + use_leap=False, casename=None, rst_filename=None): from mpi4py import MPI comm = MPI.COMM_WORLD @@ -189,8 +173,8 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, rank = comm.Get_rank() nparts = comm.Get_size() - logmgr = initialize_logmgr(use_logmgr, filename=(f"{casename}.sqlite"), - mode="wo", mpi_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) @@ -203,8 +187,6 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, t_start = time.time() t_shutdown = 720*60 - # ~~~~~~~~~~~~~~~~~~ - rst_path = "restart_data/" viz_path = "viz_data/" vizname = viz_path+casename @@ -213,16 +195,17 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, # default i/o frequencies ngarbage = 10 nviz = 25000 - nrestart = 25000 + nrestart = 25000 nhealth = 1 nstatus = 100 - order = 4 my_mechanism = "uiuc_7sp" - cantera_file = "adiabatic_flame_uiuc_7sp_phi1.00_p1.00_E:H0.0_Y.csv" - original_rst_file = "flame_1D_025um_C2H4_p=4_uiuc_7sp-000000-0000.pkl" -# transport = "power-law" -# transport = "mix-lewis" + + order = 4 + cantera_file = "adiabatic_flame_uiuc_7sp_phi1.00_p1.00_E:H0.0_Y.csv" + + # transport = "power-law" + # transport = "mix-lewis" transport = "mix" # default timestepping control @@ -232,7 +215,8 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, niter = 2000000 -###################################################### + use_sponge = False + reference_state = "flame_1D_025um_C2H4_p=4_uiuc_7sp-000000-0000.pkl" local_dt = False constant_cfl = True @@ -240,11 +224,11 @@ def main(actx_class, ctx_factory=cl.create_some_context, use_logmgr=True, dim = 2 -############################################################################## +# ############################################################################ def _compiled_stepper_wrapper(state, t, dt, rhs): return compiled_lsrk45_step(actx, state, t, dt, rhs) - + force_eval_stepper = True if integrator == "compiled_lsrk45": timestepper = _compiled_stepper_wrapper @@ -256,64 +240,61 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): print(f"\tnrestart = {nrestart}") print(f"\tnhealth = {nhealth}") print(f"\tnstatus = {nstatus}") - if (constant_cfl == False): - print(f"\tcurrent_dt = {current_dt}") - print(f"\tt_final = {t_final}") - else: + if constant_cfl: print(f"\tconstant_cfl = {constant_cfl}") print(f"\tcurrent_cfl = {current_cfl}") + else: + print(f"\tcurrent_dt = {current_dt}") + print(f"\tt_final = {t_final}") print(f"\tniter = {niter}") print(f"\torder = {order}") print(f"\tTime integration = {integrator}") -############################################################################## +# ############################################################################ - restart_step = None - if restart_file is None: + restart_step = 0 + if rst_filename is None: - xx = np.loadtxt('x_025um.dat') - yy = np.loadtxt('y_025um.dat') + xx = np.loadtxt("x_025um.dat") + yy = np.loadtxt("y_025um.dat") - coords = tuple((xx,yy)) + coords = tuple((xx, yy)) # noqa from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, axis_coords=coords, periodic=(False, True), - boundary_tag_to_face={ - "inlet": ["-x"], - "outlet": ["+x"]}) + boundary_tag_to_face={"inlet": ["-x"], + "outlet": ["+x"]}) local_mesh, global_nelements = ( generate_and_distribute_mesh(comm, generate_mesh)) local_nelements = local_mesh.nelements - else: # Restart - from mirgecom.restart import read_restart_data + else: + restart_file = f"{rst_filename}-{rank:04d}.pkl" restart_data = read_restart_data(actx, restart_file) restart_step = restart_data["step"] local_mesh = restart_data["local_mesh"] local_nelements = local_mesh.nelements global_nelements = restart_data["global_nelements"] - restart_order = int(restart_data["order"]) + # restart_order = int(restart_data["order"]) assert comm.Get_size() == restart_data["num_parts"] from mirgecom.discretization import create_discretization_collection - dcoll = create_discretization_collection(actx, local_mesh, order, - mpi_communicator=comm) + dcoll = create_discretization_collection(actx, local_mesh, order) nodes = actx.thaw(dcoll.nodes()) dd_vol = DD_VOLUME_ALL zeros = nodes[0]*0.0 - ones = nodes[0]*0.0 + 1.0 from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD quadrature_tag = DISCR_TAG_QUAD if use_overintegration else DISCR_TAG_BASE -########################################################################## +# ############################################################################ # {{{ Set up initial state using Cantera @@ -323,7 +304,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # import os # current_path = os.path.abspath(os.getcwd()) + "/" # mechanism_file = current_path + my_mechanism - mechanism_file = "uiuc_7sp" + mechanism_file = my_mechanism from mirgecom.mechanisms import get_mechanism_input mech_input = get_mechanism_input(mechanism_file) @@ -336,22 +317,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # Initial temperature, pressure, and mixture mole fractions are needed to # set up the initial state in Cantera. - if False: - cantera_data = np.genfromtxt(cantera_file, skip_header=1, delimiter=',') - - temp_unburned = cantera_data[0,2] - temp_burned = cantera_data[-1,2] - - y_unburned = cantera_data[0,4:] - y_burned = cantera_data[-1,4:] - - vel_unburned = cantera_data[0,1] - vel_burned = cantera_data[-1,1] - - mass_unburned = cantera_data[0,3] - mass_burned = cantera_data[-1,3] - - if True: + if rst_filename is None: temp_unburned = 300.0 cantera_soln.TP = temp_unburned, 101325.0 @@ -366,10 +332,23 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): vel_unburned = 0.5 vel_burned = vel_unburned*mass_unburned/mass_burned - # Uncomment next line to make pylint fail when it can't find cantera.one_atm - one_atm = cantera.one_atm - pres_unburned = one_atm - pres_burned = one_atm + else: + cantera_data = np.genfromtxt(cantera_file, skip_header=1, delimiter=",") + + temp_unburned = cantera_data[0, 2] + temp_burned = cantera_data[-1, 2] + + y_unburned = cantera_data[0, 4:] + y_burned = cantera_data[-1, 4:] + + vel_unburned = cantera_data[0, 1] + vel_burned = cantera_data[-1, 1] + + mass_unburned = cantera_data[0, 3] + mass_burned = cantera_data[-1, 3] + + pres_unburned = cantera.one_atm # pylint: disable=no-member + pres_burned = cantera.one_atm # pylint: disable=no-member # }}} @@ -383,11 +362,12 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): species_names = pyrometheus_mechanism.species_names - # }}} + # }}} if transport == "power-law": from mirgecom.transport import PowerLawTransport - transport_model = PowerLawTransport(lewis=np.ones((nspecies,)), beta=4.093e-7) + transport_model = PowerLawTransport(lewis=np.ones((nspecies,)), + beta=4.093e-7) if transport == "mix-lewis": from mirgecom.transport import MixtureAveragedTransport @@ -399,21 +379,19 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): gas_model = GasModel(eos=eos, transport=transport_model) - print(f"Pyrometheus mechanism species names {species_names}") - print(f"Unburned:") + print("Pyrometheus mechanism species names {species_names}") + print("Unburned:") print(f"T = {temp_unburned}") print(f"D = {mass_unburned}") print(f"Y = {y_unburned}") - print(f"U = {vel_unburned}") - print(f" ") - print(f"Burned:") + print(f"U = {vel_unburned}\n") + print("Burned:") print(f"T = {temp_burned}") print(f"D = {mass_burned}") print(f"Y = {y_burned}") - print(f"U = {vel_burned}") - print(f" ") + print(f"U = {vel_burned}\n") -############################################################################### +# ############################################################################ from mirgecom.limiter import bound_preserving_limiter @@ -421,8 +399,8 @@ def _limit_fluid_cv(cv, pressure, temperature, dd=None): # limit species spec_lim = make_obj_array([ - bound_preserving_limiter(dcoll, cv.species_mass_fractions[i], - mmin=0.0, mmax=1.0, modify_average=True, dd=dd) + bound_preserving_limiter(dcoll, cv.species_mass_fractions[i], + mmin=0.0, mmax=1.0, modify_average=True, dd=dd) for i in range(nspecies)]) # normalize to ensure sum_Yi = 1.0 @@ -447,17 +425,12 @@ def _limit_fluid_cv(cv, pressure, temperature, dd=None): def _get_fluid_state(cv, temp_seed): return make_fluid_state(cv=cv, gas_model=gas_model, - temperature_seed=temp_seed, - limiter_func=_limit_fluid_cv - ) + temperature_seed=temp_seed, limiter_func=_limit_fluid_cv) get_fluid_state = actx.compile(_get_fluid_state) - -############################################################################## - logmgr = initialize_logmgr(use_logmgr, filename=(f"{casename}.sqlite"), - mode="wo", mpi_comm=comm) - +# ############################################################################ + vis_timer = None if logmgr: logmgr_add_cl_device_info(logmgr, queue) @@ -485,9 +458,9 @@ def _get_fluid_state(cv, temp_seed): gc_timer = IntervalTimer("t_gc", "Time spent garbage collecting") logmgr.add_quantity(gc_timer) -############################################################################## +# ############################################################################ - if restart_file is None: + if rst_filename is None: current_t = 0.0 current_step = 0 @@ -498,7 +471,6 @@ def _get_fluid_state(cv, temp_seed): tseed = 1234.56789 from mirgecom.initializers import PlanarDiscontinuity - # use the burned conditions with a lower temperature bulk_init = PlanarDiscontinuity(dim=dim, disc_location=flame_start_loc, sigma=0.0005, nspecies=nspecies, temperature_right=temp_burned, temperature_left=temp_unburned, @@ -509,24 +481,19 @@ def _get_fluid_state(cv, temp_seed): if rank == 0: logging.info("Initializing soln.") - # for Discontinuity initial conditions + # for discontinuous initial conditions current_cv = bulk_init(x_vec=nodes, eos=eos, time=0.) else: - if local_dt: - current_t = restart_data["step"] + current_t = np.min(actx.to_numpy(restart_data["t"])) else: current_t = restart_data["t"] current_step = restart_step first_step = restart_step + 0 - if restart_order != order: - print('Wrong order...') - sys.exit() - else: - current_cv = restart_data["cv"] - tseed = restart_data["temperature_seed"] + current_cv = restart_data["cv"] + tseed = restart_data["temperature_seed"] if logmgr: logmgr_set_time(logmgr, current_step, current_t) @@ -534,57 +501,48 @@ def _get_fluid_state(cv, temp_seed): current_state = get_fluid_state(current_cv, tseed) current_state = force_evaluation(actx, current_state) -############################################################################## +# ############################################################################ # initialize the sponge field - sponge_x_thickness = 0.065 - sponge_amp = 10000.0 - - xMaxLoc = +0.100 - xMinLoc = -0.100 - - sponge_init = InitSponge(x_max=xMaxLoc, - x_min=xMinLoc, - x_thickness=sponge_x_thickness, - amplitude=sponge_amp) + sponge_init = InitSponge(x_max=+0.100, x_min=-0.100, x_thickness=0.065, + amplitude=10000.0) sponge_sigma = sponge_init(x_vec=nodes) - if False: - cantera_data = read_restart_data(actx, original_rst_file) + if use_sponge: + cantera_data = read_restart_data(actx, reference_state) ref_cv = force_evaluation(actx, cantera_data["cv"]) - if True: + else: ref_cv = force_evaluation(actx, bulk_init(x_vec=nodes, eos=eos, time=0.)) -################################################################# - - from mirgecom.boundary import ( - PrescribedFluidBoundary,PressureOutflowBoundary, - LinearizedOutflow2DBoundary, LinearizedInflow2DBoundary) +# ############################################################################ inflow_cv_cond = op.project(dcoll, dd_vol, dd_vol.trace("inlet"), ref_cv) + def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): return make_fluid_state(cv=inflow_cv_cond, gas_model=gas_model, temperature_seed=300.0) - inflow_bnd = PrescribedFluidBoundary( - boundary_state_func=inlet_bnd_state_func) - -# inflow_bnd = LinearizedInflow2DBoundary( -# free_stream_density=mass_unburned, free_stream_pressure=101325.0, -# free_stream_velocity=make_obj_array([vel_unburned, 0.0]), -# free_stream_species_mass_fractions=y_unburned) - + from mirgecom.boundary import ( + PrescribedFluidBoundary, LinearizedOutflow2DBoundary) + inflow_bnd = PrescribedFluidBoundary(boundary_state_func=inlet_bnd_state_func) outflow_bnd = LinearizedOutflow2DBoundary( free_stream_density=mass_burned, free_stream_pressure=101325.0, free_stream_velocity=make_obj_array([vel_burned, 0.0]), free_stream_species_mass_fractions=y_burned) -# outflow_bnd = PressureOutflowBoundary(boundary_pressure=101325.0) + + # from mirgecom.boundary import ( + # LinearizedInflow2DBoundary, PressureOutflowBoundary) + # inflow_bnd = LinearizedInflow2DBoundary( + # free_stream_density=mass_unburned, free_stream_pressure=101325.0, + # free_stream_velocity=make_obj_array([vel_unburned, 0.0]), + # free_stream_species_mass_fractions=y_unburned) + # outflow_bnd = PressureOutflowBoundary(boundary_pressure=101325.0) boundaries = {BoundaryDomainTag("inlet"): inflow_bnd, BoundaryDomainTag("outlet"): outflow_bnd} -#################################################################################### +# ############################################################################ visualizer = make_visualizer(dcoll) @@ -599,16 +557,16 @@ def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): if rank == 0: logger.info(init_message) -################################################################## +# ############################################################################ -# def get_production_rates(cv, temperature): -# return make_obj_array([eos.get_production_rates(cv, temperature)]) -# compute_production_rates = actx.compile(get_production_rates) + # def get_production_rates(cv, temperature): + # return make_obj_array([eos.get_production_rates(cv, temperature)]) + # compute_production_rates = actx.compile(get_production_rates) def my_write_viz(step, t, dt, state): y = state.cv.species_mass_fractions - gas_const = gas_model.eos.gas_const(cv=state.cv) + gas_const = gas_model.eos.gas_const(species_mass_fractions=y) # reaction_rates, = compute_production_rates(state.cv, state.temperature) viz_fields = [("CV_rho", state.cv.mass), @@ -618,7 +576,6 @@ def my_write_viz(step, t, dt, state): ("DV_T", state.temperature), ("DV_U", state.velocity[0]), ("DV_V", state.velocity[1]), - # ("reaction_rates", reaction_rates), ("sponge", sponge_sigma), ("R", gas_const), ("gamma", eos.gamma(state.cv, state.temperature)), @@ -631,17 +588,17 @@ def my_write_viz(step, t, dt, state): viz_fields.extend(("Y_"+species_names[i], y[i]) for i in range(nspecies)) # species diffusivity - viz_fields.extend( - ("diff_"+species_names[i], state.tv.species_diffusivity[i]) - for i in range(nspecies)) + # viz_fields.extend( + # ("diff_"+species_names[i], state.tv.species_diffusivity[i]) + # for i in range(nspecies)) - print('Writing solution file...') + print("Writing solution file...") write_visfile(dcoll, viz_fields, visualizer, vizname=vizname, step=step, t=t, overwrite=True, comm=comm) def my_write_restart(step, t, cv, tseed): restart_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - if restart_fname != restart_file: + if restart_fname != rst_filename: rst_data = { "local_mesh": local_mesh, "cv": cv, @@ -655,7 +612,7 @@ def my_write_restart(step, t, cv, tseed): write_restart_file(actx, rst_data, restart_fname, comm) -################################################################## +# ############################################################################ def my_health_check(cv, dv): health_error = False @@ -670,9 +627,17 @@ def my_health_check(cv, dv): health_error = True logger.info(f"{rank=}: NANs/Infs in temperature data.") + if check_range_local(dcoll, "vol", pressure, 99000., 103000.): + health_error = True + logger.info(f"{rank=}: Pressure range violation.") + + if check_range_local(dcoll, "vol", temperature, 250., 3000.): + health_error = True + logger.info(f"{rank=}: Temperature range violation.") + return health_error -############################################################################## +# ############################################################################ def my_pre_step(step, t, dt, state): @@ -690,7 +655,6 @@ def my_pre_step(step, t, dt, state): dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, gas_model, constant_cfl=constant_cfl, local_dt=local_dt) dt = force_evaluation(actx, dt) - #dt = force_evaluation(actx, actx.np.minimum(dt, current_dt)) else: if constant_cfl: dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, @@ -710,9 +674,7 @@ def my_pre_step(step, t, dt, state): if do_garbage: with gc_timer.start_sub_timer(): - from warnings import warn warn("Running gc.collect() to work around memory growth issue ") - import gc gc.collect() if do_health: @@ -767,16 +729,17 @@ def my_rhs(t, state): chem_rhs = eos.get_species_source_terms(fluid_state.cv, fluid_state.temperature) - #sponge_rhs = sponge_func(cv=fluid_state.cv, cv_ref=ref_cv, - # sigma=sponge_sigma) - rhs = ns_rhs + chem_rhs #+ sponge_rhs + rhs = ns_rhs + chem_rhs + if use_sponge: + sponge_rhs = sponge_func(cv=fluid_state.cv, cv_ref=ref_cv, + sigma=sponge_sigma) + rhs = rhs + sponge_rhs return make_obj_array([rhs, zeros]) def my_post_step(step, t, dt, state): if step == first_step + 1: with gc_timer.start_sub_timer(): - import gc gc.collect() # Freeze the objects that are still alive so they will not # be considered in future gc collections. @@ -791,7 +754,7 @@ def my_post_step(step, t, dt, state): return state, dt -############################################################################## +# ############################################################################ current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, t_final, constant_cfl) @@ -820,14 +783,14 @@ def my_post_step(step, t, dt, state): my_write_viz(step=current_step, t=current_t, dt=current_dt, state=current_state) my_write_restart(step=current_step, t=current_t, cv=current_state.cv, - temperature_seed=tseed) + tseed=tseed) if logmgr: logmgr.close() elif use_profiling: print(actx.tabulate_profiling_data()) - exit() + sys.exit() if __name__ == "__main__": @@ -844,10 +807,12 @@ def my_post_step(step, t, dt, state): parser.add_argument("-c", "--casename", type=ascii, dest="casename", nargs="?", action="store", help="simulation case name") + parser.add_argument("--leap", action="store_true", + help="use leap timestepper") parser.add_argument("--profiling", action="store_true", default=False, help="enable kernel profiling [OFF]") - parser.add_argument("--log", action="store_true", default=True, - help="enable logging profiling [ON]") +# parser.add_argument("--log", action="store_true", default=True, +# help="enable logging profiling [ON]") parser.add_argument("--overintegration", action="store_true", default=False, help="enable overintegration [OFF]") parser.add_argument("--esdg", action="store_true", @@ -867,10 +832,10 @@ def my_post_step(step, t, dt, state): else: print(f"Default casename {casename}") - restart_file = None + rst_filename = None if args.restart_file: - restart_file = (args.restart_file).replace("'", "") - print(f"Restarting from file: {restart_file}") + rst_filename = (args.restart_file).replace("'", "") + print(f"Restarting from file: {rst_filename}") input_file = None if args.input_file: @@ -881,7 +846,6 @@ def my_post_step(step, t, dt, state): print(f"Running {sys.argv[0]}\n") - from warnings import warn from mirgecom.simutil import ApplicationOptionsError if args.esdg: if not args.lazy and not args.numpy: @@ -891,10 +855,10 @@ def my_post_step(step, t, dt, state): 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) + lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) - main(actx_class, use_logmgr=args.log, casename=casename, - restart_file=restart_file, use_overintegration=args.overintegration) + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + use_overintegration=args.overintegration or args.esdg, + casename=casename, rst_filename=rst_filename) # vim: foldmethod=marker 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 From 03d4e4c58c0a22f98b5ee6e07e7d845ac4018b3d Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 07:04:35 -0600 Subject: [PATCH 03/16] Reduce max time --- examples/flame1d/flame1d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/flame1d/flame1d.py b/examples/flame1d/flame1d.py index 07dda71fc..93a493968 100644 --- a/examples/flame1d/flame1d.py +++ b/examples/flame1d/flame1d.py @@ -211,7 +211,7 @@ def main(actx_class, use_esdg=False, use_overintegration=False, # default timestepping control integrator = "compiled_lsrk45" current_dt = 2.0e-9 - t_final = 1.0 + t_final = 1.0e-8 niter = 2000000 From a725218017c75b2aaa34f8cada99d6b13ce3fccd Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 07:51:46 -0600 Subject: [PATCH 04/16] Update README and change max time in flame1d --- examples/README.md | 18 +++++++++++------- examples/flame1d/flame1d.py | 8 ++++---- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/examples/README.md b/examples/README.md index e5ab0b4cf..6265a4974 100644 --- a/examples/README.md +++ b/examples/README.md @@ -4,18 +4,22 @@ This directory has a collection of examples to demonstrate and test *MIRGE-Com* capabilities. All of the example exercise some unique feature of *MIRGE-Com*. 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 -- `vortex.py`: Isentropic vortex advection: outflow boundaries, verification +- `scalar-advdiff.py`: Scalar advection-diffusion verification case +- `pulse.py`: Acoustic pulse in a box, outflow boundaries test case +- `vortex.py`: Isentropic vortex advection: prescribed boundaries, verification - `hotplate.py`: Isothermal BC verification (prescribed exact soln) +- `sod.py`: Sod's shock tube case: Fluid test case with strong shock - `blasius.py`: Inflow, outflow and no-slip wall BC verification - `doublemach.py`: AV test case - `poiseuille.py`: Poiseuille flow verification case - `poiseuille-multispecies.py`: Poiseuille flow with passive scalars -- `scalar-advdiff.py`: Scalar advection-diffusion verification case +- `autoignition.py`: Chemistry verification case with Pyrometheus +- `flame1d.py`: Mixture EOS with Pyrometheus for a laminar 1D flame - `combozzle.py`: Prediction-relevant testing, kitchen sink, many options +- `heat-source.py`: Diffusion operator +- `orthotropic-diffusion.py`: Diffusion operator with orthotropic properties +- `thermally-coupled`: Fluid-solid interaction for heat transfer +- `ablation-workshop`: Verification case for a simplifed AW case 2.1 +- `multiple-volumes`: Demonstrate multiple non-interacting volumes. diff --git a/examples/flame1d/flame1d.py b/examples/flame1d/flame1d.py index 93a493968..86215e833 100644 --- a/examples/flame1d/flame1d.py +++ b/examples/flame1d/flame1d.py @@ -211,7 +211,7 @@ def main(actx_class, use_esdg=False, use_overintegration=False, # default timestepping control integrator = "compiled_lsrk45" current_dt = 2.0e-9 - t_final = 1.0e-8 + t_final = 2.0e-8 niter = 2000000 @@ -219,13 +219,13 @@ def main(actx_class, use_esdg=False, use_overintegration=False, reference_state = "flame_1D_025um_C2H4_p=4_uiuc_7sp-000000-0000.pkl" local_dt = False - constant_cfl = True + constant_cfl = False current_cfl = 0.4 - dim = 2 - # ############################################################################ + dim = 2 + def _compiled_stepper_wrapper(state, t, dt, rhs): return compiled_lsrk45_step(actx, state, t, dt, rhs) From 7ea4e9580188491f243b8929a4556dd5d2b8cff0 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 07:59:35 -0600 Subject: [PATCH 05/16] Move the files out of the folder --- examples/{flame1d => }/flame1d.py | 8 +++----- examples/{flame1d/x_025um.dat => flame1d_x_025um.dat} | 0 examples/{flame1d/y_025um.dat => flame1d_y_025um.dat} | 0 3 files changed, 3 insertions(+), 5 deletions(-) rename examples/{flame1d => }/flame1d.py (99%) rename examples/{flame1d/x_025um.dat => flame1d_x_025um.dat} (100%) rename examples/{flame1d/y_025um.dat => flame1d_y_025um.dat} (100%) diff --git a/examples/flame1d/flame1d.py b/examples/flame1d.py similarity index 99% rename from examples/flame1d/flame1d.py rename to examples/flame1d.py index 86215e833..088352819 100644 --- a/examples/flame1d/flame1d.py +++ b/examples/flame1d.py @@ -255,14 +255,12 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): restart_step = 0 if rst_filename is None: - xx = np.loadtxt("x_025um.dat") - yy = np.loadtxt("y_025um.dat") - - coords = tuple((xx, yy)) # noqa + xx = np.loadtxt("flame1d_x_025um.dat") + yy = np.loadtxt("flame1d_y_025um.dat") from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, - axis_coords=coords, + axis_coords=(xx, yy), periodic=(False, True), boundary_tag_to_face={"inlet": ["-x"], "outlet": ["+x"]}) diff --git a/examples/flame1d/x_025um.dat b/examples/flame1d_x_025um.dat similarity index 100% rename from examples/flame1d/x_025um.dat rename to examples/flame1d_x_025um.dat diff --git a/examples/flame1d/y_025um.dat b/examples/flame1d_y_025um.dat similarity index 100% rename from examples/flame1d/y_025um.dat rename to examples/flame1d_y_025um.dat From bd857bf9b599850524d60c95e0f0f3a2fd15d6b6 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 09:32:56 -0600 Subject: [PATCH 06/16] Include path to grid generation --- examples/flame1d.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 088352819..4f35298b0 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -255,8 +255,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): restart_step = 0 if rst_filename is None: - xx = np.loadtxt("flame1d_x_025um.dat") - yy = np.loadtxt("flame1d_y_025um.dat") + import mirgecom + path = mirgecom.__path__[0] + "/../examples/" + xx = np.loadtxt(path + "flame1d_x_025um.dat") + yy = np.loadtxt(path + "flame1d_y_025um.dat") from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, From 9e10712c13104010de608f396fffe5e0433367f5 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 10:26:53 -0600 Subject: [PATCH 07/16] Try this for the path --- examples/flame1d.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 4f35298b0..ca33979b1 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -241,11 +241,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): print(f"\tnhealth = {nhealth}") print(f"\tnstatus = {nstatus}") if constant_cfl: - print(f"\tconstant_cfl = {constant_cfl}") print(f"\tcurrent_cfl = {current_cfl}") else: print(f"\tcurrent_dt = {current_dt}") - print(f"\tt_final = {t_final}") + print(f"\tt_final = {t_final}") print(f"\tniter = {niter}") print(f"\torder = {order}") print(f"\tTime integration = {integrator}") @@ -255,8 +254,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): restart_step = 0 if rst_filename is None: - import mirgecom - path = mirgecom.__path__[0] + "/../examples/" +# import mirgecom +# path = mirgecom.__path__[0] + "/../examples/" + import os + path = os.path.abspath(os.getcwd()) + "/" xx = np.loadtxt(path + "flame1d_x_025um.dat") yy = np.loadtxt(path + "flame1d_y_025um.dat") From 1bfa68ec34140284587ffd836ac4988c47829aa1 Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 11:39:05 -0600 Subject: [PATCH 08/16] Try this for the path --- examples/README.md | 19 +++++++++++-------- examples/flame1d.py | 8 +++----- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/examples/README.md b/examples/README.md index 6265a4974..48313ecb4 100644 --- a/examples/README.md +++ b/examples/README.md @@ -4,22 +4,25 @@ This directory has a collection of examples to demonstrate and test *MIRGE-Com* capabilities. All of the example exercise some unique feature of *MIRGE-Com*. The examples and the unique features they exercise are as follows: +- `wave.py`: Linear wave equation verification case - `lump.py`: Lump advection, advection verification case - `scalar-lump.py`: Scalar component lump advection verification case - `scalar-advdiff.py`: Scalar advection-diffusion verification case -- `pulse.py`: Acoustic pulse in a box, outflow boundaries test case +- `poiseuille.py`: Poiseuille flow verification case +- `poiseuille-multispecies.py`: Poiseuille flow with passive scalars - `vortex.py`: Isentropic vortex advection: prescribed boundaries, verification +- `taylor-green`: Taylor-Green vortex verification case - `hotplate.py`: Isothermal BC verification (prescribed exact soln) -- `sod.py`: Sod's shock tube case: Fluid test case with strong shock - `blasius.py`: Inflow, outflow and no-slip wall BC verification -- `doublemach.py`: AV test case -- `poiseuille.py`: Poiseuille flow verification case -- `poiseuille-multispecies.py`: Poiseuille flow with passive scalars -- `autoignition.py`: Chemistry verification case with Pyrometheus -- `flame1d.py`: Mixture EOS with Pyrometheus for a laminar 1D flame +- `pulse.py`: Acoustic pulse in a box, wall boundary test case +- `pulse-mixture`: Acoustic pulse with mixture EOS and outflow boundaries +- `sod.py`: Sod's shock tube case: Fluid test case with strong shock +- `doublemach.py`: AV test case for shock capturing +- `autoignition.py`: Mixture EOS and chemistry verification case with Pyrometheus +- `flame1d.py`: Mixture EOS and transport model with Pyrometheus for 1D flame - `combozzle.py`: Prediction-relevant testing, kitchen sink, many options - `heat-source.py`: Diffusion operator - `orthotropic-diffusion.py`: Diffusion operator with orthotropic properties +- `multiple-volumes`: Demonstrate multiple non-interacting volumes. - `thermally-coupled`: Fluid-solid interaction for heat transfer - `ablation-workshop`: Verification case for a simplifed AW case 2.1 -- `multiple-volumes`: Demonstrate multiple non-interacting volumes. diff --git a/examples/flame1d.py b/examples/flame1d.py index ca33979b1..767fbfc11 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -254,12 +254,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): restart_step = 0 if rst_filename is None: -# import mirgecom -# path = mirgecom.__path__[0] + "/../examples/" import os - path = os.path.abspath(os.getcwd()) + "/" - xx = np.loadtxt(path + "flame1d_x_025um.dat") - yy = np.loadtxt(path + "flame1d_y_025um.dat") + path = os.path.dirname(os.path.abspath(__file__)) + xx = np.loadtxt(f"{path}/flame1d_x_025um.dat") + yy = np.loadtxt(f"{path}/flame1d_y_025um.dat") from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, From 654570cca45e7cf44d3067a15f89fa052982a43b Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 12:56:21 -0600 Subject: [PATCH 09/16] I think it is working fine now --- examples/flame1d.py | 104 ++++++++++++++++++++++---------------------- 1 file changed, 52 insertions(+), 52 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 767fbfc11..7644c2b8c 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -32,7 +32,7 @@ import cantera from pytools.obj_array import make_obj_array -from grudge.dof_desc import BoundaryDomainTag, DD_VOLUME_ALL +from grudge.dof_desc import BoundaryDomainTag from grudge.shortcuts import compiled_lsrk45_step from grudge.shortcuts import make_visualizer from grudge import op @@ -173,6 +173,9 @@ def main(actx_class, use_esdg=False, use_overintegration=False, 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) @@ -194,15 +197,14 @@ def main(actx_class, use_esdg=False, use_overintegration=False, # default i/o frequencies ngarbage = 10 - nviz = 25000 + nviz = 10000 nrestart = 25000 nhealth = 1 nstatus = 100 my_mechanism = "uiuc_7sp" - order = 4 - cantera_file = "adiabatic_flame_uiuc_7sp_phi1.00_p1.00_E:H0.0_Y.csv" + order = 2 # transport = "power-law" # transport = "mix-lewis" @@ -210,16 +212,15 @@ def main(actx_class, use_esdg=False, use_overintegration=False, # default timestepping control integrator = "compiled_lsrk45" - current_dt = 2.0e-9 - t_final = 2.0e-8 - + current_dt = 1.0e-9 + t_final = 1.0e-8 niter = 2000000 use_sponge = False - reference_state = "flame_1D_025um_C2H4_p=4_uiuc_7sp-000000-0000.pkl" + use_cantera = True local_dt = False - constant_cfl = False + constant_cfl = True current_cfl = 0.4 # ############################################################################ @@ -283,22 +284,19 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): from mirgecom.discretization import create_discretization_collection dcoll = create_discretization_collection(actx, local_mesh, order) - nodes = actx.thaw(dcoll.nodes()) - - dd_vol = DD_VOLUME_ALL - zeros = nodes[0]*0.0 - from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, DD_VOLUME_ALL quadrature_tag = DISCR_TAG_QUAD if use_overintegration else DISCR_TAG_BASE + dd_vol = DD_VOLUME_ALL # ############################################################################ # {{{ Set up initial state using Cantera # Use Cantera for initialization - print("Using Cantera", cantera.__version__) + print("\nUsing Cantera", cantera.__version__) # import os # current_path = os.path.abspath(os.getcwd()) + "/" @@ -316,7 +314,28 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # Initial temperature, pressure, and mixture mole fractions are needed to # set up the initial state in Cantera. - if rst_filename is None: + if use_cantera: + # Set up flame object + f = cantera.FreeFlame(cantera_soln, width=0.02) + + # Solve with mixture-averaged transport model + f.transport_model = 'mixture-averaged' + f.set_refine_criteria(ratio=2, slope=0.15, curve=0.15) + f.solve(loglevel=0, refine_grid=True, auto=True) + + temp_unburned = f.T[0] + temp_burned = f.T[-1] + + y_unburned = f.Y[:,0] + y_burned = f.Y[:,-1] + + vel_unburned = f.velocity[0] + vel_burned = f.velocity[-1] + + mass_unburned = f.density[0] + mass_burned = f.density[-1] + + else: temp_unburned = 300.0 cantera_soln.TP = temp_unburned, 101325.0 @@ -331,21 +350,6 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): vel_unburned = 0.5 vel_burned = vel_unburned*mass_unburned/mass_burned - else: - cantera_data = np.genfromtxt(cantera_file, skip_header=1, delimiter=",") - - temp_unburned = cantera_data[0, 2] - temp_burned = cantera_data[-1, 2] - - y_unburned = cantera_data[0, 4:] - y_burned = cantera_data[-1, 4:] - - vel_unburned = cantera_data[0, 1] - vel_burned = cantera_data[-1, 1] - - mass_unburned = cantera_data[0, 3] - mass_burned = cantera_data[-1, 3] - pres_unburned = cantera.one_atm # pylint: disable=no-member pres_burned = cantera.one_atm # pylint: disable=no-member @@ -378,7 +382,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): gas_model = GasModel(eos=eos, transport=transport_model) - print("Pyrometheus mechanism species names {species_names}") + print(f"Pyrometheus mechanism species names {species_names}") print("Unburned:") print(f"T = {temp_unburned}") print(f"D = {mass_unburned}") @@ -459,28 +463,27 @@ def _get_fluid_state(cv, temp_seed): # ############################################################################ + flame_start_loc = 0.0 + from mirgecom.initializers import PlanarDiscontinuity + bulk_init = PlanarDiscontinuity(dim=dim, disc_location=flame_start_loc, + sigma=0.0005, nspecies=nspecies, + temperature_right=temp_burned, temperature_left=temp_unburned, + pressure_right=pres_burned, pressure_left=pres_unburned, + velocity_right=make_obj_array([vel_burned, 0.0]), + velocity_left=make_obj_array([vel_unburned, 0.0]), + species_mass_right=y_burned, species_mass_left=y_unburned) + if rst_filename is None: current_t = 0.0 current_step = 0 first_step = 0 - flame_start_loc = 0.0 - tseed = 1234.56789 - from mirgecom.initializers import PlanarDiscontinuity - bulk_init = PlanarDiscontinuity(dim=dim, disc_location=flame_start_loc, - sigma=0.0005, nspecies=nspecies, - temperature_right=temp_burned, temperature_left=temp_unburned, - pressure_right=pres_burned, pressure_left=pres_unburned, - velocity_right=make_obj_array([vel_burned, 0.0]), - velocity_left=make_obj_array([vel_unburned, 0.0]), - species_mass_right=y_burned, species_mass_left=y_unburned) - if rank == 0: logging.info("Initializing soln.") - # for discontinuous initial conditions + current_cv = bulk_init(x_vec=nodes, eos=eos, time=0.) else: @@ -508,11 +511,7 @@ def _get_fluid_state(cv, temp_seed): sponge_sigma = sponge_init(x_vec=nodes) - if use_sponge: - cantera_data = read_restart_data(actx, reference_state) - ref_cv = force_evaluation(actx, cantera_data["cv"]) - else: - ref_cv = force_evaluation(actx, bulk_init(x_vec=nodes, eos=eos, time=0.)) + ref_cv = bulk_init(x_vec=nodes, eos=eos, time=0.) # ############################################################################ @@ -591,7 +590,8 @@ def my_write_viz(step, t, dt, state): # ("diff_"+species_names[i], state.tv.species_diffusivity[i]) # for i in range(nspecies)) - print("Writing solution file...") + if rank == 0: + logger.info("Writing solution file...") write_visfile(dcoll, viz_fields, visualizer, vizname=vizname, step=step, t=t, overwrite=True, comm=comm) @@ -626,11 +626,11 @@ def my_health_check(cv, dv): health_error = True logger.info(f"{rank=}: NANs/Infs in temperature data.") - if check_range_local(dcoll, "vol", pressure, 99000., 103000.): + if check_range_local(dcoll, "vol", pressure, 101250., 101500.): health_error = True logger.info(f"{rank=}: Pressure range violation.") - if check_range_local(dcoll, "vol", temperature, 250., 3000.): + if check_range_local(dcoll, "vol", temperature, 290., 2450.): health_error = True logger.info(f"{rank=}: Temperature range violation.") From 281b92d969357806d55577be0465b3c824e1d5ec Mon Sep 17 00:00:00 2001 From: Tulio Date: Thu, 1 Feb 2024 13:35:19 -0600 Subject: [PATCH 10/16] Fix lazy, change to coarser mesh --- examples/flame1d.py | 27 +- examples/flame1d_x_025um.dat | 715 ------------------ examples/flame1d_x_050um.dat | 499 ++++++++++++ ...lame1d_y_025um.dat => flame1d_y_050um.dat} | 3 +- 4 files changed, 512 insertions(+), 732 deletions(-) delete mode 100644 examples/flame1d_x_025um.dat create mode 100644 examples/flame1d_x_050um.dat rename examples/{flame1d_y_025um.dat => flame1d_y_050um.dat} (50%) diff --git a/examples/flame1d.py b/examples/flame1d.py index 7644c2b8c..002644ff1 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -40,13 +40,10 @@ from logpyle import IntervalTimer, set_dt from mirgecom.navierstokes import ns_operator, grad_cv_operator, grad_t_operator from mirgecom.simutil import ( - check_step, + check_step, check_naninf_local, check_range_local, get_sim_timestep, generate_and_distribute_mesh, write_visfile, - check_naninf_local, - check_range_local, - global_reduce ) from mirgecom.utils import force_evaluation from mirgecom.restart import write_restart_file, read_restart_data @@ -212,17 +209,16 @@ def main(actx_class, use_esdg=False, use_overintegration=False, # default timestepping control integrator = "compiled_lsrk45" + local_dt = False + constant_cfl = True + current_cfl = 0.4 current_dt = 1.0e-9 - t_final = 1.0e-8 + t_final = 2.5e-8 niter = 2000000 use_sponge = False use_cantera = True - local_dt = False - constant_cfl = True - current_cfl = 0.4 - # ############################################################################ dim = 2 @@ -257,8 +253,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): import os path = os.path.dirname(os.path.abspath(__file__)) - xx = np.loadtxt(f"{path}/flame1d_x_025um.dat") - yy = np.loadtxt(f"{path}/flame1d_y_025um.dat") + xx = np.loadtxt(f"{path}/flame1d_x_050um.dat") + yy = np.loadtxt(f"{path}/flame1d_y_050um.dat") from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, @@ -319,15 +315,15 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): f = cantera.FreeFlame(cantera_soln, width=0.02) # Solve with mixture-averaged transport model - f.transport_model = 'mixture-averaged' + f.transport_model = "mixture-averaged" f.set_refine_criteria(ratio=2, slope=0.15, curve=0.15) f.solve(loglevel=0, refine_grid=True, auto=True) temp_unburned = f.T[0] temp_burned = f.T[-1] - y_unburned = f.Y[:,0] - y_burned = f.Y[:,-1] + y_unburned = f.Y[:, 0] + y_burned = f.Y[:, -1] vel_unburned = f.velocity[0] vel_burned = f.velocity[-1] @@ -479,7 +475,7 @@ def _get_fluid_state(cv, temp_seed): current_step = 0 first_step = 0 - tseed = 1234.56789 + tseed = 1234.56789 + zeros if rank == 0: logging.info("Initializing soln.") @@ -500,6 +496,7 @@ def _get_fluid_state(cv, temp_seed): if logmgr: logmgr_set_time(logmgr, current_step, current_t) + current_cv = force_evaluation(actx, current_cv) current_state = get_fluid_state(current_cv, tseed) current_state = force_evaluation(actx, current_state) diff --git a/examples/flame1d_x_025um.dat b/examples/flame1d_x_025um.dat deleted file mode 100644 index cc7ce824b..000000000 --- a/examples/flame1d_x_025um.dat +++ /dev/null @@ -1,715 +0,0 @@ --2.041730741244119907e-01 --1.983270622123499416e-01 --1.926513224919013390e-01 --1.871408955788444461e-01 --1.817909665370416317e-01 --1.765968606712136624e-01 --1.715540394422544745e-01 --1.666580965015173887e-01 --1.619047538406075959e-01 --1.572898580533165225e-01 --1.528093767064319930e-01 --1.484593948162528376e-01 --1.442361114277293799e-01 --1.401358362932405843e-01 --1.361549866481058246e-01 --1.322900840800138311e-01 --1.285377514896332563e-01 --1.248947101397491904e-01 --1.213577767903471910e-01 --1.179238609171413682e-01 --1.145899620111162920e-01 --1.113531669567230098e-01 --1.082106474864382717e-01 --1.051596577094627977e-01 --1.021975317123992311e-01 --9.932168122981324676e-02 --9.652959338264238687e-02 --9.381882848247651008e-02 --9.118701789979119510e-02 --8.863186199427147693e-02 --8.615112810541737665e-02 --8.374264860167553171e-02 --8.140431898639219344e-02 --7.913409605893263754e-02 --7.692999612936025911e-02 --7.479009328511522503e-02 --7.271251770817829807e-02 --7.069545404124924493e-02 --6.873713980151230363e-02 --6.683586384060265229e-02 --6.498996484942823337e-02 --6.319782990654045085e-02 --6.145789306878531727e-02 --5.976863400300363310e-02 --5.812857665758452280e-02 --5.653628797271160328e-02 --5.499037662817479050e-02 --5.348949182765361243e-02 --5.203232211840975313e-02 --5.061759424535746538e-02 --4.924407203851058346e-02 --4.791055533283399698e-02 --4.661587891955575741e-02 --4.535891152802348780e-02 --4.413855483721546014e-02 --4.295374251604261451e-02 --4.180343929160296179e-02 --4.068664004457417382e-02 --3.960236893095399335e-02 --3.854967852938100126e-02 --3.752764901329071895e-02 --3.653538734718364545e-02 --3.557202650630299223e-02 --3.463672471904022387e-02 --3.372866473140646859e-02 --3.284705309292709763e-02 --3.199111946333547257e-02 --3.116011593946010558e-02 --3.035331640171702941e-02 --2.957001587963637226e-02 --2.880952993586874400e-02 --2.807119406813318382e-02 --2.735436312858409783e-02 --2.665841076008983873e-02 --2.598272884893036358e-02 --2.532672699343572684e-02 --2.468983198810112759e-02 --2.407148732272772901e-02 --2.347115269615161412e-02 --2.288830354413596765e-02 --2.232243058101398073e-02 --2.177303935468195492e-02 --2.123964981455377446e-02 --2.072179589209922901e-02 --2.021902509359967126e-02 --1.973089810476514813e-02 --1.925698840686755273e-02 --1.879688190405435361e-02 --1.835017656151726601e-02 --1.791648205419970469e-02 --1.749541942573605222e-02 --1.708662075732473792e-02 --1.668972884624579309e-02 --1.630439689374196463e-02 --1.593028820199067497e-02 --1.556707587990204499e-02 --1.521444255748590109e-02 --1.487208010853818818e-02 --1.453968938140448586e-02 --1.421697993758535769e-02 --1.390366979795513644e-02 --1.359948519637239721e-02 --1.330416034046682551e-02 --1.301743717939345471e-02 --1.273906517835134665e-02 --1.246880109966968901e-02 --1.220640879027002149e-02 --1.195165897531888857e-02 --1.170432905789060493e-02 --1.146420292446508754e-02 --1.123107075609079859e-02 --1.100472884504779918e-02 --1.078497941685071203e-02 --1.057163045743606498e-02 --1.036449554538300920e-02 --1.016339368902081930e-02 --9.968149168280829309e-03 --9.778591381154624895e-03 --9.594554694624330207e-03 --9.415878299934723680e-03 --9.242406072080737098e-03 --9.073986433387545500e-03 --8.910472221064058562e-03 --8.751720558614070958e-03 --8.597592730992724072e-03 --8.447954063399184046e-03 --8.302673803599631280e-03 --8.161625007677735660e-03 --8.024684429112787648e-03 --7.891732411088566226e-03 --7.762652781938836691e-03 --7.637332753638127747e-03 --7.515662823249089769e-03 --7.397536677240314738e-03 --7.282851098591018937e-03 --7.171505876601411145e-03 --7.063403719329947386e-03 --6.958450168580953589e-03 --6.856553517368337881e-03 --6.757624729783274012e-03 --6.661577363195833774e-03 --6.568327492722590699e-03 --6.477793637894199187e-03 --6.389896691458867538e-03 --6.304559850259516818e-03 --6.221708548124225101e-03 --6.141270390711320824e-03 --6.063175092252190770e-03 --5.987354414136530911e-03 --5.913742105286375582e-03 --5.842273844266806897e-03 --5.772887183082759874e-03 --5.705521492612811510e-03 --5.640117909632279269e-03 --5.576619285379335288e-03 --5.514970135619195676e-03 --5.455116592162749696e-03 --5.397006355797268583e-03 --5.340588650588063736e-03 --5.285814179511166139e-03 --5.232635081378255934e-03 --5.181004889016207343e-03 --5.130878488664703495e-03 --5.082212080556446944e-03 --5.034963140645518608e-03 --4.989090383450442250e-03 --4.944553725979494371e-03 --4.901314252706729171e-03 --4.859334181568122392e-03 --4.818576830948115633e-03 --4.779006587627720494e-03 --4.740588875666171959e-03 --4.703290126188940101e-03 --4.667077748055705431e-03 --4.631920099382662130e-03 --4.597786459894270643e-03 --4.564646158254326638e-03 --4.532464991306797870e-03 --4.501202516133172814e-03 --4.470810848816979786e-03 --4.441234589365540379e-03 --4.412411543636324089e-03 --4.384273970275990792e-03 --4.356750139496174540e-03 --4.329766046085505721e-03 --4.303247167161224594e-03 --4.277120192851587704e-03 --4.251314683807731754e-03 --4.225764622822316199e-03 --4.200409829529845422e-03 --4.175197198664148783e-03 --4.150081705905538793e-03 --4.125027103969023355e-03 --4.100006208953808906e-03 --4.075000657206401776e-03 --4.049999999999989372e-03 --4.024999999999989524e-03 --3.999999999999989675e-03 --3.974999999999989826e-03 --3.949999999999989977e-03 --3.924999999999990129e-03 --3.899999999999990280e-03 --3.874999999999990431e-03 --3.849999999999990583e-03 --3.824999999999990734e-03 --3.799999999999990885e-03 --3.774999999999991036e-03 --3.749999999999991188e-03 --3.724999999999991339e-03 --3.699999999999991490e-03 --3.674999999999991641e-03 --3.649999999999991793e-03 --3.624999999999991944e-03 --3.599999999999992095e-03 --3.574999999999992246e-03 --3.549999999999992398e-03 --3.524999999999992549e-03 --3.499999999999992700e-03 --3.474999999999992852e-03 --3.449999999999993003e-03 --3.424999999999993154e-03 --3.399999999999993305e-03 --3.374999999999993457e-03 --3.349999999999993608e-03 --3.324999999999993759e-03 --3.299999999999993910e-03 --3.274999999999994062e-03 --3.249999999999994213e-03 --3.224999999999994364e-03 --3.199999999999994515e-03 --3.174999999999994667e-03 --3.149999999999994818e-03 --3.124999999999994969e-03 --3.099999999999995121e-03 --3.074999999999995272e-03 --3.049999999999995423e-03 --3.024999999999995574e-03 --2.999999999999995726e-03 --2.974999999999995877e-03 --2.949999999999996028e-03 --2.924999999999996179e-03 --2.899999999999996331e-03 --2.874999999999996482e-03 --2.849999999999996633e-03 --2.824999999999996785e-03 --2.799999999999996936e-03 --2.774999999999997087e-03 --2.749999999999997238e-03 --2.724999999999997390e-03 --2.699999999999997541e-03 --2.674999999999997692e-03 --2.649999999999997843e-03 --2.624999999999997995e-03 --2.599999999999998146e-03 --2.574999999999998297e-03 --2.549999999999998448e-03 --2.524999999999998600e-03 --2.499999999999998751e-03 --2.474999999999998902e-03 --2.449999999999999054e-03 --2.424999999999999205e-03 --2.399999999999999356e-03 --2.374999999999999507e-03 --2.349999999999999659e-03 --2.324999999999999810e-03 --2.299999999999999961e-03 --2.275000000000000112e-03 --2.250000000000000264e-03 --2.225000000000000415e-03 --2.200000000000000566e-03 --2.175000000000000717e-03 --2.150000000000000869e-03 --2.125000000000001020e-03 --2.100000000000001171e-03 --2.075000000000001323e-03 --2.050000000000001474e-03 --2.025000000000001625e-03 --2.000000000000001776e-03 --1.975000000000001928e-03 --1.950000000000001862e-03 --1.925000000000001796e-03 --1.900000000000001731e-03 --1.875000000000001665e-03 --1.850000000000001600e-03 --1.825000000000001534e-03 --1.800000000000001469e-03 --1.775000000000001403e-03 --1.750000000000001337e-03 --1.725000000000001272e-03 --1.700000000000001206e-03 --1.675000000000001141e-03 --1.650000000000001075e-03 --1.625000000000001010e-03 --1.600000000000000944e-03 --1.575000000000000878e-03 --1.550000000000000813e-03 --1.525000000000000747e-03 --1.500000000000000682e-03 --1.475000000000000616e-03 --1.450000000000000551e-03 --1.425000000000000485e-03 --1.400000000000000419e-03 --1.375000000000000354e-03 --1.350000000000000288e-03 --1.325000000000000223e-03 --1.300000000000000157e-03 --1.275000000000000092e-03 --1.250000000000000026e-03 --1.224999999999999960e-03 --1.199999999999999895e-03 --1.174999999999999829e-03 --1.149999999999999764e-03 --1.124999999999999698e-03 --1.099999999999999633e-03 --1.074999999999999567e-03 --1.049999999999999501e-03 --1.024999999999999436e-03 --9.999999999999993703e-04 --9.749999999999993047e-04 --9.499999999999993476e-04 --9.249999999999993904e-04 --8.999999999999994333e-04 --8.749999999999994761e-04 --8.499999999999995190e-04 --8.249999999999995618e-04 --7.999999999999996047e-04 --7.749999999999996475e-04 --7.499999999999996904e-04 --7.249999999999997332e-04 --6.999999999999997760e-04 --6.749999999999998189e-04 --6.499999999999998617e-04 --6.249999999999999046e-04 --5.999999999999999474e-04 --5.749999999999999903e-04 --5.500000000000000331e-04 --5.250000000000000760e-04 --5.000000000000001188e-04 --4.750000000000001617e-04 --4.500000000000001503e-04 --4.250000000000001390e-04 --4.000000000000001276e-04 --3.750000000000001162e-04 --3.500000000000001049e-04 --3.250000000000000935e-04 --3.000000000000000821e-04 --2.750000000000000708e-04 --2.500000000000000594e-04 --2.250000000000000481e-04 --2.000000000000000367e-04 --1.750000000000000253e-04 --1.500000000000000140e-04 --1.250000000000000026e-04 --1.000000000000000048e-04 --7.500000000000000698e-05 --5.000000000000000240e-05 --2.500000000000000120e-05 --0.000000000000000000e+00 -2.500000000000000120e-05 -5.000000000000000240e-05 -7.500000000000000698e-05 -1.000000000000000048e-04 -1.250000000000000026e-04 -1.500000000000000140e-04 -1.750000000000000253e-04 -2.000000000000000367e-04 -2.250000000000000481e-04 -2.500000000000000594e-04 -2.750000000000000708e-04 -3.000000000000000821e-04 -3.250000000000000935e-04 -3.500000000000001049e-04 -3.750000000000001162e-04 -4.000000000000001276e-04 -4.250000000000001390e-04 -4.500000000000001503e-04 -4.750000000000001617e-04 -5.000000000000001188e-04 -5.250000000000000760e-04 -5.500000000000000331e-04 -5.749999999999999903e-04 -5.999999999999999474e-04 -6.249999999999999046e-04 -6.499999999999998617e-04 -6.749999999999998189e-04 -6.999999999999997760e-04 -7.249999999999997332e-04 -7.499999999999996904e-04 -7.749999999999996475e-04 -7.999999999999996047e-04 -8.249999999999995618e-04 -8.499999999999995190e-04 -8.749999999999994761e-04 -8.999999999999994333e-04 -9.249999999999993904e-04 -9.499999999999993476e-04 -9.749999999999993047e-04 -9.999999999999993703e-04 -1.024999999999999436e-03 -1.049999999999999501e-03 -1.074999999999999567e-03 -1.099999999999999633e-03 -1.124999999999999698e-03 -1.149999999999999764e-03 -1.174999999999999829e-03 -1.199999999999999895e-03 -1.224999999999999960e-03 -1.250000000000000026e-03 -1.275000000000000092e-03 -1.300000000000000157e-03 -1.325000000000000223e-03 -1.350000000000000288e-03 -1.375000000000000354e-03 -1.400000000000000419e-03 -1.425000000000000485e-03 -1.450000000000000551e-03 -1.475000000000000616e-03 -1.500000000000000682e-03 -1.525000000000000747e-03 -1.550000000000000813e-03 -1.575000000000000878e-03 -1.600000000000000944e-03 -1.625000000000001010e-03 -1.650000000000001075e-03 -1.675000000000001141e-03 -1.700000000000001206e-03 -1.725000000000001272e-03 -1.750000000000001337e-03 -1.775000000000001403e-03 -1.800000000000001469e-03 -1.825000000000001534e-03 -1.850000000000001600e-03 -1.875000000000001665e-03 -1.900000000000001731e-03 -1.925000000000001796e-03 -1.950000000000001862e-03 -1.975000000000001928e-03 -2.000000000000001776e-03 -2.025000000000001625e-03 -2.050000000000001474e-03 -2.075000000000001323e-03 -2.100000000000001171e-03 -2.125000000000001020e-03 -2.150000000000000869e-03 -2.175000000000000717e-03 -2.200000000000000566e-03 -2.225000000000000415e-03 -2.250000000000000264e-03 -2.275000000000000112e-03 -2.299999999999999961e-03 -2.324999999999999810e-03 -2.349999999999999659e-03 -2.374999999999999507e-03 -2.399999999999999356e-03 -2.424999999999999205e-03 -2.449999999999999054e-03 -2.474999999999998902e-03 -2.499999999999998751e-03 -2.524999999999998600e-03 -2.549999999999998448e-03 -2.574999999999998297e-03 -2.599999999999998146e-03 -2.624999999999997995e-03 -2.649999999999997843e-03 -2.674999999999997692e-03 -2.699999999999997541e-03 -2.724999999999997390e-03 -2.749999999999997238e-03 -2.774999999999997087e-03 -2.799999999999996936e-03 -2.824999999999996785e-03 -2.849999999999996633e-03 -2.874999999999996482e-03 -2.899999999999996331e-03 -2.924999999999996179e-03 -2.949999999999996028e-03 -2.974999999999995877e-03 -2.999999999999995726e-03 -3.024999999999995574e-03 -3.049999999999995423e-03 -3.074999999999995272e-03 -3.099999999999995121e-03 -3.124999999999994969e-03 -3.149999999999994818e-03 -3.174999999999994667e-03 -3.199999999999994515e-03 -3.224999999999994364e-03 -3.249999999999994213e-03 -3.274999999999994062e-03 -3.299999999999993910e-03 -3.324999999999993759e-03 -3.349999999999993608e-03 -3.374999999999993457e-03 -3.399999999999993305e-03 -3.424999999999993154e-03 -3.449999999999993003e-03 -3.474999999999992852e-03 -3.499999999999992700e-03 -3.524999999999992549e-03 -3.549999999999992398e-03 -3.574999999999992246e-03 -3.599999999999992095e-03 -3.624999999999991944e-03 -3.649999999999991793e-03 -3.674999999999991641e-03 -3.699999999999991490e-03 -3.724999999999991339e-03 -3.749999999999991188e-03 -3.774999999999991036e-03 -3.799999999999990885e-03 -3.824999999999990734e-03 -3.849999999999990583e-03 -3.874999999999990431e-03 -3.899999999999990280e-03 -3.924999999999990129e-03 -3.949999999999989977e-03 -3.974999999999989826e-03 -3.999999999999989675e-03 -4.024999999999989524e-03 -4.049999999999989372e-03 -4.075000657206401776e-03 -4.100006208953808906e-03 -4.125027103969023355e-03 -4.150081705905538793e-03 -4.175197198664148783e-03 -4.200409829529845422e-03 -4.225764622822316199e-03 -4.251314683807731754e-03 -4.277120192851587704e-03 -4.303247167161224594e-03 -4.329766046085505721e-03 -4.356750139496174540e-03 -4.384273970275990792e-03 -4.412411543636324089e-03 -4.441234589365540379e-03 -4.470810848816979786e-03 -4.501202516133172814e-03 -4.532464991306797870e-03 -4.564646158254326638e-03 -4.597786459894270643e-03 -4.631920099382662130e-03 -4.667077748055705431e-03 -4.703290126188940101e-03 -4.740588875666171959e-03 -4.779006587627720494e-03 -4.818576830948115633e-03 -4.859334181568122392e-03 -4.901314252706729171e-03 -4.944553725979494371e-03 -4.989090383450442250e-03 -5.034963140645518608e-03 -5.082212080556446944e-03 -5.130878488664703495e-03 -5.181004889016207343e-03 -5.232635081378255934e-03 -5.285814179511166139e-03 -5.340588650588063736e-03 -5.397006355797268583e-03 -5.455116592162749696e-03 -5.514970135619195676e-03 -5.576619285379335288e-03 -5.640117909632279269e-03 -5.705521492612811510e-03 -5.772887183082759874e-03 -5.842273844266806897e-03 -5.913742105286375582e-03 -5.987354414136530911e-03 -6.063175092252190770e-03 -6.141270390711320824e-03 -6.221708548124225101e-03 -6.304559850259516818e-03 -6.389896691458867538e-03 -6.477793637894199187e-03 -6.568327492722590699e-03 -6.661577363195833774e-03 -6.757624729783274012e-03 -6.856553517368337881e-03 -6.958450168580953589e-03 -7.063403719329947386e-03 -7.171505876601411145e-03 -7.282851098591018937e-03 -7.397536677240314738e-03 -7.515662823249089769e-03 -7.637332753638127747e-03 -7.762652781938836691e-03 -7.891732411088566226e-03 -8.024684429112787648e-03 -8.161625007677735660e-03 -8.302673803599631280e-03 -8.447954063399184046e-03 -8.597592730992724072e-03 -8.751720558614070958e-03 -8.910472221064058562e-03 -9.073986433387545500e-03 -9.242406072080737098e-03 -9.415878299934723680e-03 -9.594554694624330207e-03 -9.778591381154624895e-03 -9.968149168280829309e-03 -1.016339368902081930e-02 -1.036449554538300920e-02 -1.057163045743606498e-02 -1.078497941685071203e-02 -1.100472884504779918e-02 -1.123107075609079859e-02 -1.146420292446508754e-02 -1.170432905789060493e-02 -1.195165897531888857e-02 -1.220640879027002149e-02 -1.246880109966968901e-02 -1.273906517835134665e-02 -1.301743717939345471e-02 -1.330416034046682551e-02 -1.359948519637239721e-02 -1.390366979795513644e-02 -1.421697993758535769e-02 -1.453968938140448586e-02 -1.487208010853818818e-02 -1.521444255748590109e-02 -1.556707587990204499e-02 -1.593028820199067497e-02 -1.630439689374196463e-02 -1.668972884624579309e-02 -1.708662075732473792e-02 -1.749541942573605222e-02 -1.791648205419970469e-02 -1.835017656151726601e-02 -1.879688190405435361e-02 -1.925698840686755273e-02 -1.973089810476514813e-02 -2.021902509359967126e-02 -2.072179589209922901e-02 -2.123964981455377446e-02 -2.177303935468195492e-02 -2.232243058101398073e-02 -2.288830354413596765e-02 -2.347115269615161412e-02 -2.407148732272772901e-02 -2.468983198810112759e-02 -2.532672699343572684e-02 -2.598272884893036358e-02 -2.665841076008983873e-02 -2.735436312858409783e-02 -2.807119406813318382e-02 -2.880952993586874400e-02 -2.957001587963637226e-02 -3.035331640171702941e-02 -3.116011593946010558e-02 -3.199111946333547257e-02 -3.284705309292709763e-02 -3.372866473140646859e-02 -3.463672471904022387e-02 -3.557202650630299223e-02 -3.653538734718364545e-02 -3.752764901329071895e-02 -3.854967852938100126e-02 -3.960236893095399335e-02 -4.068664004457417382e-02 -4.180343929160296179e-02 -4.295374251604261451e-02 -4.413855483721546014e-02 -4.535891152802348780e-02 -4.661587891955575741e-02 -4.791055533283399698e-02 -4.924407203851058346e-02 -5.061759424535746538e-02 -5.203232211840975313e-02 -5.348949182765361243e-02 -5.499037662817479050e-02 -5.653628797271160328e-02 -5.812857665758452280e-02 -5.976863400300363310e-02 -6.145789306878531727e-02 -6.319782990654045085e-02 -6.498996484942823337e-02 -6.683586384060265229e-02 -6.873713980151230363e-02 -7.069545404124924493e-02 -7.271251770817829807e-02 -7.479009328511522503e-02 -7.692999612936025911e-02 -7.913409605893263754e-02 -8.140431898639219344e-02 -8.374264860167553171e-02 -8.615112810541737665e-02 -8.863186199427147693e-02 -9.118701789979119510e-02 -9.381882848247651008e-02 -9.652959338264238687e-02 -9.932168122981324676e-02 -1.021975317123992311e-01 -1.051596577094627977e-01 -1.082106474864382717e-01 -1.113531669567230098e-01 -1.145899620111162920e-01 -1.179238609171413682e-01 -1.213577767903471910e-01 -1.248947101397491904e-01 -1.285377514896332563e-01 -1.322900840800138311e-01 -1.361549866481058246e-01 -1.401358362932405843e-01 -1.442361114277293799e-01 -1.484593948162528376e-01 -1.528093767064319930e-01 -1.572898580533165225e-01 -1.619047538406075959e-01 -1.666580965015173887e-01 -1.715540394422544745e-01 -1.765968606712136624e-01 -1.817909665370416317e-01 -1.871408955788444461e-01 -1.926513224919013390e-01 -1.983270622123499416e-01 -2.041730741244119907e-01 diff --git a/examples/flame1d_x_050um.dat b/examples/flame1d_x_050um.dat new file mode 100644 index 000000000..8e0d09eb5 --- /dev/null +++ b/examples/flame1d_x_050um.dat @@ -0,0 +1,499 @@ +-2.059645113861788002e-01 +-2.000430594485161362e-01 +-1.942940769847659732e-01 +-1.887125406121930027e-01 +-1.832935732601804202e-01 +-1.780324399087118981e-01 +-1.729245434509754653e-01 +-1.679654206764740754e-01 +-1.631507383711329195e-01 +-1.584762895309958786e-01 +-1.539379896862026476e-01 +-1.495318733320344651e-01 +-1.452540904639100239e-01 +-1.411009032133037633e-01 +-1.370686825816471999e-01 +-1.331539052693592740e-01 +-1.293531505972350670e-01 +-1.256630975175028198e-01 +-1.220805217119375202e-01 +-1.186022927744954797e-01 +-1.152253714760080594e-01 +-1.119468071085445438e-01 +-1.087637349071236548e-01 +-1.056733735465208435e-01 +-1.026730227109841354e-01 +-9.976006073473490943e-02 +-9.693194231119196813e-02 +-9.418619626891726326e-02 +-9.152042341233987410e-02 +-8.893229442537153862e-02 +-8.641954783608189750e-02 +-8.397998804065506806e-02 +-8.161148338490085974e-02 +-7.931196430164434708e-02 +-7.707942150236618550e-02 +-7.491190422151360051e-02 +-7.280751851194798563e-02 +-7.076442559003962018e-02 +-6.878084022896353844e-02 +-6.685502919879258987e-02 +-6.498530975202468363e-02 +-6.317004815322088684e-02 +-6.140765825146962675e-02 +-5.969660009442957077e-02 +-5.803537858274019495e-02 +-5.642254216362429381e-02 +-5.485668156254089323e-02 +-5.333642855178030717e-02 +-5.186045475492537249e-02 +-5.042747048613417660e-02 +-4.903622362323010325e-02 +-4.768549851361449793e-02 +-4.637411491204594671e-02 +-4.510092694935803614e-02 +-4.386482213121443291e-02 +-4.266472036602646928e-02 +-4.149957302118378877e-02 +-4.036836200677342001e-02 +-3.927009888598665538e-02 +-3.820382401143639778e-02 +-3.716860568663032244e-02 +-3.616353935186714108e-02 +-3.518774679384463477e-02 +-3.424037537828880651e-02 +-3.332059730493363431e-02 +-3.242760888420045950e-02 +-3.156062983494495322e-02 +-3.071890260265805472e-02 +-2.990169169752514636e-02 +-2.910828305176504233e-02 +-2.833798339568727176e-02 +-2.759011965192244634e-02 +-2.686403834729640325e-02 +-2.615910504183422552e-02 +-2.547470377439521941e-02 +-2.481023652445443761e-02 +-2.416512268956047538e-02 +-2.353879857801293926e-02 +-2.293071691631630091e-02 +-2.234034637097975784e-02 +-2.176717108424525149e-02 +-2.121069022333796350e-02 +-2.067041754284545257e-02 +-2.014588095984301552e-02 +-1.963662214139404866e-02 +-1.914219610406495478e-02 +-1.866217082510467085e-02 +-1.819612686494905679e-02 +-1.774365700072030524e-02 +-1.730436587040112828e-02 +-1.687786962737279975e-02 +-1.646379560501520095e-02 +-1.606178199107578461e-02 +-1.567147751152295210e-02 +-1.529254112360758054e-02 +-1.492464171786450113e-02 +-1.456745782879355008e-02 +-1.422067735396738419e-02 +-1.388399728132062166e-02 +-1.355712342438201642e-02 +-1.323977016521832133e-02 +-1.293166020486521961e-02 +-1.263252432102725600e-02 +-1.234210113283505834e-02 +-1.206013687245428324e-02 +-1.178638516334673536e-02 +-1.152060680498989340e-02 +-1.126256956386674501e-02 +-1.101204797054330074e-02 +-1.076882312265646166e-02 +-1.053268249364011297e-02 +-1.030341974702229935e-02 +-1.008083455613121779e-02 +-9.864732429052497711e-03 +-9.654924538684809057e-03 +-9.451227557745306082e-03 +-9.253463498580740895e-03 +-9.061459557644270177e-03 +-8.875047964502065512e-03 +-8.694065835237789475e-03 +-8.518355030126842017e-03 +-8.347762015456020035e-03 +-8.182137729367842949e-03 +-8.021337451612330777e-03 +-7.865220677092416013e-03 +-7.713650993092499575e-03 +-7.566495960082871182e-03 +-7.423626995995853510e-03 +-7.284919263872535744e-03 +-7.150251562781936333e-03 +-7.019506221917277099e-03 +-6.892568997776831602e-03 +-6.769328974339505582e-03 +-6.649678466147926695e-03 +-6.533512924214355221e-03 +-6.420730844667198525e-03 +-6.311233680058308572e-03 +-6.204925753253561112e-03 +-6.101714173831476393e-03 +-6.001508756916831312e-03 +-5.904221944378340581e-03 +-5.809768728321553118e-03 +-5.718066576810109172e-03 +-5.629035361750455227e-03 +-5.542597288877004584e-03 +-5.458676829776566881e-03 +-5.377200655892647001e-03 +-5.298097574451947692e-03 +-5.221298466257093783e-03 +-5.146736225291216084e-03 +-5.074345700081626422e-03 +-5.004063636771345289e-03 +-4.935828623848742373e-03 +-4.869581038486991802e-03 +-4.805262994446456853e-03 +-4.742818291494481538e-03 +-4.682192366298389141e-03 +-4.623332244748784570e-03 +-4.566178026132012083e-03 +-4.510631069006149942e-03 +-4.456538895530216295e-03 +-4.403699218111612607e-03 +-4.351875361754715434e-03 +-4.300817375196196342e-03 +-4.250284855199419967e-03 +-4.200068512177228090e-03 +-4.150007631655095244e-03 +-4.100000000000002948e-03 +-4.050000000000003250e-03 +-4.000000000000003553e-03 +-3.950000000000003855e-03 +-3.900000000000003724e-03 +-3.850000000000003593e-03 +-3.800000000000003462e-03 +-3.750000000000003331e-03 +-3.700000000000003200e-03 +-3.650000000000003068e-03 +-3.600000000000002937e-03 +-3.550000000000002806e-03 +-3.500000000000002675e-03 +-3.450000000000002544e-03 +-3.400000000000002413e-03 +-3.350000000000002282e-03 +-3.300000000000002150e-03 +-3.250000000000002019e-03 +-3.200000000000001888e-03 +-3.150000000000001757e-03 +-3.100000000000001626e-03 +-3.050000000000001495e-03 +-3.000000000000001363e-03 +-2.950000000000001232e-03 +-2.900000000000001101e-03 +-2.850000000000000970e-03 +-2.800000000000000839e-03 +-2.750000000000000708e-03 +-2.700000000000000577e-03 +-2.650000000000000445e-03 +-2.600000000000000314e-03 +-2.550000000000000183e-03 +-2.500000000000000052e-03 +-2.449999999999999921e-03 +-2.399999999999999790e-03 +-2.349999999999999659e-03 +-2.299999999999999527e-03 +-2.249999999999999396e-03 +-2.199999999999999265e-03 +-2.149999999999999134e-03 +-2.099999999999999003e-03 +-2.049999999999998872e-03 +-1.999999999999998741e-03 +-1.949999999999998609e-03 +-1.899999999999998695e-03 +-1.849999999999998781e-03 +-1.799999999999998867e-03 +-1.749999999999998952e-03 +-1.699999999999999038e-03 +-1.649999999999999124e-03 +-1.599999999999999209e-03 +-1.549999999999999295e-03 +-1.499999999999999381e-03 +-1.449999999999999466e-03 +-1.399999999999999552e-03 +-1.349999999999999638e-03 +-1.299999999999999723e-03 +-1.249999999999999809e-03 +-1.199999999999999895e-03 +-1.149999999999999981e-03 +-1.100000000000000066e-03 +-1.050000000000000152e-03 +-1.000000000000000238e-03 +-9.500000000000003234e-04 +-9.000000000000003006e-04 +-8.500000000000002779e-04 +-8.000000000000002552e-04 +-7.500000000000002325e-04 +-7.000000000000002097e-04 +-6.500000000000001870e-04 +-6.000000000000001643e-04 +-5.500000000000001416e-04 +-5.000000000000001188e-04 +-4.500000000000000961e-04 +-4.000000000000000734e-04 +-3.500000000000000507e-04 +-3.000000000000000279e-04 +-2.500000000000000052e-04 +-2.000000000000000096e-04 +-1.500000000000000140e-04 +-1.000000000000000048e-04 +-5.000000000000000240e-05 +-0.000000000000000000e+00 +5.000000000000000240e-05 +1.000000000000000048e-04 +1.500000000000000140e-04 +2.000000000000000096e-04 +2.500000000000000052e-04 +3.000000000000000279e-04 +3.500000000000000507e-04 +4.000000000000000734e-04 +4.500000000000000961e-04 +5.000000000000001188e-04 +5.500000000000001416e-04 +6.000000000000001643e-04 +6.500000000000001870e-04 +7.000000000000002097e-04 +7.500000000000002325e-04 +8.000000000000002552e-04 +8.500000000000002779e-04 +9.000000000000003006e-04 +9.500000000000003234e-04 +1.000000000000000238e-03 +1.050000000000000152e-03 +1.100000000000000066e-03 +1.149999999999999981e-03 +1.199999999999999895e-03 +1.249999999999999809e-03 +1.299999999999999723e-03 +1.349999999999999638e-03 +1.399999999999999552e-03 +1.449999999999999466e-03 +1.499999999999999381e-03 +1.549999999999999295e-03 +1.599999999999999209e-03 +1.649999999999999124e-03 +1.699999999999999038e-03 +1.749999999999998952e-03 +1.799999999999998867e-03 +1.849999999999998781e-03 +1.899999999999998695e-03 +1.949999999999998609e-03 +1.999999999999998741e-03 +2.049999999999998872e-03 +2.099999999999999003e-03 +2.149999999999999134e-03 +2.199999999999999265e-03 +2.249999999999999396e-03 +2.299999999999999527e-03 +2.349999999999999659e-03 +2.399999999999999790e-03 +2.449999999999999921e-03 +2.500000000000000052e-03 +2.550000000000000183e-03 +2.600000000000000314e-03 +2.650000000000000445e-03 +2.700000000000000577e-03 +2.750000000000000708e-03 +2.800000000000000839e-03 +2.850000000000000970e-03 +2.900000000000001101e-03 +2.950000000000001232e-03 +3.000000000000001363e-03 +3.050000000000001495e-03 +3.100000000000001626e-03 +3.150000000000001757e-03 +3.200000000000001888e-03 +3.250000000000002019e-03 +3.300000000000002150e-03 +3.350000000000002282e-03 +3.400000000000002413e-03 +3.450000000000002544e-03 +3.500000000000002675e-03 +3.550000000000002806e-03 +3.600000000000002937e-03 +3.650000000000003068e-03 +3.700000000000003200e-03 +3.750000000000003331e-03 +3.800000000000003462e-03 +3.850000000000003593e-03 +3.900000000000003724e-03 +3.950000000000003855e-03 +4.000000000000003553e-03 +4.050000000000003250e-03 +4.100000000000002948e-03 +4.150007631655095244e-03 +4.200068512177228090e-03 +4.250284855199419967e-03 +4.300817375196196342e-03 +4.351875361754715434e-03 +4.403699218111612607e-03 +4.456538895530216295e-03 +4.510631069006149942e-03 +4.566178026132012083e-03 +4.623332244748784570e-03 +4.682192366298389141e-03 +4.742818291494481538e-03 +4.805262994446456853e-03 +4.869581038486991802e-03 +4.935828623848742373e-03 +5.004063636771345289e-03 +5.074345700081626422e-03 +5.146736225291216084e-03 +5.221298466257093783e-03 +5.298097574451947692e-03 +5.377200655892647001e-03 +5.458676829776566881e-03 +5.542597288877004584e-03 +5.629035361750455227e-03 +5.718066576810109172e-03 +5.809768728321553118e-03 +5.904221944378340581e-03 +6.001508756916831312e-03 +6.101714173831476393e-03 +6.204925753253561112e-03 +6.311233680058308572e-03 +6.420730844667198525e-03 +6.533512924214355221e-03 +6.649678466147926695e-03 +6.769328974339505582e-03 +6.892568997776831602e-03 +7.019506221917277099e-03 +7.150251562781936333e-03 +7.284919263872535744e-03 +7.423626995995853510e-03 +7.566495960082871182e-03 +7.713650993092499575e-03 +7.865220677092416013e-03 +8.021337451612330777e-03 +8.182137729367842949e-03 +8.347762015456020035e-03 +8.518355030126842017e-03 +8.694065835237789475e-03 +8.875047964502065512e-03 +9.061459557644270177e-03 +9.253463498580740895e-03 +9.451227557745306082e-03 +9.654924538684809057e-03 +9.864732429052497711e-03 +1.008083455613121779e-02 +1.030341974702229935e-02 +1.053268249364011297e-02 +1.076882312265646166e-02 +1.101204797054330074e-02 +1.126256956386674501e-02 +1.152060680498989340e-02 +1.178638516334673536e-02 +1.206013687245428324e-02 +1.234210113283505834e-02 +1.263252432102725600e-02 +1.293166020486521961e-02 +1.323977016521832133e-02 +1.355712342438201642e-02 +1.388399728132062166e-02 +1.422067735396738419e-02 +1.456745782879355008e-02 +1.492464171786450113e-02 +1.529254112360758054e-02 +1.567147751152295210e-02 +1.606178199107578461e-02 +1.646379560501520095e-02 +1.687786962737279975e-02 +1.730436587040112828e-02 +1.774365700072030524e-02 +1.819612686494905679e-02 +1.866217082510467085e-02 +1.914219610406495478e-02 +1.963662214139404866e-02 +2.014588095984301552e-02 +2.067041754284545257e-02 +2.121069022333796350e-02 +2.176717108424525149e-02 +2.234034637097975784e-02 +2.293071691631630091e-02 +2.353879857801293926e-02 +2.416512268956047538e-02 +2.481023652445443761e-02 +2.547470377439521941e-02 +2.615910504183422552e-02 +2.686403834729640325e-02 +2.759011965192244634e-02 +2.833798339568727176e-02 +2.910828305176504233e-02 +2.990169169752514636e-02 +3.071890260265805472e-02 +3.156062983494495322e-02 +3.242760888420045950e-02 +3.332059730493363431e-02 +3.424037537828880651e-02 +3.518774679384463477e-02 +3.616353935186714108e-02 +3.716860568663032244e-02 +3.820382401143639778e-02 +3.927009888598665538e-02 +4.036836200677342001e-02 +4.149957302118378877e-02 +4.266472036602646928e-02 +4.386482213121443291e-02 +4.510092694935803614e-02 +4.637411491204594671e-02 +4.768549851361449793e-02 +4.903622362323010325e-02 +5.042747048613417660e-02 +5.186045475492537249e-02 +5.333642855178030717e-02 +5.485668156254089323e-02 +5.642254216362429381e-02 +5.803537858274019495e-02 +5.969660009442957077e-02 +6.140765825146962675e-02 +6.317004815322088684e-02 +6.498530975202468363e-02 +6.685502919879258987e-02 +6.878084022896353844e-02 +7.076442559003962018e-02 +7.280751851194798563e-02 +7.491190422151360051e-02 +7.707942150236618550e-02 +7.931196430164434708e-02 +8.161148338490085974e-02 +8.397998804065506806e-02 +8.641954783608189750e-02 +8.893229442537153862e-02 +9.152042341233987410e-02 +9.418619626891726326e-02 +9.693194231119196813e-02 +9.976006073473490943e-02 +1.026730227109841354e-01 +1.056733735465208435e-01 +1.087637349071236548e-01 +1.119468071085445438e-01 +1.152253714760080594e-01 +1.186022927744954797e-01 +1.220805217119375202e-01 +1.256630975175028198e-01 +1.293531505972350670e-01 +1.331539052693592740e-01 +1.370686825816471999e-01 +1.411009032133037633e-01 +1.452540904639100239e-01 +1.495318733320344651e-01 +1.539379896862026476e-01 +1.584762895309958786e-01 +1.631507383711329195e-01 +1.679654206764740754e-01 +1.729245434509754653e-01 +1.780324399087118981e-01 +1.832935732601804202e-01 +1.887125406121930027e-01 +1.942940769847659732e-01 +2.000430594485161362e-01 +2.059645113861788002e-01 diff --git a/examples/flame1d_y_025um.dat b/examples/flame1d_y_050um.dat similarity index 50% rename from examples/flame1d_y_025um.dat rename to examples/flame1d_y_050um.dat index 28c76dc59..df6d1f5d0 100644 --- a/examples/flame1d_y_025um.dat +++ b/examples/flame1d_y_050um.dat @@ -1,4 +1,3 @@ 0.000000000000000000e+00 -2.500000000000000120e-05 5.000000000000000240e-05 -7.500000000000000698e-05 +1.000000000000000698e-04 From c5f279f9c7c2b4c9c5949cb00945d4efb0a820a8 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 9 Feb 2024 10:09:16 -0600 Subject: [PATCH 11/16] Fix example --- examples/flame1d.py | 84 ++++++++++++++++++++++++++++----------------- 1 file changed, 53 insertions(+), 31 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 002644ff1..8027c048b 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -428,6 +428,15 @@ def _get_fluid_state(cv, temp_seed): get_fluid_state = actx.compile(_get_fluid_state) + def get_temperature_update(cv, temperature): + y = cv.species_mass_fractions + e = eos.internal_energy(cv) / cv.mass + return make_obj_array( + [pyrometheus_mechanism.get_temperature_update_energy(e, temperature, y)] + ) + + compute_temperature_update = actx.compile(get_temperature_update) + # ############################################################################ vis_timer = None @@ -512,27 +521,27 @@ def _get_fluid_state(cv, temp_seed): # ############################################################################ - inflow_cv_cond = op.project(dcoll, dd_vol, dd_vol.trace("inlet"), ref_cv) + # inflow_cv_cond = op.project(dcoll, dd_vol, dd_vol.trace("inlet"), ref_cv) - def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): - return make_fluid_state(cv=inflow_cv_cond, gas_model=gas_model, - temperature_seed=300.0) - - from mirgecom.boundary import ( - PrescribedFluidBoundary, LinearizedOutflow2DBoundary) - inflow_bnd = PrescribedFluidBoundary(boundary_state_func=inlet_bnd_state_func) - outflow_bnd = LinearizedOutflow2DBoundary( - free_stream_density=mass_burned, free_stream_pressure=101325.0, - free_stream_velocity=make_obj_array([vel_burned, 0.0]), - free_stream_species_mass_fractions=y_burned) + # def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): + # return make_fluid_state(cv=inflow_cv_cond, gas_model=gas_model, + # temperature_seed=300.0) # from mirgecom.boundary import ( - # LinearizedInflow2DBoundary, PressureOutflowBoundary) - # inflow_bnd = LinearizedInflow2DBoundary( - # free_stream_density=mass_unburned, free_stream_pressure=101325.0, - # free_stream_velocity=make_obj_array([vel_unburned, 0.0]), - # free_stream_species_mass_fractions=y_unburned) - # outflow_bnd = PressureOutflowBoundary(boundary_pressure=101325.0) + # PrescribedFluidBoundary, LinearizedOutflow2DBoundary) + # inflow_bnd = PrescribedFluidBoundary(boundary_state_func=inlet_bnd_state_func) + # outflow_bnd = LinearizedOutflow2DBoundary( + # free_stream_density=mass_burned, free_stream_pressure=101325.0, + # free_stream_velocity=make_obj_array([vel_burned, 0.0]), + # free_stream_species_mass_fractions=y_burned) + + from mirgecom.boundary import ( + LinearizedInflow2DBoundary, PressureOutflowBoundary) + inflow_bnd = LinearizedInflow2DBoundary( + free_stream_density=mass_unburned, free_stream_pressure=101325.0, + free_stream_velocity=make_obj_array([vel_unburned, 0.0]), + free_stream_species_mass_fractions=y_unburned) + outflow_bnd = PressureOutflowBoundary(boundary_pressure=101325.0) boundaries = {BoundaryDomainTag("inlet"): inflow_bnd, BoundaryDomainTag("outlet"): outflow_bnd} @@ -554,29 +563,29 @@ def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): # ############################################################################ - # def get_production_rates(cv, temperature): - # return make_obj_array([eos.get_production_rates(cv, temperature)]) - # compute_production_rates = actx.compile(get_production_rates) + def get_production_rates(cv, temperature): + return make_obj_array([eos.get_production_rates(cv, temperature)]) + compute_production_rates = actx.compile(get_production_rates) def my_write_viz(step, t, dt, state): y = state.cv.species_mass_fractions - gas_const = gas_model.eos.gas_const(species_mass_fractions=y) + # gas_const = gas_model.eos.gas_const(species_mass_fractions=y) + # gamma = eos.gamma(state.cv, state.temperature) - # reaction_rates, = compute_production_rates(state.cv, state.temperature) + reaction_rates, = compute_production_rates(state.cv, state.temperature) viz_fields = [("CV_rho", state.cv.mass), ("CV_rhoU", state.cv.momentum), ("CV_rhoE", state.cv.energy), ("DV_P", state.pressure), ("DV_T", state.temperature), - ("DV_U", state.velocity[0]), - ("DV_V", state.velocity[1]), - ("sponge", sponge_sigma), - ("R", gas_const), - ("gamma", eos.gamma(state.cv, state.temperature)), - ("dt", dt), - ("mu", state.tv.viscosity), - ("kappa", state.tv.thermal_conductivity), + ("reaction_rates", reaction_rates), + # ("sponge", sponge_sigma), + # ("R", gas_const), + # ("gamma", gamma), + # ("dt", dt), + # ("mu", state.tv.viscosity), + # ("kappa", state.tv.thermal_conductivity), ] # species mass fractions @@ -631,6 +640,19 @@ def my_health_check(cv, dv): health_error = True logger.info(f"{rank=}: Temperature range violation.") + # temperature_update is the next temperature update in the + # `get_temperature` Newton solve. The relative size of this + # update is used to gauge convergence of the current temperature + # after a fixed number of Newton iters. + # Note: The local max jig below works around a very long compile + # in lazy mode. + temp_update, = compute_temperature_update(cv, temperature) + temp_resid = force_evaluation(actx, temp_update) / temperature + temp_resid = (actx.to_numpy(op.nodal_max_loc(dcoll, "vol", temp_resid))) + if temp_resid > 1e-8: + health_error = True + logger.info(f"{rank=}: Temperature is not converged {temp_resid=}.") + return health_error # ############################################################################ From 58229eefc645963ac009e4f73d7801eac4d8e4d6 Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 9 Feb 2024 10:12:31 -0600 Subject: [PATCH 12/16] flake8 --- examples/flame1d.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 8027c048b..919f94b3b 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -283,9 +283,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): nodes = actx.thaw(dcoll.nodes()) zeros = nodes[0]*0.0 - from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD, DD_VOLUME_ALL + from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD quadrature_tag = DISCR_TAG_QUAD if use_overintegration else DISCR_TAG_BASE - dd_vol = DD_VOLUME_ALL # ############################################################################ @@ -521,6 +520,9 @@ def get_temperature_update(cv, temperature): # ############################################################################ + # from grudge.dof_desc import DD_VOLUME_ALL + # dd_vol = DD_VOLUME_ALL + # inflow_cv_cond = op.project(dcoll, dd_vol, dd_vol.trace("inlet"), ref_cv) # def inlet_bnd_state_func(dcoll, dd_bdry, gas_model, state_minus, **kwargs): From 91f453c4a20b05265cf88e577cc911d597e2bd6d Mon Sep 17 00:00:00 2001 From: Tulio Date: Fri, 9 Feb 2024 10:23:59 -0600 Subject: [PATCH 13/16] Cosmetics --- examples/flame1d.py | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 919f94b3b..c4423a072 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -199,7 +199,7 @@ def main(actx_class, use_esdg=False, use_overintegration=False, nhealth = 1 nstatus = 100 - my_mechanism = "uiuc_7sp" + mechanism_file = "uiuc_7sp" order = 2 @@ -274,8 +274,9 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): local_mesh = restart_data["local_mesh"] local_nelements = local_mesh.nelements global_nelements = restart_data["global_nelements"] - # restart_order = int(restart_data["order"]) + restart_order = int(restart_data["order"]) + assert restart_order == order assert comm.Get_size() == restart_data["num_parts"] from mirgecom.discretization import create_discretization_collection @@ -293,11 +294,6 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # Use Cantera for initialization print("\nUsing Cantera", cantera.__version__) - # import os - # current_path = os.path.abspath(os.getcwd()) + "/" - # mechanism_file = current_path + my_mechanism - mechanism_file = my_mechanism - from mirgecom.mechanisms import get_mechanism_input mech_input = get_mechanism_input(mechanism_file) @@ -790,8 +786,7 @@ def my_post_step(step, t, dt, state): istep=current_step, dt=current_dt, t=current_t, t_final=t_final, max_steps=niter, local_dt=local_dt, force_eval=force_eval_stepper, - state=make_obj_array([current_state.cv, tseed]), - # compile_rhs=False + state=make_obj_array([current_state.cv, tseed]) ) current_cv, tseed = stepper_state current_state = make_fluid_state(current_cv, gas_model, tseed) From 5d1ad8d3f8a4d1f3a68d55c75c137f5cf764e984 Mon Sep 17 00:00:00 2001 From: Tulio Date: Mon, 19 Feb 2024 18:07:50 -0600 Subject: [PATCH 14/16] Add tpe; fix get_sim_step --- examples/flame1d.py | 105 +++++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 50 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index c4423a072..4274366a1 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -161,7 +161,7 @@ def __call__(self, x_vec): @mpi_entry_point -def main(actx_class, use_esdg=False, use_overintegration=False, +def main(actx_class, use_esdg=False, use_tpe=False, use_overintegration=False, use_leap=False, casename=None, rst_filename=None): from mpi4py import MPI @@ -199,7 +199,7 @@ def main(actx_class, use_esdg=False, use_overintegration=False, nhealth = 1 nstatus = 100 - mechanism_file = "uiuc_7sp" + mechanism_file = "wang99_reduced" order = 2 @@ -217,7 +217,10 @@ def main(actx_class, use_esdg=False, use_overintegration=False, niter = 2000000 use_sponge = False - use_cantera = True + + # use Cantera's 1D flame solution to prescribe the BC and an + # approximated initial condition with hyperbolic tangent profile + use_flame_from_cantera = True # ############################################################################ @@ -256,12 +259,16 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): xx = np.loadtxt(f"{path}/flame1d_x_050um.dat") yy = np.loadtxt(f"{path}/flame1d_y_050um.dat") + from meshmode.mesh import TensorProductElementGroup + grp_cls = TensorProductElementGroup if use_tpe else None + from meshmode.mesh.generation import generate_box_mesh generate_mesh = partial(generate_box_mesh, axis_coords=(xx, yy), periodic=(False, True), boundary_tag_to_face={"inlet": ["-x"], - "outlet": ["+x"]}) + "outlet": ["+x"]}, + group_cls=grp_cls) local_mesh, global_nelements = ( generate_and_distribute_mesh(comm, generate_mesh)) @@ -280,7 +287,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): assert comm.Get_size() == restart_data["num_parts"] from mirgecom.discretization import create_discretization_collection - dcoll = create_discretization_collection(actx, local_mesh, order) + dcoll = create_discretization_collection(actx, local_mesh, order, + tensor_product_elements=use_tpe) nodes = actx.thaw(dcoll.nodes()) zeros = nodes[0]*0.0 @@ -305,7 +313,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # Initial temperature, pressure, and mixture mole fractions are needed to # set up the initial state in Cantera. - if use_cantera: + if use_flame_from_cantera: # Set up flame object f = cantera.FreeFlame(cantera_soln, width=0.02) @@ -350,11 +358,10 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # Import Pyrometheus EOS from mirgecom.thermochemistry import get_pyrometheus_wrapper_class_from_cantera - pyrometheus_mechanism = \ - get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) - eos = PyrometheusMixture(pyrometheus_mechanism, temperature_guess=1350.0) + pyro_mech = get_pyrometheus_wrapper_class_from_cantera(cantera_soln)(actx.np) + eos = PyrometheusMixture(pyro_mech, temperature_guess=1234.56789) - species_names = pyrometheus_mechanism.species_names + species_names = pyro_mech.species_names # }}} @@ -365,11 +372,11 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): if transport == "mix-lewis": from mirgecom.transport import MixtureAveragedTransport - transport_model = MixtureAveragedTransport(pyrometheus_mechanism, + transport_model = MixtureAveragedTransport(pyro_mech, lewis=np.ones(nspecies,)) if transport == "mix": from mirgecom.transport import MixtureAveragedTransport - transport_model = MixtureAveragedTransport(pyrometheus_mechanism) + transport_model = MixtureAveragedTransport(pyro_mech) gas_model = GasModel(eos=eos, transport=transport_model) @@ -427,40 +434,11 @@ def get_temperature_update(cv, temperature): y = cv.species_mass_fractions e = eos.internal_energy(cv) / cv.mass return make_obj_array( - [pyrometheus_mechanism.get_temperature_update_energy(e, temperature, y)] + [pyro_mech.get_temperature_update_energy(e, temperature, y)] ) compute_temperature_update = actx.compile(get_temperature_update) -# ############################################################################ - - vis_timer = None - if logmgr: - logmgr_add_cl_device_info(logmgr, queue) - logmgr_add_device_memory_usage(logmgr, queue) - - logmgr.add_watches([ - ("step.max", "step = {value}, "), - ("dt.max", "dt: {value:1.6e} s, "), - ("t_sim.max", "sim time: {value:1.6e} s, "), - ("t_step.max", "------- step walltime: {value:6g} s\n") - ]) - - try: - logmgr.add_watches(["memory_usage_python.max", - "memory_usage_gpu.max"]) - except KeyError: - pass - - if use_profiling: - logmgr.add_watches(["pyopencl_array_time.max"]) - - vis_timer = IntervalTimer("t_vis", "Time spent visualizing") - logmgr.add_quantity(vis_timer) - - gc_timer = IntervalTimer("t_gc", "Time spent garbage collecting") - logmgr.add_quantity(gc_timer) - # ############################################################################ flame_start_loc = 0.0 @@ -474,7 +452,6 @@ def get_temperature_update(cv, temperature): species_mass_right=y_burned, species_mass_left=y_unburned) if rst_filename is None: - current_t = 0.0 current_step = 0 first_step = 0 @@ -497,15 +474,40 @@ def get_temperature_update(cv, temperature): current_cv = restart_data["cv"] tseed = restart_data["temperature_seed"] - if logmgr: - logmgr_set_time(logmgr, current_step, current_t) - + tseed = force_evaluation(actx, tseed) current_cv = force_evaluation(actx, current_cv) current_state = get_fluid_state(current_cv, tseed) - current_state = force_evaluation(actx, current_state) # ############################################################################ + vis_timer = None + if logmgr: + logmgr_add_cl_device_info(logmgr, queue) + logmgr_add_device_memory_usage(logmgr, queue) + logmgr_set_time(logmgr, current_step, current_t) + + logmgr.add_watches([ + ("step.max", "step = {value}, "), + ("dt.max", "dt: {value:1.6e} s, "), + ("t_sim.max", "sim time: {value:1.6e} s, "), + ("t_step.max", "------- step walltime: {value:6g} s\n") + ]) + + try: + logmgr.add_watches(["memory_usage_python.max", + "memory_usage_gpu.max"]) + except KeyError: + pass + + if use_profiling: + logmgr.add_watches(["pyopencl_array_time.max"]) + + vis_timer = IntervalTimer("t_vis", "Time spent visualizing") + logmgr.add_quantity(vis_timer) + + gc_timer = IntervalTimer("t_gc", "Time spent garbage collecting") + logmgr.add_quantity(gc_timer) + # initialize the sponge field sponge_init = InitSponge(x_max=+0.100, x_min=-0.100, x_thickness=0.065, amplitude=10000.0) @@ -669,7 +671,7 @@ def my_pre_step(step, t, dt, state): if local_dt: t = force_evaluation(actx, t) dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, - gas_model, constant_cfl=constant_cfl, local_dt=local_dt) + constant_cfl=constant_cfl, local_dt=local_dt) dt = force_evaluation(actx, dt) else: if constant_cfl: @@ -809,7 +811,9 @@ def my_post_step(step, t, dt, state): if __name__ == "__main__": - logging.basicConfig(format="%(message)s", level=logging.INFO) + logging.basicConfig( + format="%(asctime)s - %(levelname)s - %(name)s - %(message)s", + level=logging.INFO) import argparse parser = argparse.ArgumentParser(description="MIRGE-Com 1D Flame Driver") @@ -836,6 +840,7 @@ def my_post_step(step, t, dt, state): help="enable lazy evaluation [OFF]") parser.add_argument("--numpy", action="store_true", help="use numpy-based eager actx.") + parser.add_argument("--tpe", action="store_true") args = parser.parse_args() @@ -872,7 +877,7 @@ def my_post_step(step, t, dt, state): actx_class = get_reasonable_array_context_class( lazy=args.lazy, distributed=True, profiling=args.profiling, numpy=args.numpy) - main(actx_class, use_leap=args.leap, use_esdg=args.esdg, + main(actx_class, use_leap=args.leap, use_esdg=args.esdg, use_tpe=args.tpe, use_overintegration=args.overintegration or args.esdg, casename=casename, rst_filename=rst_filename) From 64ac9f87963913e045d530110da000145dd3a82e Mon Sep 17 00:00:00 2001 From: Tulio Date: Tue, 20 Feb 2024 16:38:20 -0600 Subject: [PATCH 15/16] Revert the mechanism to 7sp; Cosmetics --- examples/flame1d.py | 65 +++++++++++++++++---------------------------- 1 file changed, 24 insertions(+), 41 deletions(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 4274366a1..4de5f004c 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -42,7 +42,7 @@ from mirgecom.simutil import ( check_step, check_naninf_local, check_range_local, get_sim_timestep, - generate_and_distribute_mesh, + distribute_mesh, write_visfile, ) from mirgecom.utils import force_evaluation @@ -199,7 +199,7 @@ def main(actx_class, use_esdg=False, use_tpe=False, use_overintegration=False, nhealth = 1 nstatus = 100 - mechanism_file = "wang99_reduced" + mechanism_file = "uiuc_7sp" order = 2 @@ -209,11 +209,10 @@ def main(actx_class, use_esdg=False, use_tpe=False, use_overintegration=False, # default timestepping control integrator = "compiled_lsrk45" - local_dt = False - constant_cfl = True + constant_cfl = False current_cfl = 0.4 current_dt = 1.0e-9 - t_final = 2.5e-8 + t_final = 1.0e-8 niter = 2000000 use_sponge = False @@ -270,8 +269,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): "outlet": ["+x"]}, group_cls=grp_cls) - local_mesh, global_nelements = ( - generate_and_distribute_mesh(comm, generate_mesh)) + local_mesh, global_nelements = distribute_mesh(comm, generate_mesh) local_nelements = local_mesh.nelements else: @@ -290,7 +288,7 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): dcoll = create_discretization_collection(actx, local_mesh, order, tensor_product_elements=use_tpe) nodes = actx.thaw(dcoll.nodes()) - zeros = nodes[0]*0.0 + zeros = actx.np.zeros_like(nodes[0]) from grudge.dof_desc import DISCR_TAG_BASE, DISCR_TAG_QUAD quadrature_tag = DISCR_TAG_QUAD if use_overintegration else DISCR_TAG_BASE @@ -405,7 +403,7 @@ def _limit_fluid_cv(cv, pressure, temperature, dd=None): for i in range(nspecies)]) # normalize to ensure sum_Yi = 1.0 - aux = cv.mass*0.0 + aux = actx.np.zeros_like(cv.mass) for i in range(0, nspecies): aux = aux + spec_lim[i] spec_lim = spec_lim/aux @@ -464,10 +462,7 @@ def get_temperature_update(cv, temperature): current_cv = bulk_init(x_vec=nodes, eos=eos, time=0.) else: - if local_dt: - current_t = np.min(actx.to_numpy(restart_data["t"])) - else: - current_t = restart_data["t"] + current_t = restart_data["t"] current_step = restart_step first_step = restart_step + 0 @@ -563,9 +558,9 @@ def get_temperature_update(cv, temperature): # ############################################################################ - def get_production_rates(cv, temperature): - return make_obj_array([eos.get_production_rates(cv, temperature)]) - compute_production_rates = actx.compile(get_production_rates) + # def get_production_rates(cv, temperature): + # return make_obj_array([eos.get_production_rates(cv, temperature)]) + # compute_production_rates = actx.compile(get_production_rates) def my_write_viz(step, t, dt, state): @@ -573,13 +568,13 @@ def my_write_viz(step, t, dt, state): # gas_const = gas_model.eos.gas_const(species_mass_fractions=y) # gamma = eos.gamma(state.cv, state.temperature) - reaction_rates, = compute_production_rates(state.cv, state.temperature) + # reaction_rates, = compute_production_rates(state.cv, state.temperature) viz_fields = [("CV_rho", state.cv.mass), ("CV_rhoU", state.cv.momentum), ("CV_rhoE", state.cv.energy), ("DV_P", state.pressure), ("DV_T", state.temperature), - ("reaction_rates", reaction_rates), + # ("reaction_rates", reaction_rates), # ("sponge", sponge_sigma), # ("R", gas_const), # ("gamma", gamma), @@ -668,15 +663,9 @@ def my_pre_step(step, t, dt, state): fluid_state = get_fluid_state(cv, tseed) - if local_dt: - t = force_evaluation(actx, t) + if constant_cfl: dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, - constant_cfl=constant_cfl, local_dt=local_dt) - dt = force_evaluation(actx, dt) - else: - if constant_cfl: - dt = get_sim_timestep(dcoll, fluid_state, t, dt, current_cfl, - t_final, constant_cfl, local_dt) + t_final, constant_cfl) try: do_viz = check_step(step=step, interval=nviz) @@ -686,9 +675,7 @@ def my_pre_step(step, t, dt, state): t_elapsed = time.time() - t_start if t_shutdown - t_elapsed < 300.0: - my_write_restart(step=step, t=t, cv=fluid_state.cv, - tseed=tseed) - sys.exit() + my_write_restart(step=step, t=t, cv=fluid_state.cv, tseed=tseed) if do_garbage: with gc_timer.start_sub_timer(): @@ -704,8 +691,7 @@ def my_pre_step(step, t, dt, state): raise MyRuntimeError("Failed simulation health check.") if do_restart: - my_write_restart(step=step, t=t, cv=fluid_state.cv, - tseed=tseed) + my_write_restart(step=step, t=t, cv=fluid_state.cv, tseed=tseed) if do_viz: my_write_viz(step=step, t=t, dt=dt, state=fluid_state) @@ -765,31 +751,28 @@ def my_post_step(step, t, dt, state): "future GC collections") gc.freeze() - min_dt = np.min(actx.to_numpy(dt)) if local_dt else dt if logmgr: - set_dt(logmgr, min_dt) + set_dt(logmgr, dt) logmgr.tick_after() return state, dt # ############################################################################ - current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, - current_cfl, t_final, constant_cfl) - current_dt = force_evaluation(actx, current_dt) + if constant_cfl: + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, + current_cfl, t_final, constant_cfl) if rank == 0: logging.info("Stepping.") - (current_step, current_t, stepper_state) = \ + current_step, current_t, stepper_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, post_step_callback=my_post_step, istep=current_step, dt=current_dt, t=current_t, - t_final=t_final, max_steps=niter, local_dt=local_dt, - force_eval=force_eval_stepper, - state=make_obj_array([current_state.cv, tseed]) - ) + t_final=t_final, force_eval=force_eval_stepper, + state=make_obj_array([current_state.cv, tseed])) current_cv, tseed = stepper_state current_state = make_fluid_state(current_cv, gas_model, tseed) From 633c7182950c825e985f1610474a6dbe3c169684 Mon Sep 17 00:00:00 2001 From: Tulio Date: Mon, 26 Feb 2024 14:57:03 -0600 Subject: [PATCH 16/16] print -> logging --- examples/flame1d.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/flame1d.py b/examples/flame1d.py index 4de5f004c..c840edd12 100644 --- a/examples/flame1d.py +++ b/examples/flame1d.py @@ -298,7 +298,8 @@ def _compiled_stepper_wrapper(state, t, dt, rhs): # {{{ Set up initial state using Cantera # Use Cantera for initialization - print("\nUsing Cantera", cantera.__version__) + if rank == 0: + logging.info("\nUsing Cantera " + cantera.__version__) from mirgecom.mechanisms import get_mechanism_input mech_input = get_mechanism_input(mechanism_file)