Skip to content

simplified the initial conditions interface #39

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
zingale merged 1 commit into from
Aug 18, 2024
Merged

simplified the initial conditions interface #39

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
38 changes: 19 additions & 19 deletions examples/euler.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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"
]
},
{
Expand All @@ -228,7 +228,7 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 11,
"id": "017682fe-b0d8-48a9-88d5-19b2d96398d0",
"metadata": {},
"outputs": [],
Expand All @@ -241,7 +241,7 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 12,
"id": "64489150-696b-4b9b-a879-3d70a5cedb18",
"metadata": {},
"outputs": [
Expand Down
20 changes: 9 additions & 11 deletions examples/hse.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,9 @@
"cell_type": "code",
"execution_count": 5,
"id": "6935138d-852d-4d56-ac43-4960cd4f4b6b",
"metadata": {},
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
Expand Down Expand Up @@ -122,7 +124,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 6,
"id": "66950c33-44e4-4975-9e96-c569396b8f5c",
"metadata": {},
"outputs": [
Expand All @@ -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]) }\")"
]
},
{
Expand All @@ -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"
},
Expand All @@ -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\")"
]
}
Expand Down
11 changes: 6 additions & 5 deletions ppmpy/euler.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
110 changes: 47 additions & 63 deletions ppmpy/initial_conditions.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)