simplified the initial conditions interface (#39)

Author: zingale

examples/euler.ipynb

@@ -190,32 +190,32 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 10,
    "id": "047dc02b-b5f1-45c0-a41a-af0251834cc7",
    "metadata": {},
    "outputs": [],
    "source": [
-    "def shock_tube(g, v, U, params):\n",
+    "def shock_tube(euler):\n",
     "\n",
-    "    gamma = params[\"gamma\"]\n",
+    "    gamma = euler.params[\"gamma\"]\n",
     "    \n",
-    "    rho_l = params[\"rho_l\"]\n",
-    "    u_l = params[\"u_l\"]\n",
-    "    p_l = params[\"p_l\"]\n",
-    "    rho_r = params[\"rho_r\"]\n",
-    "    u_r = params[\"u_r\"]\n",
-    "    p_r = params[\"p_r\"]\n",
+    "    rho_l = euler.params[\"rho_l\"]\n",
+    "    u_l = euler.params[\"u_l\"]\n",
+    "    p_l = euler.params[\"p_l\"]\n",
+    "    rho_r = euler.params[\"rho_r\"]\n",
+    "    u_r = euler.params[\"u_r\"]\n",
+    "    p_r = euler.params[\"p_r\"]\n",
     "\n",
-    "    idx_l = g.x < 0.5\n",
-    "    idx_r = g.x >= 0.5\n",
+    "    idx_l = euler.grid.x < 0.5\n",
+    "    idx_r = euler.grid.x >= 0.5\n",
     "\n",
-    "    U[idx_l, v.urho] = rho_l\n",
-    "    U[idx_l, v.umx] = rho_l * u_l\n",
-    "    U[idx_l, v.uener] = p_l/(gamma - 1.0) + 0.5 * rho_l * u_l**2\n",
+    "    euler.U[idx_l, euler.v.urho] = rho_l\n",
+    "    euler.U[idx_l, euler.v.umx] = rho_l * u_l\n",
+    "    euler.U[idx_l, euler.v.uener] = p_l/(gamma - 1.0) + 0.5 * rho_l * u_l**2\n",
     "\n",
-    "    U[idx_r, v.urho] = rho_r\n",
-    "    U[idx_r, v.umx] = rho_r * u_r\n",
-    "    U[idx_r, v.uener] = p_r/(gamma - 1.0) + 0.5 * rho_r * u_r**2"
+    "    euler.U[idx_r, euler.v.urho] = rho_r\n",
+    "    euler.U[idx_r, euler.v.umx] = rho_r * u_r\n",
+    "    euler.U[idx_r, euler.v.uener] = p_r/(gamma - 1.0) + 0.5 * rho_r * u_r**2"
    ]
   },
   {
@@ -228,7 +228,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 11,
    "id": "017682fe-b0d8-48a9-88d5-19b2d96398d0",
    "metadata": {},
    "outputs": [],
@@ -241,7 +241,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 12,
    "id": "64489150-696b-4b9b-a879-3d70a5cedb18",
    "metadata": {},
    "outputs": [

examples/hse.ipynb

@@ -92,7 +92,9 @@
    "cell_type": "code",
    "execution_count": 5,
    "id": "6935138d-852d-4d56-ac43-4960cd4f4b6b",
-   "metadata": {},
+   "metadata": {
+    "scrolled": true
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -122,7 +124,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "id": "66950c33-44e4-4975-9e96-c569396b8f5c",
    "metadata": {},
    "outputs": [
@@ -140,9 +142,7 @@
    ],
    "source": [
     "for s in simulations:\n",
-    "    init = s.grid.scratch_array(nc=3)\n",
-    "    hse(s.grid, s.v, init, s.params)\n",
-    "    print(f\"{s.grid.nx:3d} : {s.grid.norm(init[:, ivar] - s.U[:, ivar]) }\")"
+    "    print(f\"{s.grid.nx:3d} : {s.grid.norm(s.U_init[:, ivar] - s.U[:, ivar]) }\")"
    ]
   },
   {
@@ -155,17 +155,17 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 7,
    "id": "d4e9df43-d976-4e00-b99b-213f52b60839",
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "[<matplotlib.lines.Line2D at 0x7ff5d0149880>]"
+       "[<matplotlib.lines.Line2D at 0x7fe7d0d9b800>]"
       ]
      },
-     "execution_count": 10,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     },
@@ -182,11 +182,9 @@
    ],
    "source": [
     "s = simulations[0]\n",
-    "init = s.grid.scratch_array(nc=3)\n",
-    "hse(s.grid, s.v, init, s.params)\n",
     "\n",
     "fig, ax = plt.subplots()\n",
-    "ax.plot(s.grid.x[s.grid.lo:s.grid.hi+1], init[s.grid.lo:s.grid.hi+1, ivar], marker=\"o\")\n",
+    "ax.plot(s.grid.x[s.grid.lo:s.grid.hi+1], s.U_init[s.grid.lo:s.grid.hi+1, ivar], marker=\"o\")\n",
     "ax.plot(s.grid.x[s.grid.lo:s.grid.hi+1], s.U[s.grid.lo:s.grid.hi+1, ivar], marker=\"x\")"
    ]
   }

ppmpy/euler.py

@@ -56,10 +56,8 @@ class Euler:
         are: "reflect", "outflow", "periodic"
     init_cond : function
         the function to call to initialize the conserved state.
-        This has the signature `init_cond(grid, v, gamma, U, params)`
-        where `grid` is a `FVGrid`, `v` is a `FluidVars`, and `params`
-        is a `dict` with any optional parameters needed for
-        initialization.
+        This has the signature `init_cond(euler)`
+        where `euler` is an `Euler` object.
     grav_func : function
         the function to call to compute the gravitational acceleration.
         This has the signature: `g = grav_func(grid, rho, params)`, where
@@ -130,12 +128,15 @@ def __init__(self, nx, C, *,
         self.g_parabola = None
 
         # initialize
-        init_cond(self.grid, self.v, self.U, self.params)
+        init_cond(self)
         for n in range(self.v.nvar):
             self.grid.ghost_fill(self.U[:, n],
                                  bc_left_type=self.bcs_left[n],
                                  bc_right_type=self.bcs_right[n])
 
+        # save ICs for later diagnostics
+        self.U_init = self.U.copy()
+
         self.use_hse_reconstruction = use_hse_reconstruction
         self.use_flattening = use_flattening
         self.use_limiting = use_limiting

ppmpy/initial_conditions.py

@@ -4,27 +4,21 @@
 import numpy as np
 
 
-def sod(g, v, U, params):  # pylint: disable=W0613
+def sod(euler):
     """Initial conditions for the classic Sod shock tube problem
 
     Parameters
     ----------
-    grid : FVGrid
-        the grid object
-    v : FluidVars
-        the fluid variables object
-    U : ndarray
-        the conserved state array
-    params : dict
-        a dictionary of parameters.
-        We expect gamma to be provided here
+    euler : Euler
+        the Euler simulation object
 
     Returns
     -------
     None
     """
 
-    gamma = params["gamma"]
+    gamma = euler.params["gamma"]
+    g = euler.grid
 
     # setup initial conditions -- this is Sod's problem
     rho_l = 1.0
@@ -37,82 +31,72 @@ def sod(g, v, U, params):  # pylint: disable=W0613
     idx_l = g.x < 0.5
     idx_r = g.x >= 0.5
 
-    U[idx_l, v.urho] = rho_l
-    U[idx_l, v.umx] = rho_l * u_l
-    U[idx_l, v.uener] = p_l/(gamma - 1.0) + 0.5 * rho_l * u_l**2
+    euler.U[idx_l, euler.v.urho] = rho_l
+    euler.U[idx_l, euler.v.umx] = rho_l * u_l
+    euler.U[idx_l, euler.v.uener] = p_l/(gamma - 1.0) + 0.5 * rho_l * u_l**2
 
-    U[idx_r, v.urho] = rho_r
-    U[idx_r, v.umx] = rho_r * u_r
-    U[idx_r, v.uener] = p_r/(gamma - 1.0) + 0.5 * rho_r * u_r**2
+    euler.U[idx_r, euler.v.urho] = rho_r
+    euler.U[idx_r, euler.v.umx] = rho_r * u_r
+    euler.U[idx_r, euler.v.uener] = p_r/(gamma - 1.0) + 0.5 * rho_r * u_r**2
 
 
-def acoustic_pulse(g, v, U, params):  # pylint: disable=W0613
+def acoustic_pulse(euler):
     """The acoustic pulse problem from McCorquodale & Colella 2011
 
     Parameters
     ----------
-    grid : FVGrid
-        the grid object
-    v : FluidVars
-        the fluid variables object
-    U : ndarray
-        the conserved state array
-    params : dict
-        a dictionary of parameters (not used)
-        We expect gamma to be provided here
+    euler : Euler
+        the Euler simulation object
 
     Returns
     -------
     None
     """
 
-    gamma = params["gamma"]
+    gamma = euler.params["gamma"]
 
-    xcenter = 0.5 * (g.xmin + g.xmax)
+    xcenter = 0.5 * (euler.grid.xmin + euler.grid.xmax)
 
     rho0 = 1.4
     delta_rho = 0.14
 
-    r = np.abs(g.x - xcenter)
+    r = np.abs(euler.grid.x - xcenter)
 
     rho = np.where(r <= 0.5, rho0 + delta_rho * np.exp(-16.0 * r**2) * np.cos(np.pi * r)**6, rho0)
     p = (rho / rho0)**gamma
     u = 0.0
 
-    U[:, v.urho] = rho
-    U[:, v.umx] = rho * u
-    U[:, v.uener] = p / (gamma - 1.0) + 0.5 * rho * u**2
+    euler.U[:, euler.v.urho] = rho
+    euler.U[:, euler.v.umx] = rho * u
+    euler.U[:, euler.v.uener] = p / (gamma - 1.0) + 0.5 * rho * u**2
 
 
-def hse(grid, v, U, params):
+def hse(euler):
     """An isothermal hydrostatic atmosphere.
+
+    The following parameters should be set in the Euler initialization:
+
+    *  `base_density` : the density at the lower boundary
+    * `base_pressure` : the pressure at the lower boundary
+    * `g_const` : the gravitational acceleration
+    * `verbose` : output the HSE error
+
     Parameters
     ----------
-    grid : FVGrid
-        the grid object
-    v : FluidVars
-        the fluid variables object
-    U : ndarray
-        the conserved state array
-    params : dict
-        a dictionary of parameters
-
-        *  `base_density` :  the density at the lower boundary
-        * `base_pressure` :  the pressure at the lower boundary
-        * `g_const` : the gravitational acceleration
-        * `gamma` : the ratio of specific heats
+    euler : Euler
+        the Euler simulation object
 
     Returns
     -------
     None
     """
 
-    rho_base = params["base_density"]
-    pres_base = params["base_pressure"]
-    g = params["g_const"]
-    gamma = params["gamma"]
+    rho_base = euler.params["base_density"]
+    pres_base = euler.params["base_pressure"]
+    g = euler.params["g_const"]
+    gamma = euler.params["gamma"]
 
-    verbose = params.get("verbose", False)
+    verbose = euler.params.get("verbose", False)
 
     # we will assume we are isothermal and constant composition.  In that case,
     # p/rho = constant
@@ -122,34 +106,34 @@ def hse(grid, v, U, params):
     # p_{i+1} = p_i + dx / 2 (rho_i + rho_{i+1} g
     # but we can write p_{i+1} = A rho_{i+1} and solve for rho_{i+1}
 
-    p = grid.scratch_array()
-    rho = grid.scratch_array()
+    p = euler.grid.scratch_array()
+    rho = euler.grid.scratch_array()
 
     # we want the base conditions to be at the lower boundary.  We will
     # set the conditions in the first zone center from the analytic expression:
     # P = P_base e^{-z/H}
 
     H = pres_base / rho_base / np.abs(g)
 
-    p[grid.lo] = pres_base * np.exp(-grid.x[grid.lo] / H)
-    rho[grid.lo] = rho_base * np.exp(-grid.x[grid.lo] / H)
+    p[euler.grid.lo] = pres_base * np.exp(-euler.grid.x[euler.grid.lo] / H)
+    rho[euler.grid.lo] = rho_base * np.exp(-euler.grid.x[euler.grid.lo] / H)
 
-    for i in range(grid.lo+1, grid.hi+1):
-        rho[i] = (p[i-1] + 0.5 * grid.dx * rho[i-1] * g) / (A - 0.5 * grid.dx * g)
+    for i in range(euler.grid.lo+1, euler.grid.hi+1):
+        rho[i] = (p[i-1] + 0.5 * euler.grid.dx * rho[i-1] * g) / (A - 0.5 * euler.grid.dx * g)
         p[i] = A * rho[i]
 
     # now check HSE
     if verbose:
         max_err = 0.0
-        for i in range(grid.lo+1, grid.hi+1):
-            dpdr = (p[i] - p[i-1])/grid.dx
+        for i in range(euler.grid.lo+1, euler.grid.hi+1):
+            dpdr = (p[i] - p[i-1]) / euler.grid.dx
             rhog = 0.5 * (rho[i] + rho[i-1]) * g
             err = np.abs(dpdr - rhog) / np.abs(rhog)
             max_err = max(max_err, err)
 
         print(f"max err = {max_err}")
 
     # now fill the conserved variables
-    U[:, v.urho] = rho[:]
-    U[:, v.umx] = 0.0
-    U[:, v.uener] = p[:] / (gamma - 1.0)
+    euler.U[:, euler.v.urho] = rho[:]
+    euler.U[:, euler.v.umx] = 0.0
+    euler.U[:, euler.v.uener] = p[:] / (gamma - 1.0)