Add numpy option (#39)

* Add numpy option

* Spell profile correctly

* Temporarily use future mirgecom@production.

* Revert back to mirgecom@production

* Clean up messaging and arg parsing a bit.

* fix numpy usage, add numpy test, test accuracy

* libegl1

* add run_numpy.sh, set casenames in run*.sh

---------

Co-authored-by: Matthias Diener <mdiener@illinois.edu>
Co-authored-by: Matthias Diener <matthias.diener@gmail.com>

.github/workflows/ci.yml

@@ -39,7 +39,7 @@ jobs:
         - uses: actions/checkout@v3
         - name: Install emirge
           run: |
-            [[ $(uname) == Linux ]] && sudo apt-get update && sudo apt-get install -y openmpi-bin libopenmpi-dev
+            [[ $(uname) == Linux ]] && sudo apt-get update && sudo apt-get install -y openmpi-bin libopenmpi-dev bsdmainutils libgl1-mesa-glx libegl1
             [[ $(uname) == Darwin ]] && brew update && brew install mpich
             ./buildMirge.sh
             cd emirge
@@ -66,12 +66,34 @@ jobs:
             cd smoke_test
             ./run_restart.sh
 
+        - name: Numpy smoke test
+          run: |
+            source emirge/config/activate_env.sh
+            cd smoke_test
+            ./run_numpy.sh
+
         - name: Lazy smoke test
           run: |
             source emirge/config/activate_env.sh
             cd smoke_test
             ./run_lazy.sh
 
+        - name: Compare smoke_test accuracy
+          run: |
+            source emirge/config/activate_env.sh
+            cd smoke_test/viz_data
+
+            TOL_LAZY=1e-2
+            TOL_NUMPY=1e-2
+
+            set -x
+            for vizfile in $(ls *-eager-*.vtu); do
+                lazy_vizfile=$(echo ${vizfile/eager/lazy})
+                python ../../emirge/mirgecom/bin/mirgecompare.py --tolerance $TOL_LAZY ${vizfile} ${lazy_vizfile}
+                numpy_vizfile=$(echo ${vizfile/eager/numpy})
+                python ../../emirge/mirgecom/bin/mirgecompare.py --tolerance $TOL_NUMPY ${vizfile} ${numpy_vizfile}
+            done
+
         - name: Lazy smoke test without wall
           run: |
             source emirge/config/activate_env.sh

driver.py

@@ -28,6 +28,8 @@
                         help="enable lazy evaluation [OFF]")
     parser.add_argument("--overintegration", action="store_true",
         help="use overintegration in the RHS computations")
+    parser.add_argument("--numpy", action="store_true",
+        help="use numpy-based eager actx.")
 
     args = parser.parse_args()
 
@@ -39,21 +41,13 @@
     else:
         print(f"Default casename {casename}")
 
-    from warnings import warn
-    if args.esdg:
-        if not args.lazy:
-            warn("ESDG requires lazy-evaluation, enabling --lazy.")
-        if not args.overintegration:
-            warn("ESDG requires overintegration, enabling --overintegration.")
-
-    lazy = args.lazy or args.esdg
-    if args.profile:
-        if lazy:
-            raise ValueError("Can't use lazy and profiling together.")
+    from mirgecom.simutil import ApplicationOptionsError
+    if args.esdg and not (args.lazy or args.numpy):
+        raise ApplicationOptionsError("ESDG requires lazy or numpy context.")
 
     from mirgecom.array_context import get_reasonable_array_context_class
     actx_class = get_reasonable_array_context_class(
-        lazy=lazy, distributed=True, profiling=args.profile)
+        lazy=args.lazy, distributed=True, profiling=args.profile, numpy=args.numpy)
 
     restart_filename = None
     if args.restart_file:

smoke_test/run.sh

@@ -1,3 +1,3 @@
 #!/bin/bash
-mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log
+mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --casename=prediction-eager
 #mpirun -n 1 python -u -O -m mpi4py driver.py -i run_params.yaml --log

smoke_test/run_lazy.sh

@@ -1,2 +1,2 @@
 #!/bin/bash
-mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --lazy
+mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --lazy --casename=prediction-lazy

smoke_test/run_numpy.sh

@@ -0,0 +1,2 @@
+#!/bin/bash
+mpirun -n 2 python -u -O -m mpi4py driver.py -i run_params.yaml --log --numpy --casename=prediction-numpy

smoke_test/run_restart.sh

@@ -1,2 +1,2 @@
 #!/bin/bash
-mpirun -n 2 python -u -O -m mpi4py driver.py -i restart_params.yaml -r restart_data/prediction-000000010 -t restart_data/prediction-000000000 --log
+mpirun -n 2 python -u -O -m mpi4py driver.py -i restart_params.yaml -r restart_data/prediction-eager-000000010 -t restart_data/prediction-eager-000000000 --log

y3prediction/prediction.py

@@ -1355,8 +1355,12 @@ def my_partitioner(mesh, tag_to_elements, num_ranks):
     else:
         logger.warning("No target file specied, using restart as target")
 
+    disc_msg = f"Making {dim}D order {order} discretization"
+    if use_overintegration:
+        disc_msg = disc_msg + f" with quadrature order {quadrature_order}"
+    disc_msg = disc_msg + "."
     if rank == 0:
-        logger.info("Making discretization")
+        logger.info(disc_msg)
 
     dcoll = create_discretization_collection(
         actx,
@@ -1413,7 +1417,7 @@ def my_partitioner(mesh, tag_to_elements, num_ranks):
         char_length_wall = force_evaluation(actx,
             characteristic_lengthscales(actx, dcoll, dd=dd_vol_wall))
         xpos_wall = wall_nodes[0]
-        char_length_wall = char_length_wall + actx.zeros_like(xpos_wall)
+        char_length_wall = char_length_wall + actx.np.zeros_like(xpos_wall)
         """
         smoothness_diffusivity_wall = \
             smooth_char_length_alpha*char_length_wall**2/current_dt
@@ -1810,8 +1814,8 @@ def update_smoothness(state, time):
 
             cv = cv + 0.*fluid_rhs
 
-            wall_mass_rhs = actx.zeros_like(wv.mass)
-            wall_ox_mass_rhs = actx.zeros_like(wv.mass)
+            wall_mass_rhs = actx.np.zeros_like(wv.mass)
+            wall_ox_mass_rhs = actx.np.zeros_like(wv.mass)
             wall_rhs = wall_time_scale * WallVars(
                 mass=wall_mass_rhs,
                 energy=wall_energy_rhs,
@@ -1978,7 +1982,7 @@ def get_production_rates(cv, temperature):
             restart_wv = WallVars(
                 mass=wall_mass,
                 energy=wall_mass * wall_cp * temp_wall,
-                ox_mass=actx.zeros_like(wall_mass))
+                ox_mass=actx.np.zeros_like(wall_mass))
 
     if use_wall:
         restart_wv = force_evaluation(actx, restart_wv)
@@ -2381,7 +2385,7 @@ def _sponge_sigma(sponge_field, x_vec):
 
     get_sponge_sigma = actx.compile(_sponge_sigma)
 
-    sponge_sigma = actx.zeros_like(restart_cv.mass)
+    sponge_sigma = actx.np.zeros_like(restart_cv.mass)
     sponge_sigma = get_sponge_sigma(sponge_sigma, fluid_nodes)
 
     def _sponge_source(cv):
@@ -2437,12 +2441,18 @@ def _sponge_source(cv):
 
             logmgr.add_watches([
                 ("memory_usage_hwm.max",
-                 "| \t memory hwm: {value:7g} Mb\n"),
-                ("memory_usage_mempool_managed.max",
-                 "| \t mempool total: {value:7g} Mb\n"),
-                ("memory_usage_mempool_active.max",
-                 "| \t mempool active: {value:7g} Mb")
-            ])
+                 "| \t memory hwm: {value:7g} Mb\n")])
+
+            from mirgecom.array_context import actx_class_is_numpy
+
+            if not actx_class_is_numpy(actx_class):
+                # numpy has no CL mempool
+                logmgr.add_watches([
+                    ("memory_usage_mempool_managed.max",
+                    "| \t mempool total: {value:7g} Mb\n"),
+                    ("memory_usage_mempool_active.max",
+                    "| \t mempool active: {value:7g} Mb")
+                ])
 
         if use_profiling:
             logmgr.add_watches(["pyopencl_array_time.max"])
@@ -2610,8 +2620,8 @@ def compute_viz_fields_coupled(fluid_state, wv, wdv, time):
 
         cv = cv + 0.*fluid_rhs
 
-        wall_mass_rhs = actx.zeros_like(wv.mass)
-        wall_ox_mass_rhs = actx.zeros_like(wv.mass)
+        wall_mass_rhs = actx.np.zeros_like(wv.mass)
+        wall_ox_mass_rhs = actx.np.zeros_like(wv.mass)
         wall_rhs = wall_time_scale * WallVars(
             mass=wall_mass_rhs,
             energy=wall_energy_rhs,
@@ -3533,17 +3543,17 @@ def unfiltered_rhs(t, state):
 
         if use_wall:
             # wall mass loss
-            wall_mass_rhs = actx.zeros_like(wv.mass)
+            wall_mass_rhs = actx.np.zeros_like(wv.mass)
             if use_wall_mass:
                 wall_mass_rhs = -wall_model.mass_loss_rate(
                     mass=wv.mass, ox_mass=wv.ox_mass,
                     temperature=wdv.temperature)
 
             # wall oxygen diffusion
-            wall_ox_mass_rhs = actx.zeros_like(wv.mass)
+            wall_ox_mass_rhs = actx.np.zeros_like(wv.mass)
             if use_wall_ox:
                 if nspecies == 0:
-                    fluid_ox_mass = actx.zeros_like(cv.mass)
+                    fluid_ox_mass = actx.np.zeros_like(cv.mass)
                 elif nspecies > 3:
                     fluid_ox_mass = cv.species_mass[i_ox]
                 else:

y3prediction/utils.py

@@ -336,7 +336,7 @@ def __call__(self, sponge_field, x_vec, *, time=0.0):
         xpos = x_vec[0]
         ypos = x_vec[1]
         actx = xpos.array_context
-        zeros = actx.zeros_like(xpos)
+        zeros = actx.np.zeros_like(xpos)
         x0 = zeros + self._x0
 
         if abs(self._direction) == 1: