Merge branch 'main' into incompressible_simplify

.github/workflows/regtest.yml

@@ -44,5 +44,8 @@ jobs:
         run: pip install .
 
       - name: Run tests via test.py
-        run: ./pyro/test.py --nproc 0
+        run: |
+          # Numpy 1.24 gets slightly diff answers: https://github.com/numpy/numpy/issues/23289
+          export NPY_DISABLE_CPU_FEATURES="AVX512_SKX"
+          ./pyro/test.py --nproc 0
 

pyproject.toml

@@ -46,3 +46,6 @@ max-line-length = 132
 [tool.codespell]
 skip = ".git,*docs/build"
 ignore-words-list = "pres"
+
+[tool.isort]
+known_first_party = ["pyro"]

pyro/advection/advective_fluxes.py

@@ -1,7 +1,5 @@
-from pyro.mesh import reconstruction
+def unsplit_fluxes(my_data, rp, dt, scalar_name, interface):
 
-
-def unsplit_fluxes(my_data, rp, dt, scalar_name):
     r"""
     Construct the fluxes through the interfaces for the linear advection
     equation:
@@ -53,45 +51,11 @@ def unsplit_fluxes(my_data, rp, dt, scalar_name):
     """
 
     myg = my_data.grid
-
     a = my_data.get_var(scalar_name)
 
-    # get the advection velocities
-    u = rp.get_param("advection.u")
-    v = rp.get_param("advection.v")
-
-    cx = u*dt/myg.dx
-    cy = v*dt/myg.dy
-
-    # --------------------------------------------------------------------------
-    # monotonized central differences
-    # --------------------------------------------------------------------------
-
-    limiter = rp.get_param("advection.limiter")
-
-    ldelta_ax = reconstruction.limit(a, myg, 1, limiter)
-    ldelta_ay = reconstruction.limit(a, myg, 2, limiter)
+    #Compute the velocities on each direction and the interface states a_{i+1/2}^{n+1/2}.
 
-    a_x = myg.scratch_array()
-
-    # upwind
-    if u < 0:
-        # a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
-        a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cx)*ldelta_ax.v(buf=1)
-    else:
-        # a_x[i,j] = a[i-1,j] + 0.5*(1.0 - cx)*ldelta_a[i-1,j]
-        a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*(1.0 - cx)*ldelta_ax.ip(-1, buf=1)
-
-    # y-direction
-    a_y = myg.scratch_array()
-
-    # upwind
-    if v < 0:
-        # a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
-        a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cy)*ldelta_ay.v(buf=1)
-    else:
-        # a_y[i,j] = a[i,j-1] + 0.5*(1.0 - cy)*ldelta_a[i,j-1]
-        a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*(1.0 - cy)*ldelta_ay.jp(-1, buf=1)
+    u, v, a_x, a_y = interface(a, myg, rp, dt)
 
     # compute the transverse flux differences.  The flux is just (u a)
     # HOTF

pyro/advection/interface.py

@@ -0,0 +1,43 @@
+from pyro.mesh import reconstruction
+
+
+def linear_interface(a, myg, rp, dt):
+
+    # get the advection velocities
+    u = rp.get_param("advection.u")
+    v = rp.get_param("advection.v")
+
+    cx = u*dt/myg.dx
+    cy = v*dt/myg.dy
+
+    # --------------------------------------------------------------------------
+    # monotonized central differences
+    # --------------------------------------------------------------------------
+
+    limiter = rp.get_param("advection.limiter")
+
+    ldelta_ax = reconstruction.limit(a, myg, 1, limiter)
+    ldelta_ay = reconstruction.limit(a, myg, 2, limiter)
+
+    a_x = myg.scratch_array()
+
+    # upwind
+    if u < 0:
+        # a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
+        a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cx)*ldelta_ax.v(buf=1)
+    else:
+        # a_x[i,j] = a[i-1,j] + 0.5*(1.0 - cx)*ldelta_a[i-1,j]
+        a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*(1.0 - cx)*ldelta_ax.ip(-1, buf=1)
+
+    # y-direction
+    a_y = myg.scratch_array()
+
+    # upwind
+    if v < 0:
+        # a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
+        a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cy)*ldelta_ay.v(buf=1)
+    else:
+        # a_y[i,j] = a[i,j-1] + 0.5*(1.0 - cy)*ldelta_a[i,j-1]
+        a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*(1.0 - cy)*ldelta_ay.jp(-1, buf=1)
+
+    return u, v, a_x, a_y

pyro/advection/simulation.py

@@ -5,6 +5,7 @@
 import numpy as np
 
 import pyro.advection.advective_fluxes as flx
+from pyro.advection.interface import linear_interface
 from pyro.mesh import patch
 from pyro.particles import particles
 from pyro.simulation_null import NullSimulation, bc_setup, grid_setup
@@ -66,7 +67,7 @@ def evolve(self):
         dtdx = self.dt/self.cc_data.grid.dx
         dtdy = self.dt/self.cc_data.grid.dy
 
-        flux_x, flux_y = flx.unsplit_fluxes(self.cc_data, self.rp, self.dt, "density")
+        flux_x, flux_y = flx.unsplit_fluxes(self.cc_data, self.rp, self.dt, "density", linear_interface)
 
         """
         do the differencing for the fluxes now.  Here, we use slices so we