From 49cdba2fbc5fca4a44cfd25248d565f3cae6005b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 11:56:25 -0500 Subject: [PATCH 1/7] move the 4th order reconstruction into mesh/ this is consistent with the other reconstruction routines --- pyro/advection_fv4/fluxes.py | 6 +++--- pyro/compressible_fv4/fluxes.py | 5 ++--- pyro/{advection_fv4/interface.py => mesh/fourth_order.py} | 4 +++- 3 files changed, 8 insertions(+), 7 deletions(-) rename pyro/{advection_fv4/interface.py => mesh/fourth_order.py} (98%) diff --git a/pyro/advection_fv4/fluxes.py b/pyro/advection_fv4/fluxes.py index 8b902cf1f..a23b66dce 100644 --- a/pyro/advection_fv4/fluxes.py +++ b/pyro/advection_fv4/fluxes.py @@ -1,5 +1,5 @@ import pyro.mesh.array_indexer as ai -from pyro.advection_fv4 import interface +from pyro.mesh import fourth_order def fluxes(my_data, rp): @@ -72,13 +72,13 @@ def fluxes(my_data, rp): 1./12.*(a.jp(-2, buf=1) + a.jp(1, buf=1)) else: - a_l, a_r = interface.states(a, myg.ng, 1) + a_l, a_r = fourth_order.states(a, myg.ng, 1) if u > 0: a_x = ai.ArrayIndexer(d=a_l, grid=myg) else: a_x = ai.ArrayIndexer(d=a_r, grid=myg) - a_l, a_r = interface.states(a, myg.ng, 2) + a_l, a_r = fourth_order.states(a, myg.ng, 2) if v > 0: a_y = ai.ArrayIndexer(d=a_l, grid=myg) else: diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index 25a1439fa..6b3d67528 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -4,9 +4,8 @@ import pyro.compressible as comp import pyro.mesh.array_indexer as ai -from pyro.advection_fv4 import interface from pyro.compressible import riemann -from pyro.mesh import reconstruction +from pyro.mesh import fourth_order, reconstruction def flux_cons(ivars, idir, gamma, q): @@ -104,7 +103,7 @@ def fluxes(myd, rp, ivars): else: for n in range(ivars.nq): - q_l[:, :, n], q_r[:, :, n] = interface.states(q_avg[:, :, n], myg.ng, idir) + q_l[:, :, n], q_r[:, :, n] = fourth_order.states(q_avg[:, :, n], myg.ng, idir) # apply flattening for n in range(ivars.nq): diff --git a/pyro/advection_fv4/interface.py b/pyro/mesh/fourth_order.py similarity index 98% rename from pyro/advection_fv4/interface.py rename to pyro/mesh/fourth_order.py index c9f355077..feaff0838 100644 --- a/pyro/advection_fv4/interface.py +++ b/pyro/mesh/fourth_order.py @@ -1,3 +1,5 @@ +"""Reconstruction routines for 4th order finite-volume methods""" + import numpy as np from numba import njit @@ -84,7 +86,7 @@ def states(a, ng, idir): # this lives on the interface d3a[i, j] = d2ac[i, j] - d2ac[i - 1, j] - # this is a look over cell centers, affecting + # this is a loop over cell centers, affecting # i-1/2,R and i+1/2,L for i in range(ilo - 1, ihi + 1): for j in range(jlo - 1, jhi + 1): From ce67a64016499428709e493996b8397951b7268e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 12:22:38 -0500 Subject: [PATCH 2/7] a few more fixes --- examples/mesh/bc_demo.py | 3 ++- pyro/mesh/__init__.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/mesh/bc_demo.py b/examples/mesh/bc_demo.py index 7ce80fe72..6c5dc5524 100644 --- a/examples/mesh/bc_demo.py +++ b/examples/mesh/bc_demo.py @@ -1,9 +1,10 @@ # test the boundary fill routine +import numpy as np + import pyro.mesh.boundary as bnd import pyro.mesh.patch as patch -import numpy as np def doit(): diff --git a/pyro/mesh/__init__.py b/pyro/mesh/__init__.py index 5f6f0e1ae..2a8a3c0e8 100644 --- a/pyro/mesh/__init__.py +++ b/pyro/mesh/__init__.py @@ -3,7 +3,7 @@ necessary to work with finite-volume data. """ -__all__ = ['patch', 'integration', 'reconstruction'] +__all__ = ['patch', 'fourth_order', 'integration', 'reconstruction'] from .array_indexer import ArrayIndexer, ArrayIndexerFC from .boundary import BC, bc_is_solid, define_bc From e061603ca40c28788e354239e7385164ce1acf4b Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 12:56:02 -0500 Subject: [PATCH 3/7] implement one-sided stencils for 4th order --- pyro/mesh/fourth_order.py | 129 +++++++++++++++++--------------------- 1 file changed, 57 insertions(+), 72 deletions(-) diff --git a/pyro/mesh/fourth_order.py b/pyro/mesh/fourth_order.py index feaff0838..799e89dd4 100644 --- a/pyro/mesh/fourth_order.py +++ b/pyro/mesh/fourth_order.py @@ -5,7 +5,7 @@ @njit(cache=True) -def states(a, ng, idir): +def states(a, ng, idir, *, lo_bc_symmetry=False, hi_bc_symmetry=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 +61,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] @@ -150,8 +176,34 @@ def states(a, ng, idir): 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] @@ -234,70 +286,3 @@ def states(a, ng, idir): 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] - - return al, ar From 0a47c60d74ff32f4e9a14080885c5ce0de4792bf Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 14:23:09 -0500 Subject: [PATCH 4/7] numba doesn't like * in arg lists --- pyro/mesh/fourth_order.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyro/mesh/fourth_order.py b/pyro/mesh/fourth_order.py index 799e89dd4..8369ff6c7 100644 --- a/pyro/mesh/fourth_order.py +++ b/pyro/mesh/fourth_order.py @@ -5,7 +5,7 @@ @njit(cache=True) -def states(a, ng, idir, *, lo_bc_symmetry=False, hi_bc_symmetry=False): +def states(a, ng, idir, lo_bc_symmetry=False, hi_bc_symmetry=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. From f49880033d3e016056c5306634f116286324095e Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 14:46:23 -0500 Subject: [PATCH 5/7] pass in solid to fourth order --- pyro/compressible_fv4/fluxes.py | 17 +++++++++++++---- pyro/compressible_fv4/simulation.py | 2 +- 2 files changed, 14 insertions(+), 5 deletions(-) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index 6b3d67528..4fda568a2 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 @@ -88,6 +88,13 @@ def fluxes(myd, rp, ivars): q_l = myg.scratch_array(nvar=ivars.nq) q_r = myg.scratch_array(nvar=ivars.nq) + 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): @@ -103,7 +110,9 @@ 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) + 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) # apply flattening for n in range(ivars.nq): @@ -121,8 +130,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 df662ae4a..930ee3ab6 100644 --- a/pyro/compressible_fv4/simulation.py +++ b/pyro/compressible_fv4/simulation.py @@ -41,7 +41,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)[:, :] = \ From 6f32c16b5d3cad13fc2c58e347e804864441236d Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 14:51:51 -0500 Subject: [PATCH 6/7] silence pylint --- pyro/compressible_fv4/fluxes.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index 4fda568a2..e607b2acb 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -88,6 +88,8 @@ def fluxes(myd, rp, ivars, solid): 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 From d77de091b5572575e4dfff3ba6aaedc179610750 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Fri, 22 Nov 2024 18:13:36 -0500 Subject: [PATCH 7/7] more BCs --- pyro/compressible_fv4/fluxes.py | 11 ++++++++++- pyro/mesh/fourth_order.py | 28 +++++++++++++++++++++++++++- 2 files changed, 37 insertions(+), 2 deletions(-) diff --git a/pyro/compressible_fv4/fluxes.py b/pyro/compressible_fv4/fluxes.py index e607b2acb..af26062c6 100644 --- a/pyro/compressible_fv4/fluxes.py +++ b/pyro/compressible_fv4/fluxes.py @@ -112,9 +112,18 @@ def fluxes(myd, rp, ivars, solid): else: for n in range(ivars.nq): + + # 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) + hi_bc_symmetry=hi_bc_symmetry, + is_normal_vel=is_normal_vel) # apply flattening for n in range(ivars.nq): diff --git a/pyro/mesh/fourth_order.py b/pyro/mesh/fourth_order.py index 8369ff6c7..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, lo_bc_symmetry=False, hi_bc_symmetry=False): +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. @@ -170,6 +172,18 @@ def states(a, ng, idir, lo_bc_symmetry=False, hi_bc_symmetry=False): 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): @@ -285,4 +299,16 @@ def states(a, ng, idir, lo_bc_symmetry=False, hi_bc_symmetry=False): if abs(dafp[i, j]) >= 2.0 * abs(dafm[i, j]): al[i, j + 1] = a[i, j] + 2.0 * dafm[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