oblig2 in the top directory of your INF5620
repo on GitHub or
Bitbucket to hold the various files of the project. Include a README
file with a short overview of the different files and the full names
of the students who did the project.The project addresses the two-dimensional, standard, linear wave equation, with damping, $$ \begin{equation} \frac{\partial^2 u}{\partial t^2} + b\frac{\partial u}{\partial t} = \frac{\partial}{\partial x}\left( q (x,y) \frac{\partial u}{\partial x}\right) + \frac{\partial}{\partial y}\left( q (x,y) \frac{\partial u}{\partial y}\right) + f(x,y,t), \end{equation} $$ with boundary condition $$ \begin{equation} \frac{\partial u}{\partial n} = 0, \end{equation} $$ in a rectangular spatial domain \( \Omega = [0,L_x]\times [0,L_y] \). The initial conditions are
$$ \begin{align} u(x,y,0)&=I(x,y),\\ u_t(x,y,0)&=V(x,y). \end{align} $$
Derive the discrete set of equations to be implemented in a program:
Compute the truncation error of the scheme at an arbitrary interior mesh point. (It easier to start with \( q=\hbox{const} \) and then generalize to variable \( q \)).
Implement the numerical scheme in a program. You may use wave2D_u0.py as a starting point (this program solves the 2D wave equation with constant wave velocity and \( u=0 \) on the boundary). You will need to include both scalar (pointwise) computation of the scheme for debugging and reference as well as a vectorized version for speed.
Construct a test case with constant solution (not 0 or 1). Make a corresponding nose test.
Invent five types of possible bugs in the implementation of the mathematical formulas. See how many of them that lead to a wrong non-constant solution.
Assume an exact solution on the form \( \uex(x,y,t)=X(x)Y(y)T(t) \) where \( X \), \( Y \), and \( T \) are polynomials of degree three or less. Construct \( X \) and \( Y \) such that the normal derivative vanishes at the four boundaries. Fit a corresponding source term \( f(x,y,t) \) in the wave equation.
It would be great if this exact solution also were an exact solution of the discrete equations when \( q \) is constant and \( b=0 \), which is often the case with lower-order polynomials because the truncation errors involve higher-order derivatives. The exact cubic solution fits the discrete equations at all inner mesh points, but not on the boundary. Add a term in the boundary equation for testing such that the exact solution also fulfills the boundary equation. Implement a corresponding nose test.
Hint. We outline some ideas in 1D. For constant \( q \) and \( b=0 \) we have the scheme \( [D_tD_t u - qD_xD_x = f]^n_{i} \) at inner points. Because \( [D_xD_x x^3]_i = 6x_i \), \( [D_xD_x x^2]_i = 2 \), \( [D_xD_x x]_i = 0 \) all are exact, the suggested \( u \) also fits the discrete equation. On the boundary we get a modified scheme, which in operator notation can be written as
$$ [D_tD_t u= qD_xD_x u + q\frac{2}{\Delta x}D_{2x} u + f]^n_i,\quad i=0,$$ and
$$ [D_tD_t u= qD_xD_x - q\frac{2}{\Delta x}D_{2x} u + f]^n_i,\quad i=N_x\tp$$ We have the results \( [D_{2x} x^3]_i = 3x_i^2 + \Delta x^2 \), and \( [D_{2x} x^2]_i = 2x_i \), \( [D_{2x} x]_i=1 \). Consider now \( [D_tD_t u= qD_xD_x u - q\frac{2}{\Delta x}D_{2x}u + f]^n_i \). Inserting \( u=X(x)T(t) \) requires the same \( f \) as for the PDE, but with an additional term \( T(t)2q\Delta x \) because of the \( D_{2x} \) operator acting on a cubic polynomial in \( x \).
This test requires the \( f \) that fits the PDE to be modified on the boundary. A possible implementation is to modify the array of \( f \) values at the boundary mesh points directly, or perform tests on the coordinates if a pointwise evaluation of \( f \) is requested:
def f(x, y, t):
if isinstance(x, np.ndarray) and isinstance(y, np.ndarray):
# Array evaluation
f_a = .... # evaluate the f that fits the PDE
# Modify boundary values
f_a[0,:] = ... # x=0
f_a[-1,:] = ... # x=Lx
f_a[:,0] = ... # y=0
f_a[:,-1] = ... # x=Ly
else:
# Assume pointwise evaluation
tol = 1E-14 # tolerance for float comparison
f_v = ... # evaluate the f that fits the PDE
# Modify boundary values
if abs(x) < tol:
f_v = ... # x=0
if abs(x-Lx) < tol:
f_v = ... # x=Lx
if abs(y) < tol:
f_v = ... # y=0
if abs(y-Ly) < tol:
f_v = ... # y=Ly
The program wave1D_dn_vc.py has a pulse function for simulating the propagation of a plug wave,
where \( I(x) \) is constant in some region of the domain and zero elsewhere.
With unit Courant number, the plug is split into two identical waves,
moving in opposite direction, exactly one cell per time
step. The discrete solution is then equal to the exact solution.
Set \( b=0 \) and \( q \) to a constant. Test the 2D program using a one-dimensional plug wave in \( x \) direction with \( c\Delta t/\Delta x=1 \) (the plug is constant in \( y \) direction and hence compatible with the \( \partial/\partial y=0 \) boundary condition). Also propagate a one-dimensional plug wave in the \( y \) direction with \( c\Delta t/\Delta y=1 \). Both test cases are essentially 1D test cases, and the results should as in the 1D case. Implement a corresponding nose test.
With no damping and constant wave velocity, our wave equation problem without any source term admits a standing wave solution:
$$ \begin{equation} \uex(x,y,t)=A\cos(k_xx)\cos(k_yy)\cos (\omega t),\quad k_x = \frac{m_x\pi}{L_x}, \ k_y = \frac{m_y\pi}{L_y}, \label{wave:app:exer:standing:waves} \end{equation} $$ for arbitrary amplitude \( A \), arbitrary integers \( m_x \) and \( m_y \), and a suitable choice of \( \omega \). This solution can be used to test the convergence rate of the numerical method.
Use the truncation error analysis to set up a model for the error in terms of the discretization parameters. Note that the truncation error is not the true error \( e^n_{i,j}=\uex(x_i,y_j,t_n) - u^n_{i,j} \), but we assume that \( e^n_{i,j} \) depends on the discretization parameters in the same way. We also assume that a norm of \( e^n_{i,j} \), e.g.,
$$ E = ||e^n_{i,j}||_{\ell^\infty} = \max_{i}\max_{j}\max_{t} |e^n_{i,j}|,$$ also depends on the discretization parameters in the same way. Introduce a common discretization parameter \( h \) such that \( \Delta t \), \( \Delta x \), and \( \Delta y \) are proportional to \( h \). Show that \( E=Ch^2 \) for some \( C \) is the expected behavior of the error. Perform experiments with decreasing \( h \), compute \( E \), and verify that \( E/h^2 \) is approximately constant.
Try to find an analytical solution of damped waves using an ansatz of the type
$$ \begin{equation} \uex(x,y,t)=\left(A\cos (\omega t) + B\sin (\omega t)\right) e^{-ct}\cos(k_xx)\cos(k_yym_yy\pi/L_y),\quad k_x = \frac{m_x\pi}{L_x}, \ k_y = \frac{m_y\pi}{L_y}\tp \label{wave:app:exer:standing:waves:damped} \end{equation} $$ That is, find \( A \), \( B \), \( \omega \), and \( c \) such that \eqref{wave:app:exer:standing:waves:damped} solves the PDE with constant \( q \), no source term, and initial condition \( u_t(x,y,0)=0 \) (as for the undamped standing waves). Make a corresponding convergence test.
Hint. The algebra can quickly be quite involved. Getting an overview of the algebra in a 1D version of this problem might be helpful. Start with relating \( A \) and \( B \) through the initial conditions (\( u=A\cos k_xx\cos k_yy \) and \( u_t=0 \) as implied by \eqref{wave:app:exer:standing:waves}) and eliminate \( B \). After having inserted \( u \) in the PDE, two equations for \( \omega \) and \( c \) arise from factoring the sine and cosine terms in time. One equation can be solved for \( \omega = \sqrt{k_x^2q + k_y^2q - c^2} \), while the other can be solved for \( c=b/2 \) by inserting the found \( \omega \) expression.
Choose some \( q(x,y)\neq 0 \) and find \( f(x,y,t) \) such that ({wave:app:exer:standing:waves:damped}) is a solution to the general 2D wave equation problem with damping and variable wave velocity. Find corresponding \( I \) and \( V \), and make a convergence test that recovers the expected convergence rate. Make a corresponding nose test.
Hint.
You may explore sympy for automating the analytical work:
>>> from sympy import *
>>> x = Symbol('x')
>>> q=x**2
>>> u=sin(x)
>>> r = diff(q*diff(u, x), x) # Derivative: (q*u_x)_x
>>> simplify(r)
x*(-x*sin(x) + 2*cos(x))
The purpose of this part is to explore what happens to a wave that enters a medium with different wave velocity. A particular physical interpretation can be wave propagation of a tsunami over a subsea hill. The unknown \( u(x,y,t) \) is then the elevation of the ocean surface, and the boundary condition \( \partial u/\partial n=0 \) means that the waves are perfectly reflected, because of a steep hill at the shore, or the condition expresses symmetry in the solution. The wave velocity is in this case given by \( q = gH(x,y) \), where \( g \) is the acceleration of gravity and \( H(x,y) \) is the stillwater depth.
It can be wise to do Problem 21, because a 1D program corresponding to the present 2D problem allows much faster experimentation with parameters and effects.
The initial surface is taken as a smooth Gaussian function
$$ \begin{equation} I(x;I_0,I_a,I_m,I_s) = I_0 + I_a\exp{\left(-\left(\frac{x-I_m}{I_s}\right)^2\right)}, \end{equation} $$ with \( I_m=0 \) reflecting the location of the peak of \( I(x) \) and \( I_s \) being a measure of the width of the function \( I(x) \) (\( I_s \) is \( \sqrt{2} \) times the standard deviation of the familiar normal distribution curve).
Three different bottom shapes can be investigated. A 2D Gaussian hill can be modeled by
$$ \begin{equation} B(x;B_0,B_a,B_{mx}, B_{my} ,B_s, b) = B_0 + B_a\exp{\left(-\left(\frac{x-B_{mx}}{B_s}\right)^2 -\left(\frac{y-B_{my}}{bB_s}\right)^2\right)}, \label{wave:app:exer:tsunami2D:hill:Gauss} \end{equation} $$ where \( b \) is a scaling parameter: \( b=1 \) gives a circular Gaussian function with circular contour lines, while \( b\neq 1 \) gives an elliptic shape with elliptic contour lines.
A less smooth hill is modeled by the "cosine hat" function
$$ \begin{equation} B(x; B_0, B_a, B_{mx}, B_{my}, B_s) = B_0 + B_a\cos{\left( \pi\frac{x-B_{mx}}{2B_s}\right)} \cos{\left( \pi\frac{y-B_{my}}{2B_s}\right)}, \label{wave:app:exer:tsunami2D:hill:cohat} \end{equation} $$ when \( 0 \leq \sqrt{x^2+y^2} \leq B_s \) and \( B=B_0 \) outside this circle.
A more dramatic hill shape is a box:
$$ \begin{equation} B(x; B_0, B_a, B_m, B_s, b) = B_0 + B_a \label{wave:app:exer:tsunami2D:hill:box} \end{equation} $$ for \( x \) and \( y \) inside a rectangle $$ B_{mx}-B_s \leq x \leq B_{mx} + B_s,\quad B_{my}-bB_s \leq y \leq B_{my} + bB_s, $$ and \( B=B_0 \) outside this rectangle. The \( b \) parameter controls the rectangular shape of the cross section of the box.
Note that the initial condition and the listed bottom shapes are symmetric around the line \( y=B_{my} \). We therefore expect the surface elevation also to be symmetric with respect to this line. This means that we can halve the computational domain by working with \( [0,L_x]\times [0, B_{my}] \). Along the upper boundary, \( y=B_{my} \), we must impose the symmetry condition \( \partial \eta/\partial n=0 \). Such a symmetry condition (\( -\eta_x=0 \)) is also needed at the \( x=0 \) boundary because the initial condition has a symmetry here. At the lower boundary \( y=0 \) we also set a Neumann condition (which becomes \( -\eta_y=0 \)). The wave motion is to be simulated until the wave hits the reflecting boundaries where \( \partial\eta/\partial n =\eta_x =0 \).
Investigate how different hill shapes, different sizes of the water gap above the hill, and different resolutions \( \Delta x = \Delta y = h \) and \( \Delta t \) influence the numerical quality of the solution. One anticipates that the less smooth hill shapes will introduce more numerical noise. Presenting the results as movies of the surface elevation is effective.
The groups can select between one or more of the following tasks.
Harmonic means are often used if the coefficient \( q \) is non-smooth or discontinuous. Investigate if harmonic averaging of \( q \) works better than the arithmetic averging for the box-shaped subsea hill. The effect might not be big unless the water gap at the top of the hill is small. It can be wise to test the effect of harmonic averaging in 1D first.
Remark. With a small gap between the obstruction and the free surface, and with abrupt changes in the bottom shape, the model PDE does not necessarily describe the wave motion in an accurate or qualitatively correct way.
Create some fancy 3D visualization of the water waves and the subsea hill. Try to make the hill transparent. Suitable tools are
Implement an open boundary condition at \( x=L_x \), \( u_t + \sqrt{q}u_x=0 \), as suggested in Problem 20. This condition is only correct in 1D, but might work satisfactorily in 2D if the wave is approximately one-dimensional when it hits the boundary. See how well this condition works in letting the tsunami pass out of the domain. The distance from the subsea hill (which disturbes the wave) and the outlet boundary \( x=L_x \) is an important parameter.
Instead of using a condition \( u_t + \sqrt{q}u_x=0 \), which is exact only for plane waves propagating in \( x \) direction, one can add an artificial domain \( [L_x,L_x+\delta]\times [0,L_y] \) where waves are sufficiently damped and absorbed. The goal of an open boundary condition is to avoid waves being reflected back into the domain. Turn on the damping parameter \( b \) in \( [L_x,L_x+\delta]\times [0,L_y] \), and test if it is wise to vary \( b \), say in a linear or exponential fashion to have a smooth transition from \( b=0 \) in the physical domain and to some significant (efficient) \( b \) value towards the artificial boundary \( x=L_x+\delta \).
Extend the program with compiled loops using one or more of the following techniques:
f2pyf2pyscipy.weavesudo apt-get install fenics on
Ubuntu will install Instant) and it is described in the
FEniCS book.
Make a parallel version of the program using different approaches: