Skip to content

Commit

Permalink
update artificial viscosity to be in sync with castro
Browse files Browse the repository at this point in the history
  • Loading branch information
zhichen3 committed Aug 29, 2024
1 parent d045795 commit af27beb
Show file tree
Hide file tree
Showing 3 changed files with 70 additions and 18 deletions.
80 changes: 65 additions & 15 deletions pyro/compressible/interface.py
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,8 @@ def states(idir, ng, dx, dt,


@njit(cache=True)
def artificial_viscosity(ng, dx, dy,
def artificial_viscosity(ng, dx, dy, Lx, Ly,
xmin, ymin, coord_type,
cvisc, u, v):
r"""
Compute the artificial viscosity. Here, we compute edge-centered
Expand Down Expand Up @@ -250,6 +251,10 @@ def artificial_viscosity(ng, dx, dy,
The number of ghost cells
dx, dy : float
Cell spacings
xmin, ymin : float
Min physical x, y boundary
Lx, Ly: ndarray
Cell size in x, y direction
cvisc : float
viscosity parameter
u, v : ndarray
Expand All @@ -265,6 +270,7 @@ def artificial_viscosity(ng, dx, dy,

avisco_x = np.zeros((qx, qy))
avisco_y = np.zeros((qx, qy))
divU = np.zeros((qx, qy))

nx = qx - 2 * ng
ny = qy - 2 * ng
Expand All @@ -273,25 +279,69 @@ def artificial_viscosity(ng, dx, dy,
jlo = ng
jhi = ng + ny

# Let's first compute divergence at the vertex
# First compute the left and right x-velocities by
# averaging over the y-interface.
# As well as the top and bottom y-velocities by
# averaging over the x-interface.
# Then a simple difference is done between the right and left,
# and top and bottom to get the divergence at the vertex.

for i in range(ilo - 1, ihi + 1):
for j in range(jlo - 1, jhi + 1):

# start by computing the divergence on the x-interface. The
# x-difference is simply the difference of the cell-centered
# x-velocities on either side of the x-interface. For the
# y-difference, first average the four cells to the node on
# each end of the edge, and: difference these to find the
# edge centered y difference.
divU_x = (u[i, j] - u[i - 1, j]) / dx + \
0.25 * (v[i, j + 1] + v[i - 1, j + 1] -
v[i, j - 1] - v[i - 1, j - 1]) / dy
# For Cartesian2d:
if coord_type == 0:
# Find the average right and left u velocity
ur = 0.5 * (u[i, j] + u[i, j - 1])
ul = 0.5 * (u[i - 1, j] + u[i - 1, j - 1])

# Find the average top and bottom v velocity
vt = 0.5 * (v[i, j] + v[i - 1, j])
vb = 0.5 * (v[i, j - 1] + v[i - 1, j - 1])

# Finite difference to get ux and vy
ux = (ur - ul) / dx
vy = (vt - vb) / dy

# Find div(U)_{i-1/2, j-1/2}
divU[i, j] = ux + vy

# For SphericalPolar:
else:
# cell-centered r-coord of right, left cell and face-centered r
rr = (i + 1 - ng) * dx + xmin
rl = (i - ng) * dx + xmin
rc = 0.5 * (rr + rl)

# cell-centered sin(theta) of top, bot cell and face-centered
sint = np.sin((j + 1 - ng) * dy + ymin)
sinb = np.sin((j - ng) * dy + ymin)
sinc = np.sin((j + 0.5 - ng) * dy + ymin)

# Find the average right and left u velocity
ur = 0.5 * (u[i, j] + u[i, j - 1])
ul = 0.5 * (u[i - 1, j] + u[i - 1, j - 1])

# Find the average top and bottom v velocity
vt = 0.5 * (v[i, j] + v[i - 1, j])
vb = 0.5 * (v[i, j - 1] + v[i - 1, j - 1])

# Finite difference to get ux and vy
ux = (ur*rr - ul*rl) / (rc * dx)
vy = (sint*vt - sinb*vb) / (rc * sinc * dy)

# Find div(U)_{i-1/2, j-1/2}
divU[i, j] = ux + vy

avisco_x[i, j] = cvisc * max(-divU_x * dx, 0.0)
# Compute divergence at the face by averaging over divergence at vertex
for i in range(ilo, ihi):
for j in range(jlo, jhi):

# now the y-interface value
divU_y = 0.25 * (u[i + 1, j] + u[i + 1, j - 1] - u[i - 1, j] - u[i - 1, j - 1]) / dx + \
(v[i, j] - v[i, j - 1]) / dy
divU_x = 0.5 * (divU[i, j] + divU[i, j + 1])
divU_y = 0.5 * (divU[i, j] + divU[i + 1, j])

avisco_y[i, j] = cvisc * max(-divU_y * dy, 0.0)
avisco_x[i, j] = cvisc * max(-divU_x * Lx[i, j], 0.0)
avisco_y[i, j] = cvisc * max(-divU_y * Ly[i, j], 0.0)

return avisco_x, avisco_y
3 changes: 2 additions & 1 deletion pyro/compressible/unsplit_fluxes.py
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,8 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt):
# =========================================================================
cvisc = rp.get_param("compressible.cvisc")

_ax, _ay = ifc.artificial_viscosity(myg.ng, myg.dx, myg.dy,
_ax, _ay = ifc.artificial_viscosity(myg.ng, myg.dx, myg.dy, myg.Lx, myg.Ly,
myg.xmin, myg.ymin, myg.coord_type,
cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng))

avisco_x = ai.ArrayIndexer(d=_ax, grid=myg)
Expand Down
5 changes: 3 additions & 2 deletions pyro/compressible_rk/fluxes.py
Original file line number Diff line number Diff line change
Expand Up @@ -194,8 +194,9 @@ def fluxes(my_data, rp, ivars, solid, tc):
# =========================================================================
cvisc = rp.get_param("compressible.cvisc")

_ax, _ay = interface.artificial_viscosity(myg.ng, myg.dx, myg.dy,
cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng))
_ax, _ay = interface.artificial_viscosity(myg.ng, myg.dx, myg.dy, myg.Lx, myg.Ly,
myg.xmin, myg.ymin, myg.coord_type,
cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng))

avisco_x = ai.ArrayIndexer(d=_ax, grid=myg)
avisco_y = ai.ArrayIndexer(d=_ay, grid=myg)
Expand Down

0 comments on commit af27beb

Please sign in to comment.