<jupyter_start><jupyter_text>Week 8 worksheet: Numerical solution of hyperbolic PDEs$$\newcommand{\vect}[1]{\bm 1}\newcommand{\grad}{\nabla}\newcommand{\pderiv}[2]{\frac{\partial 1}{\partial 2}}\newcommand{\pdderiv}[2]{\frac{\partial^2 1}{\partial 2^2}}\newcommand{\deriv}[2]{\frac{\mathrm{d} 1}{\mathrm{d} 2}}\newcommand{\Deriv}[3]{\frac{\mathrm{d}^3 1}{\mathrm{d} 2^3}}$$In this laboratory we will look at hyperbolic problems, solving the one dimensional, linear advection equation$$\pderiv{u}{t}+a\pderiv{u}{x}=0,$$where $u$ is a scalar quantity being advected at a uniform speed, $a$, over a periodic domain, $x\in[0,1]$.The initial conditions are $$u(x,0) = \begin{cases}\sin^4 2\pi x & 0 \le x \le \frac12 \\1 & \frac46 \le x \le \frac56 \\0 & \text{otherwise}\end{cases}$$If the advection speed is an integer then every second the wave form will be back at the starting position.  This means we can make a direct comparison between the initial conditions and the results from the solver.Your are going to develop a simple one dimensional solver using Python to solve the problem.  In this question you are going to write the key parts of the driver routines needed to create a grid, calculate the time step, set up the initial conditions, solve the equations using the first order upwind scheme and plot the results.  The code starts with the following imports and definitions:<jupyter_code>import numpy as np
from matplotlib import pyplot as plt

N = 500 # Number of grid points being used
AdvectionSpeed = 1.0 # Advection speed (celarity) of the wave
CourantNumber = 0.95 # Courant number
%run scripts/create_widgets.py W07<jupyter_output><empty_output><jupyter_text>*How it works: You will see cells located below each exercise, each containing a command starting with `%run scripts/show_solutions.py`. You don't need to run those yourself; the command above runs a script which automatically runs these specific cells for you. The commands in each of these cells each create the button for the corresponding exercise. The Python code to achieve this is contained in `scripts/show_solutions.py`, and relies on [IPython widgets](https://ipywidgets.readthedocs.io/en/latest/examples/Widget%20Basics.html) --- feel free to take a look at the code if you are curious.*<jupyter_code>%%javascript
MathJax.Hub.Config({
    TeX: { equationNumbers: { autoNumber: "AMS" } }
});<jupyter_output><empty_output><jupyter_text>Exercise 1 - Grid generation Part a)Write a function _GridGen(N)_ which creates an array $x$ which contains a uniformly spaced grid of $N$ grid points over the interval $[0,1]$.  The function must return both the mesh spacing $\Delta x$ and the $N$ ordinates $x_i$.You should fill in the following template.<jupyter_code>def GridGen(N):
    '''generate a 1D grid from over [0.0,1.0] with N grid points.'''
   
    return dx, x_ord
%run scripts/show_solutions.py W08_ex1_parta<jupyter_output><empty_output><jupyter_text>Write a function _InitialConditions(c)_ which returns an array $u$  so that$$u = \begin{cases}\sin^4 2\pi x & 0 \le x \le \frac12 \\1 & \frac46 \le x \le \frac56 \\0 & \text{otherwise}\end{cases}$$You should fill in the following template. def InitialConditions(x_ord):    ''' set up the intial conditions so that u(x)=sin^4(2 Pi x) for x    in [0.0,0.5] and u(x) = 1 for x in [4/6,5/6] and u(x)=0    otherwise.  These are the initial conditions for the classical    test problem.'''    return u<jupyter_code>%run scripts/show_solutions.py W08_ex1_partb<jupyter_output><empty_output><jupyter_text>Write a function _TimeStep($\Delta x$)_ which computes the stable time step, $$\Delta t = \nu \frac{\Delta x}{|a|}$$ where $\nu$ is the Courant number, and $a$ is the advection speed.  Remember that the advection speed, $a$, is given by the Python constant AdvectionSpeed and the Courant number, $\nu$ by the Python constant CourantNumber.<jupyter_code>%run scripts/show_solutions.py W08_ex1_partc<jupyter_output><empty_output><jupyter_text>Part b)The following Python function is the driver for the 1D solver.  It is passed the parameters t\_stop the time at which the solution is required, u a vector containing the initial conditions,  dx the mesh spacing $\Delta x$ and method, a function with parameters (u,dx,dt) which performs **one** time step with an advection solver and which returns both the updated solution $u^{n+1}$ and the name of the method used.<jupyter_code>def OneDSolver(t_stop, u, dx, method):
    '''One dimensional advection solver.  This calculates the time
    step and then performs the appropriate number of time steps to
    reach t_stop using the method passed in. '''

    # calculate the time step
    dt = TimeStep(dx)
    MaxStep = int(t_stop / dt)  + 1 # maximum number of time steps

    print('1D solver - performing {} time steps, dt = {:.4e}'
              .format(MaxStep, dt))
    time = 0.0
    for i in range(MaxStep):
        dt = min(dt, t_stop - time)

        # There is no need to implement periodic boundary conditions
        # here as these will be dealt with within the method using the
        # np.roll() function wich rotates the contents of an array.
        
        # call the solver
        u, scheme = method(u,dx,dt)

        # update time
        time = time + dt
        if (time >= t_stop):
            print('Done. {} seconds simulated using the {} scheme.'
                      .format(time, scheme)) 
            break
    return u, scheme<jupyter_output><empty_output><jupyter_text>Implement the first order upwind scheme $$u^{n+1}_i = u^n_i - \frac{a\Delta t}{\Delta x}\left(u^n_i - u^n_{i-1}\right)$$You should use the following template and should make use of the np.roll() function to obtain the values of $u^n_{i-1}$.  Remember that the advection speed, $a$, is given by the Python constant AdvectionSpeed.<jupyter_code>def FOU(u,dx,dt):
    ''' this is the classical First Order Upwind scheme, this is first
    order in time and space and has a stability condition that
    nu<=1.

    The periodic boundary conditions are handled using np.roll() which
    shifts an array left or right rolling round the end values.'''

    return u, 'FOU'
%run scripts/show_solutions.py W08_ex1_partd<jupyter_output><empty_output><jupyter_text>The final step is to write the main program this needs to set up the gird and the initial conditions, save a copy of the initial conditions, call the solver and plot the results.  It can be coded as follows:<jupyter_code>'''Main Program'''
# set up the grid and the initial conditions
dx,x = GridGen(N)
u = InitialConditions(x)
u0 = u.copy()

# call the solver
tOut=5
u, scheme = OneDSolver(tOut, u, dx, FOU)

# Plot the results
fig = plt.figure()
plt.plot(x,u0,'-r')
if (N>120):
    plt.plot(x,u,'-b')
else:
    plt.plot(x,u,'ob')
plt.title('Linear advection {} scheme, N={}, t={}'.format(scheme,N,tOut))
plt.show()<jupyter_output><empty_output><jupyter_text>Use the solver to explore what happens if the solution is evolved for 0.5, 1.0, 10.0 and 50.0 seconds on a grid with $N=500$ points. Briefly comment on what you observe.<jupyter_code>%run scripts/show_solutions.py W08_ex1_parte<jupyter_output><empty_output><jupyter_text>Exercise 2Implement the Lax-Friedrichs scheme, $$u^{n+1}_i = \frac{u^n_{i+1}+u^n_{i-1}}{2} - \frac{a\Delta t}{2\Delta x}\left(u^n_{i+1}-u^n_{i-1}\right).$$ You are reminded that the np.roll() function should be used to obtain the values of $u^n_{i-1}$ and $u^n_{i+1}$. Remember that the advection speed, $a$, is given by the Python constant AdvectionSpeed.  You should complete the following Python template<jupyter_code>def LaxFriedrichs(u,dx,dt):
    '''Classical Lax-Friedrichs scheme'''

    return u, 'Lax-Friedrichs'<jupyter_output><empty_output><jupyter_text>Calculate the solution at 1.0, 10.0 and 50.0 seconds on a grid with $N=500$ points. Briefly comment on what you observe.<jupyter_code>%run scripts/show_solutions.py W08_ex2_parta<jupyter_output><empty_output><jupyter_text>Calculate the solution at 1.0 on a grid with $N=50$ points. What do you observe?<jupyter_code>%run scripts/show_solutions.py W08_ex2_partb<jupyter_output><empty_output><jupyter_text>Exercise 3Implement the Lax-Wendroff scheme, $$u^{n+1}_i = u^n_i-\frac{a\Delta t}{2\Delta x}\left(u^n_{i+1}+u^n_{i-1}\right) - \frac{a^2\Delta t^2}{2\Delta x^2}\left(u^n_{i+1}-2u^n_i+u^n_{i-1}\right)$$ You are reminded that the np.roll() function should be used to obtain the values of $u^n_{i-1}$ and $u^n_{i+1}$. Remember that the advection speed, $a$, is given by the Python constant AdvectionSpeed.  You should complete the following Python template<jupyter_code>def LaxWendroff(u,dx,dt):
    '''Classical, 2nd order, Lax-Wendroff scheme'''
      
    return u, 'Lax-Wendroff'<jupyter_output><empty_output><jupyter_text>Calculate the solution at 1.0, 10.0 and 50.0 seconds on a grid with $N=500$ points. Briefly comment on what you observe.<jupyter_code>%run scripts/show_solutions.py W08_ex3_parta<jupyter_output><empty_output><jupyter_text>Exercise 4Implement the MacCormack predictor-correctorscheme, $$\begin{align*}u^p_i &= u ^n_i - \frac{a\Delta t}{\Delta x}\left(u^n_i - u^n_{i-1}\right) \\u^{n+1}_i &= \frac{u^p_i+u^n_i}{2} - \frac{a\Delta t}{2\Delta x}\left(u^p_{i+1} - u^p_i\right)\end{align*}$$ You are reminded that the np.roll() function should be used to obtain the values of $u^n_{i-1}$ and $u^n_{i+1}$. Remember that the advection speed, $a$, is given by the Python constant AdvectionSpeed.  You should complete the following Python template<jupyter_code>def MacCormack(u,dx,dt):
    '''2nd order MacCormack predictor corrector method which is
    equivalent to the classical Lax-Wendroff scheme.'''

    # Predictor
    up = 
    # Corrector
    u =
    
    return u, 'MacCormack'<jupyter_output><empty_output><jupyter_text>Check that the solution at 10.0 seconds on a grid with $N=500$ points is equivalent to the solution obtained with the Lax-Wendroff Scheme.<jupyter_code>%run scripts/show_solutions.py W08_ex4_parta
### Part b)

Using either the MacCormack or the Lax-Wendroff scheme see what happens when $a=-1$, 
compare the solution at $t=10.0$ seconds with the case where $a=1$

%run scripts/show_solutions.py W08_ex4_partb<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11313190573970862, 'fraction_numerical': 0.020338320132981324, 'mean_word_length': 4.510775862068965, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 67, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '18779', 'license': 'CC0-1.0', 'n_tokens_mistral': 3442, 'n_tokens_neox': 3288, 'n_words': 1179, 'path': 'laboratories/W08_PDE3_hyperbolic_PDE_Numerical.ipynb', 'repo_name': 'dmingram66/PDE3', 'size': '10227'}