Finite difference methods for waves on a string
The complete initial-boundary value problem
Input data in the problem
Demo of a vibrating string (\( C=0.8 \))
Demo of a vibrating string (\( C=1.0012 \))
Step 1: Discretizing the domain
The discrete solution
Step 2: Fulfilling the equation at the mesh points
Step 3: Replacing derivatives by finite differences
Step 3: Algebraic version of the PDE
Step 3: Algebraic version of the initial conditions
Step 4: Formulating a recursive algorithm
The Courant number
The finite difference stencil
The stencil for the first time level
The algorithm
Moving finite difference stencil
Sketch of an implementation (1)
PDE solvers should save memory
Sketch of an implementation (2)
Verification
A slightly generalized model problem
Discrete model for the generalized model problem
Modified equation for the first time level
Using an analytical solution of physical significance
Manufactured solution: principles
Manufactured solution: example
Testing a manufactured solution
Constructing an exact solution of the discrete equations
Analytical work with the PDE problem
Analytical work with the discrete equations (1)
Analytical work with the discrete equations (1)
Testing with the exact discrete solution
Implementation
The algorithm
What do to with the solution?
Making a solver function
Verification: exact quadratic solution
Visualization: animating \( u(x,t) \)
Making movie files
Running a case
Implementation of the case
Resulting movie for \( C=0.8 \)
The benefits of scaling
Vectorization
Operations on slices of arrays
Test the understanding
Vectorization of finite difference schemes (1)
Vectorization of finite difference schemes (2)
Vectorized implementation in the solver function
Verification of the vectorized version
Efficiency measurements
Generalization: reflecting boundaries
Neumann boundary condition
Discretization of derivatives at the boundary (1)
Discretization of derivatives at the boundary (2)
Visualization of modified boundary stencil
Implementation of Neumann conditions
Moving finite difference stencil
Index set notation
Index set notation in code
Index sets in action (1)
Index sets in action (2)
Alternative implementation via ghost cells
Implementation of ghost cells (1)
Implementation of ghost cells (2)
Generalization: variable wave velocity
The model PDE with a variable coefficient
Discretizing the variable coefficient (1)
Discretizing the variable coefficient (2)
Discretizing the variable coefficient (3)
Computing the coefficient between mesh points
Discretization of variable-coefficient wave equation in operator notation
Neumann condition and a variable coefficient
Implementation of variable coefficients
A more general model PDE with variable coefficients
Generalization: damping
Building a general 1D wave equation solver
Collection of initial conditions
Finite difference methods for 2D and 3D wave equations
Examples on wave equations written out in 2D/3D
Boundary and initial conditions
Example: 2D propagation of Gaussian function
Mesh
Discretization
Special stencil for the first time step
Variable coefficients (1)
Variable coefficients (2)
Neumann boundary condition in 2D
Implementation of 2D/3D problems
Algorithm
Scalar computations: mesh
Scalar computations: arrays
Scalar computations: initial condition
Scalar computations: primary stencil
Vectorized computations: mesh coordinates
Vectorized computations: stencil
Verification: quadratic solution (1)
Verification: quadratic solution (2)
Migrating loops to Cython
Declaring variables and annotating the code
Cython version of the functions
Visual inspection of the C translation
Building the extension module
Calling the Cython function from Python
Migrating loops to Fortran
The Fortran subroutine
Building the Fortran module with f2py
How to avoid array copying
Efficiency of translating to Fortran
Migrating loops to C via Cython
The C code
The Cython interface file
Building the extension module
Migrating loops to C via f2py
The C code and the Fortran interface file
Building the extension module
Migrating loops to C++ via f2py
Analysis of the difference equations
Properties of the solution of the wave equation
Demo of the splitting of \( I(x) \) into two waves
Effect of variable wave velocity
What happens here?
Representation of waves as sum of sine/cosine waves
Analysis of the finite difference scheme
Preliminary results
Numerical wave propagation (1)
Numerical wave propagation (2)
Numerical wave propagation (3)
Why \( C\leq 1 \) is a stability criterion
Numerical dispersion relation
The special case \( C=1 \)
Computing the error in wave velocity
Visualizing the error in wave velocity
Taylor expanding the error in wave velocity
Example on effect of wrong wave velocity (1)
Example on effect of wrong wave velocity (1)
Extending the analysis to 2D (and 3D)
Discrete wave components in 2D
Stability criterion in 2D
Stability criterion in 3D
Numerical dispersion relation in 2D (1)
Numerical dispersion relation in 2D (2)
Numerical dispersion relation in 2D (3)
Waves on a string can be modeled by the wave equation $$ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2} $$
\( u(x,t) \) is the displacement of the string
Ooops!
Mesh in time: $$ \begin{equation} 0 = t_0 < t_1 < t_2 < \cdots < t_{N_t-1} < t_{N_t} = T \end{equation} $$
Mesh in space: $$ \begin{equation} 0 = x_0 < x_1 < x_2 < \cdots < x_{N_x-1} < x_{N_x} = L \end{equation} $$
Uniform mesh with constant mesh spacings \( \Delta t \) and \( \Delta x \): $$ \begin{equation} x_i = i\Delta x,\ i=0,\ldots,N_x,\quad t_i = n\Delta t,\ n=0,\ldots,N_t \end{equation} $$
Let the PDE be satisfied at all interior mesh points: $$ \begin{equation} \frac{\partial^2}{\partial t^2} u(x_i, t_n) = c^2\frac{\partial^2}{\partial x^2} u(x_i, t_n), \label{wave:pde1:step2} \end{equation} $$ for \( i=1,\ldots,N_x-1 \) and \( n=1,\ldots,N_t-1 \).
For \( n=0 \) we have the initial conditions \( u=I(x) \) and \( u_t=0 \), and at the boundaries \( i=0,N_x \) we have the boundary condition \( u=0 \).
Widely used finite difference formula for the second-order derivative: $$ \frac{\partial^2}{\partial t^2}u(x_i,t_n)\approx \frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\Delta t^2}= [D_tD_t u]^n_i$$
and $$ \frac{\partial^2}{\partial x^2}u(x_i,t_n)\approx \frac{u_{i+1}^{n} - 2u_i^n + u^{n}_{i-1}}{\Delta x^2} = [D_xD_x u]^n_i $$
Replace derivatives by differences: $$ \begin{equation} \frac{u_i^{n+1} - 2u_i^n + u^{n-1}_i}{\Delta t^2} = c^2\frac{u_{i+1}^{n} - 2u_i^n + u^{n}_{i-1}}{\Delta x^2}, \label{wave:pde1:step3b} \end{equation} $$
In operator notation: $$ \begin{equation} [D_tD_t u = c^2 D_xD_x]^{n}_i \label{wave:pde1:step3a} \end{equation} $$
The other initial condition \( u(x,0)=I(x) \) can be computed by $$ u_i^0 = I(x_i),\quad i=0,\ldots,N_x$$
There is only one parameter, \( C \), in the discrete model: \( C \) lumps mesh parameters \( \Delta t \) and \( \Delta x \) with the only physical parameter, the wave velocity \( c \). The value \( C \) and the smoothness of \( I(x) \) govern the quality of the numerical solution.
Insert in stencil \( [D_tD_tu = c^2D_xD_x]^0_i \) to get $$ \begin{equation} u_i^1 = u^0_i - \half C^2\left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\right) \label{wave:pde1:step4:1} \end{equation} $$
web page or a movie file.
u[i] stores \( u^{n+1}_i \)u_1[i] stores \( u^n_i \)u_2[i] stores \( u^{n-1}_i \)
u is the unknown to be computed (a spatial mesh
function), u_k is the computed spatial mesh function k
time steps back in time.
The algorithm only needs to access the three most recent time levels, so we need only three arrays for \( u_i^{n+1} \), \( u_i^n \), and \( u_i^{n-1} \), \( i=0,\ldots,N_x \). Storing all the solutions in a two-dimensional array of size \( (N_x+1)\times (N_t+1) \) would be possible in this simple one-dimensional PDE problem, but not in large 2D problems and not even in small 3D problems.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 | # Given mesh points as arrays x and t (x[i], t[n])
dx = x[1] - x[0]
dt = t[1] - t[0]
C = c*dt/dx # Courant number
Nt = len(t)-1
C2 = C**2 # Help variable in the scheme
# Set initial condition u(x,0) = I(x)
for i in range(0, Nx+1):
u_1[i] = I(x[i])
# Apply special formula for first step, incorporating du/dt=0
for i in range(1, Nx):
u[i] = u_1[i] - 0.5*C**2(u_1[i+1] - 2*u_1[i] + u_1[i-1])
u[0] = 0; u[Nx] = 0 # Enforce boundary conditions
# Switch variables before next step
u_2[:], u_1[:] = u_1, u
for n in range(1, Nt):
# Update all inner mesh points at time t[n+1]
for i in range(1, Nx):
u[i] = 2u_1[i] - u_2[i] - \
C**2(u_1[i+1] - 2*u_1[i] + u_1[i-1])
# Insert boundary conditions
u[0] = 0; u[Nx] = 0
# Switch variables before next step
u_2[:], u_1[:] = u_1, u
|
Add source term \( f \) and nonzero initial condition \( u_t(x,0) \): $$ \begin{align} u_{tt} &= c^2 u_{xx} + f(x,t), \label{wave:pde2}\\ u(x,0) &= I(x), \quad &x\in [0,L] \label{wave:pde2:ic:u}\\ u_t(x,0) &= V(x), \quad &x\in [0,L] \label{wave:pde2:ic:ut}\\ u(0,t) & = 0, \quad & t>0, \label{wave:pde2:bc:0}\\ u(L,t) & = 0, \quad &t>0 \label{wave:pde2:bc:L} \end{align} $$
Writing out and solving for the unknown \( u^{n+1}_i \): $$ \begin{equation} u^{n+1}_i = -u^{n-1}_i + 2u^n_i + C^2 (u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}) + \Delta t^2 f^n_i \label{wave:pde2:step3b} \end{equation} $$
Centered difference for \( u_t(x,0) = V(x) \): $$ [D_{2t}u = V]^0_i\quad\Rightarrow\quad u^{-1}_i = u^{1}_i - 2\Delta t V_i,$$
Inserting this in the stencil \eqref{wave:pde2:step3b} for \( n=0 \) leads to $$ \begin{equation} u^{1}_i = u^0_i - \Delta t V_i + {\half} C^2 \left(u^{n}_{i+1}-2u^{n}_{i} + u^{n}_{i-1}\right) + \half\Delta t^2 f^n_i \label{wave:pde2:step3c} \end{equation} $$
PDE \( u_{tt}=c^2u_{xx}+f \): $$ -x(L-x)\sin t = -2\sin t + f\quad\Rightarrow f = (2 - x(L-x))\sin t$$
Implied initial conditions: $$ \begin{align*} u(x,0) &= I(x) = 0\\ u_t(x,0) &= V(x) = - x(L-x) \end{align*} $$
Boundary conditions: $$ u(x,0) = u(x,L) = 0 $$
Here, choose \( \uex \) such that \( \uex(x,0)=\uex(L,0)=0 \): $$ \uex (x,t) = x(L-x)(1+{\half}t), $$
Insert in the PDE and find \( f \): $$ f(x,t)=2(1+t)c^2$$
Initial conditions: $$ I(x) = x(L-x),\quad V(x)={\half}x(L-x) $$
We want to show that \( \uex \) also solves the discrete equations!
Useful preliminary result: $$ \begin{align} \lbrack D_tD_t t^2\rbrack^n &= \frac{t_{n+1}^2 - 2t_n^2 + t_{n-1}^2}{\Delta t^2} = (n+1)^2 -n^2 + (n-1)^2 = 2\\ \lbrack D_tD_t t\rbrack^n &= \frac{t_{n+1} - 2t_n + t_{n-1}}{\Delta t^2} = \frac{((n+1) -n + (n-1))\Delta t}{\Delta t^2} = 0 \end{align} $$
Hence, $$ [D_tD_t \uex]^n_i = x_i(L-x_i)[D_tD_t (1+{\half}t)]^n = x_i(L-x_i){\half}[D_tD_t t]^n = 0$$
Now, \( f^n_i = 2(1+{\half}t_n)c^2 \) and we get $$ [D_tD_t \uex - c^2D_xD_x\uex - f]^n_i = 0 - c^2(-1)2(1 + {\half}t_n + 2(1+{\half}t_n)c^2 = 0$$
Moreover, \( \uex(x_i,0)=I(x_i) \), \( \partial \uex/\partial t = V(x_i) \) at \( t=0 \), and \( \uex(x_0,t)=\uex(x_{N_x},0)=0 \). Also the modified scheme for the first time step is fulfilled by \( \uex(x_i,t_n) \).
1 2 3 4 | def user_action(u, x, t, n):
# u[i] at spatial mesh points x[i] at time t[n]
# plot u
# or store u
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 | from numpy import *
def solver(I, V, f, c, L, dt, C, T, user_action=None):
"""Solve u_tt=c^2*u_xx + f on (0,L)x(0,T]."""
Nt = int(round(T/dt))
t = linspace(0, Nt*dt, Nt+1) # Mesh points in time
dx = dt*c/float(C)
Nx = int(round(L/dx))
x = linspace(0, L, Nx+1) # Mesh points in space
C2 = C**2 # Help variable in the scheme
if f is None or f == 0 :
f = lambda x, t: 0
if V is None or V == 0:
V = lambda x: 0
u = zeros(Nx+1) # Solution array at new time level
u_1 = zeros(Nx+1) # Solution at 1 time level back
u_2 = zeros(Nx+1) # Solution at 2 time levels back
import time; t0 = time.clock() # for measuring CPU time
# Load initial condition into u_1
for i in range(0,Nx+1):
u_1[i] = I(x[i])
if user_action is not None:
user_action(u_1, x, t, 0)
# Special formula for first time step
n = 0
for i in range(1, Nx):
u[i] = u_1[i] + dt*V(x[i]) + \
0.5*C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \
0.5*dt**2*f(x[i], t[n])
u[0] = 0; u[Nx] = 0
if user_action is not None:
user_action(u, x, t, 1)
# Switch variables before next step
u_2[:] = u_1; u_1[:] = u
for n in range(1, Nt):
# Update all inner points at time t[n+1]
for i in range(1, Nx):
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1]) + \
dt**2*f(x[i], t[n])
# Insert boundary conditions
u[0] = 0; u[Nx] = 0
if user_action is not None:
if user_action(u, x, t, n+1):
break
# Switch variables before next step
u_2[:] = u_1; u_1[:] = u
cpu_time = t0 - time.clock()
return u, x, t, cpu_time
|
Exact solution of the PDE problem and the discrete equations: \( \uex (x,t) = x(L-x)(1+{\half}t) \)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | import nose.tools as nt
def test_quadratic():
"""Check that u(x,t)=x(L-x)(1+t/2) is exactly reproduced."""
def u_exact(x, t):
return x*(L-x)*(1 + 0.5*t)
def I(x):
return u_exact(x, 0)
def V(x):
return 0.5*u_exact(x, 0)
def f(x, t):
return 2*(1 + 0.5*t)*c**2
L = 2.5
c = 1.5
C = 0.75
Nx = 3 # Very coarse mesh for this exact test
dt = C*(L/Nx)/c
T = 18
u, x, t, cpu = solver(I, V, f, c, L, dt, C, T)
u_e = u_exact(x, t[-1])
diff = abs(u - u_e).max()
nt.assert_almost_equal(diff, 0, places=14)
|
Make a viz function for animating the curve, with plotting
in a user_action function plot_u:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 | def viz(I, V, f, c, L, dt, C, T, umin, umax, animate=True):
"""Run solver and visualize u at each time level."""
import scitools.std as plt
import time, glob, os
def plot_u(u, x, t, n):
"""user_action function for solver."""
plt.plot(x, u, 'r-',
xlabel='x', ylabel='u',
axis=[0, L, umin, umax],
title='t=%f' % t[n], show=True)
# Let the initial condition stay on the screen for 2
# seconds, else insert a pause of 0.2 s between each plot
time.sleep(2) if t[n] == 0 else time.sleep(0.2)
plt.savefig('frame_%04d.png' % n) # for movie making
# Clean up old movie frames
for filename in glob.glob('frame_*.png'):
os.remove(filename)
user_action = plot_u if animate else None
u, x, t, cpu = solver(I, V, f, c, L, dt, C, T, user_action)
# Make movie files
fps = 4 # Frames per second
plt.movie('frame_*.png', encoder='html', fps=fps,
output_file='movie.html')
codec2ext = dict(flv='flv', libx264='mp4', libvpx='webm',
libtheora='ogg')
filespec = 'frame_%04d.png'
movie_program = 'avconv' # or 'ffmpeg'
for codec in codec2ext:
ext = codec2ext[codec]
cmd = '%(movie_program)s -r %(fps)d -i %(filespec)s '\
'-vcodec %(codec)s movie.%(ext)s' % vars()
os.system(cmd)
|
Note: plot_u is function inside function and remembers the
local variables in viz (known as a closure).
'something_%04d.png' % frame_counter
1 2 3 4 5 6 | Terminal> scitools movie encoder=html output_file=movie.html \
fps=4 frame_*.png # web page with a player
Terminal> avconv -r 4 -i frame_%04d.png -c:v flv movie.flv
Terminal> avconv -r 4 -i frame_%04d.png -c:v libtheora movie.ogg
Terminal> avconv -r 4 -i frame_%04d.png -c:v libx264 movie.mp4
Terminal> avconv -r 4 -i frame_%04d.png -c:v libpvx movie.webm
|
%04d) is essential for correct sequence of
frames in something_*.png (Unix alphanumeric sort)frame_*.png files before making a new movie
Appropriate data:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | def guitar(C):
"""Triangular wave (pulled guitar string)."""
L = 0.75
x0 = 0.8*L
a = 0.005
freq = 440
wavelength = 2*L
c = freq*wavelength
omega = 2*pi*freq
num_periods = 1
T = 2*pi/omega*num_periods
# Choose dt the same as the stability limit for Nx=50
dt = L/50./c
def I(x):
return a*x/x0 if x < x0 else a/(L-x0)*(L-x)
umin = -1.2*a; umax = -umin
cpu = viz(I, 0, 0, c, L, dt, C, T, umin, umax, animate=True)
|
Program: wave1D_u0.py.
Insert this in the PDE (with \( f=0 \)) and dropping bars $$ u_{tt} = u_{xx}$$
Initial condition: set \( a=1 \), \( L=1 \), and \( x_0\in [0,1] \) in \eqref{wave:pde1:guitar:I}.
In the code: set a=c=L=1, x0=0.8, and there is no need to calculate with
wavelengths and frequencies to estimate \( c \)!
Just one challenge: determine the period of the waves and an appropriate end time (see the text for details).
numpy) code instead of explicit loops
1 2 3 | n = u.size
for i in range(0, n-1):
d[i] = u[i+1] - u[i]
|
numpy code: u[1:n] - u[0:n-1] or just u[1:] - u[:-1]
Newcomers to vectorization are encouraged to choose
a small array u, say with five elements,
and simulate with pen and paper
both the loop version and the vectorized version.
Finite difference schemes basically contains differences between array elements with shifted indices. Consider the updating formula
1 2 | for i in range(1, n-1):
u2[i] = u[i-1] - 2*u[i] + u[i+1]
|
The vectorization consists of replacing the loop by arithmetics on
slices of arrays of length n-2:
1 2 | u2 = u[:-2] - 2*u[1:-1] + u[2:]
u2 = u[0:n-2] - 2*u[1:n-1] + u[2:n] # alternative
|
Note: u2 gets length n-2.
If u2 is already an array of length n, do update on "inner" elements
1 2 | u2[1:-1] = u[:-2] - 2*u[1:-1] + u[2:]
u2[1:n-1] = u[0:n-2] - 2*u[1:n-1] + u[2:n] # alternative
|
Include a function evaluation too:
1 2 3 4 5 6 7 8 9 | def f(x):
return x**2 + 1
# Scalar version
for i in range(1, n-1):
u2[i] = u[i-1] - 2*u[i] + u[i+1] + f(x[i])
# Vectorized version
u2[1:-1] = u[:-2] - 2*u[1:-1] + u[2:] + f(x[1:-1])
|
Scalar loop:
1 2 3 | for i in range(1, Nx):
u[i] = 2*u_1[i] - u_2[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
|
Vectorized loop:
1 2 | u[1:-1] = - u_2[1:-1] + 2*u_1[1:-1] + \
C2*(u_1[:-2] - 2*u_1[1:-1] + u_1[2:])
|
or
1 2 | u[1:Nx] = 2*u_1[1:Nx]- u_2[1:Nx] + \
C2*(u_1[0:Nx-1] - 2*u_1[1:Nx] + u_1[2:Nx+1])
|
Program: wave1D_u0v.py
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 | def test_quadratic():
"""
Check the scalar and vectorized versions work for
a quadratic u(x,t)=x(L-x)(1+t/2) that is exactly reproduced.
"""
# The following function must work for x as array or scalar
u_exact = lambda x, t: x*(L - x)*(1 + 0.5*t)
I = lambda x: u_exact(x, 0)
V = lambda x: 0.5*u_exact(x, 0)
# f is a scalar (zeros_like(x) works for scalar x too)
f = lambda x, t: zeros_like(x) + 2*c**2*(1 + 0.5*t)
L = 2.5
c = 1.5
C = 0.75
Nx = 3 # Very coarse mesh for this exact test
dt = C*(L/Nx)/c
T = 18
def assert_no_error(u, x, t, n):
u_e = u_exact(x, t[n])
diff = abs(u - u_e).max()
nt.assert_almost_equal(diff, 0, places=13)
solver(I, V, f, c, L, dt, C, T,
user_action=assert_no_error, version='scalar')
solver(I, V, f, c, L, dt, C, T,
user_action=assert_no_error, version='vectorized')
|
Note:
wave1D_u0v.py for \( N_x \) as 50,100,200,400,800
and measuring the CPU time
For a 1D domain \( [0,L] \): $$ \left.\frac{\partial}{\partial n}\right\vert_{x=L} = \frac{\partial}{\partial x},\quad \left.\frac{\partial}{\partial n}\right\vert_{x=0} = - \frac{\partial}{\partial x} $$
Boundary condition terminology:
Discrete equation for computing \( u^3_0 \) in terms of \( u^2_0 \), \( u^1_0 \), and \( u^2_1 \):
Animation in a web page or a movie file.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | i = 0
ip1 = i+1
im1 = ip1 # i-1 -> i+1
u[i] = u_1[i] + C2*(u_1[im1] - 2*u_1[i] + u_1[ip1])
i = Nx
im1 = i-1
ip1 = im1 # i+1 -> i-1
u[i] = u_1[i] + C2*(u_1[im1] - 2*u_1[i] + u_1[ip1])
# Or just one loop over all points
for i in range(0, Nx+1):
ip1 = i+1 if i < Nx else i-1
im1 = i-1 if i > 0 else i+1
u[i] = u_1[i] + C2*(u_1[im1] - 2*u_1[i] + u_1[ip1])
|
Program wave1D_dn0.py
web page or a movie file.
| Notation | Python |
|---|---|
| \( \Ix \) | Ix |
| \( \setb{\Ix} \) | Ix[0] |
| \( \sete{\Ix} \) | Ix[-1] |
| \( \setl{\Ix} \) | Ix[1:] |
| \( \setr{\Ix} \) | Ix[:-1] |
| \( \seti{\Ix} \) | Ix[1:-1] |
Index sets for a problem in the \( x,t \) plane: $$ \begin{equation} \Ix = \{0,\ldots,N_x\},\quad \It = \{0,\ldots,N_t\}, \end{equation} $$
Implemented in Python as
1 2 | Ix = range(0, Nx+1)
It = range(0, Nt+1)
|
A finite difference scheme can with the index set notation be specified as $$ \begin{align*} u^{n+1}_i &= -u^{n-1}_i + 2u^n_i + C^2 \left(u^{n}_{i+1}-2u^{n}_{i}+u^{n}_{i-1}\right), \quad i\in\seti{\Ix},\ n\in\seti{\It}\\ u_i &= 0, \quad i=\setb{\Ix},\ n\in\seti{\It}\\ u_i &= 0, \quad i=\sete{\Ix},\ n\in\seti{\It} \end{align*} $$
Corresponding implementation:
1 2 3 4 5 6 | for n in It[1:-1]:
for i in Ix[1:-1]:
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
i = Ix[0]; u[i] = 0
i = Ix[-1]; u[i] = 0
|
Program wave1D_dn.py
Add ghost points:
1 2 3 4 5 | u = zeros(Nx+3)
u_1 = zeros(Nx+3)
u_2 = zeros(Nx+3)
x = linspace(0, L, Nx+1) # Mesh points without ghost points
|
u[-1] will always mean the last element in u0,..,Nx+2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | u = zeros(Nx+3)
Ix = range(1, u.shape[0]-1)
# Boundary values: u[Ix[0]], u[Ix[-1]]
# Set initial conditions
for i in Ix:
u_1[i] = I(x[i-Ix[0]]) # Note i-Ix[0]
# Loop over all physical mesh points
for i in Ix:
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(u_1[i-1] - 2*u_1[i] + u_1[i+1])
# Update ghost values
i = Ix[0] # x=0 boundary
u[i-1] = u[i+1]
i = Ix[-1] # x=L boundary
u[i-1] = u[i+1]
|
Program: wave1D_dn0_ghost.py.
Heterogeneous media: varying \( c=c(x) \)
This equation sampled at a mesh point \( (x_i,t_n) \): $$ \frac{\partial^2 }{\partial t^2} u(x_i,t_n) = \frac{\partial}{\partial x}\left( q(x_i) \frac{\partial}{\partial x} u(x_i,t_n)\right) + f(x_i,t_n), $$
The principal idea is to first discretize the outer derivative.
Define $$ \phi = q(x) \frac{\partial u}{\partial x} $$
Then use a centered derivative around \( x=x_i \) for the derivative of \( \phi \): $$ \left[\frac{\partial\phi}{\partial x}\right]^n_i \approx \frac{\phi_{i+\half} - \phi_{i-\half}}{\Delta x} = [D_x\phi]^n_i $$
Then discretize the inner operators: $$ \phi_{i+\half} = q_{i+\half} \left[\frac{\partial u}{\partial x}\right]^n_{i+\half} \approx q_{i+\half} \frac{u^n_{i+1} - u^n_{i}}{\Delta x} = [q D_x u]_{i+\half}^n $$
Similarly, $$ \phi_{i-\half} = q_{i-\half} \left[\frac{\partial u}{\partial x}\right]^n_{i-\half} \approx q_{i-\half} \frac{u^n_{i} - u^n_{i-1}}{\Delta x} = [q D_x u]_{i-\half}^n $$
These intermediate results are now combined to $$ \begin{equation} \left[ \frac{\partial}{\partial x}\left( q(x) \frac{\partial u}{\partial x}\right)\right]^n_i \approx \frac{1}{\Delta x^2} \left( q_{i+\half} \left({u^n_{i+1} - u^n_{i}}\right) - q_{i-\half} \left({u^n_{i} - u^n_{i-1}}\right)\right) \label{wave:pde2:var:c:formula} \end{equation} $$
In operator notation: $$ \begin{equation} \left[ \frac{\partial}{\partial x}\left( q(x) \frac{\partial u}{\partial x}\right)\right]^n_i \approx [D_xq D_x u]^n_i \label{wave:pde2:var:c:formula:op} \end{equation} $$
Many are tempted to use the chain rule on the term \( \frac{\partial}{\partial x}\left( q(x) \frac{\partial u}{\partial x}\right) \), but this is not a good idea!
The arithmetic mean in \eqref{wave:pde2:var:c:mean:arithmetic} is by far the most used averaging technique.
We clearly see the type of finite differences and averaging!
Write out and solve wrt \( u_i^{n+1} \): $$ \begin{align} u^{n+1}_i &= - u_i^{n-1} + 2u_i^n + \left(\frac{\Delta x}{\Delta t}\right)^2\times \nonumber\\ &\quad \left( \half(q_{i} + q_{i+1})(u_{i+1}^n - u_{i}^n) - \half(q_{i} + q_{i-1})(u_{i}^n - u_{i-1}^n)\right) + \nonumber\\ & \quad \Delta t^2 f^n_i \label{wave:pde2:var:c:scheme:impl} \end{align} $$
Consider \( \partial u/\partial x=0 \) at \( x=L=N_x\Delta x \): $$ \frac{u_{i+1}^{n} - u_{i-1}^n}{2\Delta x} = 0\quad u_{i+1}^n = u_{i-1}^n, \quad i=N_x $$
Insert \( u_{i+1}^n=u_{i-1}^n \) in the stencil \eqref{wave:pde2:var:c:scheme:impl} for \( i=N_x \) and obtain $$ u^{n+1}_i \approx - u_i^{n-1} + 2u_i^n + \left(\frac{\Delta x}{\Delta t}\right)^2 2q_{i}(u_{i-1}^n - u_{i}^n) + \Delta t^2 f^n_i $$
(We have used \( q_{i+\half} + q_{i-\half}\approx 2q_i \).)
Alternative: assume \( dq/dx=0 \) (simpler).
Assume c[i] holds \( c_i \) the spatial mesh points
1 2 3 4 5 | for i in range(1, Nx):
u[i] = - u_2[i] + 2*u_1[i] + \
C2*(0.5*(q[i] + q[i+1])*(u_1[i+1] - u_1[i]) - \
0.5*(q[i] + q[i-1])*(u_1[i] - u_1[i-1])) + \
dt2*f(x[i], t[n])
|
Here: C2=(dt/dx)**2
Vectorized version:
1 2 3 4 | u[1:-1] = - u_2[1:-1] + 2*u_1[1:-1] + \
C2*(0.5*(q[1:-1] + q[2:])*(u_1[2:] - u_1[1:-1]) -
0.5*(q[1:-1] + q[:-2])*(u_1[1:-1] - u_1[:-2])) + \
dt2*f(x[1:-1], t[n])
|
Neumann condition \( u_x=0 \): same ideas as in 1D (modified stencil or ghost cells).
A natural scheme is $$ \begin{equation} [\varrho D_tD_t u = D_x\overline{q}^xD_x u + f]^n_i \end{equation} $$
Or $$ \begin{equation} [D_tD_t u = \varrho^{-1}D_x\overline{q}^xD_x u + f]^n_i \end{equation} $$
No need to average \( \varrho \), just sample at \( i \)
Why do waves die out?
\( b \geq 0 \): prescribed damping coefficient.
Discretization via centered differences to ensure \( \Oof{\Delta t^2} \) error: $$ \begin{equation} [D_tD_t u + bD_{2t}u = c^2D_xD_x u + f]^n_i \label{wave:pde3:fd} \end{equation} $$
Need special formula for \( u^1_i \) + special stencil (or ghost cells) for Neumann conditions.
The program wave1D_dn_vc.py solves a fairly general 1D wave equation: $$ \begin{align} u_t &= (c^2(x)u_x)_x + f(x,t),\quad &x\in (0,L),\ t\in (0,T] \label{wave:pde2:software:ueq}\\ u(x,0) &= I(x),\quad &x\in [0,L] \label{wave:pde2:software:u0}\\ u_t(x,0) &= V(t),\quad &x\in [0,L] \label{wave:pde2:software:ut0}\\ u(0,t) &= U_0(t)\hbox{ or } u_x(0,t)=0,\quad &t\in (0,T] \label{wave:pde2:software:bc0}\\ u(L,t) &= U_L(t)\hbox{ or } u_x(L,t)=0,\quad &t\in (0,T] \label{wave:pde2:software:bcL} \end{align} $$
Can be adapted to many needs.
The function pulse in wave1D_dn_vc.py offers four
initial conditions:
gaussian)
1 2 | >>> import wave1D_dn_vc as w
>>> w.pulse(loc='left', pulse_tp='cosinehat', Nx=50, every_frame=10)
|
Constant wave velocity \( c \): $$ \begin{equation} \frac{\partial^2 u}{\partial t^2} = c^2\nabla^2 u\hbox{ for }\xpoint\in\Omega\subset\Real^d,\ t\in (0,T] \label{wave:2D3D:model1} \end{equation} $$
Variable wave velocity: $$ \begin{equation} \varrho\frac{\partial^2 u}{\partial t^2} = \nabla\cdot (q\nabla u) + f\hbox{ for }\xpoint\in\Omega\subset\Real^d,\ t\in (0,T] \label{wave:2D3D:model2} \end{equation} $$
3D, constant \( c \): $$ \begin{equation*} \nabla^2 u = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \end{equation*} $$
2D, variable \( c \): $$ \begin{equation} \varrho(x,y) \frac{\partial^2 u}{\partial t^2} = \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} $$
Compact notation: $$ \begin{align} u_{tt} &= c^2(u_{xx} + u_{yy} + u_{zz}) + f, \label{wave:2D3D:model1:v2}\\ \varrho u_{tt} &= (q u_x)_x + (q u_z)_z + (q u_z)_z + f \label{wave:2D3D:model2:v2} \end{align} $$
We need one boundary condition at each point on \( \partial\Omega \):
\( u^{n-1}_{i,j} \) and \( u^n_{i,j} \) are known, solve for \( u^{n+1}_{i,j} \): $$ u^{n+1}_{i,j} = 2u^n_{i,j} + u^{n-1}_{i,j} + c^2\Delta t^2[D_xD_x u + D_yD_y u]^n_{i,j}$$
3D wave equation: $$ \varrho u_{tt} = (qu_x)_x + (qu_y)_y + (qu_z)_z + f(x,y,z,t) $$
Just apply the 1D discretization for each term: $$ \begin{equation} [\varrho D_tD_t u = (D_x\overline{q}^x D_x u + D_y\overline{q}^y D_yu + D_z\overline{q}^z D_z u) + f]^n_{i,j,k} \end{equation} $$
Need special formula for \( u^1_{i,j,k} \) (use \( [D_{2t}u=V]^0 \) and stencil for \( n=0 \)).
Written out: $$ \begin{align*} u^{n+1}_{i,j,k} &= - u^{n-1}_{i,j,k} + 2u^{n}_{i,j,k} + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \half(q_{i,j,k} + q_{i+1,j,k})(u^{n}_{i+1,j,k} - u^{n}_{i,j,k}) - \\ &\qquad\quad \half(q_{i-1,j,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i-1,j,k})) + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \half(q_{i,j,k} + q_{i,j+1,k})(u^{n}_{i,j+1,k} - u^{n}_{i,j,k}) - \\ &\qquad\quad\half(q_{i,j-1,k} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j-1,k})) + \\ &= \frac{1}{\varrho_{i,j,k}}\frac{1}{\Delta x^2} ( \half(q_{i,j,k} + q_{i,j,k+1})(u^{n}_{i,j,k+1} - u^{n}_{i,j,k}) -\\ &\qquad\quad \half(q_{i,j,k-1} + q_{i,j,k})(u^{n}_{i,j,k} - u^{n}_{i,j,k-1})) + \\ + &\qquad \Delta t^2 f^n_{i,j,k} \end{align*} $$
Use ideas from 1D! Example: \( \frac{\partial u}{\partial n} \) at \( y=0 \), \( \frac{\partial u}{\partial n} = -\frac{\partial u}{\partial y} \)
Boundary condition discretization: $$ [-D_{2y} u = 0]^n_{i,0}\quad\Rightarrow\quad \frac{u^n_{i,1}-u^n_{i,-1}}{2\Delta y} = 0,\ i\in\Ix $$
Insert \( u^n_{i,-1}=u^n_{i,1} \) in the stencil for \( u^{n+1}_{i,j=0} \) to obtain a modified stencil on the boundary.
Pattern: use interior stencil also on the bundary, but replace \( j-1 \) by \( j+1 \)
Alternative: use ghost cells and ghost values
\( \Omega = [0,L_x]\times [0,L_y] \)
Discretization: $$ [D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j}, $$
Program: wave2D_u0.py
1 2 | def solver(I, V, f, c, Lx, Ly, Nx, Ny, dt, T,
user_action=None, version='scalar'):
|
Mesh:
1 2 3 4 5 6 7 8 | x = linspace(0, Lx, Nx+1) # mesh points in x dir
y = linspace(0, Ly, Ny+1) # mesh points in y dir
dx = x[1] - x[0]
dy = y[1] - y[0]
Nt = int(round(T/float(dt)))
t = linspace(0, N*dt, N+1) # mesh points in time
Cx2 = (c*dt/dx)**2; Cy2 = (c*dt/dy)**2 # help variables
dt2 = dt**2
|
Store \( u^{n+1}_{i,j} \), \( u^{n}_{i,j} \), and \( u^{n-1}_{i,j} \) in three two-dimensional arrays:
1 2 3 | u = zeros((Nx+1,Ny+1)) # solution array
u_1 = zeros((Nx+1,Ny+1)) # solution at t-dt
u_2 = zeros((Nx+1,Ny+1)) # solution at t-2*dt
|
\( u^{n+1}_{i,j} \) corresponds to u[i,j], etc.
1 2 3 4 5 6 7 8 9 10 | Ix = range(0, u.shape[0])
Iy = range(0, u.shape[1])
It = range(0, t.shape[0])
for i in Ix:
for j in Iy:
u_1[i,j] = I(x[i], y[j])
if user_action is not None:
user_action(u_1, x, xv, y, yv, t, 0)
|
Arguments xv and yv: for vectorized computations
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 | def advance_scalar(u, u_1, u_2, f, x, y, t, n, Cx2, Cy2, dt2,
V=None, step1=False):
Ix = range(0, u.shape[0]); Iy = range(0, u.shape[1])
if step1:
dt = sqrt(dt2) # save
Cx2 = 0.5*Cx2; Cy2 = 0.5*Cy2; dt2 = 0.5*dt2 # redefine
D1 = 1; D2 = 0
else:
D1 = 2; D2 = 1
for i in Ix[1:-1]:
for j in Iy[1:-1]:
u_xx = u_1[i-1,j] - 2*u_1[i,j] + u_1[i+1,j]
u_yy = u_1[i,j-1] - 2*u_1[i,j] + u_1[i,j+1]
u[i,j] = D1*u_1[i,j] - D2*u_2[i,j] + \
Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n])
if step1:
u[i,j] += dt*V(x[i], y[j])
# Boundary condition u=0
j = Iy[0]
for i in Ix: u[i,j] = 0
j = Iy[-1]
for i in Ix: u[i,j] = 0
i = Ix[0]
for j in Iy: u[i,j] = 0
i = Ix[-1]
for j in Iy: u[i,j] = 0
return u
|
D1 and D2: allow advance_scalar to be used also for \( u^1_{i,j} \):
1 2 | u = advance_scalar(u, u_1, u_2, f, x, y, t,
n, 0.5*Cx2, 0.5*Cy2, 0.5*dt2, D1=1, D2=0)
|
Mesh with \( 30\times 30 \) cells: vectorization reduces the CPU time by a factor of 70 (!).
Need special coordinate arrays xv and yv such that \( I(x,y) \)
and \( f(x,y,t) \) can be vectorized:
1 2 3 4 5 6 | from numpy import newaxis
xv = x[:,newaxis]
yv = y[newaxis,:]
u_1[:,:] = I(xv, yv)
f_a[:,:] = f(xv, yv, t)
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | def advance_vectorized(u, u_1, u_2, f_a, Cx2, Cy2, dt2,
V=None, step1=False):
if step1:
dt = sqrt(dt2) # save
Cx2 = 0.5*Cx2; Cy2 = 0.5*Cy2; dt2 = 0.5*dt2 # redefine
D1 = 1; D2 = 0
else:
D1 = 2; D2 = 1
u_xx = u_1[:-2,1:-1] - 2*u_1[1:-1,1:-1] + u_1[2:,1:-1]
u_yy = u_1[1:-1,:-2] - 2*u_1[1:-1,1:-1] + u_1[1:-1,2:]
u[1:-1,1:-1] = D1*u_1[1:-1,1:-1] - D2*u_2[1:-1,1:-1] + \
Cx2*u_xx + Cy2*u_yy + dt2*f_a[1:-1,1:-1]
if step1:
u[1:-1,1:-1] += dt*V[1:-1, 1:-1]
# Boundary condition u=0
j = 0
u[:,j] = 0
j = u.shape[1]-1
u[:,j] = 0
i = 0
u[i,:] = 0
i = u.shape[0]-1
u[i,:] = 0
return u
|
Manufactured solution: $$ \begin{equation} \uex(x,y,t) = x(L_x-x)y(L_y-y)(1+{\half}t) \label{wave2D3D:impl:verify:quadratic} \end{equation} $$
Requires \( f=2c^2(1+{\half}t)(y(L_y-y) + x(L_x-x)) \).
This \( \uex \) is ideal because it also solves the discrete equations!
Pure Python code:
1 2 3 4 5 6 7 8 9 | def advance_scalar(u, u_1, u_2, f, x, y, t,
n, Cx2, Cy2, dt2, D1=2, D2=1):
Ix = range(0, u.shape[0]); Iy = range(0, u.shape[1])
for i in Ix[1:-1]:
for j in Iy[1:-1]:
u_xx = u_1[i-1,j] - 2*u_1[i,j] + u_1[i+1,j]
u_yy = u_1[i,j-1] - 2*u_1[i,j] + u_1[i,j+1]
u[i,j] = D1*u_1[i,j] - D2*u_2[i,j] + \
Cx2*u_xx + Cy2*u_yy + dt2*f(x[i], y[j], t[n])
|
.pyx extension.function(a, b) \( \rightarrow \) cpdef function(int a, double b)v = 1.2 \( \rightarrow \) cdef double v = 1.2np.ndarray[np.float64_t, ndim=2, mode='c'] u
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | import numpy as np
cimport numpy as np
cimport cython
ctypedef np.float64_t DT # data type
@cython.boundscheck(False) # turn off array bounds check
@cython.wraparound(False) # turn off negative indices (u[-1,-1])
cpdef advance(
np.ndarray[DT, ndim=2, mode='c'] u,
np.ndarray[DT, ndim=2, mode='c'] u_1,
np.ndarray[DT, ndim=2, mode='c'] u_2,
np.ndarray[DT, ndim=2, mode='c'] f,
double Cx2, double Cy2, double dt2):
cdef int Nx, Ny, i, j
cdef double u_xx, u_yy
Nx = u.shape[0]-1
Ny = u.shape[1]-1
for i in xrange(1, Nx):
for j in xrange(1, Ny):
u_xx = u_1[i-1,j] - 2*u_1[i,j] + u_1[i+1,j]
u_yy = u_1[i,j-1] - 2*u_1[i,j] + u_1[i,j+1]
u[i,j] = 2*u_1[i,j] - u_2[i,j] + \
Cx2*u_xx + Cy2*u_yy + dt2*f[i,j]
|
Note: from now in we skip the code for setting boundary values
See how effective Cython can translate this code to C:
1 | Terminal> cython -a wave2D_u0_loop_cy.pyx
|
Load wave2D_u0_loop_cy.html in a browser (white: pure C, yellow: still Python):
Can click on wave2D_u0_loop_cy.c to see the generated C code...
.so file) that can be
loaded as a standard Python modulesetup.py script to build the extension module
1 2 3 4 5 6 7 8 9 10 | from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
cymodule = 'wave2D_u0_loop_cy'
setup(
name=cymodule
ext_modules=[Extension(cymodule, [cymodule + '.pyx'],)],
cmdclass={'build_ext': build_ext},
)
|
1 | Terminal> python setup.py build_ext --inplace
|
1 2 3 4 5 6 | import wave2D_u0_loop_cy
advance = wave2D_u0_loop_cy.advance
...
for n in It[1:-1: # time loop
f_a[:,:] = f(xv, yv, t[n]) # precompute, size as u
u = advance(u, u_1, u_2, f_a, x, y, t, Cx2, Cy2, dt2)
|
Efficiency:
numpy: 5.5numpy: 6advance function in pure Fortranf2py to generate C code for calling Fortran from Python
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 | subroutine advance(u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny)
integer Nx, Ny
real*8 u(0:Nx,0:Ny), u_1(0:Nx,0:Ny), u_2(0:Nx,0:Ny)
real*8 f(0:Nx, 0:Ny), Cx2, Cy2, dt2
integer i, j
Cf2py intent(in, out) u
C Scheme at interior points
do j = 1, Ny-1
do i = 1, Nx-1
u(i,j) = 2*u_1(i,j) - u_2(i,j) +
& Cx2*(u_1(i-1,j) - 2*u_1(i,j) + u_1(i+1,j)) +
& Cy2*(u_1(i,j-1) - 2*u_1(i,j) + u_1(i,j+1)) +
& dt2*f(i,j)
end do
end do
|
Note: Cf2py comment declares u as input argument and return value
back to Python
1 2 3 4 | Terminal> f2py -m wave2D_u0_loop_f77 -h wave2D_u0_loop_f77.pyf \
--overwrite-signature wave2D_u0_loop_f77.f
Terminal> f2py -c wave2D_u0_loop_f77.pyf --build-dir build_f77 \
-DF2PY_REPORT_ON_ARRAY_COPY=1 wave2D_u0_loop_f77.f
|
f2py changes the argument list (!)
1 2 3 4 5 6 | >>> import wave2D_u0_loop_f77
>>> print wave2D_u0_loop_f77.__doc__
This module 'wave2D_u0_loop_f77' is auto-generated with f2py....
Functions:
u = advance(u,u_1,u_2,f,cx2,cy2,dt2,
nx=(shape(u,0)-1),ny=(shape(u,1)-1))
|
f2py!f2py takes a copy of a numpy (C) array and transposes it
when calling Fortrannumpy arrays with Fortran storage
1 2 3 4 | order = 'Fortran' if version == 'f77' else 'C'
u = zeros((Nx+1,Ny+1), order=order)
u_1 = zeros((Nx+1,Ny+1), order=order)
u_2 = zeros((Nx+1,Ny+1), order=order)
|
Option -DF2PY_REPORT_ON_ARRAY_COPY=1 makes f2py write out
array copying:
1 2 | Terminal> f2py -c wave2D_u0_loop_f77.pyf --build-dir build_f77 \
-DF2PY_REPORT_ON_ARRAY_COPY=1 wave2D_u0_loop_f77.f
|
numpy codeadvance function in pure Cnumpy arrays transferred to C are one-dimensional in C[i,j] indices to single indices
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | /* Translate (i,j) index to single array index */
#define idx(i,j) (i)*(Ny+1) + j
void advance(double* u, double* u_1, double* u_2, double* f,
double Cx2, double Cy2, double dt2,
int Nx, int Ny)
{
int i, j;
/* Scheme at interior points */
for (i=1; i<=Nx-1; i++) {
for (j=1; j<=Ny-1; j++) {
u[idx(i,j)] = 2*u_1[idx(i,j)] - u_2[idx(i,j)] +
Cx2*(u_1[idx(i-1,j)] - 2*u_1[idx(i,j)] + u_1[idx(i+1,j)]) +
Cy2*(u_1[idx(i,j-1)] - 2*u_1[idx(i,j)] + u_1[idx(i,j+1)]) +
dt2*f[idx(i,j)];
}
}
}
}
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | import numpy as np
cimport numpy as np
cimport cython
cdef extern from "wave2D_u0_loop_c.h":
void advance(double* u, double* u_1, double* u_2, double* f,
double Cx2, double Cy2, double dt2,
int Nx, int Ny)
@cython.boundscheck(False)
@cython.wraparound(False)
def advance_cwrap(
np.ndarray[double, ndim=2, mode='c'] u,
np.ndarray[double, ndim=2, mode='c'] u_1,
np.ndarray[double, ndim=2, mode='c'] u_2,
np.ndarray[double, ndim=2, mode='c'] f,
double Cx2, double Cy2, double dt2):
advance(&u[0,0], &u_1[0,0], &u_2[0,0], &f[0,0],
Cx2, Cy2, dt2,
u.shape[0]-1, u.shape[1]-1)
return u
|
Compile and link the extension module with a setup.py file:
1 2 3 4 5 6 7 8 9 10 11 12 13 | from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
sources = ['wave2D_u0_loop_c.c', 'wave2D_u0_loop_c_cy.pyx']
module = 'wave2D_u0_loop_c_cy'
setup(
name=module,
ext_modules=[Extension(module, sources,
libraries=[], # C libs to link with
)],
cmdclass={'build_ext': build_ext},
)
|
1 | Terminal> python setup.py build_ext --inplace
|
In Python:
1 2 3 4 5 | import wave2D_u0_loop_c_cy
advance = wave2D_u0_loop_c_cy.advance_cwrap
...
f_a[:,:] = f(xv, yv, t[n])
u = advance(u, u_1, u_2, f_a, Cx2, Cy2, dt2)
|
advance function in pure Cf2py to generate C code for calling C from Pythonadvance as beforeadvance functionf2py generate the Fortran 90 moduleintent(c)):
1 2 3 4 5 6 7 8 9 | subroutine advance(u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny)
Cf2py intent(c) advance
integer Nx, Ny, N
real*8 u(0:Nx,0:Ny), u_1(0:Nx,0:Ny), u_2(0:Nx,0:Ny)
real*8 f(0:Nx, 0:Ny), Cx2, Cy2, dt2
Cf2py intent(in, out) u
Cf2py intent(c) u, u_1, u_2, f, Cx2, Cy2, dt2, Nx, Ny
return
end
|
Generate Fortran 90 module (wave2D_u0_loop_c_f2py.pyf):
1 2 3 | Terminal> f2py -m wave2D_u0_loop_c_f2py \
-h wave2D_u0_loop_c_f2py.pyf --overwrite-signature \
wave2D_u0_loop_c_f2py_signature.f
|
The compile and build step must list the C files:
1 2 3 | Terminal> f2py -c wave2D_u0_loop_c_f2py.pyf \
--build-dir tmp_build_c \
-DF2PY_REPORT_ON_ARRAY_COPY=1 wave2D_u0_loop_c.c
|
numpy C arrays to C++ array classesf2pynumpy arrays to C API and let C translate numpy arrays
into C++ array classes
Solutions: $$ u(x,t) = g_R(x-ct) + g_L(x+ct) $$
If \( u(x,0)=I(x) \) and \( u_t(x,0)=0 \): $$ u(x,t) = \half I(x-ct) + \half I(x+ct) $$
Two waves: one traveling to the right and one to the left
A wave propagates perfectly (\( C=1 \)) and hits a medium with 1/4 of the wave velocity. A part of the wave is reflected and the rest is transmitted.
We have just changed the initial condition...
Build \( I(x) \) of wave components \( e^{ikx} = \cos kx + i\sin kx \): $$ I(x) \approx \sum_{k\in K} b_k e^{ikx} $$
Our interest: one component \( e^{i(kx -\omega t)} \), \( \omega = kc \)
A similar discrete \( u^n_q = e^{i(kx_q - \tilde\omega t_n)} \) solves $$ [D_tD_t u = c^2 D_xD_x u]^n_q $$
Note: different frequency \( \tilde\omega\neq\omega \)
By \( \omega\rightarrow k \), \( t\rightarrow x \), \( n\rightarrow q \)) it follows that $$ [D_xD_x e^{ikx}]_q = -\frac{4}{\Delta x^2}\sin^2\left( \frac{k\Delta x}{2}\right)e^{ikq\Delta x} $$
Inserting a basic wave component \( u=e^{i(kx_q-\tilde\omega t_n)} \) in the scheme requires computation of $$ \begin{align*} \lbrack D_tD_t e^{ikx}e^{-i\tilde\omega t}\rbrack^n_q &= \lbrack D_tD_t e^{-i\tilde\omega t}\rbrack^ne^{ikq\Delta x}\nonumber\\ &= -\frac{4}{\Delta t^2}\sin^2\left( \frac{\tilde\omega\Delta t}{2}\right)e^{-i\tilde\omega n\Delta t}e^{ikq\Delta x}\\ \lbrack D_xD_x e^{ikx}e^{-i\tilde\omega t}\rbrack^n_q &= \lbrack D_xD_x e^{ikx}\rbrack_q e^{-i\tilde\omega n\Delta t}\nonumber\\ &= -\frac{4}{\Delta x^2}\sin^2\left( \frac{k\Delta x}{2}\right)e^{ikq\Delta x}e^{-i\tilde\omega n\Delta t} \end{align*} $$
The complete scheme, $$ \lbrack D_tD_t e^{ikx}e^{-i\tilde\omega t} = c^2D_xD_x e^{ikx}e^{-i\tilde\omega t}\rbrack^n_q $$
leads to an equation for \( \tilde\omega \): $$ \sin^2\left(\frac{\tilde\omega\Delta t}{2}\right) = C^2\sin^2\left(\frac{k\Delta x}{2}\right), $$ where \( C = \frac{c\Delta t}{\Delta x} \) is the Courant number
Taking the square root: $$ \sin\left(\frac{\tilde\omega\Delta t}{2}\right) = C\sin\left(\frac{k\Delta x}{2}\right) $$
Assume \( C>1 \). Then $$ \underbrace{\sin\left(\frac{\tilde\omega\Delta t}{2}\right)}{>1} = C\sin\left(\frac{k\Delta x}{2}\right) $$
1 2 | def r(C, p):
return 2/(C*p)*asin(C*sin(p))
|
Note: the shortest waves have the largest error, and short waves move too slowly.
For small \( p \), Taylor expand \( \tilde\omega \) as polynomial in \( p \):
1 2 3 4 5 6 7 8 9 10 11 12 13 14 | >>> C, p = symbols('C p')
>>> rs = r(C, p).series(p, 0, 7)
>>> print rs
1 - p**2/6 + p**4/120 - p**6/5040 + C**2*p**2/6 -
C**2*p**4/12 + 13*C**2*p**6/720 + 3*C**4*p**4/40 -
C**4*p**6/16 + 5*C**6*p**6/112 + O(p**7)
>>> # Factorize each term and drop the remainder O(...) term
>>> rs_factored = [factor(term) for term in rs.lseries(p)]
>>> rs_factored = sum(rs_factored)
>>> print rs_factored
p**6*(C - 1)*(C + 1)*(225*C**4 - 90*C**2 + 1)/5040 +
p**4*(C - 1)*(C + 1)*(3*C - 1)*(3*C + 1)/120 +
p**2*(C - 1)*(C + 1)/6 + 1
|
Leading error term is \( \frac{1}{6}(C^2-1)p^2 \) or $$ \frac{1}{6}\left(\frac{k\Delta x}{2}\right)^2(C^2-1) = \frac{k^2}{24}\left( c^2\Delta t^2 - \Delta x^2\right) = \Oof{\Delta t^2, \Delta x^2} $$
Smooth wave, few short waves (large \( k \)) in \( I(x) \):
Not so smooth wave, significant short waves (large \( k \)) in \( I(x) \):
is a typically solution of $$ u_{tt} = c^2(u_{xx} + u_{yy}) $$
Can build solutions by adding complex Fourier components of the form $$ e^{i(k_xx + k_yy - \omega t)} $$
This equation admits a Fourier component $$ u^n_{q,r} = e^{i(k_x q\Delta x + k_y r\Delta y - \tilde\omega n\Delta t)} $$
Inserting the expression and using formulas from the 1D analysis: $$ \sin^2\left(\frac{\tilde\omega\Delta t}{2}\right) = C_x^2\sin^2 p_x + C_y^2\sin^2 p_y $$
where $$ C_x = \frac{c^2\Delta t^2}{\Delta x^2},\quad C_y = \frac{c^2\Delta t^2}{\Delta y^2}, \quad p_x = \frac{k_x\Delta x}{2},\quad p_y = \frac{k_y\Delta y}{2} $$
Rreal-valued \( \tilde\omega \) requires $$ C_x^2 + C_y^2 \leq 1 $$
or $$ \Delta t \leq \frac{1}{c} \left( \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}\right)^{-\halfi} $$
For \( c^2=c^2(\xpoint) \) we must use the worst-case value \( \bar c = \sqrt{\max_{\xpoint\in\Omega} c^2(\xpoint)} \) and a safety factor \( \beta\leq 1 \): $$ \Delta t \leq \beta \frac{1}{\bar c} \left( \frac{1}{\Delta x^2} + \frac{1}{\Delta y^2} + \frac{1}{\Delta z^2}\right)^{-\halfi} $$
For visualization, introduce \( \theta \): $$ k_x = k\sin\theta,\quad k_y=k\cos\theta, \quad p_x=\half kh\cos\theta,\quad p_y=\half kh\sin\theta$$
Also: \( \Delta x=\Delta y=h \). Then \( C_x=C_y=c\Delta t/h\equiv C \).
Now \( \tilde\omega \) depends on
Can make color contour plots of \( 1-\tilde c/c \) in polar coordinates with \( \theta \) as the angular coordinate and \( kh \) as the radial coordinate.
