Copyright (C) 2010-2020 Luke Olson
Copyright (C) 2020 Andreas Kloeckner
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
def gaussian(x):
u = sp.exp(-100 * (x - 0.25)**2)
return u
def step(x):
u = np.zeros(x.shape)
for j in range(len(x)):
if (x[j] >= 0.6) and (x[j] <= 0.8):
u[j] = 1.0
return u
def g1(x):
return 1+gaussian(x)
def g2(x):
return 1+gaussian(x) + step(x)
g = g1
nx = 164
x = np.linspace(0, 1, nx, endpoint=False)
dx = x[1] - x[0]
xx = np.linspace(0, 1, 1000, endpoint=False)
lmbda = 0.95
nt = 250
print('tsteps = %d' % nt)
print(' dx = %g' % dx)
print('lambda = %g' % lmbda)
J = np.arange(0, nx) # all vertices
Jm1 = np.roll(J, 1)
Jp1 = np.roll(J, -1)
plt.plot(x, g(x))
tsteps = 250 dx = 0.00609756 lambda = 0.95
[<matplotlib.lines.Line2D at 0x7f1ec5d4ccd0>]
Plot the solution:
if 1:
# Burgers
def f(u):
return u**2/2
def fprime(u):
return u
else:
# advection
def f(u):
return u
def fprime(u):
return 1+0*u
steps_per_frame = 2
Implement rhs
for a Lax-Friedrichs flux:
Recall (local) Lax-Friedrichs: $$ f^{\ast} (u^{_-}, u^+) = \frac{f (u^-) + f (u^+)}{2} - \frac{\alpha}{2} (u^+ - u^-) \quad\text{with}\quad \alpha = \max \left( |f' (u^-)|, |f' (u^+)| \right).$$ Recall FV: $$ \bar{u}_{j,\ell+1} = \bar{u}_{j,\ell} - \frac{h_t}{h_x} (f (u_{j + 1 / 2,\ell}) - f (u_{j - 1 / 2,\ell})) . $$
#clear
def rhs(u):
uplus = u[Jp1]
uminus = u[J]
alpha = np.maximum(np.abs(fprime(uplus)), np.abs(fprime(uminus)))
# right-looking, between J and Jp1
fluxes = (
(f(uplus)+f(uminus))/2
- alpha/2*(uplus-uminus)
)
return - 1/dx*(fluxes[J]-fluxes[Jm1])
u = g(x)
def timestepper(n):
for i in range(steps_per_frame):
dt = dx*lmbda/np.max(np.abs(u))
u[:] = u + dt*rhs(u)
line.set_data(x, u)
return line
fig = plt.figure(figsize=(5,5))
line, = plt.plot(x, u)
ani = animation.FuncAnimation(
fig, timestepper,
frames=nt//steps_per_frame,
interval=30)
html = HTML(ani.to_jshtml())
plt.clf()
html
<Figure size 360x360 with 0 Axes>
First, need a second-order time integrator:
def rk2_step(dt, u, rhs):
k1 = rhs(u)
k2 = rhs(u+dt*k1)
return u+0.5*dt*(k1+k2)
Now upgrade the reconstruction to second order.
NOTE: It's OK to assume (here!) that the wind blows from the right to simplify upwinding.
#clear
def rhs(u):
# right-looking, between J and Jp1
fluxes = f(u[J] + 1/2*(u[J]-u[Jm1]))
return - 1/dx*(fluxes[J]-fluxes[Jm1])
u = g(x)
def timestepper(n):
# to simplify upwinding
assert np.min(u) >= 0
for i in range(steps_per_frame):
dt = 0.7*dx*lmbda/np.max(np.abs(u))
u[:] = rk2_step(dt, u, rhs)
line.set_data(x, u)
return line
fig = plt.figure(figsize=(5,5))
line, = plt.plot(x, u)
ani = animation.FuncAnimation(
fig, timestepper,
frames=nt//steps_per_frame,
interval=30)
html = HTML(ani.to_jshtml())
plt.clf()
html
<Figure size 360x360 with 0 Axes>