From af27beb60bd21d15fb55356376ab8fea8e7c3614 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 29 Aug 2024 14:29:45 -0400 Subject: [PATCH] update artificial viscosity to be in sync with castro --- pyro/compressible/interface.py | 80 +++++++++++++++++++++++------ pyro/compressible/unsplit_fluxes.py | 3 +- pyro/compressible_rk/fluxes.py | 5 +- 3 files changed, 70 insertions(+), 18 deletions(-) diff --git a/pyro/compressible/interface.py b/pyro/compressible/interface.py index c2aac43d6..b02663750 100644 --- a/pyro/compressible/interface.py +++ b/pyro/compressible/interface.py @@ -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 @@ -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 @@ -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 @@ -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 diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index e17b10693..065e8c08e 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -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) diff --git a/pyro/compressible_rk/fluxes.py b/pyro/compressible_rk/fluxes.py index 48a71acd6..5a6c46ef1 100644 --- a/pyro/compressible_rk/fluxes.py +++ b/pyro/compressible_rk/fluxes.py @@ -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)