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

update artificial viscosity to be in sync with castro #211

Merged
merged 2 commits into from
Aug 29, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
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
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
Loading