diff --git a/doc/index.rst b/doc/index.rst index dd00a3115..7456de10c 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -15,12 +15,17 @@ Here’s an example, to give you an impression: import numpy as np import pyopencl as cl + + from arraycontext import PyOpenCLArrayContext + from pytools.obj_array import flat_obj_array - from grudge.eager import EagerDGDiscretization + + from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer + import grudge.op as op + from mirgecom.wave import wave_operator from mirgecom.integrators import rk4_step - from meshmode.array_context import PyOpenCLArrayContext cl_ctx = cl.create_some_context() queue = cl.CommandQueue(cl_ctx) @@ -36,15 +41,15 @@ Here’s an example, to give you an impression: print("%d elements" % mesh.nelements) - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) fields = flat_obj_array( - [discr.zeros(actx)], - [discr.zeros(actx) for i in range(discr.dim)] + [dcoll.zeros(actx)], + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) - vis = make_visualizer(discr, order + 3) + vis = make_visualizer(dcoll, order + 3) def rhs(t, w): - return wave_operator(discr, c=1, w=w) + return wave_operator(dcoll, c=1, w=w) t = 0 t_final = 3 @@ -53,7 +58,7 @@ Here’s an example, to give you an impression: while t < t_final: fields = rk4_step(fields, t, dt, rhs) if istep % 10 == 0: - print(istep, t, discr.norm(fields[0], np.inf)) + print(istep, t, op.norm(dcoll, fields[0], np.inf)) vis.write_vtk_file("wave-eager-%04d.vtu" % istep, [("u", fields[0]), ("v", fields[1:]), ]) diff --git a/examples/autoignition-mpi.py b/examples/autoignition-mpi.py index a1af62fa5..443c63c44 100644 --- a/examples/autoignition-mpi.py +++ b/examples/autoignition-mpi.py @@ -29,13 +29,12 @@ import pyopencl.tools as cl_tools from functools import partial -from meshmode.array_context import PyOpenCLArrayContext +from arraycontext import thaw, PyOpenCLArrayContext -from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization -from grudge.shortcuts import make_visualizer +from grudge.discretization import DiscretizationCollection +from grudge.shortcuts import make_visualizer from logpyle import IntervalTimer, set_dt from mirgecom.euler import extract_vars_for_logging, units_for_logging @@ -166,17 +165,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, generate_mesh) local_nelements = local_mesh.nelements - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -273,10 +272,12 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, current_state = restart_data["state"] else: rst_state = restart_data["state"] - old_discr = EagerDGDiscretization(actx, local_mesh, order=rst_order, - mpi_communicator=comm) + old_discr = DiscretizationCollection( + actx, local_mesh, order=rst_order, + mpi_communicator=comm + ) from meshmode.discretization.connection import make_same_mesh_connection - connection = make_same_mesh_connection(actx, discr.discr_from_dd("vol"), + connection = make_same_mesh_connection(actx, dcoll.discr_from_dd("vol"), old_discr.discr_from_dd("vol")) current_state = connection(rst_state) else: @@ -292,7 +293,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=False, # }}} - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ init_message = make_init_message(dim=dim, order=order, @@ -333,7 +334,7 @@ def my_write_viz(step, t, dt, state, ts_field=None, dv=None, ("dv", dv), ("production_rates", production_rates), ("dt" if constant_cfl else "cfl", ts_field)] - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(step, t, state): @@ -357,12 +358,12 @@ def my_write_restart(step, t, state): def my_health_check(dv): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local - if check_naninf_local(discr, "vol", dv.pressure) \ - or check_range_local(discr, "vol", dv.pressure, 1e5, 2.4e5): + if check_naninf_local(dcoll, "vol", dv.pressure) \ + or check_range_local(dcoll, "vol", dv.pressure, 1e5, 2.4e5): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") - if check_range_local(discr, "vol", dv.temperature, 1.498e3, 1.52e3): + if check_range_local(dcoll, "vol", dv.temperature, 1.498e3, 1.52e3): health_error = True logger.info(f"{rank=}: Invalid temperature data found.") @@ -373,15 +374,15 @@ def my_get_timestep(t, dt, state): 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) + ts_field = current_cfl * get_inviscid_timestep(dcoll, eos=eos, cv=state) from grudge.op import nodal_min - dt = nodal_min(discr, "vol", ts_field) + dt = nodal_min(dcoll, "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) + ts_field = get_inviscid_cfl(dcoll, eos=eos, dt=dt, cv=state) from grudge.op import nodal_max - cfl = nodal_max(discr, "vol", ts_field) + cfl = nodal_max(dcoll, "vol", ts_field) return ts_field, cfl, min(t_remaining, dt) @@ -442,11 +443,11 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return (euler_operator(discr, cv=state, t=t, + return (euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) + eos.get_species_source_terms(state)) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/heat-source-mpi.py b/examples/heat-source-mpi.py index 66643596d..dcdc4c1ca 100644 --- a/examples/heat-source-mpi.py +++ b/examples/heat-source-mpi.py @@ -25,21 +25,25 @@ import numpy.linalg as la # noqa import pyopencl as cl -from meshmode.array_context import thaw, PyOpenCLArrayContext +from arraycontext import thaw, PyOpenCLArrayContext from mirgecom.profiling import PyOpenCLProfilingArrayContext from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization +from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer from grudge.dof_desc import DISCR_TAG_BASE, DTAG_BOUNDARY +import grudge.op as op + from mirgecom.integrators import rk4_step from mirgecom.diffusion import ( diffusion_operator, DirichletDiffusionBoundary, NeumannDiffusionBoundary) + from mirgecom.mpi import mpi_entry_point + import pyopencl.tools as cl_tools from mirgecom.logging_quantities import (initialize_logmgr, @@ -107,7 +111,7 @@ def main(use_profiling=False, use_logmgr=False): order = 3 - discr = EagerDGDiscretization(actx, local_mesh, order=order, + dcoll = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) if dim == 2: @@ -118,14 +122,14 @@ def main(use_profiling=False, use_logmgr=False): source_width = 0.2 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) boundaries = { DTAG_BOUNDARY("dirichlet"): DirichletDiffusionBoundary(0.), DTAG_BOUNDARY("neumann"): NeumannDiffusionBoundary(0.) } - u = discr.zeros(actx) + u = dcoll.zeros(actx) if logmgr: logmgr_add_device_name(logmgr, queue) @@ -144,12 +148,12 @@ def main(use_profiling=False, use_logmgr=False): vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) def rhs(t, u): return ( diffusion_operator( - discr, quad_tag=DISCR_TAG_BASE, + dcoll, quad_tag=DISCR_TAG_BASE, alpha=1, boundaries=boundaries, u=u) + actx.np.exp(-np.dot(nodes, nodes)/source_width**2)) @@ -160,7 +164,7 @@ def rhs(t, u): logmgr.tick_before() if istep % 10 == 0: - print(istep, t, discr.norm(u)) + print(istep, t, op.norm(dcoll, u, 2)) vis.write_vtk_file("fld-heat-source-mpi-%03d-%04d.vtu" % (rank, istep), [ ("u", u) diff --git a/examples/lump-mpi.py b/examples/lump-mpi.py index 20aec5bac..b3a5c28eb 100644 --- a/examples/lump-mpi.py +++ b/examples/lump-mpi.py @@ -29,12 +29,12 @@ import pyopencl.tools as cl_tools from functools import partial -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization -from grudge.shortcuts import make_visualizer +from grudge.discretization import DiscretizationCollection +from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator from mirgecom.simutil import ( @@ -145,17 +145,17 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, local_nelements = local_mesh.nelements order = 3 - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -189,7 +189,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, # Set the current state from time 0 current_state = initializer(nodes) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ init_message = make_init_message(dim=dim, order=order, @@ -214,7 +214,7 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): ("exact", exact), ("residual", resid)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(state, step, t): @@ -235,13 +235,13 @@ def my_write_restart(state, step, t): def my_health_check(dv, state, exact): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local - if check_naninf_local(discr, "vol", dv.pressure) \ - or check_range_local(discr, "vol", dv.pressure, .9, 1.1): + if check_naninf_local(dcoll, "vol", dv.pressure) \ + or check_range_local(dcoll, "vol", dv.pressure, .9, 1.1): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) exittol = .09 if max(component_errors) > exittol: health_error = True @@ -293,7 +293,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) status_msg = ( "------- errors=" + ", ".join("%.3g" % en for en in component_errors)) @@ -307,7 +307,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -321,10 +321,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/mixture-mpi.py b/examples/mixture-mpi.py index a8485f7d4..e8efe9f1a 100644 --- a/examples/mixture-mpi.py +++ b/examples/mixture-mpi.py @@ -29,12 +29,12 @@ import pyopencl.tools as cl_tools from functools import partial -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization -from grudge.shortcuts import make_visualizer +from grudge.discretization import DiscretizationCollection +from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator from mirgecom.simutil import ( @@ -146,17 +146,17 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, local_nelements = local_mesh.nelements order = 3 - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -193,7 +193,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, massfractions=y0s, velocity=velocity) boundaries = {BTAG_ALL: PrescribedBoundary(initializer)} - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) if rst_filename: current_t = restart_data["t"] current_step = restart_data["step"] @@ -205,7 +205,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, # Set the current state from time 0 current_state = initializer(x_vec=nodes, eos=eos) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ init_message = make_init_message(dim=dim, order=order, @@ -238,7 +238,7 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): ("exact_soln", exact), ("residual", resid)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(step, t, state): @@ -259,8 +259,8 @@ def my_write_restart(step, t, state): 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) \ - or check_range_local(discr, "vol", dv.pressure, 1e5, 1.1e5): + 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.") @@ -291,7 +291,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) from mirgecom.simutil import allsync health_errors = allsync(my_health_check(dv, component_errors), comm, op=MPI.LOR) @@ -317,7 +317,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) my_write_status(component_errors) except MyRuntimeError: @@ -327,7 +327,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -341,10 +341,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/pulse-mpi.py b/examples/pulse-mpi.py index 20064717a..b95f9475f 100644 --- a/examples/pulse-mpi.py +++ b/examples/pulse-mpi.py @@ -31,10 +31,11 @@ import pyopencl as cl import pyopencl.tools as cl_tools -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization + +from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator @@ -149,17 +150,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 1 - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -195,7 +196,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, # Set the current state from time 0 current_state = acoustic_pulse(x_vec=nodes, cv=uniform_state, eos=eos) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = "pulse" eosname = eos.__class__.__name__ @@ -215,7 +216,7 @@ def my_write_viz(step, t, state, dv=None): viz_fields = [("cv", state), ("dv", dv)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(step, t, state): @@ -236,8 +237,8 @@ def my_write_restart(step, t, state): def my_health_check(pressure): health_error = False from mirgecom.simutil import check_naninf_local, check_range_local - if check_naninf_local(discr, "vol", pressure) \ - or check_range_local(discr, "vol", pressure, .8, 1.5): + if check_naninf_local(dcoll, "vol", pressure) \ + or check_range_local(dcoll, "vol", pressure, .8, 1.5): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") return health_error @@ -279,7 +280,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -293,10 +294,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/scalar-lump-mpi.py b/examples/scalar-lump-mpi.py index 897ca80c2..2aa16cced 100644 --- a/examples/scalar-lump-mpi.py +++ b/examples/scalar-lump-mpi.py @@ -30,12 +30,12 @@ from functools import partial from pytools.obj_array import make_obj_array -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization -from grudge.shortcuts import make_visualizer +from grudge.discretization import DiscretizationCollection +from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator from mirgecom.simutil import ( @@ -144,17 +144,17 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, local_nelements = local_mesh.nelements order = 3 - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -193,7 +193,7 @@ def main(ctx_factory=cl.create_some_context, use_leap=False, # Set the current state from time 0 current_state = initializer(nodes) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ init_message = make_init_message(dim=dim, order=order, @@ -224,7 +224,7 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): ("exact", exact), ("resid", resid)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(step, t, state): @@ -245,8 +245,8 @@ def my_write_restart(step, t, state): 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) \ - or check_range_local(discr, "vol", pressure, .9, 1.1): + if check_naninf_local(dcoll, "vol", pressure) \ + or check_range_local(dcoll, "vol", pressure, .9, 1.1): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") @@ -277,7 +277,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) from mirgecom.simutil import allsync health_errors = allsync( my_health_check(dv.pressure, component_errors), @@ -305,7 +305,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) my_write_status(component_errors) except MyRuntimeError: @@ -315,7 +315,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -329,10 +329,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/sod-mpi.py b/examples/sod-mpi.py index fee600ad4..996197c1c 100644 --- a/examples/sod-mpi.py +++ b/examples/sod-mpi.py @@ -29,12 +29,12 @@ import pyopencl.tools as cl_tools from functools import partial -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization -from grudge.shortcuts import make_visualizer +from grudge.discretization import DiscretizationCollection +from grudge.shortcuts import make_visualizer from mirgecom.euler import euler_operator from mirgecom.simutil import ( @@ -144,17 +144,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, generate_mesh) local_nelements = local_mesh.nelements - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -180,7 +180,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, # Set the current state from time 0 current_state = initializer(nodes) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ @@ -213,7 +213,7 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): ("exact", exact), ("residual", resid)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(state, step, t): @@ -234,8 +234,8 @@ def my_write_restart(state, step, t): 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) \ - or check_range_local(discr, "vol", pressure, .09, 1.1): + if check_naninf_local(dcoll, "vol", pressure) \ + or check_range_local(dcoll, "vol", pressure, .09, 1.1): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") @@ -266,7 +266,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) from mirgecom.simutil import allsync health_errors = allsync( my_health_check(dv.pressure, component_errors), @@ -295,7 +295,7 @@ def my_pre_step(step, t, dt, 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) + compare_fluid_solutions(dcoll, state, exact) my_write_status(component_errors) except MyRuntimeError: @@ -305,7 +305,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -319,10 +319,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/vortex-mpi.py b/examples/vortex-mpi.py index 10c7fa234..71091dcd1 100644 --- a/examples/vortex-mpi.py +++ b/examples/vortex-mpi.py @@ -29,10 +29,11 @@ import pyopencl.tools as cl_tools from functools import partial -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import thaw, PyOpenCLArrayContext + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization + +from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer from mirgecom.profiling import PyOpenCLProfilingArrayContext @@ -148,17 +149,17 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, local_nelements = local_mesh.nelements order = 3 - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, local_mesh, order=order, mpi_communicator=comm ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) vis_timer = None if logmgr: logmgr_add_device_name(logmgr, queue) logmgr_add_device_memory_usage(logmgr, queue) - logmgr_add_many_discretization_quantities(logmgr, discr, dim, + logmgr_add_many_discretization_quantities(logmgr, dcoll, dim, extract_vars_for_logging, units_for_logging) vis_timer = IntervalTimer("t_vis", "Time spent visualizing") @@ -199,7 +200,7 @@ def main(ctx_factory=cl.create_some_context, use_logmgr=True, # Set the current state from time 0 current_state = initializer(nodes) - visualizer = make_visualizer(discr) + visualizer = make_visualizer(dcoll) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ @@ -217,7 +218,7 @@ 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) + get_inviscid_cfl(dcoll, eos, current_dt, state) if rank == 0: logger.info( f"------ {cfl=}\n" @@ -236,7 +237,7 @@ def my_write_viz(step, t, state, dv=None, exact=None, resid=None): ("exact", exact), ("residual", resid)] from mirgecom.simutil import write_visfile - write_visfile(discr, viz_fields, visualizer, vizname=casename, + write_visfile(dcoll, viz_fields, visualizer, vizname=casename, step=step, t=t, overwrite=True, vis_timer=vis_timer) def my_write_restart(step, t, state): @@ -257,8 +258,8 @@ def my_write_restart(step, t, state): 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) \ - or check_range_local(discr, "vol", pressure, .2, 1.02): + if check_naninf_local(dcoll, "vol", pressure) \ + or check_range_local(dcoll, "vol", pressure, .2, 1.02): health_error = True logger.info(f"{rank=}: Invalid pressure data found.") @@ -288,7 +289,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) from mirgecom.simutil import allsync health_errors = allsync( my_health_check(dv.pressure, component_errors), @@ -307,7 +308,7 @@ def my_pre_step(step, t, dt, state): 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) + component_errors = compare_fluid_solutions(dcoll, state, exact) my_write_status(state, component_errors) if do_viz: @@ -326,7 +327,7 @@ def my_pre_step(step, t, dt, state): my_write_restart(step=step, t=t, state=state) raise - dt = get_sim_timestep(discr, state, t, dt, current_cfl, eos, t_final, + dt = get_sim_timestep(dcoll, state, t, dt, current_cfl, eos, t_final, constant_cfl) return state, dt @@ -340,10 +341,10 @@ def my_post_step(step, t, dt, state): return state, dt def my_rhs(t, state): - return euler_operator(discr, cv=state, t=t, + return euler_operator(dcoll, cv=state, t=t, boundaries=boundaries, eos=eos) - current_dt = get_sim_timestep(discr, current_state, current_t, current_dt, + current_dt = get_sim_timestep(dcoll, current_state, current_t, current_dt, current_cfl, eos, t_final, constant_cfl) current_step, current_t, current_state = \ diff --git a/examples/wave-eager-mpi.py b/examples/wave-eager-mpi.py index 98991bab8..d69a6afd2 100644 --- a/examples/wave-eager-mpi.py +++ b/examples/wave-eager-mpi.py @@ -29,14 +29,16 @@ from pytools.obj_array import flat_obj_array -from meshmode.array_context import thaw, PyOpenCLArrayContext +from arraycontext import thaw, PyOpenCLArrayContext from mirgecom.profiling import PyOpenCLProfilingArrayContext from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import EagerDGDiscretization +from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer +import grudge.op as op + from mirgecom.mpi import mpi_entry_point from mirgecom.integrators import rk4_step from mirgecom.wave import wave_operator @@ -50,16 +52,16 @@ logmgr_add_device_memory_usage) -def bump(actx, discr, t=0): +def bump(actx, dcoll, t=0): """Create a bump.""" - source_center = np.array([0.2, 0.35, 0.1])[:discr.dim] + source_center = np.array([0.2, 0.35, 0.1])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) center_dist = flat_obj_array([ nodes[i] - source_center[i] - for i in range(discr.dim) + for i in range(dcoll.dim) ]) return ( @@ -130,7 +132,7 @@ def main(snapshot_pattern="wave-eager-{step:04d}-{rank:04d}.pkl", restart_step=N order = 3 - discr = EagerDGDiscretization(actx, local_mesh, order=order, + dcoll = DiscretizationCollection(actx, local_mesh, order=order, mpi_communicator=comm) if local_mesh.dim == 2: @@ -149,8 +151,8 @@ def main(snapshot_pattern="wave-eager-{step:04d}-{rank:04d}.pkl", restart_step=N istep = 0 fields = flat_obj_array( - bump(actx, discr), - [discr.zeros(actx) for i in range(discr.dim)] + bump(actx, dcoll), + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) else: @@ -160,10 +162,12 @@ def main(snapshot_pattern="wave-eager-{step:04d}-{rank:04d}.pkl", restart_step=N restart_fields = restart_data["fields"] old_order = restart_data["order"] if old_order != order: - old_discr = EagerDGDiscretization(actx, local_mesh, order=old_order, - mpi_communicator=comm) + old_discr = DiscretizationCollection( + actx, local_mesh, order=old_order, + mpi_communicator=comm + ) from meshmode.discretization.connection import make_same_mesh_connection - connection = make_same_mesh_connection(actx, discr.discr_from_dd("vol"), + connection = make_same_mesh_connection(actx, dcoll.discr_from_dd("vol"), old_discr.discr_from_dd("vol")) fields = connection(restart_fields) else: @@ -186,10 +190,10 @@ def main(snapshot_pattern="wave-eager-{step:04d}-{rank:04d}.pkl", restart_step=N vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) def rhs(t, w): - return wave_operator(discr, c=1, w=w) + return wave_operator(dcoll, c=1, w=w) while t < t_final: if logmgr: @@ -214,7 +218,7 @@ def rhs(t, w): ) if istep % 10 == 0: - print(istep, t, discr.norm(fields[0])) + print(istep, t, op.norm(dcoll, fields[0], np.inf)) vis.write_parallel_vtk_file( comm, "fld-wave-eager-mpi-%03d-%04d.vtu" % (rank, istep), diff --git a/examples/wave-eager.py b/examples/wave-eager.py index a06728d36..39bb87ac4 100644 --- a/examples/wave-eager.py +++ b/examples/wave-eager.py @@ -30,14 +30,14 @@ from pytools.obj_array import flat_obj_array -from grudge.eager import EagerDGDiscretization +from grudge.discretization import DiscretizationCollection from grudge.shortcuts import make_visualizer +import grudge.op as op from mirgecom.wave import wave_operator from mirgecom.integrators import rk4_step -from meshmode.dof_array import thaw -from meshmode.array_context import PyOpenCLArrayContext +from arraycontext import thaw, PyOpenCLArrayContext from mirgecom.profiling import PyOpenCLProfilingArrayContext @@ -48,16 +48,16 @@ logmgr_add_device_memory_usage) -def bump(actx, discr, t=0): +def bump(actx, dcoll, t=0): """Create a bump.""" - source_center = np.array([0.2, 0.35, 0.1])[:discr.dim] + source_center = np.array([0.2, 0.35, 0.1])[:dcoll.dim] source_width = 0.05 source_omega = 3 - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) center_dist = flat_obj_array([ nodes[i] - source_center[i] - for i in range(discr.dim) + for i in range(dcoll.dim) ]) return ( @@ -106,11 +106,11 @@ def main(use_profiling=False, use_logmgr=False): print("%d elements" % mesh.nelements) - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) fields = flat_obj_array( - bump(actx, discr), - [discr.zeros(actx) for i in range(discr.dim)] + bump(actx, dcoll), + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) if logmgr: @@ -130,10 +130,10 @@ def main(use_profiling=False, use_logmgr=False): vis_timer = IntervalTimer("t_vis", "Time spent visualizing") logmgr.add_quantity(vis_timer) - vis = make_visualizer(discr) + vis = make_visualizer(dcoll) def rhs(t, w): - return wave_operator(discr, c=1, w=w) + return wave_operator(dcoll, c=1, w=w) t = 0 t_final = 3 @@ -147,7 +147,7 @@ def rhs(t, w): if istep % 10 == 0: if use_profiling: print(actx.tabulate_profiling_data()) - print(istep, t, discr.norm(fields[0], np.inf)) + print(istep, t, op.norm(dcoll, fields[0], np.inf)) vis.write_vtk_file("fld-wave-eager-%04d.vtu" % istep, [ ("u", fields[0]), diff --git a/mirgecom/boundary.py b/mirgecom/boundary.py index 19d37ca87..ec5f4c6b7 100644 --- a/mirgecom/boundary.py +++ b/mirgecom/boundary.py @@ -33,11 +33,15 @@ """ import numpy as np -from meshmode.dof_array import thaw + +from arraycontext import thaw + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from grudge.trace_pair import TracePair from mirgecom.fluid import make_conserved +import grudge.op as op + class PrescribedBoundary: """Boundary condition prescribes boundary soln with user-specified function. @@ -59,14 +63,13 @@ def __init__(self, userfunc): """ self._userfunc = userfunc - def boundary_pair(self, discr, cv, btag, **kwargs): + def boundary_pair(self, dcoll, cv, btag, **kwargs): """Get the interior and exterior solution on the boundary.""" actx = cv.array_context - boundary_discr = discr.discr_from_dd(btag) - nodes = thaw(actx, boundary_discr.nodes()) + nodes = thaw(dcoll.nodes(btag), actx) ext_soln = self._userfunc(nodes, **kwargs) - int_soln = discr.project("vol", btag, cv) + int_soln = op.project(dcoll, "vol", btag, cv) return TracePair(btag, interior=int_soln, exterior=ext_soln) @@ -76,9 +79,9 @@ class DummyBoundary: .. automethod:: boundary_pair """ - def boundary_pair(self, discr, cv, btag, **kwargs): + def boundary_pair(self, dcoll, cv, btag, **kwargs): """Get the interior and exterior solution on the boundary.""" - dir_soln = discr.project("vol", btag, cv) + dir_soln = op.project(dcoll, "vol", btag, cv) return TracePair(btag, interior=dir_soln, exterior=dir_soln) @@ -100,7 +103,7 @@ class AdiabaticSlipBoundary: .. automethod:: boundary_pair """ - def boundary_pair(self, discr, cv, btag, **kwargs): + def boundary_pair(self, dcoll, cv, btag, **kwargs): """Get the interior and exterior solution on the boundary. The exterior solution is set such that there will be vanishing @@ -112,14 +115,14 @@ def boundary_pair(self, discr, cv, btag, **kwargs): E_plus = E_minus """ # Grab some boundary-relevant data - dim = discr.dim + dim = dcoll.dim actx = cv.mass.array_context # Grab a unit normal to the boundary - nhat = thaw(actx, discr.normal(btag)) + nhat = thaw(dcoll.normal(btag), actx) # Get the interior/exterior solns - int_cv = discr.project("vol", btag, cv) + int_cv = op.project(dcoll, "vol", btag, cv) # Subtract out the 2*wall-normal component # of velocity from the velocity at the wall to diff --git a/mirgecom/diffusion.py b/mirgecom/diffusion.py index 57c0bab66..2ce76d268 100644 --- a/mirgecom/diffusion.py +++ b/mirgecom/diffusion.py @@ -35,15 +35,19 @@ import abc import numpy as np import numpy.linalg as la # noqa + +from arraycontext import thaw + from pytools.obj_array import make_obj_array, obj_array_vectorize_n_args from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from meshmode.dof_array import thaw + from grudge.dof_desc import DOFDesc, as_dofdesc -from grudge.eager import interior_trace_pair, cross_rank_trace_pairs from grudge.trace_pair import TracePair +import grudge.op as op + -def gradient_flux(discr, quad_tag, u_tpair): +def gradient_flux(dcoll, quad_tag, u_tpair): r"""Compute the numerical flux for $\nabla u$.""" actx = u_tpair.int.array_context @@ -51,19 +55,19 @@ def gradient_flux(discr, quad_tag, u_tpair): dd_quad = dd.with_discr_tag(quad_tag) dd_allfaces_quad = dd_quad.with_dtag("all_faces") - normal_quad = thaw(actx, discr.normal(dd_quad)) + normal_quad = thaw(dcoll.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(dcoll, dd, dd_quad, a) def flux(u, normal): return -u * normal - return discr.project(dd_quad, dd_allfaces_quad, flux( + return op.project(dcoll, dd_quad, dd_allfaces_quad, flux( to_quad(u_tpair.avg), normal_quad)) -def diffusion_flux(discr, quad_tag, alpha_tpair, grad_u_tpair): +def diffusion_flux(dcoll, quad_tag, alpha_tpair, grad_u_tpair): r"""Compute the numerical flux for $\nabla \cdot (\alpha \nabla u)$.""" actx = grad_u_tpair.int[0].array_context @@ -71,10 +75,10 @@ def diffusion_flux(discr, quad_tag, alpha_tpair, grad_u_tpair): dd_quad = dd.with_discr_tag(quad_tag) dd_allfaces_quad = dd_quad.with_dtag("all_faces") - normal_quad = thaw(actx, discr.normal(dd_quad)) + normal_quad = thaw(dcoll.normal(dd_quad), actx) def to_quad(a): - return discr.project(dd, dd_quad, a) + return op.project(dcoll, dd, dd_quad, a) def flux(alpha, grad_u, normal): return -alpha * np.dot(grad_u, normal) @@ -86,7 +90,7 @@ def flux(alpha, grad_u, normal): to_quad(alpha_tpair.ext), to_quad(grad_u_tpair.ext), normal_quad) ) - return discr.project(dd_quad, dd_allfaces_quad, flux_tpair.avg) + return op.project(dcoll, dd_quad, dd_allfaces_quad, flux_tpair.avg) class DiffusionBoundary(metaclass=abc.ABCMeta): @@ -98,12 +102,12 @@ class DiffusionBoundary(metaclass=abc.ABCMeta): """ @abc.abstractmethod - def get_gradient_flux(self, discr, quad_tag, dd, alpha, u): + def get_gradient_flux(self, dcoll, quad_tag, dd, alpha, u): """Compute the flux for grad(u) on the boundary corresponding to *dd*.""" raise NotImplementedError @abc.abstractmethod - def get_diffusion_flux(self, discr, quad_tag, dd, alpha, grad_u): + def get_diffusion_flux(self, dcoll, quad_tag, dd, alpha, grad_u): """Compute the flux for diff(u) on the boundary corresponding to *dd*.""" raise NotImplementedError @@ -136,17 +140,17 @@ def __init__(self, value): """ self.value = value - def get_gradient_flux(self, discr, quad_tag, dd, alpha, u): # noqa: D102 - u_int = discr.project("vol", dd, u) + def get_gradient_flux(self, dcoll, quad_tag, dd, alpha, u): # noqa: D102 + u_int = op.project(dcoll, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=2*self.value-u_int) - return gradient_flux(discr, quad_tag, u_tpair) + return gradient_flux(dcoll, quad_tag, u_tpair) - def get_diffusion_flux(self, discr, quad_tag, dd, alpha, grad_u): # noqa: D102 - alpha_int = discr.project("vol", dd, alpha) + def get_diffusion_flux(self, dcoll, quad_tag, dd, alpha, grad_u): # noqa: D102 + alpha_int = op.project(dcoll, "vol", dd, alpha) alpha_tpair = TracePair(dd, interior=alpha_int, exterior=alpha_int) - grad_u_int = discr.project("vol", dd, grad_u) + grad_u_int = op.project(dcoll, "vol", dd, grad_u) grad_u_tpair = TracePair(dd, interior=grad_u_int, exterior=grad_u_int) - return diffusion_flux(discr, quad_tag, alpha_tpair, grad_u_tpair) + return diffusion_flux(dcoll, quad_tag, alpha_tpair, grad_u_tpair) class NeumannDiffusionBoundary(DiffusionBoundary): @@ -185,12 +189,12 @@ def __init__(self, value): """ self.value = value - def get_gradient_flux(self, discr, quad_tag, dd, alpha, u): # noqa: D102 - u_int = discr.project("vol", dd, u) + def get_gradient_flux(self, dcoll, quad_tag, dd, alpha, u): # noqa: D102 + u_int = op.project(dcoll, "vol", dd, u) u_tpair = TracePair(dd, interior=u_int, exterior=u_int) - return gradient_flux(discr, quad_tag, u_tpair) + return gradient_flux(dcoll, quad_tag, u_tpair) - def get_diffusion_flux(self, discr, quad_tag, dd, alpha, grad_u): # noqa: D102 + def get_diffusion_flux(self, dcoll, quad_tag, dd, alpha, grad_u): # noqa: D102 dd_quad = dd.with_discr_tag(quad_tag) dd_allfaces_quad = dd_quad.with_dtag("all_faces") # Compute the flux directly instead of constructing an external grad_u value @@ -198,13 +202,13 @@ def get_diffusion_flux(self, discr, quad_tag, dd, alpha, grad_u): # noqa: D102 # spatially-varying alpha case (the other approach would result in a # grad_u_tpair that lives in the quadrature discretization; diffusion_flux # would need to be modified to accept such values). - alpha_int_quad = discr.project("vol", dd_quad, alpha) - value_quad = discr.project(dd, dd_quad, self.value) + alpha_int_quad = op.project(dcoll, "vol", dd_quad, alpha) + value_quad = op.project(dcoll, dd, dd_quad, self.value) flux_quad = -alpha_int_quad*value_quad - return discr.project(dd_quad, dd_allfaces_quad, flux_quad) + return op.project(dcoll, dd_quad, dd_allfaces_quad, flux_quad) -def diffusion_operator(discr, quad_tag, alpha, boundaries, u, return_grad_u=False): +def diffusion_operator(dcoll, quad_tag, alpha, boundaries, u, return_grad_u=False): r""" Compute the diffusion operator. @@ -216,10 +220,11 @@ def diffusion_operator(discr, quad_tag, alpha, boundaries, u, return_grad_u=Fals Parameters ---------- - discr: grudge.eager.EagerDGDiscretization - the discretization to use + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. quad_tag: - quadrature tag indicating which discretization in *discr* to use for + quadrature tag indicating which discretization in *dcoll* to use for overintegration alpha: numbers.Number or meshmode.dof_array.DOFArray the diffusivity value(s) @@ -245,7 +250,7 @@ def diffusion_operator(discr, quad_tag, alpha, boundaries, u, return_grad_u=Fals if len(boundaries) != len(u): raise TypeError("boundaries must be the same length as u") return obj_array_vectorize_n_args(lambda boundaries_i, u_i: - diffusion_operator(discr, quad_tag, alpha, boundaries_i, u_i, + diffusion_operator(dcoll, quad_tag, alpha, boundaries_i, u_i, return_grad_u=return_grad_u), make_obj_array(boundaries), u) for btag, bdry in boundaries.items(): @@ -256,39 +261,43 @@ def diffusion_operator(discr, quad_tag, alpha, boundaries, u, return_grad_u=Fals dd_quad = DOFDesc("vol", quad_tag) dd_allfaces_quad = DOFDesc("all_faces", quad_tag) - grad_u = discr.inverse_mass( - discr.weak_grad(-u) + grad_u = op.inverse_mass( + dcoll, + op.weak_local_grad(dcoll, -u) - # noqa: W504 - discr.face_mass( + op.face_mass( + dcoll, dd_allfaces_quad, - gradient_flux(discr, quad_tag, interior_trace_pair(discr, u)) + gradient_flux(dcoll, quad_tag, op.interior_trace_pair(dcoll, u)) + sum( - bdry.get_gradient_flux(discr, quad_tag, as_dofdesc(btag), alpha, u) + bdry.get_gradient_flux(dcoll, quad_tag, as_dofdesc(btag), alpha, u) for btag, bdry in boundaries.items()) + sum( - gradient_flux(discr, quad_tag, u_tpair) - for u_tpair in cross_rank_trace_pairs(discr, u)) + gradient_flux(dcoll, quad_tag, u_tpair) + for u_tpair in op.cross_rank_trace_pairs(dcoll, u)) ) ) - alpha_quad = discr.project("vol", dd_quad, alpha) - grad_u_quad = discr.project("vol", dd_quad, grad_u) + alpha_quad = op.project(dcoll, "vol", dd_quad, alpha) + grad_u_quad = op.project(dcoll, "vol", dd_quad, grad_u) - diff_u = discr.inverse_mass( - discr.weak_div(dd_quad, -alpha_quad*grad_u_quad) + diff_u = op.inverse_mass( + dcoll, + op.weak_local_div(dcoll, dd_quad, -alpha_quad*grad_u_quad) - # noqa: W504 - discr.face_mass( + op.face_mass( + dcoll, dd_allfaces_quad, - diffusion_flux(discr, quad_tag, interior_trace_pair(discr, alpha), - interior_trace_pair(discr, grad_u)) + diffusion_flux(dcoll, quad_tag, op.interior_trace_pair(dcoll, alpha), + op.interior_trace_pair(dcoll, grad_u)) + sum( - bdry.get_diffusion_flux(discr, quad_tag, as_dofdesc(btag), alpha, + bdry.get_diffusion_flux(dcoll, quad_tag, as_dofdesc(btag), alpha, grad_u) for btag, bdry in boundaries.items()) + sum( - diffusion_flux(discr, quad_tag, alpha_tpair, grad_u_tpair) + diffusion_flux(dcoll, quad_tag, alpha_tpair, grad_u_tpair) for alpha_tpair, grad_u_tpair in zip( - cross_rank_trace_pairs(discr, alpha), - cross_rank_trace_pairs(discr, grad_u))) + op.cross_rank_trace_pairs(dcoll, alpha), + op.cross_rank_trace_pairs(dcoll, grad_u))) ) ) diff --git a/mirgecom/euler.py b/mirgecom/euler.py index 7f2caab08..6b904578d 100644 --- a/mirgecom/euler.py +++ b/mirgecom/euler.py @@ -53,29 +53,31 @@ """ import numpy as np # noqa -from meshmode.dof_array import thaw -from grudge.symbolic.primitives import TracePair -from grudge.eager import ( - interior_trace_pair, - cross_rank_trace_pairs -) + +from arraycontext import thaw + +from grudge.trace_pair import TracePair +import grudge.op as op + from mirgecom.fluid import ( compute_wavespeed, split_conserved, ) +from mirgecom.flux import lfr_flux +from mirgecom.inviscid import inviscid_flux -from mirgecom.inviscid import ( - inviscid_flux -) from functools import partial -from mirgecom.flux import lfr_flux -def _facial_flux(discr, eos, cv_tpair, local=False): +def _facial_flux(dcoll, eos, cv_tpair, local=False): """Return the flux across a face given the solution on both sides *cv_tpair*. Parameters ---------- + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. + eos: mirgecom.eos.GasEOS Implementing the pressure and temperature functions for returning pressure and temperature as a function of the state q. @@ -92,23 +94,23 @@ def _facial_flux(discr, eos, cv_tpair, local=False): actx = cv_tpair.int.array_context dim = cv_tpair.int.dim - euler_flux = partial(inviscid_flux, discr, eos) + euler_flux = partial(inviscid_flux, dcoll, eos) lam = actx.np.maximum( compute_wavespeed(dim, eos, cv_tpair.int), compute_wavespeed(dim, eos, cv_tpair.ext) ) - normal = thaw(actx, discr.normal(cv_tpair.dd)) + normal = thaw(dcoll.normal(cv_tpair.dd), actx) # todo: user-supplied flux routine flux_weak = lfr_flux(cv_tpair=cv_tpair, flux_func=euler_flux, normal=normal, lam=lam) if local is False: - return discr.project(cv_tpair.dd, "all_faces", flux_weak) + return op.project(dcoll, cv_tpair.dd, "all_faces", flux_weak) return flux_weak -def euler_operator(discr, eos, boundaries, cv, t=0.0): +def euler_operator(dcoll, eos, boundaries, cv, t=0.0): r"""Compute RHS of the Euler flow equations. Returns @@ -123,57 +125,63 @@ def euler_operator(discr, eos, boundaries, cv, t=0.0): Parameters ---------- - cv: :class:`mirgecom.fluid.ConservedVars` - Fluid conserved state object with the conserved variables. + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. + + eos: mirgecom.eos.GasEOS + Implementing the pressure and temperature functions for + returning pressure and temperature as a function of the state. boundaries Dictionary of boundary functions, one for each valid btag + cv: :class:`mirgecom.fluid.ConservedVars` + Fluid conserved state object with the conserved variables. + t Time - eos: mirgecom.eos.GasEOS - Implementing the pressure and temperature functions for - returning pressure and temperature as a function of the state q. - Returns ------- numpy.ndarray Agglomerated object array of DOF arrays representing the RHS of the Euler flow equations. """ - vol_weak = discr.weak_div(inviscid_flux(discr=discr, eos=eos, cv=cv).join()) + vol_weak = op.weak_local_div( + dcoll, inviscid_flux(dcoll=dcoll, eos=eos, cv=cv).join() + ) boundary_flux = ( - _facial_flux(discr=discr, eos=eos, cv_tpair=interior_trace_pair(discr, cv)) - + sum( + sum( _facial_flux( - discr, eos=eos, + dcoll, eos=eos, cv_tpair=TracePair( - part_pair.dd, - interior=split_conserved(discr.dim, part_pair.int), - exterior=split_conserved(discr.dim, part_pair.ext))) - for part_pair in cross_rank_trace_pairs(discr, cv.join())) + qt_pair.dd, + interior=split_conserved(dcoll.dim, qt_pair.int), + exterior=split_conserved(dcoll.dim, qt_pair.ext))) + for qt_pair in op.interior_trace_pairs(dcoll, cv.join())) + sum( _facial_flux( - discr=discr, eos=eos, + dcoll=dcoll, eos=eos, cv_tpair=boundaries[btag].boundary_pair( - discr, eos=eos, btag=btag, t=t, cv=cv) + dcoll, eos=eos, btag=btag, t=t, cv=cv) ) for btag in boundaries) ).join() return split_conserved( - discr.dim, discr.inverse_mass(vol_weak - discr.face_mass(boundary_flux)) + dcoll.dim, + op.inverse_mass(dcoll, vol_weak - op.face_mass(dcoll, boundary_flux)) ) -def inviscid_operator(discr, eos, boundaries, q, t=0.0): +def inviscid_operator(dcoll, eos, boundaries, q, t=0.0): """Interface :function:`euler_operator` with backwards-compatible API.""" from warnings import warn warn("Do not call inviscid_operator; it is now called euler_operator. This" "function will disappear August 1, 2021", DeprecationWarning, stacklevel=2) - return euler_operator(discr, eos, boundaries, split_conserved(discr.dim, q), t) + return euler_operator(dcoll, eos, boundaries, split_conserved(dcoll.dim, q), t) # By default, run unitless diff --git a/mirgecom/filter.py b/mirgecom/filter.py index 2737970b9..291f77902 100644 --- a/mirgecom/filter.py +++ b/mirgecom/filter.py @@ -94,7 +94,7 @@ def make_spectral_filter(actx, group, cutoff, mode_response_function): lambda grp: grp.discretization_key() ) def _spectral_filter_scaling(group): - mode_ids = group.mode_ids() + mode_ids = group.basis_obj().mode_ids order = group.order nmodes = len(mode_ids) diff --git a/mirgecom/fluid.py b/mirgecom/fluid.py index 21b2f5d21..0b27fad7b 100644 --- a/mirgecom/fluid.py +++ b/mirgecom/fluid.py @@ -328,7 +328,7 @@ def make_conserved(dim, mass, energy, momentum, species_mass=None): ) -def velocity_gradient(discr, cv, grad_cv): +def velocity_gradient(dcoll, cv, grad_cv): r""" Compute the gradient of fluid velocity. @@ -355,8 +355,9 @@ def velocity_gradient(discr, cv, grad_cv): Parameters ---------- - discr: grudge.eager.EagerDGDiscretization - the discretization to use + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. cv: ConservedVars the fluid conserved variables grad_cv: ConservedVars @@ -378,7 +379,7 @@ def velocity_gradient(discr, cv, grad_cv): for i in range(cv.dim)]) -def species_mass_fraction_gradient(discr, cv, grad_cv): +def species_mass_fraction_gradient(dcoll, cv, grad_cv): r""" Compute the gradient of species mass fractions. @@ -393,8 +394,9 @@ def species_mass_fraction_gradient(discr, cv, grad_cv): Parameters ---------- - discr: grudge.eager.EagerDGDiscretization - the discretization to use + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. cv: ConservedVars the fluid conserved variables grad_cv: ConservedVars diff --git a/mirgecom/initializers.py b/mirgecom/initializers.py index 1174a2de4..a60d3e89e 100644 --- a/mirgecom/initializers.py +++ b/mirgecom/initializers.py @@ -38,8 +38,9 @@ """ import numpy as np + +from arraycontext import thaw from pytools.obj_array import make_obj_array -from meshmode.dof_array import thaw from mirgecom.eos import IdealSingleGas from mirgecom.fluid import make_conserved @@ -373,7 +374,7 @@ def __call__(self, x_vec, *, t=0, eos=None): return make_conserved(dim=self._dim, mass=mass, energy=energy, momentum=mom) - def exact_rhs(self, discr, cv, t=0.0): + def exact_rhs(self, dcoll, cv, t=0.0): """ Create the RHS for the lump-of-mass solution at time *t*, locations *x_vec*. @@ -382,14 +383,17 @@ def exact_rhs(self, discr, cv, t=0.0): Parameters ---------- - q + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. + cv State array which expects at least the canonical conserved quantities (mass, energy, momentum) for the fluid at each point. t: float Time at which RHS is desired """ actx = cv.array_context - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) lump_loc = self._center + t * self._velocity # coordinates relative to lump center rel_center = make_obj_array( @@ -550,7 +554,7 @@ def __call__(self, x_vec, *, t=0, eos=None): return make_conserved(dim=self._dim, mass=mass, energy=energy, momentum=mom, species_mass=species_mass) - def exact_rhs(self, discr, cv, t=0.0): + def exact_rhs(self, dcoll, cv, t=0.0): """ Create a RHS for multi-component lump soln at time *t*, locations *x_vec*. @@ -559,14 +563,17 @@ def exact_rhs(self, discr, cv, t=0.0): Parameters ---------- - q + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. + cv State array which expects at least the canonical conserved quantities (mass, energy, momentum) for the fluid at each point. t: float Time at which RHS is desired """ actx = cv.array_context - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) loc_update = t * self._velocity mass = 0 * nodes[0] + self._rho0 @@ -747,20 +754,23 @@ def __call__(self, x_vec, *, t=0, eos=None): return make_conserved(dim=self._dim, mass=mass, energy=energy, momentum=mom, species_mass=species_mass) - def exact_rhs(self, discr, cv, t=0.0): + def exact_rhs(self, dcoll, cv, t=0.0): """ Create the RHS for the uniform solution. (Hint - it should be all zero). Parameters ---------- - q + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. + cv State array which expects at least the canonical conserved quantities (mass, energy, momentum) for the fluid at each point. (unused) t: float Time at which RHS is desired (unused) """ actx = cv.array_context - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) mass = nodes[0].copy() mass[:] = 1.0 massrhs = 0.0 * mass diff --git a/mirgecom/inviscid.py b/mirgecom/inviscid.py index e3fe348d2..1d4082e0f 100644 --- a/mirgecom/inviscid.py +++ b/mirgecom/inviscid.py @@ -40,7 +40,7 @@ from mirgecom.fluid import make_conserved -def inviscid_flux(discr, eos, cv): +def inviscid_flux(dcoll, eos, cv): r"""Compute the inviscid flux vectors from fluid conserved vars *cv*. The inviscid fluxes are @@ -66,15 +66,16 @@ def inviscid_flux(discr, eos, cv): (mom / cv.mass) * cv.species_mass.reshape(-1, 1))) -def get_inviscid_timestep(discr, eos, cv): +def get_inviscid_timestep(dcoll, 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 + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. eos: mirgecom.eos.GasEOS Implementing the pressure and temperature functions for returning pressure and temperature as a function of the state q. @@ -88,18 +89,19 @@ def get_inviscid_timestep(discr, eos, cv): 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) + characteristic_lengthscales(cv.array_context, dcoll) + / compute_wavespeed(dcoll, eos, cv) ) -def get_inviscid_cfl(discr, eos, dt, cv): +def get_inviscid_cfl(dcoll, eos, dt, cv): """Return node-local CFL based on current state and timestep. Parameters ---------- - discr: :class:`grudge.eager.EagerDGDiscretization` - the discretization to use + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. eos: mirgecom.eos.GasEOS Implementing the pressure and temperature functions for returning pressure and temperature as a function of the state q. @@ -113,4 +115,4 @@ def get_inviscid_cfl(discr, eos, dt, cv): :class:`meshmode.dof_array.DOFArray` The CFL at each node. """ - return dt / get_inviscid_timestep(discr, eos=eos, cv=cv) + return dt / get_inviscid_timestep(dcoll, eos=eos, cv=cv) diff --git a/mirgecom/io.py b/mirgecom/io.py index e6c7eb817..97862df81 100644 --- a/mirgecom/io.py +++ b/mirgecom/io.py @@ -31,6 +31,8 @@ from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa +import grudge.op as op + def make_init_message(*, dim, order, dt, t_final, nstatus, nviz, cfl, constant_cfl, @@ -51,12 +53,16 @@ def make_init_message(*, dim, order, dt, t_final, ) -def make_status_message(*, discr, t, step, dt, cfl, dependent_vars): +def make_status_message(*, dcoll, t, step, dt, cfl, dependent_vars): r"""Make simulation status and health message.""" dv = dependent_vars - from functools import partial - _min = partial(discr.nodal_min, "vol") - _max = partial(discr.nodal_max, "vol") + + def _min(u): + return op.nodal_min(dcoll, "vol", u) + + def _max(u): + return op.nodal_max(dcoll, "vol", u) + statusmsg = ( f"Status: {step=} {t=}\n" f"------- P({_min(dv.pressure):.3g}, {_max(dv.pressure):.3g})\n" diff --git a/mirgecom/logging_quantities.py b/mirgecom/logging_quantities.py index 5e94e1965..b7e8113e8 100644 --- a/mirgecom/logging_quantities.py +++ b/mirgecom/logging_quantities.py @@ -43,7 +43,10 @@ MultiPostLogQuantity, add_run_info, add_general_quantities, add_simulation_quantities) from meshmode.array_context import PyOpenCLArrayContext -from meshmode.discretization import Discretization + +from grudge.discretization import DiscretizationCollection +import grudge.op as grudge_op + import pyopencl as cl from typing import Optional, Callable @@ -98,21 +101,21 @@ def logmgr_add_device_memory_usage(logmgr: LogManager, queue: cl.CommandQueue): logmgr.add_quantity(DeviceMemoryUsage()) -def logmgr_add_many_discretization_quantities(logmgr: LogManager, discr, dim, +def logmgr_add_many_discretization_quantities(logmgr: LogManager, dcoll, dim, extract_vars_for_logging, units_for_logging): """Add default discretization quantities to the logmgr.""" for op in ["min", "max", "L2_norm"]: for quantity in ["pressure", "temperature"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + dcoll, quantity, op, extract_vars_for_logging, units_for_logging)) for quantity in ["mass", "energy"]: logmgr.add_quantity(DiscretizationBasedQuantity( - discr, quantity, op, extract_vars_for_logging, units_for_logging)) + dcoll, quantity, op, extract_vars_for_logging, units_for_logging)) for d in range(dim): logmgr.add_quantity(DiscretizationBasedQuantity( - discr, "momentum", op, extract_vars_for_logging, units_for_logging, + dcoll, "momentum", op, extract_vars_for_logging, units_for_logging, axis=d)) @@ -237,7 +240,7 @@ class DiscretizationBasedQuantity(PostLogQuantity, StateConsumer): Possible rank aggregation operations (``op``) are: min, max, L2_norm. """ - def __init__(self, discr: Discretization, quantity: str, op: str, + def __init__(self, dcoll: DiscretizationCollection, quantity: str, op: str, extract_vars_for_logging, units_logging, name: str = None, axis: Optional[int] = None): unit = units_logging(quantity) @@ -248,21 +251,21 @@ def __init__(self, discr: Discretization, quantity: str, op: str, LogQuantity.__init__(self, name, unit) StateConsumer.__init__(self, extract_vars_for_logging) - self.discr = discr + self.dcoll = dcoll self.quantity = quantity self.axis = axis - from functools import partial - if op == "min": - self._discr_reduction = partial(self.discr.nodal_min, "vol") + self._nodal_reduction = lambda u: grudge_op.nodal_min(self.dcoll, + "vol", u) self.rank_aggr = min elif op == "max": - self._discr_reduction = partial(self.discr.nodal_max, "vol") + self._nodal_reduction = lambda u: grudge_op.nodal_max(self.dcoll, + "vol", u) self.rank_aggr = max elif op == "L2_norm": - self._discr_reduction = partial(self.discr.norm, p=2) + self._nodal_reduction = lambda u: grudge_op.norm(self.dcoll, u, 2) self.rank_aggr = max else: raise ValueError(f"unknown operation {op}") @@ -282,7 +285,7 @@ def __call__(self): if self.axis is not None: # e.g. momentum quantity = quantity[self.axis] - return self._discr_reduction(quantity) + return self._nodal_reduction(quantity) # }}} diff --git a/mirgecom/simutil.py b/mirgecom/simutil.py index 3b0a03c1e..d31a53742 100644 --- a/mirgecom/simutil.py +++ b/mirgecom/simutil.py @@ -74,7 +74,7 @@ def check_step(step, interval): return False -def get_sim_timestep(discr, state, t, dt, cfl, eos, +def get_sim_timestep(dcoll, state, t, dt, cfl, eos, t_final, constant_cfl=False): """Return the maximum stable timestep for a typical fluid simulation. @@ -95,8 +95,9 @@ def get_sim_timestep(discr, state, t, dt, cfl, eos, Parameters ---------- - discr - Grudge discretization or discretization collection? + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. state: :class:`~mirgecom.fluid.ConservedVars` The fluid state. t: float @@ -123,13 +124,13 @@ def get_sim_timestep(discr, state, t, dt, cfl, eos, 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) + dcoll, "vol", + get_inviscid_timestep(dcoll=dcoll, eos=eos, cv=state) ) return min(t_remaining, mydt) -def write_visfile(discr, io_fields, visualizer, vizname, +def write_visfile(dcoll, io_fields, visualizer, vizname, step=0, t=0, overwrite=False, vis_timer=None): """Write VTK output for the fields specified in *io_fields*. @@ -144,7 +145,7 @@ def write_visfile(discr, io_fields, visualizer, vizname, from contextlib import nullcontext from mirgecom.io import make_rank_fname, make_par_fname - comm = discr.mpi_communicator + comm = dcoll.mpi_communicator rank = 0 if comm is not None: rank = comm.Get_rank() @@ -176,24 +177,24 @@ def allsync(local_values, comm=None, op=None): return comm.allreduce(local_values, op=op) -def check_range_local(discr, dd, field, min_value, max_value): +def check_range_local(dcoll, dd, field, min_value, max_value): """Check for any negative values.""" return ( - op.nodal_min_loc(discr, dd, field) < min_value - or op.nodal_max_loc(discr, dd, field) > max_value + op.nodal_min_loc(dcoll, dd, field) < min_value + or op.nodal_max_loc(dcoll, dd, field) > max_value ) -def check_naninf_local(discr, dd, field): +def check_naninf_local(dcoll, dd, field): """Check for any NANs or Infs in the field.""" - s = op.nodal_sum_loc(discr, dd, field) + s = op.nodal_sum_loc(dcoll, dd, field) return np.isnan(s) or (s == np.inf) -def compare_fluid_solutions(discr, red_state, blue_state): +def compare_fluid_solutions(dcoll, red_state, blue_state): """Return inf norm of (*red_state* - *blue_state*) for each component.""" resid = red_state - blue_state - return [discr.norm(v, np.inf) for v in resid.join()] + return [op.norm(dcoll, v, np.inf) for v in resid.join()] def generate_and_distribute_mesh(comm, generate_mesh): diff --git a/mirgecom/wave.py b/mirgecom/wave.py index 9254b6a50..8781ac0ec 100644 --- a/mirgecom/wave.py +++ b/mirgecom/wave.py @@ -29,20 +29,23 @@ import numpy as np import numpy.linalg as la # noqa + +from arraycontext import thaw + from pytools.obj_array import flat_obj_array from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from meshmode.dof_array import thaw + from grudge.trace_pair import TracePair -from grudge.eager import interior_trace_pair, cross_rank_trace_pairs +import grudge.op as op -def _flux(discr, c, w_tpair): +def _flux(dcoll, c, w_tpair): u = w_tpair[0] v = w_tpair[1:] actx = w_tpair.int[0].array_context - normal = thaw(actx, discr.normal(w_tpair.dd)) + normal = thaw(dcoll.normal(w_tpair.dd), actx) flux_weak = flat_obj_array( np.dot(v.avg, normal), @@ -55,16 +58,17 @@ def _flux(discr, c, w_tpair): 0.5*normal*np.dot(normal, v.ext-v.int), ) - return discr.project(w_tpair.dd, "all_faces", c*flux_weak) + return op.project(dcoll, w_tpair.dd, "all_faces", c*flux_weak) -def wave_operator(discr, c, w): +def wave_operator(dcoll, c, w): """Compute the RHS of the wave equation. Parameters ---------- - discr: grudge.eager.EagerDGDiscretization - the discretization to use + dcoll: :class:`grudge.discretization.DiscretizationCollection` + An object containing connections and mappings to different + discretizations over the mesh. c: float the (constant) wave speed w: numpy.ndarray @@ -78,25 +82,27 @@ def wave_operator(discr, c, w): u = w[0] v = w[1:] - dir_u = discr.project("vol", BTAG_ALL, u) - dir_v = discr.project("vol", BTAG_ALL, v) + dir_u = op.project(dcoll, "vol", BTAG_ALL, u) + dir_v = op.project(dcoll, "vol", BTAG_ALL, v) dir_bval = flat_obj_array(dir_u, dir_v) dir_bc = flat_obj_array(-dir_u, dir_v) return ( - discr.inverse_mass( + op.inverse_mass( + dcoll, flat_obj_array( - -c*discr.weak_div(v), - -c*discr.weak_grad(u) + -c*op.weak_local_div(dcoll, v), + -c*op.weak_local_grad(dcoll, u) ) + # noqa: W504 - discr.face_mass( - _flux(discr, c=c, w_tpair=interior_trace_pair(discr, w)) - + _flux(discr, c=c, + op.face_mass( + dcoll, + _flux(dcoll, c=c, w_tpair=op.interior_trace_pair(dcoll, w)) + + _flux(dcoll, c=c, w_tpair=TracePair(BTAG_ALL, interior=dir_bval, exterior=dir_bc)) + sum( - _flux(discr, c=c, w_tpair=tpair) - for tpair in cross_rank_trace_pairs(discr, w)) + _flux(dcoll, c=c, w_tpair=tpair) + for tpair in op.cross_rank_trace_pairs(dcoll, w)) ) ) ) diff --git a/test/test_bc.py b/test/test_bc.py index 47806879c..130d91a6a 100644 --- a/test/test_bc.py +++ b/test/test_bc.py @@ -29,15 +29,19 @@ import logging import pytest -from meshmode.dof_array import thaw +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from mirgecom.initializers import Lump from mirgecom.boundary import AdiabaticSlipBoundary from mirgecom.eos import IdealSingleGas -from grudge.eager import EagerDGDiscretization -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) + +from grudge.discretization import DiscretizationCollection +import grudge.op as op logger = logging.getLogger(__name__) @@ -63,11 +67,14 @@ def test_slipwall_identity(actx_factory, dim): ) order = 3 - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) eos = IdealSingleGas() orig = np.zeros(shape=(dim,)) - nhat = thaw(actx, discr.normal(BTAG_ALL)) + nhat = thaw(dcoll.normal(BTAG_ALL), actx) + + def bnd_norm(u): + return op.norm(dcoll, u, p=np.inf, dd=BTAG_ALL) logger.info(f"Number of {dim}d elems: {mesh.nelements}") @@ -81,10 +88,8 @@ def test_slipwall_identity(actx_factory, dim): wall = AdiabaticSlipBoundary() uniform_state = initializer(nodes) - from functools import partial - bnd_norm = partial(discr.norm, p=np.inf, dd=BTAG_ALL) - bnd_pair = wall.boundary_pair(discr, btag=BTAG_ALL, + bnd_pair = wall.boundary_pair(dcoll, btag=BTAG_ALL, eos=eos, cv=uniform_state) # check that mass and energy are preserved @@ -128,13 +133,13 @@ def test_slipwall_flux(actx_factory, dim, order): a=(-0.5,) * dim, b=(0.5,) * dim, nelements_per_axis=(nel_1d,) * dim ) - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) - nhat = thaw(actx, discr.normal(BTAG_ALL)) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) + nhat = thaw(dcoll.normal(BTAG_ALL), actx) h = 1.0 / nel_1d - from functools import partial - bnd_norm = partial(discr.norm, p=np.inf, dd=BTAG_ALL) + def bnd_norm(u): + return op.norm(dcoll, u, p=np.inf, dd=BTAG_ALL) logger.info(f"Number of {dim}d elems: {mesh.nelements}") # for velocities in each direction @@ -148,7 +153,7 @@ def test_slipwall_flux(actx_factory, dim, order): from mirgecom.initializers import Uniform initializer = Uniform(dim=dim, velocity=vel) uniform_state = initializer(nodes) - bnd_pair = wall.boundary_pair(discr, btag=BTAG_ALL, + bnd_pair = wall.boundary_pair(dcoll, btag=BTAG_ALL, eos=eos, cv=uniform_state) # Check the total velocity component normal @@ -158,7 +163,7 @@ def test_slipwall_flux(actx_factory, dim, order): err_max = max(err_max, bnd_norm(np.dot(avg_state.momentum, nhat))) from mirgecom.euler import _facial_flux - bnd_flux = _facial_flux(discr, eos, cv_tpair=bnd_pair, local=True) + bnd_flux = _facial_flux(dcoll, eos, cv_tpair=bnd_pair, local=True) err_max = max(err_max, bnd_norm(bnd_flux.mass), bnd_norm(bnd_flux.energy)) diff --git a/test/test_diffusion.py b/test/test_diffusion.py index 2103badc3..fecb8d981 100644 --- a/test/test_diffusion.py +++ b/test/test_diffusion.py @@ -27,16 +27,22 @@ import pymbolic as pmbl import pymbolic.primitives as prim import mirgecom.symbolic as sym + +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from mirgecom.diffusion import ( diffusion_operator, DirichletDiffusionBoundary, NeumannDiffusionBoundary) -from meshmode.dof_array import thaw, DOFArray -from grudge.dof_desc import DTAG_BOUNDARY, DISCR_TAG_BASE, DISCR_TAG_QUAD -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) +from meshmode.dof_array import DOFArray + +from grudge.dof_desc import DTAG_BOUNDARY, DISCR_TAG_BASE, DISCR_TAG_QUAD +import grudge.op as op import pytest @@ -111,7 +117,7 @@ def get_mesh(n): sym_u *= sym_sin(sym_coords[i]) sym_u *= sym_cos(sym_coords[dim-1]) - def get_boundaries(discr, actx, t): + def get_boundaries(dcoll, actx, t): boundaries = {} for i in range(dim-1): @@ -145,8 +151,8 @@ def get_mesh(n): sym_u *= sym_sin(sym_coords[i]) sym_u *= sym_cos(sym_coords[dim-1]) - def get_boundaries(discr, actx, t): - nodes = thaw(actx, discr.nodes()) + def get_boundaries(dcoll, actx, t): + nodes = thaw(dcoll.nodes(), actx) def sym_eval(expr): return sym.EvaluationMapper({"x": nodes, "t": t})(expr) @@ -159,15 +165,15 @@ def sym_eval(expr): for i in range(dim-1): lower_btag = DTAG_BOUNDARY("-"+str(i)) upper_btag = DTAG_BOUNDARY("+"+str(i)) - upper_grad_u = discr.project("vol", upper_btag, exact_grad_u) - normal = thaw(actx, discr.normal(upper_btag)) + upper_grad_u = op.project(dcoll, "vol", upper_btag, exact_grad_u) + normal = thaw(dcoll.normal(upper_btag), actx) upper_grad_u_dot_n = np.dot(upper_grad_u, normal) boundaries[lower_btag] = NeumannDiffusionBoundary(0.) boundaries[upper_btag] = NeumannDiffusionBoundary(upper_grad_u_dot_n) lower_btag = DTAG_BOUNDARY("-"+str(dim-1)) upper_btag = DTAG_BOUNDARY("+"+str(dim-1)) - upper_u = discr.project("vol", upper_btag, exact_u) + upper_u = op.project(dcoll, "vol", upper_btag, exact_u) boundaries[lower_btag] = DirichletDiffusionBoundary(0.) boundaries[upper_btag] = DirichletDiffusionBoundary(upper_u) @@ -201,7 +207,7 @@ def get_mesh(n): sym_u *= sym_sin(sym_coords[i]) sym_u *= sym_cos(sym_coords[dim-1]) - def get_boundaries(discr, actx, t): + def get_boundaries(dcoll, actx, t): boundaries = {} for i in range(dim-1): @@ -260,11 +266,11 @@ def test_diffusion_accuracy(actx_factory, problem, nsteps, dt, scales, order, for n in scales: mesh = p.get_mesh(n) - from grudge.eager import EagerDGDiscretization + from grudge.discretization import DiscretizationCollection from meshmode.discretization.poly_element import \ QuadratureSimplexGroupFactory, \ PolynomialWarpAndBlendGroupFactory - discr = EagerDGDiscretization( + dcoll = DiscretizationCollection( actx, mesh, discr_tag_to_group_factory={ DISCR_TAG_BASE: PolynomialWarpAndBlendGroupFactory(order), @@ -272,7 +278,7 @@ def test_diffusion_accuracy(actx_factory, problem, nsteps, dt, scales, order, } ) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) def sym_eval(expr, t): return sym.EvaluationMapper({"x": nodes, "t": t})(expr) @@ -285,8 +291,8 @@ def sym_eval(expr, t): discr_tag = DISCR_TAG_BASE def get_rhs(t, u): - return (diffusion_operator(discr, quad_tag=discr_tag, alpha=alpha, - boundaries=p.get_boundaries(discr, actx, t), u=u) + return (diffusion_operator(dcoll, quad_tag=discr_tag, alpha=alpha, + boundaries=p.get_boundaries(dcoll, actx, t), u=u) + sym_eval(sym_f, t)) t = 0. @@ -302,13 +308,13 @@ def get_rhs(t, u): expected_u = sym_eval(p.sym_u, t) rel_linf_err = ( - discr.norm(u - expected_u, np.inf) - / discr.norm(expected_u, np.inf)) + op.norm(dcoll, u - expected_u, np.inf) + / op.norm(dcoll, expected_u, np.inf)) eoc_rec.add_data_point(1./n, rel_linf_err) if visualize: from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, discr.order+3) + vis = make_visualizer(dcoll, dcoll.order+3) vis.write_vtk_file("diffusion_accuracy_{order}_{n}.vtu".format( order=order, n=n), [ ("u", u), @@ -335,10 +341,10 @@ def test_diffusion_discontinuous_alpha(actx_factory, order, visualize=False): mesh = get_box_mesh(1, -1, 1, n) - from grudge.eager import EagerDGDiscretization - discr = EagerDGDiscretization(actx, mesh, order=order) + from grudge.discretization import DiscretizationCollection + dcoll = DiscretizationCollection(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) # Set up a 1D heat equation interface problem, apply the diffusion operator to # the exact steady state solution, and check that it's zero @@ -371,13 +377,13 @@ def test_diffusion_discontinuous_alpha(actx_factory, order, visualize=False): def get_rhs(t, u): return diffusion_operator( - discr, quad_tag=DISCR_TAG_BASE, alpha=alpha, boundaries=boundaries, u=u) + dcoll, quad_tag=DISCR_TAG_BASE, alpha=alpha, boundaries=boundaries, u=u) rhs = get_rhs(0, u_steady) if visualize: from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, discr.order+3) + vis = make_visualizer(dcoll, dcoll.order+3) vis.write_vtk_file("diffusion_discontinuous_alpha_rhs_{order}.vtu" .format(order=order), [ ("alpha", alpha), @@ -385,7 +391,7 @@ def get_rhs(t, u): ("rhs", rhs), ]) - linf_err = discr.norm(rhs, np.inf) + linf_err = op.norm(dcoll, rhs, np.inf) assert(linf_err < 1e-11) # Now check stability @@ -415,7 +421,7 @@ def get_rhs(t, u): ("u_steady", u_steady), ]) - linf_diff = discr.norm(u - u_steady, np.inf) + linf_diff = op.norm(dcoll, u - u_steady, np.inf) assert linf_diff < 0.1 @@ -451,21 +457,21 @@ def test_diffusion_compare_to_nodal_dg(actx_factory, problem, order, t = 1.23456789 - from grudge.eager import EagerDGDiscretization - discr_mirgecom = EagerDGDiscretization(actx, mesh, order=order) - nodes_mirgecom = thaw(actx, discr_mirgecom.nodes()) + from grudge.discretization import DiscretizationCollection + dcoll_mirgecom = DiscretizationCollection(actx, mesh, order=order) + nodes_mirgecom = thaw(dcoll_mirgecom.nodes(), actx) def sym_eval_mirgecom(expr): return sym.EvaluationMapper({"x": nodes_mirgecom, "t": t})(expr) u_mirgecom = sym_eval_mirgecom(p.sym_u) - diffusion_u_mirgecom = diffusion_operator(discr_mirgecom, - quad_tag=DISCR_TAG_BASE, alpha=discr_mirgecom.zeros(actx)+1., - boundaries=p.get_boundaries(discr_mirgecom, actx, t), u=u_mirgecom) + diffusion_u_mirgecom = diffusion_operator(dcoll_mirgecom, + quad_tag=DISCR_TAG_BASE, alpha=dcoll_mirgecom.zeros(actx)+1., + boundaries=p.get_boundaries(dcoll_mirgecom, actx, t), u=u_mirgecom) discr_ndg = ndgctx.get_discr(actx) - nodes_ndg = thaw(actx, discr_ndg.nodes()) + nodes_ndg = thaw(discr_ndg.nodes(), actx) def sym_eval_ndg(expr): return sym.EvaluationMapper({"x": nodes_ndg, "t": t})(expr) @@ -513,10 +519,10 @@ def test_diffusion_obj_array_vectorize(actx_factory): mesh = p.get_mesh(n) - from grudge.eager import EagerDGDiscretization - discr = EagerDGDiscretization(actx, mesh, order=4) + from grudge.discretization import DiscretizationCollection + dcoll = DiscretizationCollection(actx, mesh, order=4) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) t = 1.23456789 @@ -528,24 +534,24 @@ def sym_eval(expr): u1 = sym_eval(sym_u1) u2 = sym_eval(sym_u2) - boundaries = p.get_boundaries(discr, actx, t) + boundaries = p.get_boundaries(dcoll, actx, t) - diffusion_u1 = diffusion_operator(discr, quad_tag=DISCR_TAG_BASE, alpha=alpha, + diffusion_u1 = diffusion_operator(dcoll, quad_tag=DISCR_TAG_BASE, alpha=alpha, boundaries=boundaries, u=u1) assert isinstance(diffusion_u1, DOFArray) expected_diffusion_u1 = sym_eval(sym_diffusion_u1) rel_linf_err = ( - discr.norm(diffusion_u1 - expected_diffusion_u1, np.inf) - / discr.norm(expected_diffusion_u1, np.inf)) + op.norm(dcoll, diffusion_u1 - expected_diffusion_u1, np.inf) + / op.norm(dcoll, expected_diffusion_u1, np.inf)) assert rel_linf_err < 1.e-5 boundaries_vector = [boundaries, boundaries] u_vector = make_obj_array([u1, u2]) diffusion_u_vector = diffusion_operator( - discr, quad_tag=DISCR_TAG_BASE, alpha=alpha, + dcoll, quad_tag=DISCR_TAG_BASE, alpha=alpha, boundaries=boundaries_vector, u=u_vector ) @@ -557,8 +563,8 @@ def sym_eval(expr): sym_eval(sym_diffusion_u2) ]) rel_linf_err = ( - discr.norm(diffusion_u_vector - expected_diffusion_u_vector, np.inf) - / discr.norm(expected_diffusion_u_vector, np.inf)) + op.norm(dcoll, diffusion_u_vector - expected_diffusion_u_vector, np.inf) + / op.norm(dcoll, expected_diffusion_u_vector, np.inf)) assert rel_linf_err < 1.e-5 diff --git a/test/test_eos.py b/test/test_eos.py index 80c182052..855331ed5 100644 --- a/test/test_eos.py +++ b/test/test_eos.py @@ -27,18 +27,17 @@ import logging import numpy as np import numpy.linalg as la # noqa -import pyopencl as cl -import pyopencl.clrandom -import pyopencl.clmath import pytest + +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from pytools.obj_array import make_obj_array from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from meshmode.dof_array import thaw -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) import cantera import pyrometheus as pyro @@ -47,12 +46,14 @@ Vortex2D, Lump, MixtureInitializer ) -from grudge.eager import EagerDGDiscretization from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests, ) from mirgecom.mechanisms import get_mechanism_cti +from grudge.discretization import DiscretizationCollection +import grudge.op as op + logger = logging.getLogger(__name__) @@ -60,7 +61,7 @@ [("uiuc", 1e-12), ("sanDiego", 1e-8)]) @pytest.mark.parametrize("y0", [0, 1]) -def test_pyrometheus_mechanisms(ctx_factory, mechname, rate_tol, y0): +def test_pyrometheus_mechanisms(actx_factory, mechname, rate_tol, y0): """Test known pyrometheus mechanisms. This test reproduces a pyrometheus-native test in the MIRGE context. @@ -68,9 +69,7 @@ def test_pyrometheus_mechanisms(ctx_factory, mechname, rate_tol, y0): Tests that the Pyrometheus mechanism code gets the same thermo properties as the corresponding mechanism in Cantera. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() dim = 1 nel_1d = 2 @@ -85,7 +84,7 @@ def test_pyrometheus_mechanisms(ctx_factory, mechname, rate_tol, y0): logger.info(f"Number of elements {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) # Pyrometheus initialization mech_cti = get_mechanism_cti(mechname) @@ -120,7 +119,7 @@ def test_pyrometheus_mechanisms(ctx_factory, mechname, rate_tol, y0): can_r = cantera_soln.net_rates_of_progress can_omega = cantera_soln.net_production_rates - ones = discr.zeros(actx) + 1.0 + ones = dcoll.zeros(actx) + 1.0 tin = can_t * ones pin = can_p * ones yin = make_obj_array([can_y[i] * ones for i in range(nspecies)]) @@ -149,33 +148,31 @@ def test_pyrometheus_mechanisms(ctx_factory, mechname, rate_tol, y0): print(f"can_omega = {can_omega}") print(f"prom_omega = {prom_omega}") - assert discr.norm((prom_c - can_c) / can_c, np.inf) < 1e-14 - assert discr.norm((prom_t - can_t) / can_t, np.inf) < 1e-14 - assert discr.norm((prom_rho - can_rho) / can_rho, np.inf) < 1e-14 - assert discr.norm((prom_p - can_p) / can_p, np.inf) < 1e-14 - assert discr.norm((prom_e - can_e) / can_e, np.inf) < 1e-6 - assert discr.norm((prom_k - can_k) / can_k, np.inf) < 1e-10 + assert op.norm(dcoll, (prom_c - can_c) / can_c, np.inf) < 1e-14 + assert op.norm(dcoll, (prom_t - can_t) / can_t, np.inf) < 1e-14 + assert op.norm(dcoll, (prom_rho - can_rho) / can_rho, np.inf) < 1e-14 + assert op.norm(dcoll, (prom_p - can_p) / can_p, np.inf) < 1e-14 + assert op.norm(dcoll, (prom_e - can_e) / can_e, np.inf) < 1e-6 + assert op.norm(dcoll, (prom_k - can_k) / can_k, np.inf) < 1e-10 # Pyro chem test comparisons for i, rate in enumerate(can_r): - assert discr.norm((prom_r[i] - rate), np.inf) < rate_tol + assert op.norm(dcoll, (prom_r[i] - rate), np.inf) < rate_tol for i, rate in enumerate(can_omega): - assert discr.norm((prom_omega[i] - rate), np.inf) < rate_tol + assert op.norm(dcoll, (prom_omega[i] - rate), np.inf) < rate_tol @pytest.mark.parametrize("mechname", ["uiuc", "sanDiego"]) @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("y0", [0, 1]) @pytest.mark.parametrize("vel", [0.0, 1.0]) -def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): +def test_pyrometheus_eos(actx_factory, mechname, dim, y0, vel): """Test PyrometheusMixture EOS for all available mechanisms. Tests that the PyrometheusMixture EOS gets the same thermo properties (p, T, e) as the Pyrometheus-native mechanism code. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 4 @@ -189,8 +186,8 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): logger.info(f"Number of elements {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Pyrometheus initialization mech_cti = get_mechanism_cti(mechname) @@ -215,7 +212,7 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): print(f"Testing {mechname}(t,P) = ({tempin}, {pressin})") - ones = discr.zeros(actx) + 1.0 + ones = dcoll.zeros(actx) + 1.0 tin = tempin * ones pin = pressin * ones yin = y0s * ones @@ -246,17 +243,17 @@ def test_pyrometheus_eos(ctx_factory, mechname, dim, y0, vel): print(f"pyro_eos.e = {internal_energy}") tol = 1e-14 - assert discr.norm((cv.mass - pyro_rho) / pyro_rho, np.inf) < tol - assert discr.norm((temperature - pyro_t) / pyro_t, np.inf) < tol - assert discr.norm((internal_energy - pyro_e) / pyro_e, np.inf) < tol - assert discr.norm((p - pyro_p) / pyro_p, np.inf) < tol + assert op.norm(dcoll, (cv.mass - pyro_rho) / pyro_rho, np.inf) < tol + assert op.norm(dcoll, (temperature - pyro_t) / pyro_t, np.inf) < tol + assert op.norm(dcoll, (internal_energy - pyro_e) / pyro_e, np.inf) < tol + assert op.norm(dcoll, (p - pyro_p) / pyro_p, np.inf) < tol @pytest.mark.parametrize(("mechname", "rate_tol"), [("uiuc", 1e-12), ("sanDiego", 1e-8)]) @pytest.mark.parametrize("y0", [0, 1]) -def test_pyrometheus_kinetics(ctx_factory, mechname, rate_tol, y0): +def test_pyrometheus_kinetics(actx_factory, mechname, rate_tol, y0): """Test known pyrometheus reaction mechanisms. This test reproduces a pyrometheus-native test in the MIRGE context. @@ -266,9 +263,7 @@ def test_pyrometheus_kinetics(ctx_factory, mechname, rate_tol, y0): are integrated in time and verified against a homogeneous reactor in Cantera. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() dim = 1 nel_1d = 4 @@ -283,8 +278,8 @@ def test_pyrometheus_kinetics(ctx_factory, mechname, rate_tol, y0): logger.info(f"Number of elements {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - ones = discr.zeros(actx) + 1.0 + dcoll = DiscretizationCollection(actx, mesh, order=order) + ones = dcoll.zeros(actx) + 1.0 # Pyrometheus initialization mech_cti = get_mechanism_cti(mechname) @@ -346,34 +341,32 @@ def test_pyrometheus_kinetics(ctx_factory, mechname, rate_tol, y0): # Print print(f"can_r = {can_r}") print(f"pyro_r = {pyro_r}") - abs_diff = discr.norm(pyro_r - can_r, np.inf) + abs_diff = op.norm(dcoll, pyro_r - can_r, np.inf) if abs_diff > 1e-14: min_r = (np.abs(can_r)).min() if min_r > 0: - assert discr.norm((pyro_r - can_r) / can_r, np.inf) < rate_tol + assert op.norm(dcoll, (pyro_r - can_r) / can_r, np.inf) < rate_tol else: - assert discr.norm(pyro_r, np.inf) < rate_tol + assert op.norm(dcoll, pyro_r, np.inf) < rate_tol print(f"can_omega = {can_omega}") print(f"pyro_omega = {pyro_omega}") for i, omega in enumerate(can_omega): omin = np.abs(omega).min() if omin > 1e-12: - assert discr.norm((pyro_omega[i] - omega) / omega, np.inf) < 1e-8 + assert op.norm(dcoll, (pyro_omega[i] - omega) / omega, np.inf) < 1e-8 else: - assert discr.norm(pyro_omega[i], np.inf) < 1e-12 + assert op.norm(dcoll, pyro_omega[i], np.inf) < 1e-12 @pytest.mark.parametrize("dim", [1, 2, 3]) -def test_idealsingle_lump(ctx_factory, dim): +def test_idealsingle_lump(actx_factory, dim): """Test IdealSingle EOS with mass lump. Tests that the IdealSingleGas EOS returns the correct (uniform) pressure for the Lump solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 4 @@ -386,8 +379,8 @@ def test_idealsingle_lump(ctx_factory, dim): order = 3 logger.info(f"Number of elements {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Init soln with Vortex center = np.zeros(shape=(dim,)) @@ -399,14 +392,14 @@ def test_idealsingle_lump(ctx_factory, dim): p = eos.pressure(cv) exp_p = 1.0 - errmax = discr.norm(p - exp_p, np.inf) + errmax = op.norm(dcoll, p - exp_p, np.inf) exp_ke = 0.5 * cv.mass ke = eos.kinetic_energy(cv) - kerr = discr.norm(ke - exp_ke, np.inf) + kerr = op.norm(dcoll, ke - exp_ke, np.inf) te = eos.total_energy(cv, p) - terr = discr.norm(te - cv.energy, np.inf) + terr = op.norm(dcoll, te - cv.energy, np.inf) logger.info(f"lump_soln = {cv}") logger.info(f"pressure = {p}") @@ -416,15 +409,13 @@ def test_idealsingle_lump(ctx_factory, dim): assert terr < 1e-15 -def test_idealsingle_vortex(ctx_factory): +def test_idealsingle_vortex(actx_factory): r"""Test EOS with isentropic vortex. Tests that the IdealSingleGas EOS returns the correct pressure (p) for the Vortex2D solution field (i.e. $p = \rho^{\gamma}$). """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() dim = 2 nel_1d = 4 @@ -438,8 +429,8 @@ def test_idealsingle_vortex(ctx_factory): order = 3 logger.info(f"Number of elements {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) eos = IdealSingleGas() # Init soln with Vortex vortex = Vortex2D() @@ -448,14 +439,14 @@ def test_idealsingle_vortex(ctx_factory): gamma = eos.gamma() p = eos.pressure(cv) exp_p = cv.mass ** gamma - errmax = discr.norm(p - exp_p, np.inf) + errmax = op.norm(dcoll, p - exp_p, np.inf) exp_ke = 0.5 * np.dot(cv.momentum, cv.momentum) / cv.mass ke = eos.kinetic_energy(cv) - kerr = discr.norm(ke - exp_ke, np.inf) + kerr = op.norm(dcoll, ke - exp_ke, np.inf) te = eos.total_energy(cv, p) - terr = discr.norm(te - cv.energy, np.inf) + terr = op.norm(dcoll, te - cv.energy, np.inf) logger.info(f"vortex_soln = {cv}") logger.info(f"pressure = {p}") diff --git a/test/test_euler.py b/test/test_euler.py index 243917c28..0d4a8009a 100644 --- a/test/test_euler.py +++ b/test/test_euler.py @@ -33,15 +33,19 @@ import math from functools import partial +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from pytools.obj_array import ( flat_obj_array, make_obj_array, ) -from meshmode.dof_array import thaw from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa -from grudge.eager import interior_trace_pair -from grudge.symbolic.primitives import TracePair + from mirgecom.euler import euler_operator from mirgecom.fluid import ( split_conserved, @@ -51,11 +55,10 @@ from mirgecom.initializers import Vortex2D, Lump, MulticomponentLump from mirgecom.boundary import PrescribedBoundary, DummyBoundary from mirgecom.eos import IdealSingleGas -from grudge.eager import EagerDGDiscretization -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) +from grudge.discretization import DiscretizationCollection +from grudge.symbolic.primitives import TracePair +import grudge.op as op from grudge.shortcuts import make_visualizer from mirgecom.inviscid import ( @@ -90,7 +93,7 @@ def test_inviscid_flux(actx_factory, nspecies, dim): ) order = 3 - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) eos = IdealSingleGas() logger.info(f"Number of {dim}d elems: {mesh.nelements}") @@ -100,7 +103,7 @@ def rand(): return DOFArray( actx, tuple(actx.from_numpy(np.random.rand(grp.nelements, grp.nunit_dofs)) - for grp in discr.discr_from_dd("vol").groups) + for grp in dcoll.discr_from_dd("vol").groups) ) mass = rand() @@ -136,7 +139,7 @@ def rand(): # }}} - flux = inviscid_flux(discr, eos, cv) + flux = inviscid_flux(dcoll, eos, cv) flux_resid = flux - expected_flux for i in range(numeq, dim): @@ -172,7 +175,7 @@ def test_inviscid_flux_components(actx_factory, dim): ) order = 3 - discr = EagerDGDiscretization(actx, mesh, order=order) + dcoll = DiscretizationCollection(actx, mesh, order=order) eos = IdealSingleGas() logger.info(f"Number of {dim}d elems: {mesh.nelements}") @@ -183,24 +186,24 @@ def test_inviscid_flux_components(actx_factory, dim): # the expected values (and p0 within tolerance) # === with V = 0, fixed P = p0 tolerance = 1e-15 - nodes = thaw(actx, discr.nodes()) - mass = discr.zeros(actx) + np.dot(nodes, nodes) + 1.0 - mom = make_obj_array([discr.zeros(actx) for _ in range(dim)]) - p_exact = discr.zeros(actx) + p0 + nodes = thaw(dcoll.nodes(), actx) + mass = dcoll.zeros(actx) + np.dot(nodes, nodes) + 1.0 + mom = make_obj_array([dcoll.zeros(actx) for _ in range(dim)]) + p_exact = dcoll.zeros(actx) + p0 energy = p_exact / 0.4 + 0.5 * np.dot(mom, mom) / mass cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) p = eos.pressure(cv) - flux = inviscid_flux(discr, eos, cv) - assert discr.norm(p - p_exact, np.inf) < tolerance + flux = inviscid_flux(dcoll, eos, cv) + assert op.norm(dcoll, p - p_exact, np.inf) < tolerance logger.info(f"{dim}d flux = {flux}") # for velocity zero, these components should be == zero - assert discr.norm(flux.mass, 2) == 0.0 - assert discr.norm(flux.energy, 2) == 0.0 + assert op.norm(dcoll, flux.mass, 2) == 0.0 + assert op.norm(dcoll, flux.energy, 2) == 0.0 # The momentum diagonal should be p # Off-diagonal should be identically 0 - assert discr.norm(flux.momentum - p0*np.identity(dim), np.inf) < tolerance + assert op.norm(dcoll, flux.momentum - p0*np.identity(dim), np.inf) < tolerance @pytest.mark.parametrize(("dim", "livedim"), [ @@ -231,34 +234,34 @@ def test_inviscid_mom_flux_components(actx_factory, dim, livedim): ) order = 3 - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) tolerance = 1e-15 for livedim in range(dim): - mass = discr.zeros(actx) + 1.0 + np.dot(nodes, nodes) - mom = make_obj_array([discr.zeros(actx) for _ in range(dim)]) + mass = dcoll.zeros(actx) + 1.0 + np.dot(nodes, nodes) + mom = make_obj_array([dcoll.zeros(actx) for _ in range(dim)]) mom[livedim] = mass - p_exact = discr.zeros(actx) + p0 + p_exact = dcoll.zeros(actx) + p0 energy = ( p_exact / (eos.gamma() - 1.0) + 0.5 * np.dot(mom, mom) / mass ) cv = make_conserved(dim, mass=mass, energy=energy, momentum=mom) p = eos.pressure(cv) - assert discr.norm(p - p_exact, np.inf) < tolerance - flux = inviscid_flux(discr, eos, cv) + assert op.norm(dcoll, p - p_exact, np.inf) < tolerance + flux = inviscid_flux(dcoll, eos, cv) logger.info(f"{dim}d flux = {flux}") vel_exact = mom / mass # first two components should be nonzero in livedim only - assert discr.norm(flux.mass - mom, np.inf) == 0 + assert op.norm(dcoll, flux.mass - mom, np.inf) == 0 eflux_exact = (energy + p_exact)*vel_exact - assert discr.norm(flux.energy - eflux_exact, np.inf) == 0 + assert op.norm(dcoll, flux.energy - eflux_exact, np.inf) == 0 logger.info("Testing momentum") xpmomflux = mass*np.outer(vel_exact, vel_exact) + p_exact*np.identity(dim) - assert discr.norm(flux.momentum - xpmomflux, np.inf) < tolerance + assert op.norm(dcoll, flux.momentum - xpmomflux, np.inf) < tolerance @pytest.mark.parametrize("nspecies", [0, 10]) @@ -291,14 +294,14 @@ def test_facial_flux(actx_factory, nspecies, order, dim): logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - zeros = discr.zeros(actx) + dcoll = DiscretizationCollection(actx, mesh, order=order) + zeros = dcoll.zeros(actx) ones = zeros + 1.0 - mass_input = discr.zeros(actx) + 1.0 - energy_input = discr.zeros(actx) + 2.5 + mass_input = dcoll.zeros(actx) + 1.0 + energy_input = dcoll.zeros(actx) + 2.5 mom_input = flat_obj_array( - [discr.zeros(actx) for i in range(discr.dim)] + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) mass_frac_input = flat_obj_array( [ones / ((i + 1) * 10) for i in range(nspecies)] @@ -313,11 +316,13 @@ def test_facial_flux(actx_factory, nspecies, order, dim): # Check the boundary facial fluxes as called on an interior boundary interior_face_flux = _facial_flux( - discr, eos=IdealSingleGas(), cv_tpair=interior_trace_pair(discr, cv)) + dcoll, eos=IdealSingleGas(), + cv_tpair=op.interior_trace_pair(dcoll, cv) + ) def inf_norm(data): if len(data) > 0: - return discr.norm(data, np.inf, dd="all_faces") + return op.norm(dcoll, data, np.inf, dd="all_faces") else: return 0.0 @@ -337,19 +342,19 @@ def inf_norm(data): # (Explanation courtesy of Mike Campbell, # https://github.com/illinois-ceesd/mirgecom/pull/44#discussion_r463304292) - nhat = thaw(actx, discr.normal("int_faces")) - mom_flux_exact = discr.project("int_faces", "all_faces", p0*nhat) + nhat = thaw(dcoll.normal("int_faces"), actx) + mom_flux_exact = op.project(dcoll, "int_faces", "all_faces", p0*nhat) print(f"{mom_flux_exact=}") print(f"{interior_face_flux.momentum=}") momerr = inf_norm(interior_face_flux.momentum - mom_flux_exact) assert momerr < tolerance eoc_rec0.add_data_point(1.0 / nel_1d, momerr) - # Check the boundary facial fluxes as called on a domain boundary - dir_mass = discr.project("vol", BTAG_ALL, mass_input) - dir_e = discr.project("vol", BTAG_ALL, energy_input) - dir_mom = discr.project("vol", BTAG_ALL, mom_input) - dir_mf = discr.project("vol", BTAG_ALL, species_mass_input) + # Check the boundary facial fluxes as called on a boundary + dir_mass = op.project(dcoll, "vol", BTAG_ALL, mass_input) + dir_e = op.project(dcoll, "vol", BTAG_ALL, energy_input) + dir_mom = op.project(dcoll, "vol", BTAG_ALL, mom_input) + dir_mf = op.project(dcoll, "vol", BTAG_ALL, species_mass_input) dir_bval = make_conserved(dim, mass=dir_mass, energy=dir_e, momentum=dir_mom, species_mass=dir_mf) @@ -357,7 +362,7 @@ def inf_norm(data): momentum=dir_mom, species_mass=dir_mf) boundary_flux = _facial_flux( - discr, eos=IdealSingleGas(), + dcoll, eos=IdealSingleGas(), cv_tpair=TracePair(BTAG_ALL, interior=dir_bval, exterior=dir_bc) ) @@ -365,8 +370,8 @@ def inf_norm(data): assert inf_norm(boundary_flux.energy) < tolerance assert inf_norm(boundary_flux.species_mass) < tolerance - nhat = thaw(actx, discr.normal(BTAG_ALL)) - mom_flux_exact = discr.project(BTAG_ALL, "all_faces", p0*nhat) + nhat = thaw(dcoll.normal(BTAG_ALL), actx) + mom_flux_exact = op.project(dcoll, BTAG_ALL, "all_faces", p0*nhat) momerr = inf_norm(boundary_flux.momentum - mom_flux_exact) assert momerr < tolerance @@ -411,15 +416,15 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): f"Number of {dim}d elements: {mesh.nelements}" ) - discr = EagerDGDiscretization(actx, mesh, order=order) - zeros = discr.zeros(actx) + dcoll = DiscretizationCollection(actx, mesh, order=order) + zeros = dcoll.zeros(actx) ones = zeros + 1.0 - mass_input = discr.zeros(actx) + 1 - energy_input = discr.zeros(actx) + 2.5 + mass_input = dcoll.zeros(actx) + 1 + energy_input = dcoll.zeros(actx) + 2.5 mom_input = make_obj_array( - [discr.zeros(actx) for i in range(discr.dim)] + [dcoll.zeros(actx) for i in range(dcoll.dim)] ) mass_frac_input = flat_obj_array( @@ -433,12 +438,12 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): species_mass=species_mass_input) expected_rhs = split_conserved( - dim, make_obj_array([discr.zeros(actx) + dim, make_obj_array([dcoll.zeros(actx) for i in range(num_equations)]) ) boundaries = {BTAG_ALL: DummyBoundary()} - inviscid_rhs = euler_operator(discr, eos=IdealSingleGas(), + inviscid_rhs = euler_operator(dcoll, eos=IdealSingleGas(), boundaries=boundaries, cv=cv, t=0.0) rhs_resid = inviscid_rhs - expected_rhs @@ -459,26 +464,26 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): f"rhoy_rhs = {rhoy_rhs}\n" ) - assert discr.norm(rho_resid, np.inf) < tolerance - assert discr.norm(rhoe_resid, np.inf) < tolerance + assert op.norm(dcoll, rho_resid, np.inf) < tolerance + assert op.norm(dcoll, rhoe_resid, np.inf) < tolerance for i in range(dim): - assert discr.norm(mom_resid[i], np.inf) < tolerance + assert op.norm(dcoll, mom_resid[i], np.inf) < tolerance for i in range(nspecies): - assert discr.norm(rhoy_resid[i], np.inf) < tolerance + assert op.norm(dcoll, rhoy_resid[i], np.inf) < tolerance - err_max = discr.norm(rho_resid, np.inf) + err_max = op.norm(dcoll, rho_resid, np.inf) eoc_rec0.add_data_point(1.0 / nel_1d, err_max) # set a non-zero, but uniform velocity component for i in range(len(mom_input)): - mom_input[i] = discr.zeros(actx) + (-1.0) ** i + mom_input[i] = dcoll.zeros(actx) + (-1.0) ** i cv = make_conserved( dim, mass=mass_input, energy=energy_input, momentum=mom_input, species_mass=species_mass_input) boundaries = {BTAG_ALL: DummyBoundary()} - inviscid_rhs = euler_operator(discr, eos=IdealSingleGas(), + inviscid_rhs = euler_operator(dcoll, eos=IdealSingleGas(), boundaries=boundaries, cv=cv, t=0.0) rhs_resid = inviscid_rhs - expected_rhs @@ -487,15 +492,15 @@ def test_uniform_rhs(actx_factory, nspecies, dim, order): mom_resid = rhs_resid.momentum rhoy_resid = rhs_resid.species_mass - assert discr.norm(rho_resid, np.inf) < tolerance - assert discr.norm(rhoe_resid, np.inf) < tolerance + assert op.norm(dcoll, rho_resid, np.inf) < tolerance + assert op.norm(dcoll, rhoe_resid, np.inf) < tolerance for i in range(dim): - assert discr.norm(mom_resid[i], np.inf) < tolerance + assert op.norm(dcoll, mom_resid[i], np.inf) < tolerance for i in range(nspecies): - assert discr.norm(rhoy_resid[i], np.inf) < tolerance + assert op.norm(dcoll, rhoy_resid[i], np.inf) < tolerance - err_max = discr.norm(rho_resid, np.inf) + err_max = op.norm(dcoll, rho_resid, np.inf) eoc_rec1.add_data_point(1.0 / nel_1d, err_max) logger.info( @@ -538,8 +543,8 @@ def test_vortex_rhs(actx_factory, order): f"Number of {dim}d elements: {mesh.nelements}" ) - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Init soln with Vortex and expected RHS = 0 vortex = Vortex2D(center=[0, 0], velocity=[0, 0]) @@ -547,10 +552,10 @@ def test_vortex_rhs(actx_factory, order): boundaries = {BTAG_ALL: PrescribedBoundary(vortex)} inviscid_rhs = euler_operator( - discr, eos=IdealSingleGas(), boundaries=boundaries, + dcoll, eos=IdealSingleGas(), boundaries=boundaries, cv=vortex_soln, t=0.0) - err_max = discr.norm(inviscid_rhs.join(), np.inf) + err_max = op.norm(dcoll, inviscid_rhs.join(), np.inf) eoc_rec.add_data_point(1.0 / nel_1d, err_max) logger.info( @@ -591,8 +596,8 @@ def test_lump_rhs(actx_factory, dim, order): logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Init soln with Lump and expected RHS = 0 center = np.zeros(shape=(dim,)) @@ -601,10 +606,10 @@ def test_lump_rhs(actx_factory, dim, order): lump_soln = lump(nodes) boundaries = {BTAG_ALL: PrescribedBoundary(lump)} inviscid_rhs = euler_operator( - discr, eos=IdealSingleGas(), boundaries=boundaries, cv=lump_soln, t=0.0) - expected_rhs = lump.exact_rhs(discr, cv=lump_soln, t=0) + dcoll, eos=IdealSingleGas(), boundaries=boundaries, cv=lump_soln, t=0.0) + expected_rhs = lump.exact_rhs(dcoll, cv=lump_soln, t=0) - err_max = discr.norm((inviscid_rhs-expected_rhs).join(), np.inf) + err_max = op.norm(dcoll, (inviscid_rhs-expected_rhs).join(), np.inf) if err_max > maxxerr: maxxerr = err_max @@ -650,8 +655,8 @@ def test_multilump_rhs(actx_factory, dim, order, v0): logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) centers = make_obj_array([np.zeros(shape=(dim,)) for i in range(nspecies)]) spec_y0s = np.ones(shape=(nspecies,)) @@ -669,11 +674,11 @@ def test_multilump_rhs(actx_factory, dim, order, v0): boundaries = {BTAG_ALL: PrescribedBoundary(lump)} inviscid_rhs = euler_operator( - discr, eos=IdealSingleGas(), boundaries=boundaries, cv=lump_soln, t=0.0) - expected_rhs = lump.exact_rhs(discr, cv=lump_soln, t=0) + dcoll, eos=IdealSingleGas(), boundaries=boundaries, cv=lump_soln, t=0.0) + expected_rhs = lump.exact_rhs(dcoll, cv=lump_soln, t=0) print(f"inviscid_rhs = {inviscid_rhs}") print(f"expected_rhs = {expected_rhs}") - err_max = discr.norm((inviscid_rhs-expected_rhs).join(), np.inf) + err_max = op.norm(dcoll, (inviscid_rhs-expected_rhs).join(), np.inf) if err_max > maxxerr: maxxerr = err_max @@ -719,11 +724,11 @@ def _euler_flow_stepper(actx, parameters): dim = mesh.dim istep = 0 - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) cv = initializer(nodes) - sdt = cfl * get_inviscid_timestep(discr, eos=eos, cv=cv) + sdt = cfl * get_inviscid_timestep(dcoll, eos=eos, cv=cv) initname = initializer.__class__.__name__ eosname = eos.__class__.__name__ @@ -736,7 +741,7 @@ def _euler_flow_stepper(actx, parameters): f"EOS: {eosname}" ) - vis = make_visualizer(discr, order) + vis = make_visualizer(dcoll, order) def write_soln(state, write_status=True): dv = eos.dependent_vars(cv=state) @@ -767,7 +772,7 @@ def write_soln(state, write_status=True): return maxerr def rhs(t, q): - return euler_operator(discr, eos=eos, boundaries=boundaries, cv=q, t=t) + return euler_operator(dcoll, eos=eos, boundaries=boundaries, cv=q, t=t) filter_order = 8 eta = .5 @@ -797,20 +802,20 @@ def rhs(t, q): cv = rk4_step(cv, t, dt, rhs) cv = split_conserved( - dim, filter_modally(discr, "vol", cutoff, frfunc, cv.join()) + dim, filter_modally(dcoll, "vol", cutoff, frfunc, cv.join()) ) t += dt istep += 1 - sdt = cfl * get_inviscid_timestep(discr, eos=eos, cv=cv) + sdt = cfl * get_inviscid_timestep(dcoll, eos=eos, cv=cv) if nstepstatus > 0: logger.info("Writing final dump.") maxerr = max(write_soln(cv, False)) else: expected_result = initializer(nodes, t=t) - maxerr = discr.norm((cv - expected_result).join(), np.inf) + maxerr = op.norm(dcoll, (cv - expected_result).join(), np.inf) logger.info(f"Max Error: {maxerr}") if maxerr > exittol: diff --git a/test/test_filter.py b/test/test_filter.py index 78b7d5028..31b151106 100644 --- a/test/test_filter.py +++ b/test/test_filter.py @@ -30,17 +30,21 @@ import numpy as np from functools import partial -from meshmode.dof_array import thaw -from grudge.eager import EagerDGDiscretization -from meshmode.array_context import ( # noqa +from arraycontext import ( # noqa + thaw, pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) + as pytest_generate_tests +) + from pytools.obj_array import ( make_obj_array ) -from meshmode.dof_array import thaw # noqa + from mirgecom.filter import make_spectral_filter +from grudge.discretization import DiscretizationCollection +import grudge.op as op + @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("order", [2, 3, 4]) @@ -65,8 +69,8 @@ def test_filter_coeff(actx_factory, filter_order, order, dim): a=(-0.5,) * dim, b=(0.5,) * dim, nelements_per_axis=(nel_1d,) * dim ) - discr = EagerDGDiscretization(actx, mesh, order=order) - vol_discr = discr.discr_from_dd("vol") + dcoll = DiscretizationCollection(actx, mesh, order=order) + vol_discr = dcoll.discr_from_dd("vol") eta = .5 # just filter half the modes # number of modes see e.g.: @@ -125,7 +129,7 @@ def test_filter_coeff(actx_factory, filter_order, order, dim): frfunc = partial(xmrfunc, alpha=alpha, filter_order=filter_order) for group in vol_discr.groups: - mode_ids = group.mode_ids() + mode_ids = group.basis_obj().mode_ids filter_coeff = actx.thaw( make_spectral_filter( actx, group, cutoff=cutoff, @@ -167,8 +171,8 @@ def test_filter_function(actx_factory, dim, order, do_viz=False): a=(0.0,) * dim, b=(1.0,) * dim, nelements_per_axis=(nel_1d,) * dim ) - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # number of modes see e.g.: # JSH/TW Nodal DG Methods, Section 10.1 @@ -189,10 +193,10 @@ def test_filter_function(actx_factory, dim, order, do_viz=False): uniform_soln = initr(t=0, x_vec=nodes).join() from mirgecom.filter import filter_modally - filtered_soln = filter_modally(discr, "vol", cutoff, + filtered_soln = filter_modally(dcoll, "vol", cutoff, frfunc, uniform_soln) soln_resid = uniform_soln - filtered_soln - max_errors = [discr.norm(v, np.inf) for v in soln_resid] + max_errors = [op.norm(dcoll, v, np.inf) for v in soln_resid] tol = 1e-14 @@ -206,7 +210,7 @@ def polyfn(coeff): # , x_vec): r = nodes[0] result = 0 for n, a in enumerate(coeff): - result += a * r ** n + result = result + a * r ** n return make_obj_array([result]) # Any order {cutoff} and below fields should be unharmed @@ -214,10 +218,10 @@ def polyfn(coeff): # , x_vec): field_order = int(cutoff) coeff = [1.0 / (i + 1) for i in range(field_order + 1)] field = polyfn(coeff=coeff) - filtered_field = filter_modally(discr, "vol", cutoff, + filtered_field = filter_modally(dcoll, "vol", cutoff, frfunc, field) soln_resid = field - filtered_field - max_errors = [discr.norm(v, np.inf) for v in soln_resid] + max_errors = [op.norm(dcoll, v, np.inf) for v in soln_resid] logger.info(f"Field = {field}") logger.info(f"Filtered = {filtered_field}") logger.info(f"Max Errors (poly) = {max_errors}") @@ -228,16 +232,16 @@ def polyfn(coeff): # , x_vec): tol = 1e-1 if do_viz is True: from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, discr.order) + vis = make_visualizer(dcoll) from grudge.dof_desc import DD_VOLUME_MODAL, DD_VOLUME - modal_map = discr.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) + modal_map = dcoll.connection_from_dds(DD_VOLUME, DD_VOLUME_MODAL) for field_order in range(cutoff+1, cutoff+4): coeff = [1.0 / (i + 1) for i in range(field_order+1)] field = polyfn(coeff=coeff) - filtered_field = filter_modally(discr, "vol", cutoff, + filtered_field = filter_modally(dcoll, "vol", cutoff, frfunc, field) unfiltered_spectrum = modal_map(field) @@ -253,6 +257,6 @@ def polyfn(coeff): # , x_vec): ] vis.write_vtk_file(f"filter_test_{field_order}.vtu", io_fields) field_resid = unfiltered_spectrum - filtered_spectrum - max_errors = [discr.norm(v, np.inf) for v in field_resid] + max_errors = [op.norm(dcoll, v, np.inf) for v in field_resid] # fields should be different, but not too different assert(tol > np.max(max_errors) > threshold) diff --git a/test/test_fluid.py b/test/test_fluid.py index 682220ed9..b9cc335ab 100644 --- a/test/test_fluid.py +++ b/test/test_fluid.py @@ -31,17 +31,18 @@ import logging import pytest -from pytools.obj_array import ( - make_obj_array, - obj_array_vectorize +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests ) -from meshmode.dof_array import thaw +from pytools.obj_array import make_obj_array + from mirgecom.fluid import split_conserved, join_conserved -from grudge.eager import EagerDGDiscretization -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) + +from grudge.discretization import DiscretizationCollection +import grudge.op as op logger = logging.getLogger(__name__) @@ -64,14 +65,14 @@ def test_velocity_gradient_sanity(actx_factory, dim, mass_exp, vel_fac): ) order = 3 - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) - zeros = discr.zeros(actx) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) + zeros = dcoll.zeros(actx) ones = zeros + 1.0 mass = 1*ones for i in range(mass_exp): - mass *= (mass + i) + mass = mass * (mass + i) energy = zeros + 2.5 velocity = vel_fac * nodes mom = mass * velocity @@ -79,14 +80,14 @@ def test_velocity_gradient_sanity(actx_factory, dim, mass_exp, vel_fac): q = join_conserved(dim, mass=mass, energy=energy, momentum=mom) cv = split_conserved(dim, q) - grad_q = obj_array_vectorize(discr.grad, q) + grad_q = op.local_grad(dcoll, q) grad_cv = split_conserved(dim, grad_q) - grad_v = velocity_gradient(discr, cv, grad_cv) + grad_v = velocity_gradient(dcoll, cv, grad_cv) tol = 1e-12 exp_result = vel_fac * np.eye(dim) * ones - grad_v_err = [discr.norm(grad_v[i] - exp_result[i], np.inf) + grad_v_err = [op.norm(dcoll, grad_v[i] - exp_result[i], np.inf) for i in range(dim)] assert max(grad_v_err) < tol @@ -114,9 +115,9 @@ def test_velocity_gradient_eoc(actx_factory, dim): a=(1.0,) * dim, b=(2.0,) * dim, nelements_per_axis=(nel_1d,) * dim ) - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) - zeros = discr.zeros(actx) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) + zeros = dcoll.zeros(actx) energy = zeros + 2.5 mass = nodes[dim-1]*nodes[dim-1] @@ -126,10 +127,10 @@ def test_velocity_gradient_eoc(actx_factory, dim): q = join_conserved(dim, mass=mass, energy=energy, momentum=mom) cv = split_conserved(dim, q) - grad_q = obj_array_vectorize(discr.grad, q) + grad_q = op.local_grad(dcoll, q) grad_cv = split_conserved(dim, grad_q) - grad_v = velocity_gradient(discr, cv, grad_cv) + grad_v = velocity_gradient(dcoll, cv, grad_cv) def exact_grad_row(xdata, gdim, dim): exact_grad_row = make_obj_array([zeros for _ in range(dim)]) @@ -137,7 +138,7 @@ def exact_grad_row(xdata, gdim, dim): return exact_grad_row comp_err = make_obj_array([ - discr.norm(grad_v[i] - exact_grad_row(nodes[i], i, dim), np.inf) + op.norm(dcoll, grad_v[i] - exact_grad_row(nodes[i], i, dim), np.inf) for i in range(dim)]) err_max = comp_err.max() eoc.add_data_point(h, err_max) diff --git a/test/test_init.py b/test/test_init.py index 76802f03c..79cbdf63f 100644 --- a/test/test_init.py +++ b/test/test_init.py @@ -25,14 +25,15 @@ import logging import numpy as np import numpy.linalg as la # noqa -import pyopencl as cl -import pyopencl.clrandom -import pyopencl.clmath from pytools.obj_array import make_obj_array import pytest -from meshmode.array_context import PyOpenCLArrayContext -from meshmode.dof_array import thaw +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from meshmode.mesh import BTAG_ALL, BTAG_NONE # noqa from mirgecom.initializers import Vortex2D @@ -44,25 +45,21 @@ from mirgecom.initializers import SodShock1D from mirgecom.eos import IdealSingleGas -from grudge.eager import EagerDGDiscretization -from pyopencl.tools import ( # noqa - pytest_generate_tests_for_pyopencl as pytest_generate_tests, -) +from grudge.discretization import DiscretizationCollection +import grudge.op as op logger = logging.getLogger(__name__) @pytest.mark.parametrize("dim", [1, 2, 3]) @pytest.mark.parametrize("nspecies", [0, 10]) -def test_uniform_init(ctx_factory, dim, nspecies): +def test_uniform_init(actx_factory, dim, nspecies): """Test the uniform flow initializer. Simple test to check that uniform initializer creates the expected solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 4 from meshmode.mesh.generation import generate_regular_rect_mesh @@ -74,8 +71,8 @@ def test_uniform_init(ctx_factory, dim, nspecies): order = 3 logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) velocity = np.ones(shape=(dim,)) from mirgecom.initializers import Uniform @@ -86,7 +83,7 @@ def test_uniform_init(ctx_factory, dim, nspecies): def inf_norm(data): if len(data) > 0: - return discr.norm(data, np.inf) + return op.norm(dcoll, data, np.inf) else: return 0.0 @@ -109,14 +106,12 @@ def inf_norm(data): assert mferrmax < 1e-15 -def test_lump_init(ctx_factory): +def test_lump_init(actx_factory): """ Simple test to check that Lump initializer creates the expected solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() dim = 2 nel_1d = 4 @@ -129,8 +124,8 @@ def test_lump_init(ctx_factory): order = 3 logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Init soln with Vortex center = np.zeros(shape=(dim,)) @@ -142,7 +137,7 @@ def test_lump_init(ctx_factory): p = 0.4 * (cv.energy - 0.5 * np.dot(cv.momentum, cv.momentum) / cv.mass) exp_p = 1.0 - errmax = discr.norm(p - exp_p, np.inf) + errmax = op.norm(dcoll, p - exp_p, np.inf) logger.info(f"lump_soln = {cv}") logger.info(f"pressure = {p}") @@ -150,14 +145,12 @@ def test_lump_init(ctx_factory): assert errmax < 1e-15 -def test_vortex_init(ctx_factory): +def test_vortex_init(actx_factory): """ Simple test to check that Vortex2D initializer creates the expected solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() dim = 2 nel_1d = 4 @@ -170,8 +163,8 @@ def test_vortex_init(ctx_factory): order = 3 logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) # Init soln with Vortex vortex = Vortex2D() @@ -179,7 +172,7 @@ def test_vortex_init(ctx_factory): gamma = 1.4 p = 0.4 * (cv.energy - 0.5 * np.dot(cv.momentum, cv.momentum) / cv.mass) exp_p = cv.mass ** gamma - errmax = discr.norm(p - exp_p, np.inf) + errmax = op.norm(dcoll, p - exp_p, np.inf) logger.info(f"vortex_soln = {cv}") logger.info(f"pressure = {p}") @@ -187,14 +180,12 @@ def test_vortex_init(ctx_factory): assert errmax < 1e-15 -def test_shock_init(ctx_factory): +def test_shock_init(actx_factory): """ Simple test to check that Shock1D initializer creates the expected solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 10 dim = 2 @@ -208,8 +199,8 @@ def test_shock_init(ctx_factory): order = 3 print(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) initr = SodShock1D() cv = initr(t=0.0, x_vec=nodes) @@ -221,18 +212,16 @@ def test_shock_init(ctx_factory): eos = IdealSingleGas() p = eos.pressure(cv) - assert discr.norm(actx.np.where(nodes_x < 0.5, p-xpl, p-xpr), np.inf) < tol + assert op.norm(dcoll, actx.np.where(nodes_x < 0.5, p-xpl, p-xpr), np.inf) < tol @pytest.mark.parametrize("dim", [1, 2, 3]) -def test_uniform(ctx_factory, dim): +def test_uniform(actx_factory, dim): """ Simple test to check that Uniform initializer creates the expected solution field. """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 2 @@ -245,8 +234,8 @@ def test_uniform(ctx_factory, dim): order = 1 print(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) print(f"DIM = {dim}, {len(nodes)}") print(f"Nodes={nodes}") @@ -255,26 +244,24 @@ def test_uniform(ctx_factory, dim): cv = initr(t=0.0, x_vec=nodes) tol = 1e-15 - assert discr.norm(cv.mass - 1.0, np.inf) < tol - assert discr.norm(cv.energy - 2.5, np.inf) < tol + assert op.norm(dcoll, cv.mass - 1.0, np.inf) < tol + assert op.norm(dcoll, cv.energy - 2.5, np.inf) < tol print(f"Uniform Soln:{cv}") eos = IdealSingleGas() p = eos.pressure(cv) print(f"Press:{p}") - assert discr.norm(p - 1.0, np.inf) < tol + assert op.norm(dcoll, p - 1.0, np.inf) < tol @pytest.mark.parametrize("dim", [1, 2, 3]) -def test_pulse(ctx_factory, dim): +def test_pulse(actx_factory, dim): """ Test of Gaussian pulse generator. If it looks, walks, and quacks like a Gaussian, then ... """ - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 10 @@ -287,8 +274,8 @@ def test_pulse(ctx_factory, dim): order = 1 print(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) print(f"DIM = {dim}, {len(nodes)}") print(f"Nodes={nodes}") @@ -307,33 +294,31 @@ def test_pulse(ctx_factory, dim): print(f"exact: {pulse_check}") pulse_resid = pulse - pulse_check print(f"pulse residual: {pulse_resid}") - assert(discr.norm(pulse_resid, np.inf) < tol) + assert(op.norm(dcoll, pulse_resid, np.inf) < tol) # proper scaling with amplitude? amp = 2.0 pulse = 0 pulse = make_pulse(amp=amp, r0=r0, w=w, r=nodes) pulse_resid = pulse - (pulse_check + pulse_check) - assert(discr.norm(pulse_resid, np.inf) < tol) + assert(op.norm(dcoll, pulse_resid, np.inf) < tol) # proper scaling with r? amp = 1.0 rcheck = np.sqrt(2.0) * nodes pulse = make_pulse(amp=amp, r0=r0, w=w, r=rcheck) - assert(discr.norm(pulse - (pulse_check * pulse_check), np.inf) < tol) + assert(op.norm(dcoll, pulse - (pulse_check * pulse_check), np.inf) < tol) # proper scaling with w? w = w / np.sqrt(2.0) pulse = make_pulse(amp=amp, r0=r0, w=w, r=nodes) - assert(discr.norm(pulse - (pulse_check * pulse_check), np.inf) < tol) + assert(op.norm(dcoll, pulse - (pulse_check * pulse_check), np.inf) < tol) @pytest.mark.parametrize("dim", [1, 2, 3]) -def test_multilump(ctx_factory, dim): +def test_multilump(actx_factory, dim): """Test the multi-component lump initializer.""" - cl_ctx = ctx_factory() - queue = cl.CommandQueue(cl_ctx) - actx = PyOpenCLArrayContext(queue) + actx = actx_factory() nel_1d = 4 nspecies = 10 @@ -347,8 +332,8 @@ def test_multilump(ctx_factory, dim): order = 3 logger.info(f"Number of elements: {mesh.nelements}") - discr = EagerDGDiscretization(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + dcoll = DiscretizationCollection(actx, mesh, order=order) + nodes = thaw(dcoll.nodes(), actx) rho0 = 1.5 centers = make_obj_array([np.zeros(shape=(dim,)) for i in range(nspecies)]) @@ -372,11 +357,11 @@ def test_multilump(ctx_factory, dim): print(f"get_num_species = {numcvspec}") assert get_num_species(dim, cv.join()) == nspecies - assert discr.norm(cv.mass - rho0) == 0.0 + assert op.norm(dcoll, cv.mass - rho0) == 0.0 p = 0.4 * (cv.energy - 0.5 * np.dot(cv.momentum, cv.momentum) / cv.mass) exp_p = 1.0 - errmax = discr.norm(p - exp_p, np.inf) + errmax = op.norm(dcoll, p - exp_p, np.inf) assert len(cv.species_mass) == nspecies species_mass = cv.species_mass @@ -391,7 +376,7 @@ def test_multilump(ctx_factory, dim): print(f"exp_mass = {exp_mass}") print(f"mass_resid = {mass_resid}") - assert discr.norm(mass_resid, np.inf) == 0.0 + assert op.norm(dcoll, mass_resid, np.inf) == 0.0 logger.info(f"lump_soln = {cv}") logger.info(f"pressure = {p}") diff --git a/test/test_wave.py b/test/test_wave.py index d33c98259..aaf0cf9b5 100644 --- a/test/test_wave.py +++ b/test/test_wave.py @@ -23,16 +23,20 @@ import numpy as np import pyopencl.array as cla # noqa import pyopencl.clmath as clmath # noqa + +from arraycontext import ( # noqa + thaw, + pytest_generate_tests_for_pyopencl_array_context + as pytest_generate_tests +) + from pytools.obj_array import flat_obj_array, make_obj_array import pymbolic as pmbl import pymbolic.primitives as prim import mirgecom.symbolic as sym from mirgecom.wave import wave_operator -from meshmode.dof_array import thaw -from meshmode.array_context import ( # noqa - pytest_generate_tests_for_pyopencl_array_context - as pytest_generate_tests) +import grudge.op as op import pytest @@ -164,10 +168,10 @@ def test_wave_accuracy(actx_factory, problem, order, visualize=False): for n in [8, 10, 12] if p.dim == 3 else [8, 12, 16]: mesh = p.mesh_factory(n) - from grudge.eager import EagerDGDiscretization - discr = EagerDGDiscretization(actx, mesh, order=order) + from grudge.discretization import DiscretizationCollection + dcoll = DiscretizationCollection(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) def sym_eval(expr, t): return sym.EvaluationMapper({"c": p.c, "x": nodes, "t": t})(expr) @@ -179,19 +183,19 @@ def sym_eval(expr, t): fields = flat_obj_array(u, v) - rhs = wave_operator(discr, c=p.c, w=fields) + rhs = wave_operator(dcoll, c=p.c, w=fields) rhs[0] = rhs[0] + sym_eval(sym_f, t_check) expected_rhs = sym_eval(sym_rhs, t_check) rel_linf_err = ( - discr.norm(rhs - expected_rhs, np.inf) - / discr.norm(expected_rhs, np.inf)) + op.norm(dcoll, rhs - expected_rhs, np.inf) + / op.norm(dcoll, expected_rhs, np.inf)) eoc_rec.add_data_point(1./n, rel_linf_err) if visualize: from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, discr.order) + vis = make_visualizer(dcoll) vis.write_vtk_file("wave_accuracy_{order}_{n}.vtu".format(order=order, n=n), [ ("u", fields[0]), @@ -228,16 +232,16 @@ def test_wave_stability(actx_factory, problem, timestep_scale, order, mesh = p.mesh_factory(8) - from grudge.eager import EagerDGDiscretization - discr = EagerDGDiscretization(actx, mesh, order=order) + from grudge.discretization import DiscretizationCollection + dcoll = DiscretizationCollection(actx, mesh, order=order) - nodes = thaw(actx, discr.nodes()) + nodes = thaw(dcoll.nodes(), actx) def sym_eval(expr, t): return sym.EvaluationMapper({"c": p.c, "x": nodes, "t": t})(expr) def get_rhs(t, w): - result = wave_operator(discr, c=p.c, w=w) + result = wave_operator(dcoll, c=p.c, w=w) result[0] += sym_eval(sym_f, t) return result @@ -260,7 +264,7 @@ def get_rhs(t, w): if visualize: from grudge.shortcuts import make_visualizer - vis = make_visualizer(discr, discr.order) + vis = make_visualizer(dcoll) vis.write_vtk_file("wave_stability.vtu", [ ("u", fields[0]), @@ -269,8 +273,8 @@ def get_rhs(t, w): ("v_expected", expected_fields[1:]), ]) - err = discr.norm(fields-expected_fields, np.inf) - max_err = discr.norm(expected_fields, np.inf) + err = op.norm(dcoll, fields-expected_fields, np.inf) + max_err = op.norm(dcoll, expected_fields, np.inf) assert err < max_err