diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index d59352544..cbd55a903 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -37,7 +37,7 @@ def flux_cons(ivars, idir, gamma, q): return flux -def fluxes(myd, rp, ivars): +def fluxes(myd, rp, ivars, solid): alpha = 0.3 beta = 0.3 @@ -112,6 +112,15 @@ def fluxes(myd, rp, ivars): q_l = myg.scratch_array(nvar=ivars.nq) q_r = myg.scratch_array(nvar=ivars.nq) + lo_bc_symmetry = False + hi_bc_symmetry = False + if idir == 1: + lo_bc_symmetry = solid.xl + hi_bc_symmetry = solid.xr + elif idir == 2: + lo_bc_symmetry = solid.yl + hi_bc_symmetry = solid.yr + if nolimit: for n in range(ivars.nq): @@ -127,7 +136,18 @@ def fluxes(myd, rp, ivars): else: for n in range(ivars.nq): - q_l[:, :, n], q_r[:, :, n] = fourth_order.states(q_avg[:, :, n], myg.ng, idir) + + # for symmetry BCs, we want to odd-reflect the interface state on velocity + is_normal_vel = False + if idir == 1 and n == ivars.iu: + is_normal_vel = True + elif idir == 2 and n == ivars.iv: + is_normal_vel = True + + q_l[:, :, n], q_r[:, :, n] = fourth_order.states(q_avg[:, :, n], myg.ng, idir, + lo_bc_symmetry=lo_bc_symmetry, + hi_bc_symmetry=hi_bc_symmetry, + is_normal_vel=is_normal_vel) # apply flattening for n in range(ivars.nq): @@ -145,8 +165,8 @@ def fluxes(myd, rp, ivars): # solve the Riemann problem to find the face-average q _q = riemann.riemann_prim(idir, myg.ng, ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix, ivars.naux, - 0, 0, - gamma, q_l, q_r) + 0, 0, + gamma, q_l, q_r) q_int_avg = ai.ArrayIndexer(_q, grid=myg) diff --git a/pyro/compressible_fv4/simulation.py b/pyro/compressible_fv4/simulation.py index c008c6aa9..ba9a42003 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -45,7 +45,7 @@ def substep(self, myd): k = myg.scratch_array(nvar=self.ivars.nvar) - flux_x, flux_y = flx.fluxes(myd, self.rp, self.ivars) + flux_x, flux_y = flx.fluxes(myd, self.rp, self.ivars, self.solid) for n in range(self.ivars.nvar): k.v(n=n)[:, :] = \ diff --git a/pyro/mesh/fourth_order.py b/pyro/mesh/fourth_order.py index feaff0838..8cd1058cf 100644 --- a/pyro/mesh/fourth_order.py +++ b/pyro/mesh/fourth_order.py @@ -5,7 +5,9 @@ @njit(cache=True) -def states(a, ng, idir): +def states(a, ng, idir, + lo_bc_symmetry=False, hi_bc_symmetry=False, + is_normal_vel=False): r""" Predict the cell-centered state to the edges in one-dimension using the reconstructed, limited slopes. We use a fourth-order Godunov method. @@ -61,8 +63,34 @@ def states(a, ng, idir): for j in range(jlo - 1, jhi + 1): # interpolate to the edges - a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \ - (1.0 / 12.0) * (a[i - 2, j] + a[i + 1, j]) + + if i == ilo+1 and lo_bc_symmetry: + # use a stencil for the interface that is one zone + # from the left physical boundary, MC Eq. 22 + a_int[i, j] = (1.0 / 12.0) * (3.0 * a[i-1, j] + 13.0 * a[i, j] - + 5.0 * a[i+1, j] + a[i+2, j]) + + elif i == ilo and lo_bc_symmetry: + # use a stencil for when the interface is on the + # left physical boundary MC Eq. 21 + a_int[i, j] = (1.0 / 12.0) * (25.0 * a[i, j] - 23.0 * a[i+1, j] + + 13.0 * a[i+2, j] - 3.0 * a[i+3, j]) + + elif i == ihi and hi_bc_symmetry: + # use a stencil for the interface that is one zone + # from the right physical boundary, MC Eq. 22 + a_int[i, j] = (1.0 / 12.0) * (3.0 * a[i, j] + 13.0 * a[i-1, j] - + 5.0 * a[i-2, j] + a[i-3, j]) + + elif i == ihi+1 and hi_bc_symmetry: + # use a stencil for when the interface is on the + # right physical boundary MC Eq. 21 + a_int[i, j] = (1.0 / 12.0) * (25.0 * a[i-1, j] - 23.0 * a[i-2, j] + + 13.0 * a[i-3, j] - 3.0 * a[i-4, j]) + + else: + a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \ + (1.0 / 12.0) * (a[i - 2, j] + a[i + 1, j]) al[i, j] = a_int[i, j] ar[i, j] = a_int[i, j] @@ -144,14 +172,52 @@ def states(a, ng, idir): if abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j]): al[i + 1, j] = a[i, j] + 2.0 * dafm[i, j] + if lo_bc_symmetry: + if is_normal_vel: + al[ilo, :] = -ar[ilo, :] + else: + al[ilo, :] = ar[ilo, :] + + if hi_bc_symmetry: + if is_normal_vel: + ar[ihi+1, :] = -al[ihi+1, :] + else: + ar[ihi+1, :] = al[ihi+1, :] + elif idir == 2: for i in range(ilo - 1, ihi + 1): for j in range(jlo - 2, jhi + 3): # interpolate to the edges - a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \ - (1.0 / 12.0) * (a[i, j - 2] + a[i, j + 1]) + + if j == jlo+1 and lo_bc_symmetry: + # use a stencil for the interface that is one zone + # from the left physical boundary, MC Eq. 22 + a_int[i, j] = (1.0 / 12.0) * (3.0 * a[i, j-1] + 13.0 * a[i, j] - + 5.0 * a[i, j+1] + a[i, j+2]) + + elif j == jlo and lo_bc_symmetry: + # use a stencil for when the interface is on the + # left physical boundary MC Eq. 21 + a_int[i, j] = (1.0 / 12.0) * (25.0 * a[i, j] - 23.0 * a[i, j+1] + + 13.0 * a[i, j+2] - 3.0 * a[i, j+3]) + + elif j == jhi and hi_bc_symmetry: + # use a stencil for the interface that is one zone + # from the right physical boundary, MC Eq. 22 + a_int[i, j] = (1.0 / 12.0) * (3.0 * a[i, j] + 13.0 * a[i, j-1] - + 5.0 * a[i, j-2] + a[i, j-3]) + + elif j == jhi+1 and hi_bc_symmetry: + # use a stencil for when the interface is on the + # right physical boundary MC Eq. 21 + a_int[i, j] = (1.0 / 12.0) * (25.0 * a[i, j-1] - 23.0 * a[i, j-2] + + 13.0 * a[i, j-3] - 3.0 * a[i, j-4]) + + else: + a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \ + (1.0 / 12.0) * (a[i, j - 2] + a[i, j + 1]) al[i, j] = a_int[i, j] ar[i, j] = a_int[i, j] @@ -233,71 +299,16 @@ def states(a, ng, idir): if abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j]): al[i, j + 1] = a[i, j] + 2.0 * dafm[i, j] - return al, ar - - -@njit(cache=True) -def states_nolimit(a, qx, qy, ng, idir): - r""" - Predict the cell-centered state to the edges in one-dimension using the - reconstructed slopes (and without slope limiting). We use a fourth-order - Godunov method. - - Our convention here is that: - - ``al[i,j]`` will be :math:`al_{i-1/2,j}`, - - ``al[i+1,j]`` will be :math:`al_{i+1/2,j}`. - - Parameters - ---------- - a : ndarray - The cell-centered state. - ng : int - The number of ghost cells - idir : int - Are we predicting to the edges in the x-direction (1) or y-direction (2)? - - Returns - ------- - out : ndarray, ndarray - The state predicted to the left and right edges. - """ - - a_int = np.zeros((qx, qy)) - al = np.zeros((qx, qy)) - ar = np.zeros((qx, qy)) - - nx = qx - 2 * ng - ny = qy - 2 * ng - ilo = ng - ihi = ng + nx - jlo = ng - jhi = ng + ny - - # we need interface values on all faces of the domain - if idir == 1: - - for i in range(ilo - 2, ihi + 3): - for j in range(jlo - 1, jhi + 1): - - # interpolate to the edges - a_int[i, j] = (7.0 / 12.0) * (a[i - 1, j] + a[i, j]) - \ - (1.0 / 12.0) * (a[i - 2, j] + a[i + 1, j]) - - al[i, j] = a_int[i, j] - ar[i, j] = a_int[i, j] - - elif idir == 2: - - for i in range(ilo - 1, ihi + 1): - for j in range(jlo - 2, jhi + 3): - - # interpolate to the edges - a_int[i, j] = (7.0 / 12.0) * (a[i, j - 1] + a[i, j]) - \ - (1.0 / 12.0) * (a[i, j - 2] + a[i, j + 1]) - - al[i, j] = a_int[i, j] - ar[i, j] = a_int[i, j] + if lo_bc_symmetry: + if is_normal_vel: + al[:, jlo] = -ar[:, jlo] + else: + al[:, jlo] = ar[:, jlo] + + if hi_bc_symmetry: + if is_normal_vel: + ar[:, jhi+1] = -al[:, jhi+1] + else: + ar[:, jhi+1] = al[:, jhi+1] return al, ar