add dloga in grid (#298)

adds dloga in grid, mainly used for primitive geometric terms later

Author: zhichen3

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
@@ -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())) *