Merge branch 'main' into fix_transverse_pres_gradient

pyro/mesh/patch.py

@@ -220,6 +220,13 @@ def __init__(self, nx, ny, *, ng=1,
 
         self.Ay = self.Lx
 
+        # Spatial derivative of log(area). It's zero for Cartesian
+
+        self.dlogAx = ArrayIndexer(np.zeros_like(self.Ax),
+                                   grid=self)
+        self.dlogAy = ArrayIndexer(np.zeros_like(self.Ay),
+                                   grid=self)
+
         # Volume of the cell.
 
         self.V = ArrayIndexer(np.full((self.qx, self.qy), self.dx * self.dy),
@@ -235,7 +242,9 @@ def __str__(self):
 class SphericalPolar(Grid2d):
     """
     This class defines a spherical polar grid.
-    This is technically a 2D geometry but assumes azimuthal symmetry.
+    This is technically a 3D geometry but assumes azimuthal symmetry,
+    and zero velocity in phi-direction.
+    Hence 2D.
 
     Define:
     r = x
@@ -245,15 +254,15 @@ class SphericalPolar(Grid2d):
     def __init__(self, nx, ny, *, ng=1,
                  xmin=0.2, xmax=1.0, ymin=0.0, ymax=1.0):
 
+        super().__init__(nx, ny, ng=ng, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
+
         # Make sure theta is within [0, PI]
         assert ymin >= 0.0 and ymax <= np.pi, "y or \u03b8 should be within [0, \u03c0]."
 
         # Make sure the ghost cells doesn't extend out negative x(r)
-        assert xmin - ng*(xmax-xmin)/nx >= 0.0, \
+        assert xmin - ng*self.dx >= 0.0, \
             "xmin (r-direction), must be large enough so ghost cell doesn't have negative x."
 
-        super().__init__(nx, ny, ng=ng, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
-
         self.coord_type = 1
 
         # Length of the side along r-direction, dr
@@ -279,6 +288,12 @@ def __init__(self, nx, ny, *, ng=1,
         self.Ay = np.abs(np.pi * np.sin(self.yl2d) *
                          (self.xr2d**2 - self.xl2d**2))
 
+        # dlogAx = 1 / r^2 d( r^2 ) / dr = 2 / r
+        self.dlogAx = 2.0 / self.x2d
+
+        # dlogAy = 1 / (r sin(theta)) d( sin(theta) )/dtheta = cot(theta) / r
+        self.dlogAy = 1.0 / (np.tan(self.y2d) * self.x2d)
+
         # Returns an array of the volume of each cell.
         # dV = dL_r * dL_theta * dL_phi
         #    = (dr) * (r * dtheta) * (r * sin(theta) * dphi)

pyro/mesh/tests/test_patch.py

@@ -92,6 +92,12 @@ def test_Ax(self):
     def test_Ay(self):
         assert np.all(self.g.Ay.v() == 0.25)
 
+    def test_dlogAx(self):
+        assert np.all(self.g.dlogAx.v() == 0.0)
+
+    def test_dlogAy(self):
+        assert np.all(self.g.dlogAy.v() == 0.0)
+
     def test_V(self):
         assert np.all(self.g.V.v() == 0.1 * 0.25)
 
@@ -135,6 +141,21 @@ def test_Ay(self):
                           (self.g.xr2d.v()**2 - self.g.xl2d.v()**2))
         assert_array_equal(self.g.Ay.jp(1), area_y_r)
 
+    def test_dlogAx(self):
+        i = 4
+        r = self.g.xmin + (i + 0.5 - self.g.ng) * self.g.dx
+        dlogAx = 2.0 / r
+        assert (self.g.dlogAx[i, :] == dlogAx).all()
+
+    def test_dlogAy(self):
+        i = 4
+        j = 6
+        r = self.g.xmin + (i + 0.5 - self.g.ng) * self.g.dx
+        tan = np.tan(self.g.ymin + (j + 0.5 - self.g.ng) * self.g.dy)
+        dlogAy = 1.0 / (r * tan)
+
+        assert self.g.dlogAy[i, j] == dlogAy
+
     def test_V(self):
         volume = np.abs(-2.0 * np.pi / 3.0 *
                  (np.cos(self.g.yr2d.v()) - np.cos(self.g.yl2d.v())) *

pyro/simulation_null.py

@@ -1,4 +1,5 @@
 import h5py
+import numpy as np
 
 import pyro.mesh.boundary as bnd
 import pyro.util.profile_pyro as profile
@@ -52,6 +53,19 @@ def grid_setup(rp, ng=1):
                           ymin=ymin, ymax=ymax,
                           ng=ng)
 
+    # Make sure y-boundary condition is reflect when
+    # ymin = 0 and ymax = pi
+
+    if grid_type == "SphericalPolar":
+        if ymin <= 0.05:
+            rp.set_param("mesh.ylboundary", "reflect")
+            msg.warning("With SphericalPolar grid," +
+                        " mesh.ylboundary auto set to reflect when ymin ~ 0")
+        if abs(np.pi - ymax) <= 0.05:
+            rp.set_param("mesh.yrboundary", "reflect")
+            msg.warning("With SphericalPolar grid," +
+                        " mesh.yrboundary auto set to reflect when ymax ~ pi")
+
     return my_grid