Merge branch 'main' into spherical_compressible_prim_geom

examples/mesh/bc_demo.py

@@ -1,9 +1,10 @@
 # test the boundary fill routine
 
 
+import numpy as np
+
 import pyro.mesh.boundary as bnd
 import pyro.mesh.patch as patch
-import numpy as np
 
 
 def doit():

pyro/advection_fv4/fluxes.py

@@ -1,5 +1,5 @@
 import pyro.mesh.array_indexer as ai
-from pyro.advection_fv4 import interface
+from pyro.mesh import fourth_order
 
 
 def fluxes(my_data, rp):
@@ -72,13 +72,13 @@ def fluxes(my_data, rp):
                              1./12.*(a.jp(-2, buf=1) + a.jp(1, buf=1))
 
     else:
-        a_l, a_r = interface.states(a, myg.ng, 1)
+        a_l, a_r = fourth_order.states(a, myg.ng, 1)
         if u > 0:
             a_x = ai.ArrayIndexer(d=a_l, grid=myg)
         else:
             a_x = ai.ArrayIndexer(d=a_r, grid=myg)
 
-        a_l, a_r = interface.states(a, myg.ng, 2)
+        a_l, a_r = fourth_order.states(a, myg.ng, 2)
         if v > 0:
             a_y = ai.ArrayIndexer(d=a_l, grid=myg)
         else:

pyro/compressible/unsplit_fluxes.py

@@ -287,9 +287,6 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
     tm_source.begin()
 
     myg = my_data.grid
-
-    pres = my_data.get_var("pressure")
-
     dens_src = my_aux.get_var("dens_src")
     xmom_src = my_aux.get_var("xmom_src")
     ymom_src = my_aux.get_var("ymom_src")
@@ -303,11 +300,6 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
     ymom_src[:, :] = U_src[:, :, ivars.iymom]
     E_src[:, :] = U_src[:, :, ivars.iener]
 
-    # apply non-conservative pressure gradient
-    if myg.coord_type == 1:
-        xmom_src.v()[:, :] += -(pres.ip(1) - pres.v()) / myg.Lx.v()
-        ymom_src.v()[:, :] += -(pres.jp(1) - pres.v()) / myg.Ly.v()
-
     my_aux.fill_BC("dens_src")
     my_aux.fill_BC("xmom_src")
     my_aux.fill_BC("ymom_src")
@@ -412,25 +404,49 @@ def apply_transverse_flux(U_xl, U_xr, U_yl, U_yr,
         with source terms added.
     """
 
+    myg = my_data.grid
+
     # Use Riemann Solver to get interface flux using the left and right states
 
-    F_x = riemann.riemann_flux(1, U_xl, U_xr,
-                               my_data, rp, ivars,
-                               solid.xl, solid.xr, tc)
+    if myg.coord_type == 1:
+        # We need pressure from interface state for transverse update for
+        # SphericalPolar geometry. So we need interface conserved states.
+        F_x, U_x = riemann.riemann_flux(1, U_xl, U_xr,
+                                        my_data, rp, ivars,
+                                        solid.xl, solid.xr, tc,
+                                        return_cons=True)
 
-    F_y = riemann.riemann_flux(2, U_yl, U_yr,
-                               my_data, rp, ivars,
-                               solid.yl, solid.yr, tc)
+        F_y, U_y = riemann.riemann_flux(2, U_yl, U_yr,
+                                        my_data, rp, ivars,
+                                        solid.yl, solid.yr, tc,
+                                        return_cons=True)
 
-    # Now we update the conserved states using the transverse fluxes.
+        gamma = rp.get_param("eos.gamma")
 
-    myg = my_data.grid
+        # Find primitive variable since we need pressure in transverse update.
+        qx = comp.cons_to_prim(U_x, gamma, ivars, myg)
+        qy = comp.cons_to_prim(U_y, gamma, ivars, myg)
+
+    else:
+        # Directly calculate the interface flux using Riemann Solver
+        F_x = riemann.riemann_flux(1, U_xl, U_xr,
+                                   my_data, rp, ivars,
+                                   solid.xl, solid.xr, tc,
+                                   return_cons=False)
+
+        F_y = riemann.riemann_flux(2, U_yl, U_yr,
+                                   my_data, rp, ivars,
+                                   solid.yl, solid.yr, tc,
+                                   return_cons=False)
+
+    # Now we update the conserved states using the transverse fluxes.
 
     tm_transverse = tc.timer("transverse flux addition")
     tm_transverse.begin()
 
     b = (2, 1)
-    hdtV = 0.5*dt / myg.V
+    hdt = 0.5*dt
+    hdtV = hdt / myg.V
 
     for n in range(ivars.nvar):
 
@@ -454,6 +470,25 @@ def apply_transverse_flux(U_xl, U_xr, U_yl, U_yr,
             - hdtV.v(buf=b)*(F_x.ip(1, buf=b, n=n)*myg.Ax.ip(1, buf=b) -
                              F_x.v(buf=b, n=n)*myg.Ax.v(buf=b))
 
+    # apply non-conservative pressure gradient for momentum in spherical geometry
+    # Note that we should only apply this pressure gradient
+    # to the momentum corresponding to the transverse direction.
+    # The momentum in the normal direction already updated pressure during reconstruction.
+
+    if myg.coord_type == 1:
+
+        U_xl.v(buf=b, n=ivars.iymom)[:, :] += - hdt * (qy.ip_jp(-1, 1, buf=b, n=ivars.ip) -
+                                                       qy.ip(-1, buf=b, n=ivars.ip)) / myg.Ly.v(buf=b)
+
+        U_xr.v(buf=b, n=ivars.iymom)[:, :] += - hdt * (qy.jp(1, buf=b, n=ivars.ip) -
+                                                       qy.v(buf=b, n=ivars.ip)) / myg.Ly.v(buf=b)
+
+        U_yl.v(buf=b, n=ivars.ixmom)[:, :] += - hdt * (qx.ip_jp(1, -1, buf=b, n=ivars.ip) -
+                                                       qx.jp(-1, buf=b, n=ivars.ip)) / myg.Lx.v(buf=b)
+
+        U_yr.v(buf=b, n=ivars.ixmom)[:, :] += - hdt * (qx.ip(1, buf=b, n=ivars.ip) -
+                                                       qx.v(buf=b, n=ivars.ip)) / myg.Lx.v(buf=b)
+
     tm_transverse.end()
 
     return U_xl, U_xr, U_yl, U_yr

pyro/compressible_fv4/fluxes.py

@@ -4,9 +4,8 @@
 
 import pyro.compressible as comp
 import pyro.mesh.array_indexer as ai
-from pyro.advection_fv4 import interface
 from pyro.compressible import riemann
-from pyro.mesh import reconstruction
+from pyro.mesh import fourth_order, reconstruction
 
 
 def flux_cons(ivars, idir, gamma, q):
@@ -104,7 +103,7 @@ def fluxes(myd, rp, ivars):
 
         else:
             for n in range(ivars.nq):
-                q_l[:, :, n], q_r[:, :, n] = interface.states(q_avg[:, :, n], myg.ng, idir)
+                q_l[:, :, n], q_r[:, :, n] = fourth_order.states(q_avg[:, :, n], myg.ng, idir)
 
             # apply flattening
             for n in range(ivars.nq):

pyro/mesh/__init__.py

@@ -3,7 +3,7 @@
 necessary to work with finite-volume data.
 """
 
-__all__ = ['patch', 'integration', 'reconstruction']
+__all__ = ['patch', 'fourth_order', 'integration', 'reconstruction']
 
 from .array_indexer import ArrayIndexer, ArrayIndexerFC
 from .boundary import BC, bc_is_solid, define_bc

pyro/advection_fv4/interface.py → pyro/mesh/fourth_order.py

@@ -1,3 +1,5 @@
+"""Reconstruction routines for 4th order finite-volume methods"""
+
 import numpy as np
 from numba import njit
 
@@ -84,7 +86,7 @@ def states(a, ng, idir):
                 # this lives on the interface
                 d3a[i, j] = d2ac[i, j] - d2ac[i - 1, j]
 
-        # this is a look over cell centers, affecting
+        # this is a loop over cell centers, affecting
         # i-1/2,R and i+1/2,L
         for i in range(ilo - 1, ihi + 1):
             for j in range(jlo - 1, jhi + 1):