Skip to content

Commit

Permalink
make advection interface more readable
Browse files Browse the repository at this point in the history
  • Loading branch information
zhichen3 committed May 29, 2024
1 parent 62423e1 commit 829d164
Show file tree
Hide file tree
Showing 2 changed files with 15 additions and 23 deletions.
31 changes: 12 additions & 19 deletions pyro/advection/interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,6 @@ def linear_interface(a, myg, rp, dt):
u = rp.get_param("advection.u")
v = rp.get_param("advection.v")

cx = u*dt/myg.dx
cy = v*dt/myg.dy

# --------------------------------------------------------------------------
# monotonized central differences
# --------------------------------------------------------------------------
Expand All @@ -27,30 +24,26 @@ def linear_interface(a, myg, rp, dt):

# upwind
if u < 0:
# a_x[i,j] = a[i,j] - 0.5*(1.0 + cx)*ldelta_a[i,j]
# a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cx)*ldelta_ax.v(buf=1)
a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*myg.dx*ldelta_ax.v(buf=1) \
- 0.5*dt*ldelta_Fx.v(buf=1)
# a_x[i,j] = a[i,j] - 0.5*ldelta_ax[i,j] - 0.5*dt/dx*ldelta_Fx[i,j]
a_x.v(buf=1)[:, :] = a.v(buf=1) - 0.5*ldelta_ax.v(buf=1) \
- 0.5*dt/myg.dx*ldelta_Fx.v(buf=1)

else:
# a_x[i,j] = a[i-1,j] + 0.5*(1.0 - cx)*ldelta_a[i-1,j]
# a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*(1.0 - cx)*ldelta_ax.ip(-1, buf=1)
a_x.v(buf=1)[:, :] = a.v(buf=1) + 0.5*myg.dy*ldelta_ax.ip(-1, buf=1)
- 0.5*dt*ldelta_Fx.ip(-1, buf=1)
# a_x[i,j] = a[i-1,j] + 0.5*ldelta_ax[i-1,j] - 0.5*dt/dx*ldelta_Fx[i-1,j]
a_x.v(buf=1)[:, :] = a.ip(-1, buf=1) + 0.5*ldelta_ax.ip(-1, buf=1) \
- 0.5*dt/myg.dx*ldelta_Fx.ip(-1, buf=1)

# y-direction
a_y = myg.scratch_array()

# upwind
if v < 0:
# a_y[i,j] = a[i,j] - 0.5*(1.0 + cy)*ldelta_a[i,j]
# a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*(1.0 + cy)*ldelta_ay.v(buf=1)
a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*myg.dy*ldelta_ay.v(buf=1) \
- 0.5*dt*ldelta_Fy.v(buf=1)
# a_y[i,j] = a[i,j] - 0.5*ldelta_ay[i,j] - 0.5*dt/dy*ldelta_Fy[i,j]
a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*ldelta_ay.v(buf=1) \
- 0.5*dt/myg.dy*ldelta_Fy.v(buf=1)
else:
# a_y[i,j] = a[i,j-1] + 0.5*(1.0 - cy)*ldelta_a[i,j-1]
# a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*(1.0 - cy)*ldelta_ay.jp(-1, buf=1)
a_y.v(buf=1)[:, :] = a.v(buf=1) - 0.5*myg.dy*ldelta_ay.jp(-1, buf=1) \
- 0.5*dt*ldelta_Fy.jp(-1, buf=1)
# a_y[i,j] = a[i,j-1] + 0.5*ldelta_ay[i,j-1] - 0.5*dt/dy*ldelta_Fy[i,j-1]
a_y.v(buf=1)[:, :] = a.jp(-1, buf=1) + 0.5*ldelta_ay.jp(-1, buf=1) \
- 0.5*dt/myg.dy*ldelta_Fy.jp(-1, buf=1)

return u, v, a_x, a_y
7 changes: 3 additions & 4 deletions pyro/mesh/patch.py
Original file line number Diff line number Diff line change
Expand Up @@ -234,7 +234,7 @@ def A_y(self, j=0, buf=0):

def V(self, i=0, j=0, buf=0):
"""
Returns the area
Returns the volume
parameter:
-----------
i : shifts the array in the x-direction by index i
Expand All @@ -245,6 +245,7 @@ def V(self, i=0, j=0, buf=0):
return self.vol[self.ilo+i+buf:self.ihi+1+i+buf,
self.jlo+j+buf:self.jhi+1+j+buf]


class SphericalPolar(Grid2d):
"""
This class defines a spherical polar grid.
Expand Down Expand Up @@ -311,7 +312,7 @@ def A_y(self, j=0, buf=0):

def V(self, i=0, j=0, buf=0):
"""
Returns the area
Returns the volume
parameter:
-----------
i : shifts the array in the x-direction by index i
Expand All @@ -323,8 +324,6 @@ def V(self, i=0, j=0, buf=0):
self.jlo+j+buf:self.jhi+1+j+buf]




class CellCenterData2d:
"""
A class to define cell-centered data that lives on a grid. A
Expand Down

0 comments on commit 829d164

Please sign in to comment.