implement well-balancing with characteristic tracing (#30)

Author: zingale

ppmpy/euler.py

@@ -82,7 +82,7 @@ def __init__(self, nx, C, *,
                  gamma=1.4,
                  init_cond=None, grav_func=None,
                  params=None,
-                 use_hse_reconstruction=False,
+                 use_hse_reconstruction=False, hse_as_perturbation=False,
                  use_limiting=True, use_flattening=True):
 
         self.grid = FVGrid(nx, ng=4)
@@ -138,6 +138,8 @@ def __init__(self, nx, C, *,
         self.U_init = self.U.copy()
 
         self.use_hse_reconstruction = use_hse_reconstruction
+        self.hse_as_perturbation = hse_as_perturbation
+
         self.use_flattening = use_flattening
         self.use_limiting = use_limiting
 
@@ -198,7 +200,7 @@ def construct_parabola(self):
         for ivar in range(self.v.nvar):
             if ivar == self.v.qp and self.use_hse_reconstruction:
                 self.q_parabola.append(HSEPPMInterpolant(self.grid, q[:, ivar], q[:, self.v.qrho], g,
-                                                      limit=self.use_limiting, chi_flat=chi))
+                                                      limit=self.use_limiting, chi_flat=chi, leave_as_perturbation=self.hse_as_perturbation))
             else:
                 self.q_parabola.append(PPMInterpolant(self.grid, q[:, ivar],
                                                       limit=self.use_limiting, chi_flat=chi))
@@ -272,13 +274,23 @@ def interface_states(self):
             q_ref_m = Im[i, 0, :]
 
             # build eigensystem
-            ev, lvec, rvec = eigen(q_ref_m[self.v.qrho], q_ref_m[self.v.qu], q_ref_m[self.v.qp], self.gamma)
+            if self.use_hse_reconstruction and self.hse_as_perturbation:
+                # we need to make the reference state a pressure again
+                ev, lvec, rvec = eigen(q_ref_m[self.v.qrho],
+                                       q_ref_m[self.v.qu],
+                                       self.q_parabola[self.v.qp].p_hse_m[i] + q_ref_m[self.v.qp],
+                                       self.gamma)
+            else:
+                ev, lvec, rvec = eigen(q_ref_m[self.v.qrho],
+                                       q_ref_m[self.v.qu],
+                                       q_ref_m[self.v.qp],
+                                       self.gamma)
 
             # loop over waves and compute l . (qref - I) for each wave
             beta_xm = np.zeros(3)
             for iwave in range(3):
                 dq = q_ref_m - Im[i, iwave, :]
-                if self.grav_func is not None:
+                if self.grav_func is not None and not self.hse_as_perturbation:
                     dq[self.v.qu] -= 0.5 * self.dt * Im_g[i, iwave]
                 beta_xm[iwave] = lvec[iwave, :] @ dq
 
@@ -289,19 +301,31 @@ def interface_states(self):
                 if ev[iwave] <= 0:
                     q_right[i, :] -= beta_xm[iwave] * rvec[iwave, :]
 
+            if self.use_hse_reconstruction and self.hse_as_perturbation:
+                q_right[i, self.v.qp] += self.q_parabola[self.v.qp].p_hse_m[i]
+
             # left state on interface i+1 -- this uses the "p" reconstructuion
 
             # reference state -- fastest wave moving to the right
             q_ref_p = Ip[i, 2, :]
 
             # build eigensystem
-            ev, lvec, rvec = eigen(q_ref_p[self.v.qrho], q_ref_p[self.v.qu], q_ref_p[self.v.qp], self.gamma)
+            if self.use_hse_reconstruction and self.hse_as_perturbation:
+                ev, lvec, rvec = eigen(q_ref_p[self.v.qrho],
+                                       q_ref_p[self.v.qu],
+                                       self.q_parabola[self.v.qp].p_hse_p[i] + q_ref_p[self.v.qp],
+                                       self.gamma)
+            else:
+                ev, lvec, rvec = eigen(q_ref_p[self.v.qrho],
+                                       q_ref_p[self.v.qu],
+                                       q_ref_p[self.v.qp],
+                                       self.gamma)
 
             # loop over waves and compute l . (qref - I) for each wave
             beta_xp = np.zeros(3)
             for iwave in range(3):
                 dq = q_ref_p - Ip[i, iwave, :]
-                if self.grav_func is not None:
+                if self.grav_func is not None and not self.hse_as_perturbation:
                     dq[self.v.qu] -= 0.5 * self.dt * Ip_g[i, iwave]
                 beta_xp[iwave] = lvec[iwave, :] @ dq
 
@@ -312,6 +336,9 @@ def interface_states(self):
                 if ev[iwave] >= 0:
                     q_left[i+1, :] -= beta_xp[iwave] * rvec[iwave, :]
 
+            if self.use_hse_reconstruction and self.hse_as_perturbation:
+                q_left[i+1, self.v.qp] += self.q_parabola[self.v.qp].p_hse_p[i]
+
         # enforce reflection
         for n in range(self.v.nvar):
             if self.bcs_left[n] == "reflect":

ppmpy/reconstruction.py

@@ -277,13 +277,18 @@ class HSEPPMInterpolant(PPMInterpolant):
         the flattening coefficient.
     """
 
-    def __init__(self, grid, p, rho, g, *, limit=True, chi_flat=None):
+    def __init__(self, grid, p, rho, g, *, limit=True, chi_flat=None, leave_as_perturbation=False):
 
         super().__init__(grid, p, limit=limit, chi_flat=chi_flat)
 
         self.rho = rho
         self.g = g
 
+        self.p_hse_m = None
+        self.p_hse_p = None
+
+        self.leave_as_perturbation = leave_as_perturbation
+
     def construct_parabola(self):
         """compute the coefficients of a parabolic interpolant for the
         pressure in each zone by subtracting off HSE.  This will give
@@ -377,10 +382,19 @@ def construct_parabola(self):
             self.am[:] = (1.0 - self.chi_flat[:]) * p_hse[:] + self.chi_flat[:] * self.am[:]
             self.ap[:] = (1.0 - self.chi_flat[:]) * p_hse[:] + self.chi_flat[:] * self.ap[:]
 
-        # finally, add back in the HSE correction
-        self.ap[:] += self.a[:] + 0.5 * self.grid.dx * self.rho[:] * self.g[:]
-        self.am[:] += self.a[:] - 0.5 * self.grid.dx * self.rho[:] * self.g[:]
+        self.p_hse_p = self.a[:] + 0.5 * self.grid.dx * self.rho[:] * self.g[:]
+        self.p_hse_m = self.a[:] - 0.5 * self.grid.dx * self.rho[:] * self.g[:]
 
-        self.a6 = 6.0 * self.a - 3.0 * (self.am + self.ap)
+        if not self.leave_as_perturbation:
+
+            # finally, add back in the HSE correction
+            self.ap[:] += self.p_hse_p[:]
+            self.am[:] += self.p_hse_m[:]
+
+            self.a6 = 6.0 * self.a - 3.0 * (self.am + self.ap)
+
+        else:
+            # the cell-center state here is zero, since we subtract off pressur
+            self.a6 = - 3.0 * (self.am + self.ap)
 
         self.initialized = True