simplify incompressible solver based on burgers solver (#169)

based on the previous burgers pr.

Simplifies the get_interface_states function in incomp_interface.py from burgers/burgers_interface.py

Also, let incompressible inherit from burgers to eliminate the method_compute_timestep function.

Author: zhichen3

pyro/incompressible/incomp_interface.py

@@ -32,12 +32,29 @@ def mac_vels(grid,  dt,
 
     # get the full u and v left and right states (including transverse
     # terms) on both the x- and y-interfaces
-    # pylint: disable-next=unused-variable
-    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(grid, dt,
-                                                                          u, v,
-                                                                          ldelta_ux, ldelta_vx,
-                                                                          ldelta_uy, ldelta_vy,
-                                                                          gradp_x, gradp_y)
+
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                burgers_interface.get_interface_states(grid, dt,
+                                                       u, v,
+                                                       ldelta_ux, ldelta_vx,
+                                                       ldelta_uy, ldelta_vy)
+
+    # apply transverse terms on both x- and y- interfaces
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                burgers_interface.apply_transverse_corrections(grid, dt,
+                                                               u_xl, u_xr,
+                                                               u_yl, u_yr,
+                                                               v_xl, v_xr,
+                                                               v_yl, v_yr)
+
+    # apply pressure gradient correction terms to the interface state
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                apply_gradp_corrections(dt,
+                                        u_xl, u_xr,
+                                        u_yl, u_yr,
+                                        v_xl, v_xr,
+                                        v_yl, v_yr,
+                                        gradp_x, gradp_y)
 
     # Riemann problem -- this follows Burger's equation.  We don't use
     # any input velocity for the upwinding.  Also, we only care about
@@ -82,13 +99,29 @@ def states(grid, dt,
         x and y velocities predicted to the interfaces
     """
 
-    # get the full u and v left and right states (including transverse
-    # terms) on both the x- and y-interfaces
-    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = get_interface_states(grid, dt,
-                                                                          u, v,
-                                                                          ldelta_ux, ldelta_vx,
-                                                                          ldelta_uy, ldelta_vy,
-                                                                          gradp_x, gradp_y)
+    # get the full u and v left and right states without transverse terms and gradp
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                burgers_interface.get_interface_states(grid, dt,
+                                                       u, v,
+                                                       ldelta_ux, ldelta_vx,
+                                                       ldelta_uy, ldelta_vy)
+
+    # apply transverse terms on both x- and y- interfaces
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                burgers_interface.apply_transverse_corrections(grid, dt,
+                                                               u_xl, u_xr,
+                                                               u_yl, u_yr,
+                                                               v_xl, v_xr,
+                                                               v_yl, v_yr)
+
+    # apply pressure gradient correction terms to the interface state
+    u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr = \
+                apply_gradp_corrections(dt,
+                                        u_xl, u_xr,
+                                        u_yl, u_yr,
+                                        v_xl, v_xr,
+                                        v_yl, v_yr,
+                                        gradp_x, gradp_y)
 
     # upwind using the MAC velocity to determine which state exists on
     # the interface
@@ -100,128 +133,51 @@ def states(grid, dt,
     return u_xint, v_xint, u_yint, v_yint
 
 
-def get_interface_states(grid, dt,
-                         u, v,
-                         ldelta_ux, ldelta_vx,
-                         ldelta_uy, ldelta_vy,
-                         gradp_x, gradp_y):
+def apply_gradp_corrections(dt,
+                            u_xl, u_xr,
+                            u_yl, u_yr,
+                            v_xl, v_xr,
+                            v_yl, v_yr,
+                            gradp_x, gradp_y):
     r"""
-    Compute the unsplit predictions of u and v on both the x- and
-    y-interfaces.  This includes the transverse terms.
-
     Parameters
     ----------
     grid : Grid2d
         The grid object
     dt : float
         The timestep
-    u, v : ndarray
-        x-velocity and y-velocity
-    ldelta_ux, ldelta_uy: ndarray
-        Limited slopes of the x-velocity in the x and y directions
-    ldelta_vx, ldelta_vy: ndarray
-        Limited slopes of the y-velocity in the x and y directions
-    gradp_x, gradp_y : ndarray
-        Pressure gradients in the x and y directions
-
+    u_xl, u_xr : ndarray ndarray
+        left and right states of x-velocity in x-interface.
+    u_yl, u_yr : ndarray ndarray
+        left and right states of x-velocity in y-interface.
+    v_xl, v_xr : ndarray ndarray
+        left and right states of y-velocity in x-interface.
+    v_yl, u_yr : ndarray ndarray
+        left and right states of y-velocity in y-interface.
     Returns
     -------
     out : ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray
-        unsplit predictions of u and v on both the x- and
-        y-interfaces
+        correct the interface states of the left and right states of u and v on
+        both the x- and y-interfaces interface states with the pressure gradient
+        terms.
     """
 
-    u_xl = grid.scratch_array()
-    u_xr = grid.scratch_array()
-    u_yl = grid.scratch_array()
-    u_yr = grid.scratch_array()
-
-    v_xl = grid.scratch_array()
-    v_xr = grid.scratch_array()
-    v_yl = grid.scratch_array()
-    v_yr = grid.scratch_array()
-
-    # first predict u and v to both interfaces, considering only the normal
-    # part of the predictor.  These are the 'hat' states.
-
-    dtdx = dt / grid.dx
-    dtdy = dt / grid.dy
-
-    # u on x-edges
-    u_xl.ip(1, buf=2)[:, :] = u.v(buf=2) + \
-        0.5 * (1.0 - dtdx * u.v(buf=2)) * ldelta_ux.v(buf=2)
-    u_xr.v(buf=2)[:, :] = u.v(buf=2) - \
-        0.5 * (1.0 + dtdx * u.v(buf=2)) * ldelta_ux.v(buf=2)
-
-    # v on x-edges
-    v_xl.ip(1, buf=2)[:, :] = v.v(buf=2) + \
-        0.5 * (1.0 - dtdx * u.v(buf=2)) * ldelta_vx.v(buf=2)
-    v_xr.v(buf=2)[:, :] = v.v(buf=2) - \
-        0.5 * (1.0 + dtdx * u.v(buf=2)) * ldelta_vx.v(buf=2)
-
-    # u on y-edges
-    u_yl.jp(1, buf=2)[:, :] = u.v(buf=2) + \
-        0.5 * (1.0 - dtdy * v.v(buf=2)) * ldelta_uy.v(buf=2)
-    u_yr.v(buf=2)[:, :] = u.v(buf=2) - \
-        0.5 * (1.0 + dtdy * v.v(buf=2)) * ldelta_uy.v(buf=2)
-
-    # v on y-edges
-    v_yl.jp(1, buf=2)[:, :] = v.v(buf=2) + \
-        0.5 * (1.0 - dtdy * v.v(buf=2)) * ldelta_vy.v(buf=2)
-    v_yr.v(buf=2)[:, :] = v.v(buf=2) - \
-        0.5 * (1.0 + dtdy * v.v(buf=2)) * ldelta_vy.v(buf=2)
-
-    # now get the normal advective velocities on the interfaces by solving
-    # the Riemann problem.
-    uhat_adv = burgers_interface.riemann(grid, u_xl, u_xr)
-    vhat_adv = burgers_interface.riemann(grid, v_yl, v_yr)
-
-    # now that we have the advective velocities, upwind the left and right
-    # states using the appropriate advective velocity.
-
-    # on the x-interfaces, we upwind based on uhat_adv
-    u_xint = burgers_interface.upwind(grid, u_xl, u_xr, uhat_adv)
-    v_xint = burgers_interface.upwind(grid, v_xl, v_xr, uhat_adv)
-
-    # on the y-interfaces, we upwind based on vhat_adv
-    u_yint = burgers_interface.upwind(grid, u_yl, u_yr, vhat_adv)
-    v_yint = burgers_interface.upwind(grid, v_yl, v_yr, vhat_adv)
-
-    # at this point, these states are the `hat' states -- they only
-    # considered the normal to the interface portion of the predictor.
-
-    ubar = grid.scratch_array()
-    vbar = grid.scratch_array()
-
-    ubar.v(buf=2)[:, :] = 0.5 * (uhat_adv.v(buf=2) + uhat_adv.ip(1, buf=2))
-    vbar.v(buf=2)[:, :] = 0.5 * (vhat_adv.v(buf=2) + vhat_adv.jp(1, buf=2))
-
-    # v du/dy is the transerse term for the u states on x-interfaces
-    vu_y = grid.scratch_array()
-    vu_y.v(buf=2)[:, :] = vbar.v(buf=2) * (u_yint.jp(1, buf=2) - u_yint.v(buf=2))
-
-    u_xl.ip(1, buf=2)[:, :] += -0.5 * dtdy * vu_y.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2)
-    u_xr.v(buf=2)[:, :] += -0.5 * dtdy * vu_y.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2)
-
-    # v dv/dy is the transverse term for the v states on x-interfaces
-    vv_y = grid.scratch_array()
-    vv_y.v(buf=2)[:, :] = vbar.v(buf=2) * (v_yint.jp(1, buf=2) - v_yint.v(buf=2))
-
-    v_xl.ip(1, buf=2)[:, :] += -0.5 * dtdy * vv_y.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2)
-    v_xr.v(buf=2)[:, :] += -0.5 * dtdy * vv_y.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2)
+    # Apply pressure gradient correction terms
 
-    # u dv/dx is the transverse term for the v states on y-interfaces
-    uv_x = grid.scratch_array()
-    uv_x.v(buf=2)[:, :] = ubar.v(buf=2) * (v_xint.ip(1, buf=2) - v_xint.v(buf=2))
+    # transverse term for the u states on x-interfaces
+    u_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
+    u_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
 
-    v_yl.jp(1, buf=2)[:, :] += -0.5 * dtdx * uv_x.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2)
-    v_yr.v(buf=2)[:, :] += -0.5 * dtdx * uv_x.v(buf=2) - 0.5 * dt * gradp_y.v(buf=2)
+    # transverse term for the v states on x-interfaces
+    v_xl.ip(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
+    v_xr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
 
-    # u du/dx is the transverse term for the u states on y-interfaces
-    uu_x = grid.scratch_array()
-    uu_x.v(buf=2)[:, :] = ubar.v(buf=2) * (u_xint.ip(1, buf=2) - u_xint.v(buf=2))
+    # transverse term for the v states on y-interfaces
+    v_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
+    v_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_y.v(buf=2)
 
-    u_yl.jp(1, buf=2)[:, :] += -0.5 * dtdx * uu_x.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2)
-    u_yr.v(buf=2)[:, :] += -0.5 * dtdx * uu_x.v(buf=2) - 0.5 * dt * gradp_x.v(buf=2)
+    # transverse term for the u states on y-interfaces
+    u_yl.jp(1, buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
+    u_yr.v(buf=2)[:, :] += - 0.5 * dt * gradp_x.v(buf=2)
 
     return u_xl, u_xr, u_yl, u_yr, v_xl, v_xr, v_yl, v_yr

pyro/incompressible/simulation.py

@@ -4,14 +4,15 @@
 import numpy as np
 
 import pyro.mesh.array_indexer as ai
+from pyro.burgers import Simulation as burgers_simulation
 from pyro.incompressible import incomp_interface
 from pyro.mesh import patch, reconstruction
 from pyro.multigrid import MG
 from pyro.particles import particles
-from pyro.simulation_null import NullSimulation, bc_setup, grid_setup
+from pyro.simulation_null import bc_setup, grid_setup
 
 
-class Simulation(NullSimulation):
+class Simulation(burgers_simulation):
 
     def initialize(self):
         """
@@ -51,27 +52,6 @@ def initialize(self):
         problem = importlib.import_module(f"pyro.incompressible.problems.{self.problem_name}")
         problem.init_data(self.cc_data, self.rp)
 
-    def method_compute_timestep(self):
-        """
-        The timestep() function computes the advective timestep
-        (CFL) constraint.  The CFL constraint says that information
-        cannot propagate further than one zone per timestep.
-
-        We use the driver.cfl parameter to control what fraction of the CFL
-        step we actually take.
-        """
-
-        cfl = self.rp.get_param("driver.cfl")
-
-        u = self.cc_data.get_var("x-velocity")
-        v = self.cc_data.get_var("y-velocity")
-
-        # the timestep is min(dx/|u|, dy|v|)
-        xtmp = self.cc_data.grid.dx/(abs(u))
-        ytmp = self.cc_data.grid.dy/(abs(v))
-
-        self.dt = cfl*float(min(xtmp.min(), ytmp.min()))
-
     def preevolve(self):
         """
         preevolve is called before we being the timestepping loop.  For