Skip to content

Commit

Permalink
get rid of redundant A_r, fix A_y calcualtion for sphericalpolar and …
Browse files Browse the repository at this point in the history
…some cleanup
  • Loading branch information
zhichen3 committed Aug 1, 2024
1 parent a6ac3d6 commit 0f1ed3c
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 44 deletions.
40 changes: 11 additions & 29 deletions pyro/mesh/patch.py
Original file line number Diff line number Diff line change
Expand Up @@ -214,13 +214,11 @@ def __init__(self, nx, ny, ng=1,

# This is area of the side that is perpendicular to x.

self.Ax_l = self.Ly.copy()
self.Ax_r = self.Ly.copy()
self.Ax = self.Ly

# This is area of the side that is perpendicular to y.

self.Ay_l = self.Lx.copy()
self.Ay_r = self.Lx.copy()
self.Ay = self.Lx

# Volume of the cell.

Expand Down Expand Up @@ -272,41 +270,25 @@ def __init__(self, nx, ny, ng=1,
# dL_theta x dL_phi = r^2 * sin(theta) * dtheta * dphi

# 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.v(buf=self.ng)**2 *
(np.cos(self.yr2d.v(buf=self.ng)) -
np.cos(self.yl2d.v(buf=self.ng))))
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.v(buf=self.ng)**2 *
(np.cos(self.yr2d.v(buf=self.ng)) -
np.cos(self.yl2d.v(buf=self.ng))))
self.Ax_r = ArrayIndexer(area_x_r, grid=self)
self.Ax = np.abs(-2.0 * np.pi * self.xl2d**2 *
(np.cos(self.yr2d) - np.cos(self.yl2d)))

# Returns an array of the face area that points in the theta(y) direction.
# dL_phi x dL_r = dr * r * sin(theta) * dphi

# 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.v(buf=self.ng)) *
(self.xr2d.v(buf=self.ng)**2 -
self.xl2d.v(buf=self.ng)**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.v(buf=self.ng)) *
(self.xr2d.v(buf=self.ng)**2 -
self.xl2d.v(buf=self.ng)**2))
self.Ay_r = ArrayIndexer(area_y_r, grid=self)
# dAtheta_l = pi * sin(theta{i-1/2}) * (r{i+1/2}^2 - r{i-1/2}^2)
self.Ay = np.abs(np.pi * np.sin(self.yl2d) *
(self.xr2d**2 - self.xl2d**2))

# Returns an array of the volume of each cell.
# dV = dL_r * dL_theta * dL_phi
# = (dr) * (r * dtheta) * (r * sin(theta) * dphi)
# dV = - 2*np.pi / 3 * (cos(theta{i+1/2}) - cos(theta{i-1/2})) * (r{i+1/2}^3 - r{i-1/2}^3)

volume = np.abs(-2.0 * np.pi / 3.0 *
(np.cos(self.yr2d.v(buf=self.ng)) - np.cos(self.yl2d.v(buf=self.ng))) *
(self.xr2d.v(buf=self.ng)**3 - self.xl2d.v(buf=self.ng)**3))
self.V = ArrayIndexer(volume, grid=self)
self.V = np.abs(-2.0 * np.pi / 3.0 *
(np.cos(self.yr2d) - np.cos(self.yl2d)) *
(self.xr2d - self.xl2d) *
(self.xr2d**2 + self.xl2d**2 + self.xr2d*self.xl2d))

def __str__(self):
""" print out some basic information about the grid object """
Expand Down
28 changes: 13 additions & 15 deletions pyro/mesh/tests/test_patch.py
Original file line number Diff line number Diff line change
Expand Up @@ -87,12 +87,10 @@ def teardown_method(self):
self.g = None

def test_Ax(self):
assert np.all(self.g.Ax_l.v() == 0.1)
assert np.all(self.g.Ax_r.v() == 0.1)
assert np.all(self.g.Ax.v() == 0.1)

def test_Ay(self):
assert np.all(self.g.Ay_l.v() == 0.25)
assert np.all(self.g.Ay_r.v() == 0.25)
assert np.all(self.g.Ay.v() == 0.25)

def test_V(self):
assert np.all(self.g.V.v() == 0.1 * 0.25)
Expand Down Expand Up @@ -120,22 +118,22 @@ def teardown_method(self):
def test_Ax(self):
area_x_l = np.abs(-2.0 * np.pi * self.g.xl2d.v()**2 *
(np.cos(self.g.yr2d.v()) - np.cos(self.g.yl2d.v())))
assert_array_equal(self.g.Ax_l.v(), area_x_l)
assert_array_equal(self.g.Ax.v(), area_x_l)

area_x_r = np.abs(-2.0 * np.pi * self.g.xr2d.v()**2 *
(np.cos(self.g.yr2d.v()) - np.cos(self.g.yl2d.v())))
assert_array_equal(self.g.Ax_r.v(), area_x_r)
assert_array_equal(self.g.Ax.ip(1), area_x_r)

def test_Ay(self):
area_y_l = np.abs(0.5 * np.pi *
np.sin(self.g.yl2d.v()) *
(self.g.xr2d.v()**2 - self.g.xl2d.v()**2))
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.v()) *
(self.g.xr2d.v()**2 - self.g.xl2d.v()**2))
assert_array_equal(self.g.Ay_r.v(), area_y_r)
area_y_l = np.abs(np.pi *
np.sin(self.g.yl2d.v()) *
(self.g.xr2d.v()**2 - self.g.xl2d.v()**2))
assert_array_equal(self.g.Ay.v(), area_y_l)

area_y_r = np.abs(np.pi *
np.sin(self.g.yr2d.v()) *
(self.g.xr2d.v()**2 - self.g.xl2d.v()**2))
assert_array_equal(self.g.Ay.jp(1), area_y_r)

def test_V(self):
volume = np.abs(-2.0 * np.pi / 3.0 *
Expand Down

0 comments on commit 0f1ed3c

Please sign in to comment.