Modernize Mirgecom

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:]), ])
 

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 = \

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)

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 = \