add left and right area

python-hydro · Jul 20, 2024 · b93c187 · b93c187
1 parent 1de6124
commit b93c187
Show file tree

Hide file tree

Showing 2 changed files with 38 additions and 17 deletions.
diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
@@ -196,19 +196,18 @@ def __init__(self, nx, ny, ng=1,
         super().__init__(nx, ny, ng, xmin, xmax, ymin, ymax)
 
         # This is length of the side that is perpendicular to x.
-        # self.area_x = self.dy
 
         area_x = np.full((self.qx, self.qy), self.dy)
-        self.Ax = ArrayIndexer(area_x, grid=self)
+        self.Ax_l = ArrayIndexer(area_x, grid=self)
+        self.Ax_r = ArrayIndexer(area_x, grid=self)
 
         # This is length of the side that is perpendicular to y.
-        # self.area_y = self.dx
 
         area_y = np.full((self.qx, self.qy), self.dx)
-        self.Ay = ArrayIndexer(area_y, grid=self)
+        self.Ay_l = ArrayIndexer(area_y, grid=self)
+        self.Ay_r = ArrayIndexer(area_y, grid=self)
 
         # Volume (Area) of the cell.
-        # self.vol = self.dx * self.dy
 
         volume = np.full((self.qx, self.qy), self.dx * self.dy)
         self.V = ArrayIndexer(volume, grid=self)
@@ -244,19 +243,29 @@ def __init__(self, nx, ny, ng=1,
 
         # Returns an array of the face area that points in the r(x) direction.
         # dL_theta x dL_phi = r^2 * sin(theta) * dtheta * dphi
-        # dA_r = - r{i-1/2}^2 * 2pi * cos(theta{i+1/2}) - cos(theta{i-1/2})
 
-        area_x = np.abs(-2.0 * np.pi * self.xl2d**2 *
+        # dAr_l = - r{i-1/2}^2 * 2pi * cos(theta{i+1/2}) - cos(theta{i-1/2})
+        area_x_l = np.abs(-2.0 * np.pi * self.xl2d**2 *
                  (np.cos(self.yr2d) - np.cos(self.yl2d)))
-        self.Ax = ArrayIndexer(area_x, grid=self)
+        self.Ax_l = ArrayIndexer(area_x_l, grid=self)
+
+        # dAr_r = - r{i+1/2}^2 * 2pi * cos(theta{i+1/2}) - cos(theta{i-1/2})
+        area_x_r = np.abs(-2.0 * np.pi * self.xr2d**2 *
+                 (np.cos(self.yr2d) - np.cos(self.yl2d)))
+        self.Ax_r = ArrayIndexer(area_x_r, grid=self)
 
         # Returns an array of the face area that points in the theta(y) direction.
         # dL_phi x dL_r = dr * r * sin(theta) * dphi
-        # dA_theta = 0.5 * pi * sin(theta{i-1/2}) * (r{i+1/2}^2 - r{i-1/2}^2)
 
-        area_y = np.abs(0.5 * np.pi * np.sin(self.yl2d) *
+        # dAtheta_l = 0.5 * pi * sin(theta{i-1/2}) * (r{i+1/2}^2 - r{i-1/2}^2)
+        area_y_l = np.abs(0.5 * np.pi * np.sin(self.yl2d) *
+                 (self.xr2d**2 - self.xl2d**2))
+        self.Ay_l = ArrayIndexer(area_y_l, grid=self)
+
+        # dAtheta_r = 0.5 * pi * sin(theta{i+1/2}) * (r{i+1/2}^2 - r{i-1/2}^2)
+        area_y_r = np.abs(0.5 * np.pi * np.sin(self.yr2d) *
                  (self.xr2d**2 - self.xl2d**2))
-        self.Ay = ArrayIndexer(area_y, grid=self)
+        self.Ay_r = ArrayIndexer(area_y_r, grid=self)
 
         # Returns an array of the volume of each cell.
         # dV = dL_r * dL_theta * dL_phi

diff --git a/pyro/mesh/tests/test_patch.py b/pyro/mesh/tests/test_patch.py
@@ -87,10 +87,12 @@ def teardown_method(self):
         self.g = None
 
     def test_Ax(self):
-        assert np.all(self.g.Ax.v() == 0.1)
+        assert np.all(self.g.Ax_l.v() == 0.1)
+        assert np.all(self.g.Ax_r.v() == 0.1)
 
     def test_Ay(self):
-        assert np.all(self.g.Ay.v() == 0.25)
+        assert np.all(self.g.Ay_l.v() == 0.25)
+        assert np.all(self.g.Ay_r.v() == 0.25)
 
     def test_V(self):
         assert np.all(self.g.V.v() == 0.1 * 0.25)
@@ -121,23 +123,33 @@ def test_Ax(self):
         jlo = self.g.jlo
         jhi = self.g.jhi
 
-        area_x = np.abs(-2.0 * np.pi * self.g.xl2d[ilo:ihi+1, jlo:jhi+1]**2 *
+        area_x_l = np.abs(-2.0 * np.pi * self.g.xl2d[ilo:ihi+1, jlo:jhi+1]**2 *
                  (np.cos(self.g.yr2d[ilo:ihi+1, jlo:jhi+1]) -
                   np.cos(self.g.yl2d[ilo:ihi+1, jlo:jhi+1])))
+        assert_array_equal(self.g.Ax_l.v(), area_x_l)
 
-        assert_array_equal(self.g.Ax.v(), area_x)
+        area_x_r = np.abs(-2.0 * np.pi * self.g.xr2d[ilo:ihi+1, jlo:jhi+1]**2 *
+                 (np.cos(self.g.yr2d[ilo:ihi+1, jlo:jhi+1]) -
+                  np.cos(self.g.yl2d[ilo:ihi+1, jlo:jhi+1])))
+        assert_array_equal(self.g.Ax_r.v(), area_x_r)
 
     def test_Ay(self):
         ilo = self.g.ilo
         ihi = self.g.ihi
         jlo = self.g.jlo
         jhi = self.g.jhi
 
-        area_y = np.abs(0.5 * np.pi *
+        area_y_l = np.abs(0.5 * np.pi *
                  np.sin(self.g.yl2d[ilo:ihi+1, jlo:jhi+1]) *
                  (self.g.xr2d[ilo:ihi+1, jlo:jhi+1]**2 -
                   self.g.xl2d[ilo:ihi+1, jlo:jhi+1]**2))
-        assert_array_equal(self.g.Ay.v(), area_y)
+        assert_array_equal(self.g.Ay_l.v(), area_y_l)
+
+        area_y_r = np.abs(0.5 * np.pi *
+                 np.sin(self.g.yr2d[ilo:ihi+1, jlo:jhi+1]) *
+                 (self.g.xr2d[ilo:ihi+1, jlo:jhi+1]**2 -
+                  self.g.xl2d[ilo:ihi+1, jlo:jhi+1]**2))
+        assert_array_equal(self.g.Ay_r.v(), area_y_r)
 
     def test_V(self):
         ilo = self.g.ilo