From e0b136f6cfef263aa15e521c681b6c5487a3466b Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Thu, 9 Mar 2023 14:13:20 -0500
Subject: [PATCH 01/13] Included PolarGrid class. right now it assumes no ghost
 cells

---
 pyro/mesh/patch.py | 179 +++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 174 insertions(+), 5 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 52b3e85c1..42e33b4ff 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -34,9 +34,9 @@
 import h5py
 import numpy as np
 
-import pyro.mesh.boundary as bnd
-from pyro.mesh.array_indexer import ArrayIndexer, ArrayIndexerFC
-from pyro.util import msg
+# import pyro.mesh.boundary as bnd
+# from pyro.mesh.array_indexer import ArrayIndexer, ArrayIndexerFC
+# from pyro.util import msg
 
 
 class Grid2d:
@@ -885,5 +885,174 @@ def do_demo():
     mydata.pretty_print("a")
 
 
-if __name__ == "__main__":
-    do_demo()
+# **c checking
+class PolarGrid(Grid2d):
+    """
+    the 2-d grid class.  The grid object will contain the coordinate
+    information (at various centerings).
+
+    A basic (1-d) representation of the layout is::
+
+       |     |      |     X     |     |      |     |     X     |      |     |
+       +--*--+- // -+--*--X--*--+--*--+- // -+--*--+--*--X--*--+- // -+--*--+
+          0          ng-1    ng   ng+1         ... ng+nx-1 ng+nx      2ng+nx-1
+
+                            ilo                      ihi
+
+       |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->|
+
+    The '*' marks the data locations.
+    """
+
+    # pylint: disable=too-many-instance-attributes
+
+    def __init__(self, nx, ny, ng=1,
+                 xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0):
+
+        """
+        Create a PolarGrid object.
+
+        The only data that we require is the number of points that
+        make up the mesh in each direction.  Optionally we take the
+        extrema of the domain (default is [0,1]x[0,1]) and number of
+        ghost cells (default is 1).
+
+        Note that the Grid2d object only defines the discretization,
+        it does not know about the boundary conditions, as these can
+        vary depending on the variable.
+
+        Parameters
+        ----------
+        nx : int
+            Number of zones in the r-direction
+        ny : int
+            Number of zones in the theta-direction
+        ng : int, optional
+            Number of ghost cells
+        xmin : float, optional
+            Physical coordinate at the lower x boundary
+        xmax : float, optional
+            Physical coordinate at the upper x boundary
+        ymin : float, optional
+            Physical coordinate at the lower y boundary
+        ymax : float, optional
+            Physical coordinate at the upper y boundary
+        """
+
+        # pylint: disable=too-many-arguments
+
+        # size of grid
+        self.nx = int(nx)
+        self.ny = int(ny)
+        self.ng = int(ng)
+
+        self.qx = int(2*ng + nx)
+        self.qy = int(2*ng + ny)
+
+        # domain extrema
+        self.xmin = xmin
+        self.xmax = xmax
+
+        self.ymin = ymin
+        self.ymax = ymax
+
+        # compute the indices of the block interior (excluding guardcells)
+        self.ilo = self.ng
+        self.ihi = self.ng + self.nx-1
+
+        self.jlo = self.ng
+        self.jhi = self.ng + self.ny-1
+
+        # center of the grid (for convenience)
+        self.ic = self.ilo + self.nx//2 - 1
+        self.jc = self.jlo + self.ny//2 - 1
+
+        # define the coordinate information at the left, center, and right
+        # zone coordinates
+        self.dx = (xmax - xmin)/nx
+
+        self.xl = (np.arange(self.qx) - ng)*self.dx + xmin
+        self.xr = (np.arange(self.qx) + 1.0 - ng)*self.dx + xmin
+        self.x = 0.5*(self.xl + self.xr)
+
+        self.dy = (ymax - ymin)/ny
+
+        self.yl = (np.arange(self.qy) - ng)*self.dy + ymin
+        self.yr = (np.arange(self.qy) + 1.0 - ng)*self.dy + ymin
+        self.y = 0.5*(self.yl + self.yr)
+
+        # 2-d versions of the zone coordinates (replace with meshgrid?)
+        x2d = np.repeat(self.x, self.qy)
+        x2d.shape = (self.qx, self.qy)
+        self.x2d = x2d
+
+        y2d = np.repeat(self.y, self.qx)
+        y2d.shape = (self.qy, self.qx)
+        y2d = np.transpose(y2d)
+        self.y2d = y2d
+
+
+    def face_area(self):
+        """
+        Return an array of the face areas. 
+        The shape of the returned array is (ni, nj).
+        """
+        tr = lambda arr: arr.transpose(1, 2, 0)
+        x = self.cell_vertices()[:,0]
+        y = self.cell_vertices()[0,:]
+        r0 = x[:-1, :-1]
+        r1 = x[+1:, :-1]
+        t0 = y[:-1, :-1]
+        t1 = y[+1:, +1:]
+
+        # ** the area of the arc
+
+        area_i = np.pi * (r0 * np.sin(t0) + r0 * np.sin(t1)) * np.sqrt(np.square(r0 * np.sin(t1) - r0 * np.sin(t0)) + np.square(r0 * np.cos(t1) - r0 * np.cos(t0)))
+        area_j = np.pi * (r0 * np.sin(t0) + r1 * np.sin(t0)) * np.sqrt(np.square(r1 * np.sin(t0) - r0 * np.sin(t0)) + np.square(r1 * np.cos(t0) - r0 * np.cos(t0)))
+        return tr(np.array([area_i, area_j]))
+
+    def cell_volumes(self):
+        """
+        Return an array of the cell volume data for the given coordinate box
+        The shape of the returned array is (ni, nj).
+        """
+
+        x = self.cell_vertices()[:,0]
+        y = self.cell_vertices()[0,:]
+
+        r0 = x[:-1, :-1]
+        r1 = x[+1:, :-1]
+        t0 = y[:-1, :-1]
+        t1 = y[+1:, +1:]
+
+        return (r1 ** 3 - r0 ** 3) * (np.cos(t1) - np.cos(t0)) * (-2.0 * np.pi) / 3.0
+        # return r1
+
+    def cell_vertices(self):
+        """
+        Return the coordinates of cell vertices
+        The arrays range from 0 to 1 for now
+        """
+        # if self.ng == 0:
+        #     xv = np.linspace(0, 1, self.nx + 1)[:-1]
+        #     yv = np.linspace(0, 1, self.ny + 1)[:-1]
+        # else:
+        #     xv = np.linspace(0, 1, self.nx + 1)
+        #     yv = np.linspace(0, 1, self.ny + 1)
+
+        xv = np.linspace(0, 1, self.nx + 1)[:-1]
+        yv = np.linspace(0, 1, self.ny + 1)[:-1]
+
+     
+        tr = lambda arr: arr.transpose(1, 2, 0)
+        x, y = np.meshgrid(xv, yv, indexing="ij")
+        return tr(np.array([x, y]))
+
+   
+
+    
+
+
+
+# if __name__ == "__main__":
+#     do_demo()

From eb71efc415e53e2329961147a056cc27ff16c418 Mon Sep 17 00:00:00 2001
From: Sabina Sagynbayeva <83788641+ssagynbayeva@users.noreply.github.com>
Date: Thu, 9 Mar 2023 14:34:40 -0500
Subject: [PATCH 02/13] Update pyro/mesh/patch.py

Co-authored-by: Eric T. Johnson <yut23@users.noreply.github.com>
---
 pyro/mesh/patch.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 42e33b4ff..4dea2db2e 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -1054,5 +1054,5 @@ def cell_vertices(self):
 
 
 
-# if __name__ == "__main__":
-#     do_demo()
+if __name__ == "__main__":
+    do_demo()

From b716ca308311674793547d139d4089a19aa0e7de Mon Sep 17 00:00:00 2001
From: Sabina Sagynbayeva <83788641+ssagynbayeva@users.noreply.github.com>
Date: Thu, 9 Mar 2023 14:34:53 -0500
Subject: [PATCH 03/13] Update pyro/mesh/patch.py

Co-authored-by: Eric T. Johnson <yut23@users.noreply.github.com>
---
 pyro/mesh/patch.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 4dea2db2e..89f85b1bd 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -34,9 +34,9 @@
 import h5py
 import numpy as np
 
-# import pyro.mesh.boundary as bnd
-# from pyro.mesh.array_indexer import ArrayIndexer, ArrayIndexerFC
-# from pyro.util import msg
+import pyro.mesh.boundary as bnd
+from pyro.mesh.array_indexer import ArrayIndexer, ArrayIndexerFC
+from pyro.util import msg
 
 
 class Grid2d:

From 196e9df08990e5cde8cb294301faf278163e2d3c Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Wed, 15 Mar 2023 11:53:27 -0400
Subject: [PATCH 04/13] added A_x and A_y

---
 pyro/mesh/patch.py | 38 ++++++++++++++++++++++++++++++++++++--
 1 file changed, 36 insertions(+), 2 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 89f85b1bd..83d2d4025 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -1011,6 +1011,41 @@ def face_area(self):
         area_j = np.pi * (r0 * np.sin(t0) + r1 * np.sin(t0)) * np.sqrt(np.square(r1 * np.sin(t0) - r0 * np.sin(t0)) + np.square(r1 * np.cos(t0) - r0 * np.cos(t0)))
         return tr(np.array([area_i, area_j]))
 
+    def area_x(self):
+        """
+        Return an array of the face areas. 
+        The shape of the returned array is (ni, nj).
+        """
+        tr = lambda arr: arr.transpose(1, 2, 0)
+        x = self.cell_vertices()[:,0]
+        y = self.cell_vertices()[0,:]
+        r0 = x[:-1, :-1]
+        r1 = x[+1:, :-1]
+        t0 = y[:-1, :-1]
+        t1 = y[+1:, +1:]
+
+        area = r1 - r0
+        return area
+
+    def area_y(self):
+        """
+        Return an array of the face areas. 
+        The shape of the returned array is (ni, nj).
+        """
+        tr = lambda arr: arr.transpose(1, 2, 0)
+        x = self.cell_vertices()[:,0]
+        y = self.cell_vertices()[0,:]
+        r0 = x[:-1, :-1]
+        r1 = x[+1:, :-1]
+        t0 = y[:-1, :-1]
+        t1 = y[+1:, +1:]
+
+        # ** the area of a part of an annulus
+
+        area = 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0)
+        return area
+
+
     def cell_volumes(self):
         """
         Return an array of the cell volume data for the given coordinate box
@@ -1025,8 +1060,7 @@ def cell_volumes(self):
         t0 = y[:-1, :-1]
         t1 = y[+1:, +1:]
 
-        return (r1 ** 3 - r0 ** 3) * (np.cos(t1) - np.cos(t0)) * (-2.0 * np.pi) / 3.0
-        # return r1
+        return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0) * (r1 - r0)
 
     def cell_vertices(self):
         """

From bbe762df96d598d9c92333b0ce8904fec852d59e Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Tue, 21 Mar 2023 17:23:15 -0400
Subject: [PATCH 05/13] fixed the areas. They are now using self.xr and self.yr
 instead of cell_vertices

---
 pyro/mesh/patch.py | 50 +++++++++++++++++++++++++++-------------------
 1 file changed, 30 insertions(+), 20 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 83d2d4025..2ee33c245 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -1016,13 +1016,15 @@ def area_x(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-        tr = lambda arr: arr.transpose(1, 2, 0)
-        x = self.cell_vertices()[:,0]
-        y = self.cell_vertices()[0,:]
-        r0 = x[:-1, :-1]
-        r1 = x[+1:, :-1]
-        t0 = y[:-1, :-1]
-        t1 = y[+1:, +1:]
+        # tr = lambda arr: arr.transpose(1, 2, 0)
+        # x = self.cell_vertices()[:,0]
+        # y = self.cell_vertices()[0,:]
+        # r0 = x[:-1, :-1]
+        # r1 = x[+1:, :-1]
+        # t0 = y[:-1, :-1]
+        # t1 = y[+1:, +1:]
+        r1 = self.xl
+        r0 = self.xr
 
         area = r1 - r0
         return area
@@ -1032,13 +1034,16 @@ def area_y(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-        tr = lambda arr: arr.transpose(1, 2, 0)
-        x = self.cell_vertices()[:,0]
-        y = self.cell_vertices()[0,:]
-        r0 = x[:-1, :-1]
-        r1 = x[+1:, :-1]
-        t0 = y[:-1, :-1]
-        t1 = y[+1:, +1:]
+        # tr = lambda arr: arr.transpose(1, 2, 0)
+        # x = self.cell_vertices()[:,0]
+        # y = self.cell_vertices()[0,:]
+        # r0 = x[:-1, :-1]
+        # r1 = x[+1:, :-1]
+        # t0 = y[:-1, :-1]
+        # t1 = y[+1:, +1:]
+
+        r1, t1 = np.meshgrid(self.xl, self.yl)
+        r0, t0 = np.meshgrid(self.xr, self.yr)
 
         # ** the area of a part of an annulus
 
@@ -1052,13 +1057,18 @@ def cell_volumes(self):
         The shape of the returned array is (ni, nj).
         """
 
-        x = self.cell_vertices()[:,0]
-        y = self.cell_vertices()[0,:]
+        # x = self.cell_vertices()[:,0]
+        # y = self.cell_vertices()[0,:]
 
-        r0 = x[:-1, :-1]
-        r1 = x[+1:, :-1]
-        t0 = y[:-1, :-1]
-        t1 = y[+1:, +1:]
+        # r0 = x[:-1, :-1]
+        # r1 = x[+1:, :-1]
+        # t0 = y[:-1, :-1]
+        # t1 = y[+1:, +1:]
+
+        r1 = self.xl
+        r0 = self.xr
+        t1 = self.yl
+        t0 = self.yr
 
         return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0) * (r1 - r0)
 

From c377f3bcb7266ccc441fe36105d6f28c104b2d78 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Tue, 21 Mar 2023 17:28:05 -0400
Subject: [PATCH 06/13] removed _init_

---
 pyro/mesh/patch.py | 121 +++++++--------------------------------------
 1 file changed, 18 insertions(+), 103 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index bbad5d92f..91c96d728 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -905,110 +905,25 @@ class PolarGrid(Grid2d):
 
     # pylint: disable=too-many-instance-attributes
 
-    def __init__(self, nx, ny, ng=1,
-                 xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0):
-
-        """
-        Create a PolarGrid object.
-
-        The only data that we require is the number of points that
-        make up the mesh in each direction.  Optionally we take the
-        extrema of the domain (default is [0,1]x[0,1]) and number of
-        ghost cells (default is 1).
-
-        Note that the Grid2d object only defines the discretization,
-        it does not know about the boundary conditions, as these can
-        vary depending on the variable.
-
-        Parameters
-        ----------
-        nx : int
-            Number of zones in the r-direction
-        ny : int
-            Number of zones in the theta-direction
-        ng : int, optional
-            Number of ghost cells
-        xmin : float, optional
-            Physical coordinate at the lower x boundary
-        xmax : float, optional
-            Physical coordinate at the upper x boundary
-        ymin : float, optional
-            Physical coordinate at the lower y boundary
-        ymax : float, optional
-            Physical coordinate at the upper y boundary
-        """
-
-        # pylint: disable=too-many-arguments
-
-        # size of grid
-        self.nx = int(nx)
-        self.ny = int(ny)
-        self.ng = int(ng)
-
-        self.qx = int(2*ng + nx)
-        self.qy = int(2*ng + ny)
 
-        # domain extrema
-        self.xmin = xmin
-        self.xmax = xmax
-
-        self.ymin = ymin
-        self.ymax = ymax
-
-        # compute the indices of the block interior (excluding guardcells)
-        self.ilo = self.ng
-        self.ihi = self.ng + self.nx-1
-
-        self.jlo = self.ng
-        self.jhi = self.ng + self.ny-1
-
-        # center of the grid (for convenience)
-        self.ic = self.ilo + self.nx//2 - 1
-        self.jc = self.jlo + self.ny//2 - 1
-
-        # define the coordinate information at the left, center, and right
-        # zone coordinates
-        self.dx = (xmax - xmin)/nx
-
-        self.xl = (np.arange(self.qx) - ng)*self.dx + xmin
-        self.xr = (np.arange(self.qx) + 1.0 - ng)*self.dx + xmin
-        self.x = 0.5*(self.xl + self.xr)
-
-        self.dy = (ymax - ymin)/ny
-
-        self.yl = (np.arange(self.qy) - ng)*self.dy + ymin
-        self.yr = (np.arange(self.qy) + 1.0 - ng)*self.dy + ymin
-        self.y = 0.5*(self.yl + self.yr)
-
-        # 2-d versions of the zone coordinates (replace with meshgrid?)
-        x2d = np.repeat(self.x, self.qy)
-        x2d.shape = (self.qx, self.qy)
-        self.x2d = x2d
-
-        y2d = np.repeat(self.y, self.qx)
-        y2d.shape = (self.qy, self.qx)
-        y2d = np.transpose(y2d)
-        self.y2d = y2d
-
-
-    def face_area(self):
-        """
-        Return an array of the face areas. 
-        The shape of the returned array is (ni, nj).
-        """
-        tr = lambda arr: arr.transpose(1, 2, 0)
-        x = self.cell_vertices()[:,0]
-        y = self.cell_vertices()[0,:]
-        r0 = x[:-1, :-1]
-        r1 = x[+1:, :-1]
-        t0 = y[:-1, :-1]
-        t1 = y[+1:, +1:]
-
-        # ** the area of the arc
-
-        area_i = np.pi * (r0 * np.sin(t0) + r0 * np.sin(t1)) * np.sqrt(np.square(r0 * np.sin(t1) - r0 * np.sin(t0)) + np.square(r0 * np.cos(t1) - r0 * np.cos(t0)))
-        area_j = np.pi * (r0 * np.sin(t0) + r1 * np.sin(t0)) * np.sqrt(np.square(r1 * np.sin(t0) - r0 * np.sin(t0)) + np.square(r1 * np.cos(t0) - r0 * np.cos(t0)))
-        return tr(np.array([area_i, area_j]))
+    # def face_area(self):
+    #     """
+    #     Return an array of the face areas. 
+    #     The shape of the returned array is (ni, nj).
+    #     """
+    #     tr = lambda arr: arr.transpose(1, 2, 0)
+    #     x = self.cell_vertices()[:,0]
+    #     y = self.cell_vertices()[0,:]
+    #     r0 = x[:-1, :-1]
+    #     r1 = x[+1:, :-1]
+    #     t0 = y[:-1, :-1]
+    #     t1 = y[+1:, +1:]
+
+    #     # ** the area of the arc
+
+    #     area_i = np.pi * (r0 * np.sin(t0) + r0 * np.sin(t1)) * np.sqrt(np.square(r0 * np.sin(t1) - r0 * np.sin(t0)) + np.square(r0 * np.cos(t1) - r0 * np.cos(t0)))
+    #     area_j = np.pi * (r0 * np.sin(t0) + r1 * np.sin(t0)) * np.sqrt(np.square(r1 * np.sin(t0) - r0 * np.sin(t0)) + np.square(r1 * np.cos(t0) - r0 * np.cos(t0)))
+    #     return tr(np.array([area_i, area_j]))
 
     def area_x(self):
         """

From 7405339136a8d0fb1e4ee30dad2ea80f9f82edec Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Mon, 3 Apr 2023 16:13:52 -0400
Subject: [PATCH 07/13] changed areas and the cell volume

---
 pyro/mesh/patch.py | 40 +++++++++++-----------------------------
 1 file changed, 11 insertions(+), 29 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 91c96d728..921204111 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -930,15 +930,10 @@ def area_x(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-        # tr = lambda arr: arr.transpose(1, 2, 0)
-        # x = self.cell_vertices()[:,0]
-        # y = self.cell_vertices()[0,:]
-        # r0 = x[:-1, :-1]
-        # r1 = x[+1:, :-1]
-        # t0 = y[:-1, :-1]
-        # t1 = y[+1:, +1:]
-        r1 = self.xl
-        r0 = self.xr
+        r1 = self.xr
+        r0 = self.xl
+
+        # ** this is just \Delta r
 
         area = r1 - r0
         return area
@@ -948,20 +943,13 @@ def area_y(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-        # tr = lambda arr: arr.transpose(1, 2, 0)
-        # x = self.cell_vertices()[:,0]
-        # y = self.cell_vertices()[0,:]
-        # r0 = x[:-1, :-1]
-        # r1 = x[+1:, :-1]
-        # t0 = y[:-1, :-1]
-        # t1 = y[+1:, +1:]
 
-        r1, t1 = np.meshgrid(self.xl, self.yl)
-        r0, t0 = np.meshgrid(self.xr, self.yr)
+        t1 = self.yr
+        t0 = self.yl
 
-        # ** the area of a part of an annulus
+        # ** this is just \Delta \theta
 
-        area = 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0)
+        area = t1 - t0
         return area
 
 
@@ -971,20 +959,14 @@ def cell_volumes(self):
         The shape of the returned array is (ni, nj).
         """
 
-        # x = self.cell_vertices()[:,0]
-        # y = self.cell_vertices()[0,:]
-
-        # r0 = x[:-1, :-1]
-        # r1 = x[+1:, :-1]
-        # t0 = y[:-1, :-1]
-        # t1 = y[+1:, +1:]
-
         r1 = self.xl
         r0 = self.xr
         t1 = self.yl
         t0 = self.yr
 
-        return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0) * (r1 - r0)
+        # ** this is just the face area
+
+        return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0)
 
     def cell_vertices(self):
         """

From df648dcd784439b995d206258bf697a8945154c0 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Mon, 3 Apr 2023 16:25:43 -0400
Subject: [PATCH 08/13] fixed flake8 errors

---
 pyro/mesh/patch.py | 47 ----------------------------------------------
 1 file changed, 47 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 921204111..f7336922d 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -905,26 +905,6 @@ class PolarGrid(Grid2d):
 
     # pylint: disable=too-many-instance-attributes
 
-
-    # def face_area(self):
-    #     """
-    #     Return an array of the face areas. 
-    #     The shape of the returned array is (ni, nj).
-    #     """
-    #     tr = lambda arr: arr.transpose(1, 2, 0)
-    #     x = self.cell_vertices()[:,0]
-    #     y = self.cell_vertices()[0,:]
-    #     r0 = x[:-1, :-1]
-    #     r1 = x[+1:, :-1]
-    #     t0 = y[:-1, :-1]
-    #     t1 = y[+1:, +1:]
-
-    #     # ** the area of the arc
-
-    #     area_i = np.pi * (r0 * np.sin(t0) + r0 * np.sin(t1)) * np.sqrt(np.square(r0 * np.sin(t1) - r0 * np.sin(t0)) + np.square(r0 * np.cos(t1) - r0 * np.cos(t0)))
-    #     area_j = np.pi * (r0 * np.sin(t0) + r1 * np.sin(t0)) * np.sqrt(np.square(r1 * np.sin(t0) - r0 * np.sin(t0)) + np.square(r1 * np.cos(t0) - r0 * np.cos(t0)))
-    #     return tr(np.array([area_i, area_j]))
-
     def area_x(self):
         """
         Return an array of the face areas. 
@@ -952,7 +932,6 @@ def area_y(self):
         area = t1 - t0
         return area
 
-
     def cell_volumes(self):
         """
         Return an array of the cell volume data for the given coordinate box
@@ -968,31 +947,5 @@ def cell_volumes(self):
 
         return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0)
 
-    def cell_vertices(self):
-        """
-        Return the coordinates of cell vertices
-        The arrays range from 0 to 1 for now
-        """
-        # if self.ng == 0:
-        #     xv = np.linspace(0, 1, self.nx + 1)[:-1]
-        #     yv = np.linspace(0, 1, self.ny + 1)[:-1]
-        # else:
-        #     xv = np.linspace(0, 1, self.nx + 1)
-        #     yv = np.linspace(0, 1, self.ny + 1)
-
-        xv = np.linspace(0, 1, self.nx + 1)[:-1]
-        yv = np.linspace(0, 1, self.ny + 1)[:-1]
-
-     
-        tr = lambda arr: arr.transpose(1, 2, 0)
-        x, y = np.meshgrid(xv, yv, indexing="ij")
-        return tr(np.array([x, y]))
-
-   
-
-    
-
-
-
 if __name__ == "__main__":
     do_demo()

From 8e8e0d6eb83e42201e3f99dedcb2f503510ed06b Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Mon, 3 Apr 2023 16:29:55 -0400
Subject: [PATCH 09/13] added a line before main

---
 pyro/mesh/patch.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index f7336922d..8c6e2fa88 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -947,5 +947,6 @@ def cell_volumes(self):
 
         return 0.5 * (r1 ** 2 - r0 ** 2) * (t1 - t0)
 
+
 if __name__ == "__main__":
     do_demo()

From 55c837708b9ab38e6d602a5b53825ef464b37765 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Wed, 5 Apr 2023 13:55:27 -0400
Subject: [PATCH 10/13] recomputed areas

---
 pyro/mesh/patch.py | 23 +++++++++++++++++------
 1 file changed, 17 insertions(+), 6 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 8c6e2fa88..9c30b465c 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -900,6 +900,15 @@ class PolarGrid(Grid2d):
 
        |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->|
 
+              ____
+             /    |  
+            /     |
+           /      | 
+          /       |
+         /        |
+        /_________|
+
+
     The '*' marks the data locations.
     """
 
@@ -912,10 +921,12 @@ def area_x(self):
         """
         r1 = self.xr
         r0 = self.xl
+        t1 = self.yr
+        t0 = self.yl
 
-        # ** this is just \Delta r
+        # ** this is just 1/2*dr*d\theta
 
-        area = r1 - r0
+        area = 0.5 * (r1 - r0) * (t1 - t0)
         return area
 
     def area_y(self):
@@ -924,12 +935,12 @@ def area_y(self):
         The shape of the returned array is (ni, nj).
         """
 
-        t1 = self.yr
-        t0 = self.yl
+        r1 = self.xr
+        r0 = self.xl
 
-        # ** this is just \Delta \theta
+        # ** this is just dr
 
-        area = t1 - t0
+        area = r1 - r0
         return area
 
     def cell_volumes(self):

From 8fa6b97557ce3f0c9440236c94a173e116bd29e5 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Wed, 5 Apr 2023 14:08:24 -0400
Subject: [PATCH 11/13] ASCII picture

---
 pyro/mesh/patch.py | 22 +++++++---------------
 1 file changed, 7 insertions(+), 15 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 3021023b3..2193521ae 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -890,26 +890,18 @@ class PolarGrid(Grid2d):
     the 2-d grid class.  The grid object will contain the coordinate
     information (at various centerings).
 
-    A basic (1-d) representation of the layout is::
-
-       |     |      |     X     |     |      |     |     X     |      |     |
-       +--*--+- // -+--*--X--*--+--*--+- // -+--*--+--*--X--*--+- // -+--*--+
-          0          ng-1    ng   ng+1         ... ng+nx-1 ng+nx      2ng+nx-1
+    A basic representation of the layout is::
 
-                            ilo                      ihi
-
-       |<- ng guardcells->|<---- nx interior zones ----->|<- ng guardcells->|
-
-              ____
-             /    |  
+              *---* \theta_{i+1/2}
+             /    | 
             /     |
-           /      | 
+           /      * \theta_i 
           /       |
          /        |
-        /_________|
+        *____*____* \theta_{i-1/2}
+r_{i-1/2}   r_i   r_{i+1/2}
 
-
-    The '*' marks the data locations.
+    The '*' marks represent the vertices; i index is the data location.
     """
 
     # pylint: disable=too-many-instance-attributes

From 54ecf2d35458968a59daa75e730f6d5fac287015 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Wed, 5 Apr 2023 18:22:07 -0400
Subject: [PATCH 12/13] added unit tests for the PolarGrid()

---
 pyro/mesh/tests/test_patch.py | 76 +++++++++++++++++++++++++++++++++++
 1 file changed, 76 insertions(+)

diff --git a/pyro/mesh/tests/test_patch.py b/pyro/mesh/tests/test_patch.py
index d006f9078..46d07212a 100644
--- a/pyro/mesh/tests/test_patch.py
+++ b/pyro/mesh/tests/test_patch.py
@@ -244,3 +244,79 @@ def test_bcs():
     # top
     assert_array_equal(d[myg.ilo:myg.ihi+1, myg.jhi-1:myg.jhi+1],
                        -np.fliplr(d[myg.ilo:myg.ihi+1, myg.jhi+1:myg.jhi+3]))
+
+
+# PolarGrid tests
+class TestPolarGrid(object):
+    @classmethod
+    def setup_class(cls):
+        """ this is run once for each class before any tests """
+        pass
+
+    @classmethod
+    def teardown_class(cls):
+        """ this is run once for each class after all tests """
+        pass
+
+    def setup_method(self):
+        """ this is run before each test """
+        self.g = patch.PolarGrid(4, 6, ng=2, ymax=1.5)
+
+    def teardown_method(self):
+        """ this is run after each test """
+        self.g = None
+
+    def test_dx_dy(self):
+        assert self.g.dx == 0.25
+        assert self.g.dy == 0.25
+
+    def test_grid_coords(self):
+        assert_array_equal(self.g.x[self.g.ilo:self.g.ihi+1],
+                           np.array([0.125, 0.375, 0.625, 0.875]))
+        assert_array_equal(self.g.y[self.g.jlo:self.g.jhi+1],
+                           np.array([0.125, 0.375, 0.625, 0.875, 1.125, 1.375]))
+
+    def test_grid_2d_coords(self):
+        assert_array_equal(self.g.x, self.g.x2d[:, self.g.jc])
+        assert_array_equal(self.g.y, self.g.y2d[self.g.ic, :])
+
+    def test_scratch_array(self):
+        q = self.g.scratch_array()
+        assert q.shape == (self.g.qx, self.g.qy)
+
+    def test_coarse_like(self):
+        q = self.g.coarse_like(2)
+        assert q.qx == 2*self.g.ng + self.g.nx//2
+        assert q.qy == 2*self.g.ng + self.g.ny//2
+
+    def test_fine_like(self):
+        q = self.g.fine_like(2)
+        assert q.qx == 2*self.g.ng + 2*self.g.nx
+        assert q.qy == 2*self.g.ng + 2*self.g.ny
+
+    def test_norm(self):
+        q = self.g.scratch_array()
+        # there are 24 elements, the norm L2 norm is
+        # sqrt(dx*dy*24)
+        q.v()[:, :] = np.array([[1, 1, 1, 1, 1, 1],
+                               [1, 1, 1, 1, 1, 1],
+                               [1, 1, 1, 1, 1, 1],
+                               [1, 1, 1, 1, 1, 1]])
+
+        assert q.norm() == np.sqrt(24*self.g.dx*self.g.dy)
+
+    def test_equality(self):
+        g2 = patch.PolarGrid(2, 5, ng=1)
+        assert g2 != self.g
+
+    def test_area_x(self):
+        A = self.g.area_x()
+        assert A[0,0] == (self.g.yr[0] - self.g.yl[0]) * (self.g.xr[0] - self.g.xl[0]) * 0.5
+
+    def test_area_y(self):
+        A = self.g.area_y()
+        assert A[0] == (self.g.xr - self.g.xl)
+
+    def test_cell_volumes(self):
+        V = self.g.cell_volumes()
+        assert V[0,0] == (self.g.yr[0] - self.g.yl[0]) * (self.g.xr[0] ** 2 - self.g.xl[0] ** 2) * 0.5

From 6173abdd4fbb9b46772a4e536ae1f59245a663f9 Mon Sep 17 00:00:00 2001
From: Sabina-star <83788641+Sabina-star@users.noreply.github.com>
Date: Thu, 6 Apr 2023 15:22:36 -0400
Subject: [PATCH 13/13] changes

---
 pyro/mesh/patch.py | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/pyro/mesh/patch.py b/pyro/mesh/patch.py
index 2193521ae..5801aadab 100644
--- a/pyro/mesh/patch.py
+++ b/pyro/mesh/patch.py
@@ -911,14 +911,12 @@ def area_x(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-        r1 = self.xr
-        r0 = self.xl
-        t1 = self.yr
-        t0 = self.yl
+        r1, t1 = np.meshgrid(self.xr, self.yr)
+        r0, t0 = np.meshgrid(self.xl, self.yl)
 
-        # ** this is just 1/2*dr*d\theta
+        # ** this is just 1/2*r*d\theta
 
-        area = 0.5 * (r1 - r0) * (t1 - t0)
+        area = 0.5 * r0 * (t1 - t0)
         return area
 
     def area_y(self):
@@ -926,9 +924,8 @@ def area_y(self):
         Return an array of the face areas. 
         The shape of the returned array is (ni, nj).
         """
-
-        r1 = self.xr
-        r0 = self.xl
+        r1, t1 = np.meshgrid(self.xr, self.yr)
+        r0, t0 = np.meshgrid(self.xl, self.yl)
 
         # ** this is just dr
 
@@ -940,11 +937,8 @@ def cell_volumes(self):
         Return an array of the cell volume data for the given coordinate box
         The shape of the returned array is (ni, nj).
         """
-
-        r1 = self.xl
-        r0 = self.xr
-        t1 = self.yl
-        t0 = self.yr
+        r1, t1 = np.meshgrid(self.xr, self.yr)
+        r0, t0 = np.meshgrid(self.xl, self.yl)
 
         # ** this is just the face area