#!/usr/bin/env python # coding: utf-8 # # Experimenting with Dispersion and Dissipation # # Copyright (C) 2010-2020 Luke Olson
# Copyright (C) 2020 Andreas Kloeckner # #
# MIT License # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal # in the Software without restriction, including without limitation the rights # to use, copy, modify, merge, publish, distribute, sublicense, and/or sell # copies of the Software, and to permit persons to whom the Software is # furnished to do so, subject to the following conditions: # # 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. #
# In: import numpy as np import scipy as sp import matplotlib.pyplot as plt # ## Problem Description # Consider # $$u_t + a u_x = 0$$ # with periodic BC on the interval $[0,1]$. # In: a = 1.0 T = 4.0 / a # end time for four cycles # ## Set up the Grid # - dx will be the grid spacing in the $x$-direction # - x will be the grid coordinates # - xx will be really fine grid coordinates # In: nx = 90 k = 10 x = np.linspace(0, 2*np.pi, nx, endpoint=False) dx = x - x xx = np.linspace(0, 2*np.pi, 1000, endpoint=False) # Now define an initial condition: # In: def f(x): u = np.sin(k * x) return u # In: plt.plot(xx, f(xx), lw=3, clip_on=False) # ## Setting the Time Step # # Now we need a time step. Let # $$\Delta t = \Delta x \frac{\lambda}{a}$$ # with CFL number $\lambda$. # In: lmbda = 0.9 dt = dx * lmbda / a nt = int(T/dt) print('T = %g' % T) print('tsteps = %d' % nt) print(' dx = %g' % dx) print(' dt = %g' % dt) print('lambda = %g' % lmbda) # Make an index list, called $J$, so that we can access $J+1$ and $J-1$ easily. # In: J = np.arange(0, nx) # all vertices Jm1 = np.roll(J, 1) Jp1 = np.roll(J, -1) # ## Run and Plot Final Solution # In: plotit = True uETBS = f(x) uLW = f(x) if plotit: fig = plt.figure(figsize=(10,10)) ax = fig.add_subplot(111) ax.set_title('u vs x') for n in range(0, nt): uETBS[J] = uETBS[J] - lmbda * (uETBS[J] - uETBS[Jm1]) uLW[J] = uLW[J] - lmbda * (1.0 / 2.0) * (uLW[Jp1] - uLW[Jm1]) + (lmbda**2 / 2.0) * (uLW[Jp1] - 2 * uLW[J] + uLW[Jm1]) uex = f((xx - a * (n+1) * dt) % 1.0) if plotit: ax.plot(x, uETBS, '-', clip_on=False, label="ETBS") ax.plot(x, uLW, '-', clip_on=False, label="Lax-Wendroff") ax.plot(xx, uex, 'k-', clip_on=False, label='exact') ax.legend(frameon=False) plt.show() # In[ ]: # In[ ]: