more BCs

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):

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