Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

add support for spherical polar grid: add geometric terms (volume/area) #201

Merged
merged 25 commits into from
Aug 1, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
b74c474
initial update
zhichen3 May 2, 2024
b37fbe8
update
zhichen3 May 2, 2024
62432bc
Merge branch 'spherical_polar_grid' of github.com:zhichen3/pyro2 into…
zhichen3 May 15, 2024
62423e1
update
zhichen3 May 20, 2024
829d164
make advection interface more readable
zhichen3 May 29, 2024
4f2d04e
revert advection files
zhichen3 Jul 12, 2024
4ba2390
Merge branch 'main' into spherical_polar_grid
zhichen3 Jul 12, 2024
d501d70
fix flake8
zhichen3 Jul 12, 2024
11f26d6
update example
zhichen3 Jul 12, 2024
c9a084b
use ArrayIndexer for area and volume
zhichen3 Jul 14, 2024
8f26cb1
fix flake8
zhichen3 Jul 14, 2024
c823b37
add unit test to cartesian2d and spherical polar
zhichen3 Jul 16, 2024
c921ebb
add __str__ and update mesh-example to demonstrate sphericalPolar gri…
zhichen3 Jul 16, 2024
2174795
fix flake8
zhichen3 Jul 16, 2024
c0e9164
fix pytest
zhichen3 Jul 16, 2024
1de6124
update notebook
zhichen3 Jul 16, 2024
b93c187
add left and right area
zhichen3 Jul 20, 2024
c246e8a
Merge branch 'main' into spherical_polar_grid
zhichen3 Jul 20, 2024
1dab7ae
change x2d y2d to arrayindexer
zhichen3 Jul 20, 2024
1ebe1ad
remove unused variable
zhichen3 Jul 20, 2024
bb91e21
add differential length of each direction
zhichen3 Jul 24, 2024
43db316
change Lx and Ly to 1d array or float
zhichen3 Jul 25, 2024
0151d78
revert
zhichen3 Jul 25, 2024
a6ac3d6
add coord_type to grid class
zhichen3 Jul 31, 2024
0f1ed3c
get rid of redundant A_r, fix A_y calcualtion for sphericalpolar and …
zhichen3 Aug 1, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 9 additions & 8 deletions examples/examples.ipynb

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions pyro/_defaults
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ store_images = 0 ; store vis images to files (1=yes, 0=no)

[mesh]

grid_type = Cartesian2d ; Geometry of the Grid ('Cartesian2d' or 'SphericalPolar')
xmin = 0.0 ; domain minimum x-coordinate
xmax = 1.0 ; domain maximum x-coordinate
ymin = 0.0 ; domain minimum y-coordinate
Expand Down
4 changes: 2 additions & 2 deletions pyro/compressible/interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -201,8 +201,8 @@ def states(idir, ng, dx, dt,

# construct the states
for m in range(nvar):
sum_l = np.dot(betal, rvec[:, m])
sum_r = np.dot(betar, rvec[:, m])
sum_l = np.dot(betal, np.ascontiguousarray(rvec[:, m]))
zingale marked this conversation as resolved.
Show resolved Hide resolved
sum_r = np.dot(betar, np.ascontiguousarray(rvec[:, m]))

if idir == 1:
q_l[i + 1, j, m] = q_l[i + 1, j, m] + sum_l
Expand Down
121 changes: 118 additions & 3 deletions pyro/mesh/mesh-examples.ipynb

Large diffs are not rendered by default.

128 changes: 120 additions & 8 deletions pyro/mesh/patch.py
Original file line number Diff line number Diff line change
Expand Up @@ -133,15 +133,18 @@ def __init__(self, nx, ny, ng=1,
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
# 2-d versions of the zone coordinates
x2d, y2d = np.meshgrid(self.x, self.y, indexing='ij')
self.x2d = ArrayIndexer(d=x2d, grid=self)
self.y2d = ArrayIndexer(d=y2d, grid=self)

y2d = np.repeat(self.y, self.qx)
y2d.shape = (self.qy, self.qx)
y2d = np.transpose(y2d)
self.y2d = y2d
xl2d, yl2d = np.meshgrid(self.xl, self.yl, indexing='ij')
self.xl2d = ArrayIndexer(d=xl2d, grid=self)
self.yl2d = ArrayIndexer(d=yl2d, grid=self)

xr2d, yr2d = np.meshgrid(self.xr, self.yr, indexing='ij')
self.xr2d = ArrayIndexer(d=xr2d, grid=self)
self.yr2d = ArrayIndexer(d=yr2d, grid=self)

def scratch_array(self, nvar=1):
"""
Expand Down Expand Up @@ -186,6 +189,115 @@ def __eq__(self, other):
return result


class Cartesian2d(Grid2d):
"""
This class defines a 2D Cartesian Grid.
Define:
x = x
y = y
"""

def __init__(self, nx, ny, ng=1,
xmin=0.0, xmax=1.0, ymin=0.0, ymax=1.0):

super().__init__(nx, ny, ng, xmin, xmax, ymin, ymax)

self.coord_type = 0

# Length of the side in x- and y-direction

self.Lx = ArrayIndexer(np.full((self.qx, self.qy), self.dx),
grid=self)
self.Ly = ArrayIndexer(np.full((self.qx, self.qy), self.dy),
grid=self)

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

self.Ax = self.Ly

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

self.Ay = self.Lx

# Volume of the cell.

self.V = ArrayIndexer(np.full((self.qx, self.qy), self.dx * self.dy),
grid=self)

def __str__(self):
""" print out some basic information about the grid object """
return f"Cartesian 2D Grid: xmin = {self.xmin}, xmax = {self.xmax}, " + \
f"ymin = {self.ymin}, ymax = {self.ymax}, " + \
f"nx = {self.nx}, ny = {self.ny}, ng = {self.ng}"


class SphericalPolar(Grid2d):
"""
This class defines a spherical polar grid.
This is technically a 2D geometry but assumes azimuthal symmetry.
Define:
r = x
theta = y
"""

def __init__(self, nx, ny, ng=1,
xmin=0.2, xmax=1.0, ymin=0.0, ymax=1.0):

# 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, \
"xmin (r-direction), must be large enough so ghost cell doesn't have negative x."

super().__init__(nx, ny, ng, xmin, xmax, ymin, ymax)

self.coord_type = 1

# Length of the side along r-direction, dr

self.Lx = ArrayIndexer(np.full((self.qx, self.qy), self.dx),
grid=self)

# Length of the side along theta-direction, r*dtheta

self.Ly = ArrayIndexer(np.full((self.qx, self.qy), self.x2d*self.dy),
grid=self)

# 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

# dAr_l = - r{i-1/2}^2 * 2pi * cos(theta{i+1/2}) - cos(theta{i-1/2})
self.Ax = np.abs(-2.0 * np.pi * self.xl2d**2 *
(np.cos(self.yr2d) - np.cos(self.yl2d)))

zingale marked this conversation as resolved.
Show resolved Hide resolved
# 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 = 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)

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 """
return "Spherical Polar 2D Grid: Define x : r, y : \u03b8. " + \
f"xmin (r) = {self.xmin}, xmax= {self.xmax}, " + \
f"ymin = {self.ymin}, ymax = {self.ymax}, " + \
f"nx = {self.nx}, ny = {self.ny}, ng = {self.ng}"


class CellCenterData2d:
"""
A class to define cell-centered data that lives on a grid. A
Expand Down
74 changes: 74 additions & 0 deletions pyro/mesh/tests/test_patch.py
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,80 @@ def test_equality(self):
assert g2 != self.g


# Cartesian2d tests
class TestCartesian2d:
@classmethod
def setup_class(cls):
""" this is run once for each class before any tests """

@classmethod
def teardown_class(cls):
""" this is run once for each class after all tests """

def setup_method(self):
""" this is run before each test """
self.g = patch.Cartesian2d(4, 10, ng=2)

def teardown_method(self):
""" this is run after each test """
self.g = None

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

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

def test_V(self):
assert np.all(self.g.V.v() == 0.1 * 0.25)


# SphericalPolar Grid tests
class TestSphericalPolar:
@classmethod
def setup_class(cls):
""" this is run once for each class before any tests """

@classmethod
def teardown_class(cls):
""" this is run once for each class after all tests """

def setup_method(self):
""" this is run before each test """
self.g = patch.SphericalPolar(4, 10, xmin=1.0, xmax=2.0,
ymax=np.pi, ng=2)

def teardown_method(self):
""" this is run after each test """
self.g = None

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.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.ip(1), area_x_r)

def test_Ay(self):
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 *
(np.cos(self.g.yr2d.v()) - np.cos(self.g.yl2d.v())) *
(self.g.xr2d.v()**3 - self.g.xl2d.v()**3))
assert_array_equal(self.g.V.v(), volume)


# CellCenterData2d tests
class TestCellCenterData2d:
@classmethod
Expand Down
23 changes: 19 additions & 4 deletions pyro/simulation_null.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,11 +34,26 @@ def grid_setup(rp, ng=1):
ymax = rp.get_param("mesh.ymax")
except KeyError:
ymax = 1.0
msg.warning("mesh.ynax not set, defaulting to 1.0")
msg.warning("mesh.ymax not set, defaulting to 1.0")

try:
grid_type = rp.get_param("mesh.grid_type")
except KeyError:
grid_type = "Cartesian2d"
msg.warning("mesh.grid_type not set, defaulting to Cartesian2D")

if grid_type == "Cartesian2d":
create_grid = patch.Cartesian2d
elif grid_type == "SphericalPolar":
create_grid = patch.SphericalPolar
else:
raise ValueError("Unsupported grid type!")

my_grid = create_grid(nx, ny,
xmin=xmin, xmax=xmax,
ymin=ymin, ymax=ymax,
ng=ng)

my_grid = patch.Grid2d(nx, ny,
xmin=xmin, xmax=xmax,
ymin=ymin, ymax=ymax, ng=ng)
return my_grid


Expand Down