Start also this weeks assignments by updating your repository! Either press the sync fork button on gihub.com, or pull in the latest changes through
git pull
This should pull in a new file Wave1D.py containing a new wave equation solver. Study this solver and try to understand how it works. This week we will focus even more on boundary conditions as this is normally the most complicated part of any PDE solver.
i) The pulled in solver is implemented with Dirichlet boundary conditions. The first task is to implement Neumann boundary conditions. Neumann conditions are implemented for x=0 using
\(\frac{u^{n+1}_1 - u^{n+1}_{-1}}{2 \Delta t} = 0\)
This gives us \(u^{n+1}_{-1} = u^{n+1}_1\), which can be used in the discretized PDE for the boundary node \(i=0\). We get Eq. (2.38) for x=0 and (2.40) for x=L. These equations can most easily and elegantly be implemented simply by modifying the first and last line of the difference matrix D2.
ii) Next assignment is to implement open boundaries. This is also considered in problem 2.12 in the FDM book. The open boundaries require that at \(x=0\) and \(x = L\), we have, respectively:
\(\frac{\partial u}{\partial t} - c \frac{\partial u}{\partial x} = 0 \quad \text{and} \quad \frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0\)
This is easiest to implement explicitly using first order forward and backward stencils. For the left boundary this means
\(\frac{u_0^{n+1}-u_0^n}{\Delta t} - c \frac{u_1^n - u_0^n}{\Delta x}=0\)
which can be solved for \(u^{n+1}_0\), which is exactly the value you wish to set at the boundary. A similar expression derived with a backward stencil for space will give \(u^{n+1}_N\) for the right hand side.
This open boundary implementation described above works, but it is only first order accurate. Using second order accurate stencils for the left boundary leads to
(1) \(\frac{u_0^{n+1}-u_0^{n-1}}{2\Delta t} - c \frac{u_1^n - u_{-1}^n}{2 \Delta x}=0\)
This stencil can not be used explicitly, because it depends on \(u^{n}_{-1}\). However, you can use (1) to solve for \(u^{n}_{-1}\), and then use this to eliminate \(u^{n}_{-1}\) in the discrete PDE for \(n=0\) . This will give you an expression to use for \(u^{n+1}_0\). A similar procedure works for the right boundary.
iii) Modify the solver such that you can use any boundary condition on any side of the domain. That is, you should be able to choose Dirichlet on the left and then Neumann of the right etc.
Suggestion: Instead of just `bc` (see `apply_bcs`) use a dictionary such that {'left': 0, 'right': 1} would indicate Dirichlet on the left and Neumann on the right. The solver is already set to use
- 0 Dirichlet
- 1 Neumann
- 2 Open boundary
- 3 Periodic
You need to modify `apply_bcs` such that both sides are treated individually. You should also modify `D2`, where the Neumann condition is implemented implicitly in the differentiation matrix.
iv) Implement periodic boundary conditions. The periodic boundary condition cannot be mixed with other conditions.
v) Verify the implementations of the different boundaries using a Gaussian pulse that travels in both directions. The initial condition is
(2) \(u(x, t) = \exp(-200 (x-L/2-ct)^2)\)
The Wave1D solver is already implemented with this pulse as default initial condition. However, there are still two possibilities and you choose the behaviour with the `ic` keyword when calling the solver.
ic=0 - uses Eq. (2) for initial condition and in addition \(\frac{\partial u}{\partial t}=0 \)
This leads to two pulses that travels in both directions.
ic=1 - uses Eq. (2) for initial condition AND a second initial condition at \(u(x, -\Delta t)\). This leads to a single pulse that travels in only one direction.
For verification
- The Neumann condition should simply send the wave back exactly as it was only moving in the other direction.
- The open boundary should see the wave disappear.
All boundary conditions are tested in the function `test_pulse_bcs`. When this function passes the boundary conditions should be implemented correctly (even though the test may be rather simple).
Compare the implemented boundary conditions also visually with this online solver.