diff --git a/examples/euler.ipynb b/examples/euler.ipynb index cf86eca..73cd0f6 100644 --- a/examples/euler.ipynb +++ b/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": [ diff --git a/examples/hse.ipynb b/examples/hse.ipynb index 7c171f5..e9e62fa 100644 --- a/examples/hse.ipynb +++ b/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": [ - "[]" + "[]" ] }, - "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\")" ] } diff --git a/ppmpy/euler.py b/ppmpy/euler.py index 1bfb16a..807c132 100644 --- a/ppmpy/euler.py +++ b/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 diff --git a/ppmpy/initial_conditions.py b/ppmpy/initial_conditions.py index a95619c..1c8eb27 100644 --- a/ppmpy/initial_conditions.py +++ b/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,8 +106,8 @@ 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: @@ -131,18 +115,18 @@ def hse(grid, v, U, params): 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) @@ -150,6 +134,6 @@ def hse(grid, v, U, params): 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)