diff --git a/.ci-support/production-setup.sh b/.ci-support/production-setup.sh index df547769a..73e7c30dc 100644 --- a/.ci-support/production-setup.sh +++ b/.ci-support/production-setup.sh @@ -1,7 +1,7 @@ #!/bin/bash set -x -PRODUCTION_CHANGE_OWNER="illinois-ceesd" +PRODUCTION_CHANGE_OWNER="" PRODUCTION_CHANGE_BRANCH="" CURRENT_BRANCH=$(git rev-parse --abbrev-ref HEAD) git config user.email "stupid@dumb.com" diff --git a/examples/autoignition-mpi.py b/examples/autoignition-mpi.py index 893b24e63..965054130 100644 --- a/examples/autoignition-mpi.py +++ b/examples/autoignition-mpi.py @@ -43,12 +43,12 @@ from mirgecom.euler import euler_operator from mirgecom.simutil import ( + get_sim_timestep, generate_and_distribute_mesh, write_visfile ) from mirgecom.io import make_init_message from mirgecom.mpi import mpi_entry_point - from mirgecom.integrators import rk4_step from mirgecom.steppers import advance_state from mirgecom.boundary import AdiabaticSlipBoundary @@ -144,12 +144,11 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - restarting = rst_filename is not None - if restarting: # read the grid from restart data - rst_fname = f"{rst_filename}-{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_fname) + 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"] @@ -262,7 +261,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, my_boundary = AdiabaticSlipBoundary() boundaries = {BTAG_ALL: my_boundary} - if restarting: + if rst_filename: current_step = rst_step current_t = rst_time if logmgr: @@ -315,14 +314,23 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, f" {eq_pressure=}, {eq_temperature=}," f" {eq_density=}, {eq_mass_fractions=}") - def my_write_viz(step, t, state, dv=None, production_rates=None): + def my_write_status(dt, cfl): + status_msg = f"------ {dt=}" if constant_cfl else f"----- {cfl=}" + if rank == 0: + logger.info(status_msg) + + def my_write_viz(step, t, dt, state, ts_field=None, dv=None, + production_rates=None, cfl=None): if dv is None: dv = eos.dependent_vars(state) if production_rates is None: production_rates = eos.get_production_rates(state) + if ts_field is None: + ts_field, cfl, dt = my_get_timestep(t=t, dt=dt, state=state) viz_fields = [("cv", state), ("dv", dv), - ("production_rates", production_rates)] + ("production_rates", production_rates), + ("dt" if constant_cfl else "cfl", ts_field)] write_visfile(discr, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True) @@ -358,6 +366,23 @@ def my_health_check(dv): return health_error + def my_get_timestep(t, dt, state): + # richer interface to calculate {dt,cfl} returns node-local estimates + t_remaining = max(0, t_final - t) + if constant_cfl: + from mirgecom.inviscid import get_inviscid_timestep + ts_field = current_cfl * get_inviscid_timestep(discr, eos=eos, cv=state) + from grudge.op import nodal_min + dt = nodal_min(discr, "vol", ts_field) + cfl = current_cfl + else: + from mirgecom.inviscid import get_inviscid_cfl + ts_field = get_inviscid_cfl(discr, eos=eos, dt=dt, cv=state) + from grudge.op import nodal_max + cfl = nodal_max(discr, "vol", ts_field) + + return ts_field, cfl, min(t_remaining, dt) + def my_pre_step(step, t, dt, state): try: dv = None @@ -369,6 +394,7 @@ def my_pre_step(step, t, dt, state): 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: dv = eos.dependent_vars(state) @@ -379,6 +405,11 @@ def my_pre_step(step, t, dt, state): logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") + ts_field, cfl, dt = my_get_timestep(t=t, dt=dt, state=state) + + if do_status: + my_write_status(dt, cfl) + if do_restart: my_write_restart(step=step, t=t, state=state) @@ -386,18 +417,18 @@ def my_pre_step(step, t, dt, state): production_rates = eos.get_production_rates(state) if dv is None: dv = eos.dependent_vars(state) - my_write_viz(step=step, t=t, state=state, dv=dv, - production_rates=production_rates) + my_write_viz(step=step, t=t, dt=dt, state=state, dv=dv, + production_rates=production_rates, + ts_field=ts_field, cfl=cfl) except MyRuntimeError: if rank == 0: logger.info("Errors detected; attempting graceful exit.") - my_write_viz(step=step, t=t, state=state) + my_write_viz(step=step, t=t, dt=dt, state=state) my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -413,6 +444,9 @@ def my_rhs(t, state): boundaries=boundaries, eos=eos) + eos.get_species_source_terms(state)) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, @@ -425,8 +459,11 @@ def my_rhs(t, state): final_dv = eos.dependent_vars(current_state) final_dm = eos.get_production_rates(current_state) - my_write_viz(step=current_step, t=current_t, state=current_state, dv=final_dv, - production_rates=final_dm) + ts_field, cfl, dt = my_get_timestep(t=current_t, dt=current_dt, + state=current_state) + my_write_viz(step=current_step, t=current_t, dt=dt, state=current_state, + dv=final_dv, production_rates=final_dm, ts_field=ts_field, cfl=cfl) + my_write_status(dt=dt, cfl=cfl) my_write_restart(step=current_step, t=current_t, state=current_state) if logmgr: diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py index 75eae69c4..1cc4553c2 100644 --- a/examples/lump-mpi.py +++ b/examples/lump-mpi.py @@ -37,7 +37,10 @@ from mirgecom.euler import euler_operator -from mirgecom.simutil import generate_and_distribute_mesh +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 @@ -69,8 +72,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_leap=False, - use_profiling=False, rst_step=None, rst_name=None, - casename="lump", use_logmgr=True): + use_profiling=False, rst_filename=None, casename="lump", + use_logmgr=True): """Drive example.""" cl_ctx = ctx_factory() @@ -96,40 +99,36 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, actx = PyOpenCLArrayContext(queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) - dim = 3 - nel_1d = 16 - order = 3 + # timestepping control + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step t_final = 0.01 current_cfl = 1.0 - vel = np.zeros(shape=(dim,)) - orig = np.zeros(shape=(dim,)) - vel[:dim] = 1.0 current_dt = .001 current_t = 0 - eos = IdealSingleGas() - initializer = Lump(dim=dim, center=orig, velocity=vel) - boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + current_step = 0 constant_cfl = False + + # some i/o frequencies nstatus = 1 nhealth = 1 nrestart = 10 nviz = 1 - current_step = 0 - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step + + dim = 3 rst_path = "restart_data/" rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) + 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_fname) + 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"] @@ -137,6 +136,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, else: # generate the grid from scratch box_ll = -5.0 box_ur = 5.0 + nel_1d = 16 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) @@ -144,6 +144,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, generate_mesh) local_nelements = local_mesh.nelements + order = 3 discr = EagerDGDiscretization( actx, local_mesh, order=order, mpi_communicator=comm ) @@ -167,9 +168,17 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - if rst_step: + # soln setup, init + eos = IdealSingleGas() + vel = np.zeros(shape=(dim,)) + orig = np.zeros(shape=(dim,)) + vel[:dim] = 1.0 + initializer = Lump(dim=dim, center=orig, velocity=vel) + boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -208,17 +217,18 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): def my_write_restart(state, step, t): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "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) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "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, state, exact): health_error = False @@ -265,10 +275,6 @@ def my_pre_step(step, t, dt, state): logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed solution health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) @@ -299,8 +305,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -315,6 +322,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py index eb1ec13ea..5de084950 100644 --- a/examples/mixture-mpi.py +++ b/examples/mixture-mpi.py @@ -37,7 +37,10 @@ from mirgecom.euler import euler_operator -from mirgecom.simutil import generate_and_distribute_mesh +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 @@ -72,8 +75,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_leap=False, - use_profiling=False, rst_step=None, rst_name=None, - casename="uiuc_mixture", use_logmgr=True): + use_profiling=False, rst_filename=None, casename="uiuc_mixture", + use_logmgr=True): """Drive example.""" cl_ctx = ctx_factory() @@ -99,41 +102,40 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, actx = PyOpenCLArrayContext(queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) - dim = 3 - nel_1d = 16 - order = 3 + # timestepping control + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step t_final = 0.002 current_cfl = 1.0 - velocity = np.zeros(shape=(dim,)) - velocity[:dim] = 1.0 current_dt = .001 current_t = 0 + current_step = 0 constant_cfl = False + + # some i/o frequencies nstatus = 1 nhealth = 1 nrestart = 5 nviz = 1 - current_step = 0 - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step + dim = 3 rst_path = "restart_data/" rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) - + 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_fname) + 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["nparts"] == 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 @@ -143,6 +145,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, generate_mesh) local_nelements = local_mesh.nelements + order = 3 discr = EagerDGDiscretization( actx, local_mesh, order=order, mpi_communicator=comm ) @@ -182,14 +185,16 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, y0s[nspecies-1] = 1.0 - spec_sum # Mixture defaults to STP (p, T) = (1atm, 300K) + velocity = np.zeros(shape=(dim,)) + velocity[:dim] = 1.0 initializer = MixtureInitializer(dim=dim, nspecies=nspecies, massfractions=y0s, velocity=velocity) boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} nodes = thaw(actx, discr.nodes()) - if rst_step: + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -211,6 +216,13 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, if rank == 0: logger.info(init_message) + def my_write_status(component_errors): + status_msg = ( + "------- errors=" + + ", ".join("%.3g" % en for en in component_errors)) + if rank == 0: + logger.info(status_msg) + def my_write_viz(step, t, state, dv=None, exact=None, resid=None): viz_fields = [("cv", state)] if dv is None: @@ -229,19 +241,20 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): def my_write_restart(step, t, state): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "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(state, dv, exact): + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "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(discr, "vol", dv.pressure) \ @@ -249,8 +262,6 @@ def my_health_check(state, dv, exact): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) exittol = .09 if max(component_errors) > exittol: health_error = True @@ -263,6 +274,7 @@ def my_pre_step(step, t, dt, state): try: dv = None exact = None + component_errors = None if logmgr: logmgr.tick_before() @@ -276,18 +288,16 @@ def my_pre_step(step, t, dt, state): if do_health: dv = eos.dependent_vars(state) exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) from mirgecom.simutil import allsync - health_errors = allsync(my_health_check(state, dv, exact), comm, + health_errors = allsync(my_health_check(dv, component_errors), comm, op=MPI.LOR) if health_errors: if rank == 0: logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) @@ -301,15 +311,12 @@ def my_pre_step(step, t, dt, state): resid=resid) if do_status: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, t=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) - status_msg = ( - "------- errors=" - + ", ".join("%.3g" % en for en in component_errors)) - if rank == 0: - logger.info(status_msg) + if component_errors is None: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) + my_write_status(component_errors) except MyRuntimeError: if rank == 0: @@ -318,8 +325,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -334,6 +342,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, diff --git a/examples/pulse-mpi.py b/examples/pulse-mpi.py index 7afdda6d3..3831897b4 100644 --- a/examples/pulse-mpi.py +++ b/examples/pulse-mpi.py @@ -38,7 +38,10 @@ from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator -from mirgecom.simutil import generate_and_distribute_mesh +from mirgecom.simutil import ( + get_sim_timestep, + generate_and_distribute_mesh +) from mirgecom.io import make_init_message from mirgecom.integrators import rk4_step @@ -76,7 +79,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_logmgr=True, use_leap=False, use_profiling=False, casename="pulse", - rst_step=None, rst_name=None): + rst_filename=None): """Drive the example.""" cl_ctx = ctx_factory() @@ -102,41 +105,34 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, actx = PyOpenCLArrayContext(queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) - dim = 2 - nel_1d = 16 - order = 1 + # timestepping control + current_step = 0 + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step t_final = 0.1 current_cfl = 1.0 - vel = np.zeros(shape=(dim,)) - orig = np.zeros(shape=(dim,)) current_dt = .01 current_t = 0 - eos = IdealSingleGas() - initializer = Lump(dim=dim, center=orig, velocity=vel, rhoamp=0.0) - boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} - wall = AdiabaticSlipBoundary() - boundaries = {BTAG_ALL: wall} constant_cfl = False + + # some i/o frequencies nstatus = 1 nrestart = 5 nviz = 10 nhealth = 1 - current_step = 0 - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step + dim = 2 rst_path = "restart_data/" rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) - + 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_fname) + 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"] @@ -145,12 +141,14 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, from meshmode.mesh.generation import generate_regular_rect_mesh box_ll = -0.5 box_ur = 0.5 + nel_1d = 16 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 = 1 discr = EagerDGDiscretization( actx, local_mesh, order=order, mpi_communicator=comm ) @@ -171,12 +169,19 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, ("t_log.max", "log walltime: {value:6g} s") ]) + eos = IdealSingleGas() + vel = np.zeros(shape=(dim,)) + orig = np.zeros(shape=(dim,)) + initializer = Lump(dim=dim, center=orig, velocity=vel, rhoamp=0.0) + boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + wall = AdiabaticSlipBoundary() + boundaries = {BTAG_ALL: wall} uniform_state = initializer(nodes) acoustic_pulse = AcousticPulse(dim=dim, amplitude=1.0, width=.1, center=orig) - if rst_step: + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -210,17 +215,18 @@ def my_write_viz(step, t, state, dv=None): def my_write_restart(step, t, state): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "t": t, - "step": step, - "order": order, - "global_nelements": global_nelements, - "num_parts": num_parts - } - from mirgecom.restart import write_restart_file - write_restart_file(actx, rst_data, rst_fname, comm) + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": num_parts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) def my_health_check(pressure): health_error = False @@ -253,10 +259,6 @@ def my_pre_step(step, t, dt, state): logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) @@ -272,8 +274,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -288,6 +291,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, diff --git a/examples/run_examples.sh b/examples/run_examples.sh index c6c195a99..71a57d69d 100755 --- a/examples/run_examples.sh +++ b/examples/run_examples.sh @@ -7,6 +7,7 @@ origin=$(pwd) examples_dir=${1-$origin} declare -i exitcode=0 echo "Running examples in $examples_dir ..." +failed_examples="" for example in $examples_dir/*.py do if [[ "$example" == *"-mpi.py" ]] @@ -23,6 +24,7 @@ do else ((exitcode=exitcode+1)) echo "Example $example failed." + failed_examples="$failed_examples $example" fi done echo "Done running examples!" @@ -30,7 +32,7 @@ if [[ $exitcode -eq 0 ]] then echo "No errors." else - echo "Errors detected ($exitcode)." + echo "Errors detected ($exitcode):($failed_examples )" exit $exitcode fi #rm -f examples/*.vtu diff --git a/examples/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py index c7ce42fdc..067ff22ea 100644 --- a/examples/scalar-lump-mpi.py +++ b/examples/scalar-lump-mpi.py @@ -38,7 +38,10 @@ from mirgecom.euler import euler_operator -from mirgecom.simutil import generate_and_distribute_mesh +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 @@ -70,7 +73,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_leap=False, - use_profiling=False, rst_step=None, rst_name=None, + use_profiling=False, rst_filename=None, casename="lumpy-scalars", use_logmgr=True): """Drive example.""" cl_ctx = ctx_factory() @@ -97,50 +100,34 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, actx = PyOpenCLArrayContext(queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) - dim = 3 - nel_1d = 16 - order = 3 + # timestepping control + current_step = 0 + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step t_final = 0.005 current_cfl = 1.0 current_dt = .001 current_t = 0 constant_cfl = False + + # some i/o frequencies nstatus = 1 nrestart = 5 nviz = 1 nhealth = 1 - current_step = 0 - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step - box_ll = -5.0 - box_ur = 5.0 - - nspecies = 4 - centers = make_obj_array([np.zeros(shape=(dim,)) for i in range(nspecies)]) - spec_y0s = np.ones(shape=(nspecies,)) - spec_amplitudes = np.ones(shape=(nspecies,)) - eos = IdealSingleGas() - - velocity = np.ones(shape=(dim,)) - - initializer = MulticomponentLump(dim=dim, nspecies=nspecies, - spec_centers=centers, velocity=velocity, - spec_y0s=spec_y0s, - spec_amplitudes=spec_amplitudes) - boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + dim = 3 rst_path = "restart_data/" rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) - + 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_fname) + 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"] @@ -148,6 +135,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, else: # generate the grid from scratch box_ll = -5.0 box_ur = 5.0 + nel_1d = 16 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) @@ -155,6 +143,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, generate_mesh) local_nelements = local_mesh.nelements + order = 3 discr = EagerDGDiscretization( actx, local_mesh, order=order, mpi_communicator=comm ) @@ -178,9 +167,22 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - if rst_step: + # soln setup and init + nspecies = 4 + centers = make_obj_array([np.zeros(shape=(dim,)) for i in range(nspecies)]) + spec_y0s = np.ones(shape=(nspecies,)) + spec_amplitudes = np.ones(shape=(nspecies,)) + eos = IdealSingleGas() + velocity = np.ones(shape=(dim,)) + + initializer = MulticomponentLump(dim=dim, nspecies=nspecies, + spec_centers=centers, velocity=velocity, + spec_y0s=spec_y0s, + spec_amplitudes=spec_amplitudes) + boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -202,6 +204,12 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, if rank == 0: logger.info(init_message) + def my_write_status(component_errors): + if rank == 0: + logger.info( + "------- errors=" + + ", ".join("%.3g" % en for en in component_errors)) + def my_write_viz(step, t, state, dv=None, exact=None, resid=None): if dv is None: dv = eos.dependent_vars(state) @@ -219,19 +227,20 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): def my_write_restart(step, t, state): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "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(state, pressure, exact): + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "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(pressure, component_errors): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local if check_naninf_local(discr, "vol", pressure) \ @@ -239,8 +248,6 @@ def my_health_check(state, pressure, exact): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) exittol = .09 if max(component_errors) > exittol: health_error = True @@ -253,6 +260,7 @@ def my_pre_step(step, t, dt, state): try: dv = None exact = None + component_errors = None if logmgr: logmgr.tick_before() @@ -266,18 +274,18 @@ def my_pre_step(step, t, dt, state): if do_health: dv = eos.dependent_vars(state) exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) from mirgecom.simutil import allsync - health_errors = allsync(my_health_check(state, dv.pressure, exact), - comm, op=MPI.LOR) + health_errors = allsync( + my_health_check(dv.pressure, component_errors), + comm, op=MPI.LOR + ) if health_errors: if rank == 0: logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) @@ -291,15 +299,12 @@ def my_pre_step(step, t, dt, state): resid=resid) if do_status: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, t=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) - status_msg = ( - "------- errors=" - + ", ".join("%.3g" % en for en in component_errors)) - if rank == 0: - logger.info(status_msg) + if component_errors is None: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) + my_write_status(component_errors) except MyRuntimeError: if rank == 0: @@ -308,8 +313,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -324,6 +330,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, dt=current_dt, diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py index 9e3165f4e..2e4bb350f 100644 --- a/examples/sod-mpi.py +++ b/examples/sod-mpi.py @@ -37,7 +37,10 @@ from mirgecom.euler import euler_operator -from mirgecom.simutil import generate_and_distribute_mesh +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 @@ -70,7 +73,7 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_logmgr=True, use_leap=False, use_profiling=False, casename="sod1d", - rst_step=None, rst_name=None): + rst_filename=None): """Drive the example.""" cl_ctx = ctx_factory() @@ -97,7 +100,6 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) dim = 1 - nel_1d = 24 order = 1 # tolerate large errors; case is unstable t_final = 0.01 @@ -123,17 +125,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) - + 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_fname) + 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["nparts"] == num_parts else: # generate the grid from scratch from meshmode.mesh.generation import generate_regular_rect_mesh + nel_1d = 24 box_ll = -5.0 box_ur = 5.0 generate_mesh = partial(generate_regular_rect_mesh, a=(box_ll,)*dim, @@ -162,9 +164,9 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, ("t_log.max", "log walltime: {value:6g} s") ]) - if rst_step: + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -187,6 +189,13 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, if rank == 0: logger.info(init_message) + def my_write_status(component_errors): + if rank == 0: + logger.info( + "------- errors=" + + ", ".join("%.3g" % en for en in component_errors) + ) + def my_write_viz(step, t, state, dv=None, exact=None, resid=None): if dv is None: dv = eos.dependent_vars(state) @@ -204,19 +213,20 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): def my_write_restart(state, step, t): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "t": t, - "step": step, - "order": order, - "global_nelements": global_nelements, - "num_parts": num_parts - } - from mirgecom.restart import write_restart_file - write_restart_file(actx, rst_data, rst_fname, comm) - - def my_health_check(state, pressure, exact): + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": num_parts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(pressure, component_errors): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local if check_naninf_local(discr, "vol", pressure) \ @@ -224,8 +234,6 @@ def my_health_check(state, pressure, exact): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) exittol = .09 if max(component_errors) > exittol: health_error = True @@ -238,6 +246,7 @@ def my_pre_step(step, t, dt, state): try: dv = None exact = None + component_errors = None if logmgr: logmgr.tick_before() @@ -251,18 +260,18 @@ def my_pre_step(step, t, dt, state): if do_health: dv = eos.dependent_vars(state) exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) from mirgecom.simutil import allsync - health_errors = allsync(my_health_check(state, dv.pressure, exact), - comm, op=MPI.LOR) + health_errors = allsync( + my_health_check(dv.pressure, component_errors), + comm, op=MPI.LOR + ) if health_errors: if rank == 0: logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) @@ -276,15 +285,13 @@ def my_pre_step(step, t, dt, state): resid=resid) if do_status: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, t=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) - status_msg = ( - "------- errors=" - + ", ".join("%.3g" % en for en in component_errors)) - if rank == 0: - logger.info(status_msg) + if component_errors is None: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = \ + compare_fluid_solutions(discr, state, exact) + my_write_status(component_errors) except MyRuntimeError: if rank == 0: @@ -293,8 +300,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -309,6 +317,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, dt=current_dt, diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index 17ab5a830..356f5c99b 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -39,6 +39,7 @@ from mirgecom.euler import euler_operator from mirgecom.simutil import ( + get_sim_timestep, generate_and_distribute_mesh, check_step ) @@ -73,8 +74,8 @@ class MyRuntimeError(RuntimeError): @mpi_entry_point def main(ctx_factory=cl.create_some_context, use_logmgr=True, - use_leap=False, use_profiling=False, casename="sod1d", - rst_step=None, rst_name=None): + use_leap=False, use_profiling=False, rst_filename=None, + casename="vortex"): """Drive the example.""" cl_ctx = ctx_factory() @@ -100,32 +101,26 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, actx = PyOpenCLArrayContext(queue, allocator=cl_tools.MemoryPool(cl_tools.ImmediateAllocator(queue))) - dim = 2 - nel_1d = 16 - order = 3 + # timestepping control + current_step = 0 + if use_leap: + from leap.rk import RK4MethodBuilder + timestepper = RK4MethodBuilder("state") + else: + timestepper = rk4_step t_final = 0.01 current_cfl = 1.0 - vel = np.zeros(shape=(dim,)) - orig = np.zeros(shape=(dim,)) - vel[:dim] = 1.0 current_dt = .001 current_t = 0 - eos = IdealSingleGas() - initializer = Vortex2D(center=orig, velocity=vel) - casename = "vortex" - boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} constant_cfl = False + + # some i/o frequencies nrestart = 10 nstatus = 1 nviz = 10 nhealth = 10 - current_step = 0 - if use_leap: - from leap.rk import RK4MethodBuilder - timestepper = RK4MethodBuilder("state") - else: - timestepper = rk4_step + dim = 2 if dim != 2: raise ValueError("This example must be run with dim = 2.") @@ -133,25 +128,26 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, rst_pattern = ( rst_path + "{cname}-{step:04d}-{rank:04d}.pkl" ) - if rst_step: # read the grid from restart data - rst_fname = rst_pattern.format(cname=rst_name, step=rst_step, rank=rank) - + 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_fname) + 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["nparts"] == num_parts else: # generate the grid from scratch - from meshmode.mesh.generation import generate_regular_rect_mesh + 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 discr = EagerDGDiscretization( actx, local_mesh, order=order, mpi_communicator=comm ) @@ -185,9 +181,16 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - if rst_step: + # soln setup and init + eos = IdealSingleGas() + vel = np.zeros(shape=(dim,)) + orig = np.zeros(shape=(dim,)) + vel[:dim] = 1.0 + initializer = Vortex2D(center=orig, velocity=vel) + boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} + if rst_filename: current_t = restart_data["t"] - current_step = rst_step + current_step = restart_data["step"] current_state = restart_data["state"] if logmgr: from mirgecom.logging_quantities import logmgr_set_time @@ -210,6 +213,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, if rank == 0: logger.info(init_message) + def my_write_status(state, component_errors, cfl=None): + if cfl is None: + from mirgecom.inviscid import get_inviscid_cfl + cfl = current_cfl if constant_cfl else \ + get_inviscid_cfl(discr, eos, current_dt, state) + if rank == 0: + logger.info( + f"------ {cfl=}\n" + "------- errors=" + + ", ".join("%.3g" % en for en in component_errors)) + def my_write_viz(step, t, state, dv=None, exact=None, resid=None): if dv is None: dv = eos.dependent_vars(state) @@ -227,19 +241,20 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): def my_write_restart(step, t, state): rst_fname = rst_pattern.format(cname=casename, step=step, rank=rank) - rst_data = { - "local_mesh": local_mesh, - "state": state, - "t": t, - "step": step, - "order": order, - "global_nelements": global_nelements, - "num_parts": num_parts - } - from mirgecom.restart import write_restart_file - write_restart_file(actx, rst_data, rst_fname, comm) - - def my_health_check(state, pressure, exact): + if rst_fname != rst_filename: + rst_data = { + "local_mesh": local_mesh, + "state": state, + "t": t, + "step": step, + "order": order, + "global_nelements": global_nelements, + "num_parts": num_parts + } + from mirgecom.restart import write_restart_file + write_restart_file(actx, rst_data, rst_fname, comm) + + def my_health_check(pressure, component_errors): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local if check_naninf_local(discr, "vol", pressure) \ @@ -247,8 +262,6 @@ def my_health_check(state, pressure, exact): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) exittol = .1 if max(component_errors) > exittol: health_error = True @@ -261,6 +274,7 @@ def my_pre_step(step, t, dt, state): try: dv = None exact = None + component_errors = None if logmgr: logmgr.tick_before() @@ -273,21 +287,29 @@ def my_pre_step(step, t, dt, state): if do_health: dv = eos.dependent_vars(state) exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) from mirgecom.simutil import allsync - health_errors = allsync(my_health_check(state, dv.pressure, exact), - comm, op=MPI.LOR) + health_errors = allsync( + my_health_check(dv.pressure, component_errors), + comm, op=MPI.LOR + ) if health_errors: if rank == 0: logger.info("Fluid solution failed health check.") raise MyRuntimeError("Failed simulation health check.") - if step == rst_step: # don't do viz or restart @ restart - do_viz = False - do_restart = False - if do_restart: my_write_restart(step=step, t=t, state=state) + if do_status: + if component_errors is None: + if exact is None: + exact = initializer(x_vec=nodes, eos=eos, t=t) + from mirgecom.simutil import compare_fluid_solutions + component_errors = compare_fluid_solutions(discr, state, exact) + my_write_status(state, component_errors) + if do_viz: if dv is None: dv = eos.dependent_vars(state) @@ -297,17 +319,6 @@ def my_pre_step(step, t, dt, state): my_write_viz(step=step, t=t, state=state, dv=dv, exact=exact, resid=resid) - if do_status: - if exact is None: - exact = initializer(x_vec=nodes, eos=eos, t=t) - from mirgecom.simutil import compare_fluid_solutions - component_errors = compare_fluid_solutions(discr, state, exact) - status_msg = ( - "------- errors=" - + ", ".join("%.3g" % en for en in component_errors)) - if rank == 0: - logger.info(status_msg) - except MyRuntimeError: if rank == 0: logger.info("Errors detected; attempting graceful exit.") @@ -315,8 +326,9 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - t_remaining = max(0, t_final - t) - return state, min(dt, t_remaining) + dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + constant_cfl) + return state, dt def my_post_step(step, t, dt, state): # Logmgr needs to know about EOS, dt, dim? @@ -331,6 +343,9 @@ def my_rhs(t, state): return euler_operator(discr, cv=state, t=t, boundaries=boundaries, eos=eos) + current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_cfl, eos, t_final, constant_cfl) + current_step, current_t, current_state = \ advance_state(rhs=my_rhs, timestepper=timestepper, pre_step_callback=my_pre_step, diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index 0f67f1f60..e3fe348d2 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -66,30 +66,51 @@ def inviscid_flux(discr, eos, cv): (mom / cv.mass) * cv.species_mass.reshape(-1, 1))) -def get_inviscid_timestep(discr, eos, cfl, cv): - """Routine (will) return the (local) maximum stable inviscid timestep. - - Currently, it's a hack waiting for the geometric_factor helpers port - from grudge. +def get_inviscid_timestep(discr, eos, cv): + """Return node-local stable timestep estimate for an inviscid fluid. + + The maximum stable timestep is computed from the acoustic wavespeed. + + Parameters + ---------- + discr: grudge.eager.EagerDGDiscretization + the discretization to use + eos: mirgecom.eos.GasEOS + Implementing the pressure and temperature functions for + returning pressure and temperature as a function of the state q. + cv: :class:`~mirgecom.fluid.ConservedVars` + Fluid soluition + Returns + ------- + class:`~meshmode.dof_array.DOFArray` + The maximum stable timestep at each node. """ - dim = discr.dim - mesh = discr.mesh - order = max([grp.order for grp in discr.discr_from_dd("vol").groups]) - nelements = mesh.nelements - nel_1d = nelements ** (1.0 / (1.0 * dim)) - - # This roughly reproduces the timestep AK used in wave toy - dt = (1.0 - 0.25 * (dim - 1)) / (nel_1d * order ** 2) - return cfl * dt - -# dt_ngf = dt_non_geometric_factor(discr.mesh) -# dt_gf = dt_geometric_factor(discr.mesh) -# wavespeeds = compute_wavespeed(w,eos=eos) -# max_v = clmath.max(wavespeeds) -# return c*dt_ngf*dt_gf/max_v + from grudge.dt_utils import characteristic_lengthscales + from mirgecom.fluid import compute_wavespeed + return ( + characteristic_lengthscales(cv.array_context, discr) + / compute_wavespeed(discr, eos, cv) + ) def get_inviscid_cfl(discr, eos, dt, cv): - """Calculate and return CFL based on current state and timestep.""" - wanted_dt = get_inviscid_timestep(discr, eos=eos, cfl=1.0, cv=cv) - return dt / wanted_dt + """Return node-local CFL based on current state and timestep. + + Parameters + ---------- + discr: :class:`grudge.eager.EagerDGDiscretization` + the discretization to use + eos: mirgecom.eos.GasEOS + Implementing the pressure and temperature functions for + returning pressure and temperature as a function of the state q. + dt: float or :class:`~meshmode.dof_array.DOFArray` + A constant scalar dt or node-local dt + cv: :class:`~mirgecom.fluid.ConservedVars` + Fluid solution + + Returns + ------- + :class:`meshmode.dof_array.DOFArray` + The CFL at each node. + """ + return dt / get_inviscid_timestep(discr, eos=eos, cv=cv) diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index ac6219fa7..3b0a03c1e 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -4,7 +4,7 @@ ----------------- .. autofunction:: check_step -.. autofunction:: inviscid_sim_timestep +.. autofunction:: get_sim_timestep .. autofunction:: write_visfile .. autofunction:: allsync @@ -48,7 +48,6 @@ import logging import numpy as np -from mirgecom.inviscid import get_inviscid_timestep # bad smell? import grudge.op as op logger = logging.getLogger(__name__) @@ -75,14 +74,58 @@ def check_step(step, interval): return False -def inviscid_sim_timestep(discr, state, t, dt, cfl, eos, - t_final, constant_cfl=False): - """Return the maximum stable dt.""" +def get_sim_timestep(discr, state, t, dt, cfl, eos, + t_final, constant_cfl=False): + """Return the maximum stable timestep for a typical fluid simulation. + + This routine returns *dt*, the users defined constant timestep, or + *max_dt*, the maximum domain-wide stability-limited + timestep for a fluid simulation. It calls the collective: + :func:`~grudge.op.nodal_min` on the inside which makes it + domain-wide regardless of parallel decomposition. + + Two modes are supported: + - Constant DT mode: returns the minimum of (t_final-t, dt) + - Constant CFL mode: returns (cfl * max_dt) + + .. important:: + The current implementation is calculating an acoustic-limited + timestep and CFL for an inviscid fluid. The addition of viscous + fluxes includes modification to this routine. + + Parameters + ---------- + discr + Grudge discretization or discretization collection? + state: :class:`~mirgecom.fluid.ConservedVars` + The fluid state. + t: float + Current time + t_final: float + Final time + dt: float + The current timestep + cfl: float + The current CFL number + eos: :class:`~mirgecom.eos.GasEOS` + Gas equation-of-state supporting speed_of_sound + constant_cfl: bool + True if running constant CFL mode + + Returns + ------- + float + The maximum stable DT based on inviscid fluid acoustic wavespeed. + """ mydt = dt t_remaining = max(0, t_final - t) - if constant_cfl is True: - mydt = get_inviscid_timestep(discr=discr, cv=state, - cfl=cfl, eos=eos) + if constant_cfl: + from mirgecom.inviscid import get_inviscid_timestep + from grudge.op import nodal_min + mydt = cfl * nodal_min( + discr, "vol", + get_inviscid_timestep(discr=discr, eos=eos, cv=state) + ) return min(t_remaining, mydt) diff --git a/test/test_euler.py b/test/test_euler.py index 3a16c77f3..243917c28 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -721,8 +721,9 @@ def _euler_flow_stepper(actx, parameters): discr = EagerDGDiscretization(actx, mesh, order=order) nodes = thaw(actx, discr.nodes()) + cv = initializer(nodes) - sdt = get_inviscid_timestep(discr, eos=eos, cfl=cfl, cv=cv) + sdt = cfl * get_inviscid_timestep(discr, eos=eos, cv=cv) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ @@ -802,7 +803,7 @@ def rhs(t, q): t += dt istep += 1 - sdt = get_inviscid_timestep(discr, eos=eos, cfl=cfl, cv=cv) + sdt = cfl * get_inviscid_timestep(discr, eos=eos, cv=cv) if nstepstatus > 0: logger.info("Writing final dump.")