[](https://mybinder.org/v2/gh/schrodingersket/hedges/HEAD?urlpath=notebooks%2Fmain.ipynb)
# Hyperbolic Educational Discontinuous Galerkin Equation Solver
This repository contains a Discontinuous Galerkin Python solver for 1D Hyperbolic PDEs. It is
intended to be used primarily for education and as such prioritizes clarity in the codebase over
computational efficiency or cleverness. It is designed to be extensible and modular so that support
for new hyperbolic systems is easy to implement, and comes out of the box with an example for
the 1D shallow water equations (SWE) with variable bathymetry, prescribed inflow, and free outflow.
## Usage
Once you've installed the dependencies listed in `requirements.txt`
(`pip install -r requirements.txt`), running `main.py` will generate a full
animation of the SWE for several different variations on Gaussian inflows:
```python main.py```
Each simulation result is saved in a sub-directory of the `out` folder according to the format
`<epoch timestamp>_<six-digit random hex>`; the simulation data itself is saved in `solution.csv`,
the parameters for the simulation in `parameters.csv`, and a GIF of the solution animation is
saved at `swe_1d.gif`.
The main script (`main.py`) instantiates several instances of `simulation.SWEFlowRunner`, which is
configured specifically for running simultaneous simulations for prescribed inflow and transmissive
outflow for the Shallow Water Equations. While that class may be useful as a reference
implementation for running several simulations, the following is a full example of
running a single simulation for the purpose of experimentation/modification:
```python
#!/usr/bin/env python
# coding: utf-8
# # 1D Discontinuous Galerkin Shallow Water Solver
#
# We solve the 1D Shallow Water Equations in conservative form:
#
# \begin{align*}
# h_t + q_x &= 0 \\
# q_t + \left[ \frac{q^2}{h} + \frac{1}{2} g h^2 \right]_x &= -g h b_x - C_f \left(\frac{q}{h}\right)^2
# \end{align*}
#
# We neglect friction so that $C_f = 0$.
#
import numpy as np
import matplotlib.pyplot as plt
import hedges.bc as bc
import hedges.fluxes as fluxes
import hedges.quadrature as quadrature
import hedges.rk as rk
import hedges.swe_1d.pde as pde
# Physical parameters
#
g = 1.0 # gravity
# Domain
#
tspan = (0.0, 4.0)
xspan = (-1, 1)
# Bathymetry parameters
#
b_smoothness = 0.1
b_amplitude = 0.02
b_slope = 0.05
assert(b_smoothness > 0)
# Inflow parameters
#
inflow_amplitude = 0.05
inflow_smoothness = 1.0
inflow_peak_time = 2.0
assert(inflow_amplitude > 0)
# Initial waveheight
#
h0 = 0.2
assert(h0 > 0)
def swe_bathymetry(x):
"""
Describes bathymetry with an upslope which is perturbed by a hyperbolic tangent function.
"""
return b_slope * x + b_amplitude * np.arctan(x / b_smoothness)
def swe_bathymetry_derivative(x):
"""
Derivative of swe_bathymetry
"""
return b_slope + b_amplitude / (
b_smoothness * (1 + np.square(x / b_smoothness))
)
def q_bc(t):
"""
Describes a Gaussian inflow, where the function transitions to a constant value upon attaining
its maximum value.
:param t: Time
:return:
"""
t_np = np.array(t)
return inflow_amplitude * np.exp(
-np.square(
np.minimum(t_np, inflow_peak_time * np.ones(t_np.shape)) - inflow_peak_time
) / (2 * np.square(inflow_smoothness))
)
def initial_condition(x):
"""
Creates initial conditions for (h, uh).
:param x: Computational domain
:return:
"""
initial_height = h0 * np.ones(x.shape) - swe_bathymetry(x) # horizontal water surface
initial_flow = q_bc(0) * np.ones(x.shape) # Start with inflow BC
initial_values = np.array((
initial_height,
initial_flow,
))
# Verify consistency of initial condition
#
if not np.allclose(initial_values[1][0], q_bc(0)):
raise ValueError('Initial flow condition must match prescribed inflow.')
return initial_values
# Plot bathymetry and ICs
#
xl, xr = xspan
t0, tf = tspan
xx = np.linspace(xl, xr, num=100)
tt = np.linspace(t0, tf, num=100)
fig, (h_ax, hv_ax, q_bc_ax) = plt.subplots(3, 1)
ic = initial_condition(xx)
bb = swe_bathymetry(xx)
qq_bc = q_bc(tt)
# Plot initial wave height and bathymetry
#
h_ax.plot(xx, ic[0] + bb)
h_ax.plot(xx, bb)
h_ax.set_title('Initial wave height $h(x, 0)$')
# Plot initial flow rate
#
hv_ax.plot(xx, ic[1])
hv_ax.set_title('Initial flow rate $q(x, 0)$')
# Plot flow rate at left boundary over simulation time
#
q_bc_ax.plot(tt, qq_bc)
q_bc_ax.set_title('Boundary flow rate $q({}, t)$'.format(xl))
plt.tight_layout()
plt.show()
# Instantiate solver with bathymetry
#
solver = pde.ShallowWater1D(
b=swe_bathymetry,
b_x=swe_bathymetry_derivative,
gravity=g
)
t_interval_ms = 20
dt = t_interval_ms / 1000
surface_flux = fluxes.lax_friedrichs_flux
print('Integrating ODE system...')
solution = solver.solve(
tspan=tspan,
xspan=xspan,
cell_count=16,
polydeg=4,
initial_condition=initial_condition,
intercell_flux=surface_flux,
left_boundary_flux=pde.ShallowWater1D.bc_prescribed_inflow(
q_bc,
gravity=g,
surface_flux=surface_flux,
),
right_boundary_flux=bc.transmissive_boundary(
surface_flux=surface_flux,
direction=bc.Direction.DOWNSTREAM,
),
quad_rule=quadrature.gll,
**{
'method': rk.SSPRK33,
't_eval': np.arange(tspan[0], tspan[1], dt),
'max_step': dt, # max time step for ODE solver
'rtol': 1.0e-6,
'atol': 1.0e-6,
}
)
# Plot solution animation
#
ani, plt = solver.plot_animation(solution, frame_interval=t_interval_ms)
# Save animation to file
#
movie_name = 'swe_1d.gif'
print('Writing movie to {}...'.format(movie_name))
ani.save(movie_name, progress_callback=lambda i, n: print(
f'Saving animation frame {i + 1}/{n}'
) if i % 50 == 0 else None)
print('Animation written to {}.'.format(movie_name))
plt.show()
```
The `plot_animation` function of the `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver` base class returns
a tuple of the form
([FuncAnimation](https://matplotlib.org/3.5.1/api/_as_gen/matplotlib.animation.FuncAnimation.html),
[matplotlib](https://matplotlib.org/3.5.1/api/matplotlib_configuration_api.html#matplotlib)) which
may be used to save individual frames or modify plots as desired. Subclass implementations may
optionally override the default plotting behavior, with a reference implementation provided in
`hedges.swe_1d.ShallowWater1D`.
A Jupyter notebook is provided in this repository at `hedges/main.ipynb`, and may be viewed (and
interacted with) at
[BinderHub](https://mybinder.org/v2/gh/schrodingersket/hedges/HEAD?urlpath=notebooks%2Fmain.ipynb).
## Extending the Solver
Currently, only the 1D Shallow Water equations are implemented. However, support can easily be added
for other 1D hyperbolic PDE systems by simply subclassing `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver`
and implementing the flux `F(u)`, source `S(u)`, and flux Jacobian `F'(u)` terms for the desired
system in conservative form. See `hedges.swe_1d.ShallowWater1D` for a reference implementation.
A local Lax-Friedrichs flux is used as the approximate Riemann solver between cell
interfaces, but the `solve` method of the `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver` class accepts
a function reference which can be used to implement e.g. HLL flux (or exact Riemann solvers).
Raw data
{
"_id": null,
"home_page": "",
"name": "hedges",
"maintainer": "",
"docs_url": null,
"requires_python": ">=3.8",
"maintainer_email": "",
"keywords": "Partial differential equations,Hyperbolic,Scientific computing,Scientific machine learning,Discontinuous Galerkin",
"author": "",
"author_email": "schrodingersket <schrodingersket@gmail.com>",
"download_url": "https://files.pythonhosted.org/packages/2f/03/61669f2296e6068578e72fda29e131109a9e160ae9867ec3e13e00aa90e6/hedges-0.0.1a5.tar.gz",
"platform": null,
"description": "\n[](https://mybinder.org/v2/gh/schrodingersket/hedges/HEAD?urlpath=notebooks%2Fmain.ipynb)\n\n# Hyperbolic Educational Discontinuous Galerkin Equation Solver\n\nThis repository contains a Discontinuous Galerkin Python solver for 1D Hyperbolic PDEs. It is \nintended to be used primarily for education and as such prioritizes clarity in the codebase over\ncomputational efficiency or cleverness. It is designed to be extensible and modular so that support\nfor new hyperbolic systems is easy to implement, and comes out of the box with an example for\nthe 1D shallow water equations (SWE) with variable bathymetry, prescribed inflow, and free outflow.\n\n## Usage\n\nOnce you've installed the dependencies listed in `requirements.txt` \n(`pip install -r requirements.txt`), running `main.py` will generate a full \nanimation of the SWE for several different variations on Gaussian inflows:\n\n```python main.py```\n\nEach simulation result is saved in a sub-directory of the `out` folder according to the format\n`<epoch timestamp>_<six-digit random hex>`; the simulation data itself is saved in `solution.csv`,\nthe parameters for the simulation in `parameters.csv`, and a GIF of the solution animation is\nsaved at `swe_1d.gif`.\n\nThe main script (`main.py`) instantiates several instances of `simulation.SWEFlowRunner`, which is\nconfigured specifically for running simultaneous simulations for prescribed inflow and transmissive \noutflow for the Shallow Water Equations. While that class may be useful as a reference \nimplementation for running several simulations, the following is a full example of \nrunning a single simulation for the purpose of experimentation/modification:\n\n```python\n#!/usr/bin/env python\n# coding: utf-8\n\n# # 1D Discontinuous Galerkin Shallow Water Solver\n#\n# We solve the 1D Shallow Water Equations in conservative form:\n#\n# \\begin{align*}\n# h_t + q_x &= 0 \\\\\n# q_t + \\left[ \\frac{q^2}{h} + \\frac{1}{2} g h^2 \\right]_x &= -g h b_x - C_f \\left(\\frac{q}{h}\\right)^2\n# \\end{align*}\n#\n# We neglect friction so that $C_f = 0$.\n#\n\nimport numpy as np\nimport matplotlib.pyplot as plt\n\n\nimport hedges.bc as bc\nimport hedges.fluxes as fluxes\nimport hedges.quadrature as quadrature\nimport hedges.rk as rk\nimport hedges.swe_1d.pde as pde\n\n\n# Physical parameters\n#\ng = 1.0 # gravity\n\n# Domain\n#\ntspan = (0.0, 4.0)\nxspan = (-1, 1)\n\n# Bathymetry parameters\n#\nb_smoothness = 0.1\nb_amplitude = 0.02\nb_slope = 0.05\nassert(b_smoothness > 0)\n\n# Inflow parameters\n#\ninflow_amplitude = 0.05\ninflow_smoothness = 1.0\ninflow_peak_time = 2.0\nassert(inflow_amplitude > 0)\n\n# Initial waveheight\n#\nh0 = 0.2\nassert(h0 > 0)\n\n\ndef swe_bathymetry(x):\n \"\"\"\n Describes bathymetry with an upslope which is perturbed by a hyperbolic tangent function.\n \"\"\"\n return b_slope * x + b_amplitude * np.arctan(x / b_smoothness)\n\n\ndef swe_bathymetry_derivative(x):\n \"\"\"\n Derivative of swe_bathymetry\n \"\"\"\n return b_slope + b_amplitude / (\n b_smoothness * (1 + np.square(x / b_smoothness))\n )\n\n\ndef q_bc(t):\n \"\"\"\n Describes a Gaussian inflow, where the function transitions to a constant value upon attaining\n its maximum value.\n\n :param t: Time\n :return:\n \"\"\"\n t_np = np.array(t)\n return inflow_amplitude * np.exp(\n -np.square(\n np.minimum(t_np, inflow_peak_time * np.ones(t_np.shape)) - inflow_peak_time\n ) / (2 * np.square(inflow_smoothness))\n )\n\n\ndef initial_condition(x):\n \"\"\"\n Creates initial conditions for (h, uh).\n\n :param x: Computational domain\n :return:\n \"\"\"\n initial_height = h0 * np.ones(x.shape) - swe_bathymetry(x) # horizontal water surface\n initial_flow = q_bc(0) * np.ones(x.shape) # Start with inflow BC\n\n initial_values = np.array((\n initial_height,\n initial_flow,\n ))\n\n # Verify consistency of initial condition\n #\n if not np.allclose(initial_values[1][0], q_bc(0)):\n raise ValueError('Initial flow condition must match prescribed inflow.')\n\n return initial_values\n\n\n# Plot bathymetry and ICs\n#\nxl, xr = xspan\nt0, tf = tspan\n\nxx = np.linspace(xl, xr, num=100)\ntt = np.linspace(t0, tf, num=100)\n\nfig, (h_ax, hv_ax, q_bc_ax) = plt.subplots(3, 1)\n\nic = initial_condition(xx)\nbb = swe_bathymetry(xx)\nqq_bc = q_bc(tt)\n\n# Plot initial wave height and bathymetry\n#\nh_ax.plot(xx, ic[0] + bb)\nh_ax.plot(xx, bb)\nh_ax.set_title('Initial wave height $h(x, 0)$')\n\n# Plot initial flow rate\n#\nhv_ax.plot(xx, ic[1])\nhv_ax.set_title('Initial flow rate $q(x, 0)$')\n\n# Plot flow rate at left boundary over simulation time\n#\nq_bc_ax.plot(tt, qq_bc)\nq_bc_ax.set_title('Boundary flow rate $q({}, t)$'.format(xl))\n\nplt.tight_layout()\nplt.show()\n\n# Instantiate solver with bathymetry\n#\nsolver = pde.ShallowWater1D(\n b=swe_bathymetry,\n b_x=swe_bathymetry_derivative,\n gravity=g\n)\n\n\nt_interval_ms = 20\ndt = t_interval_ms / 1000\nsurface_flux = fluxes.lax_friedrichs_flux\nprint('Integrating ODE system...')\nsolution = solver.solve(\n tspan=tspan,\n xspan=xspan,\n cell_count=16,\n polydeg=4,\n initial_condition=initial_condition,\n intercell_flux=surface_flux,\n left_boundary_flux=pde.ShallowWater1D.bc_prescribed_inflow(\n q_bc,\n gravity=g,\n surface_flux=surface_flux,\n ),\n right_boundary_flux=bc.transmissive_boundary(\n surface_flux=surface_flux,\n direction=bc.Direction.DOWNSTREAM,\n ),\n quad_rule=quadrature.gll,\n **{\n 'method': rk.SSPRK33,\n 't_eval': np.arange(tspan[0], tspan[1], dt),\n 'max_step': dt, # max time step for ODE solver\n 'rtol': 1.0e-6,\n 'atol': 1.0e-6,\n }\n)\n\n# Plot solution animation\n#\nani, plt = solver.plot_animation(solution, frame_interval=t_interval_ms)\n\n# Save animation to file\n#\nmovie_name = 'swe_1d.gif'\nprint('Writing movie to {}...'.format(movie_name))\n\nani.save(movie_name, progress_callback=lambda i, n: print(\n f'Saving animation frame {i + 1}/{n}'\n) if i % 50 == 0 else None)\nprint('Animation written to {}.'.format(movie_name))\n\nplt.show()\n```\n\nThe `plot_animation` function of the `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver` base class returns \na tuple of the form \n([FuncAnimation](https://matplotlib.org/3.5.1/api/_as_gen/matplotlib.animation.FuncAnimation.html), \n[matplotlib](https://matplotlib.org/3.5.1/api/matplotlib_configuration_api.html#matplotlib)) which \nmay be used to save individual frames or modify plots as desired. Subclass implementations may \noptionally override the default plotting behavior, with a reference implementation provided in \n`hedges.swe_1d.ShallowWater1D`.\n\nA Jupyter notebook is provided in this repository at `hedges/main.ipynb`, and may be viewed (and\ninteracted with) at \n[BinderHub](https://mybinder.org/v2/gh/schrodingersket/hedges/HEAD?urlpath=notebooks%2Fmain.ipynb).\n\n## Extending the Solver\n\nCurrently, only the 1D Shallow Water equations are implemented. However, support can easily be added\nfor other 1D hyperbolic PDE systems by simply subclassing `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver`\nand implementing the flux `F(u)`, source `S(u)`, and flux Jacobian `F'(u)` terms for the desired \nsystem in conservative form. See `hedges.swe_1d.ShallowWater1D` for a reference implementation.\n\nA local Lax-Friedrichs flux is used as the approximate Riemann solver between cell \ninterfaces, but the `solve` method of the `hedges.hyperbolic_solver_1d.Hyperbolic1DSolver` class accepts\na function reference which can be used to implement e.g. HLL flux (or exact Riemann solvers).\n",
"bugtrack_url": null,
"license": "LGPL-2.1",
"summary": "A Discontinuous Galerkin solver oriented toward prototyping and education",
"version": "0.0.1a5",
"project_urls": {
"Bug Tracker": "https://github.com/schrodingersket/hedges/issues",
"Homepage": "https://github.com/schrodingersket/hedges"
},
"split_keywords": [
"partial differential equations",
"hyperbolic",
"scientific computing",
"scientific machine learning",
"discontinuous galerkin"
],
"urls": [
{
"comment_text": "",
"digests": {
"blake2b_256": "f1d1c51da12fef652975a435d388c4d519d060658da150c3ccf13b17f1399dd8",
"md5": "ebd4f7fddb03a45d655d02cd14277ad1",
"sha256": "47c46b80e8f0827ef769ae31ea38208ee148a7d8617e74ae65be6894b5789b6a"
},
"downloads": -1,
"filename": "hedges-0.0.1a5-py3-none-any.whl",
"has_sig": false,
"md5_digest": "ebd4f7fddb03a45d655d02cd14277ad1",
"packagetype": "bdist_wheel",
"python_version": "py3",
"requires_python": ">=3.8",
"size": 17877,
"upload_time": "2023-08-12T20:00:37",
"upload_time_iso_8601": "2023-08-12T20:00:37.198698Z",
"url": "https://files.pythonhosted.org/packages/f1/d1/c51da12fef652975a435d388c4d519d060658da150c3ccf13b17f1399dd8/hedges-0.0.1a5-py3-none-any.whl",
"yanked": false,
"yanked_reason": null
},
{
"comment_text": "",
"digests": {
"blake2b_256": "2f0361669f2296e6068578e72fda29e131109a9e160ae9867ec3e13e00aa90e6",
"md5": "2b4e66c679277f105e85aba26237b28d",
"sha256": "72b3b6e9f5961db250836c093416192d51a8c71f35f11bb2ca3aa871741a6d4b"
},
"downloads": -1,
"filename": "hedges-0.0.1a5.tar.gz",
"has_sig": false,
"md5_digest": "2b4e66c679277f105e85aba26237b28d",
"packagetype": "sdist",
"python_version": "source",
"requires_python": ">=3.8",
"size": 20348,
"upload_time": "2023-08-12T20:00:38",
"upload_time_iso_8601": "2023-08-12T20:00:38.652604Z",
"url": "https://files.pythonhosted.org/packages/2f/03/61669f2296e6068578e72fda29e131109a9e160ae9867ec3e13e00aa90e6/hedges-0.0.1a5.tar.gz",
"yanked": false,
"yanked_reason": null
}
],
"upload_time": "2023-08-12 20:00:38",
"github": true,
"gitlab": false,
"bitbucket": false,
"codeberg": false,
"github_user": "schrodingersket",
"github_project": "hedges",
"travis_ci": false,
"coveralls": false,
"github_actions": false,
"requirements": [],
"lcname": "hedges"
}
JSON
