Leapfrog Scheme for Advection Equation
Finite difference
IV
ex) advection equation
2
Keywords
•
leapfrog scheme
•
physical and computational modes
Advection equation
•
It is assumed that
u
is positive constant for simplicity
2
move eastward
2. Leapfrog method
•
Here, we consider the center method i
n time and space
•
To check the stability, it is again assumed that
•
If the
upwind
and
downwind
methods in space are used, the
center
method in time is
unstable
(not
shown)
3
2. Leapfrog method
•
If the Courant number
C
is defined as
udt
/
dx
, this equation
becomes
4
2. Leapfrog method
•
It is necessary to separate the sign of the inside of the square
root
5
(i)
(ii)
2. Leapfrog method
(i)
6
neutral
neutral
real part
imaginary part
2. Leapfrog method
•
What is
•
Because
7
CFL condition
CFL condition
2. Leapfrog method
(ii)
8
•
Because
C
sin
kdx
> 1, |
| > 1
•
Thus, if we satisfy the CFL condition, the center method in time
and space is stable for the advection equation
unstable
unstable
Summary of stability
9
Note
•
At the first timestep (
t
= 2
dt
), it is necessary to obtain the data at
t
=
dt,
as well as the initial condition (
t
= 0)
•
In general, the forward method is used to calculate the data at
t
=
dt
10
Example
11
1.
Why are the values negative?
2.
Why is the shape changed?
h
x
no numerical
diffusion
Why are the values negative?
•
Recalling that
•
=1 is a
physical mode
physical mode
, and
=
‒1 is a
computational mode
computational mode
,
which is
unrealistic
•
The forward and backward methods have only one mode, but
the center method has two modes
•
One cause of a computational mode is inconsistency between
values at the current and previous timesteps
12
Example: physical mode
•
Initial conditions:
h
(
t
= 0) = sin(2
x
/16
dx
)
h
(
t
=
dt
) = sin(2
(
x
‒
udt
)/16
dx
)
13
h
x
Example: computational mode
•
Initial conditions:
h
(
t
= 0) = sin(2
x
/16
dx
)
h
(
t
=
dt
) =
‒
sin(2
(
x
‒
udt
)/16
dx
)
14
h
x
Example: computational mode
•
Initial conditions:
h
(
t
= 0) = sin(2
x
/16
dx
)
h
(
t
=
dt
) =
‒
sin(2
(
x
‒
udt
)/16
dx
)
15
h
time
0
Asselin filter
•
Several ways are proposed to reduce the effect of the
computational mode
•
h
* is corrected values by the Asselin filter
•
is from 0.1 to 0.3, and is around 0.1
16
Asselin filter
•
This example is
= 0.4
•
If
= 0.5, the summation becomes zero
17
time
time
h
Example: Asselin filter
18
a little bit better
Why is the shape changed?
•
Wavenumber dependency of propagation speed is checked
19
Why is the shape changed?
20
where
h
0
=
h
(0). Therefore,
Why is the shape changed?
•
The propagation speed is not equal to
u
, but depends on
wavenumber
k
(i.e. wavelength
)
•
When
kdx
<< 1 (large wavelength), because sin
kdx
kdx
, the
propagation speed is equal to
u
•
However, when
kdx
=
(i.e.
= 2
dx
), the propagation speed is
zero
21
Why is the shape changed?
•
In this center difference method,
= 2
dx
is
the shortest
the shortest
wavelength
wavelength
•
This is because
h
(
x
+
dx
) –
h
(
x
–
dx
) = 0
•
Since
the propagation speed depends on wavenumber
,
the
shape of the waves is changed (
dispersive
dispersive
)
22
x
h
Advection+damping equation
23
Advection+damping equation
24
h
x
Example
25
%matplotlib inline
import matplotlib as mpl
import matplotlib.pyplot as plt
plt.contourf(h.T)
plt.colorbar()
time
x
The provided images illustrate the Leapfrog scheme applied to an advection equation, focusing on the center method in time and space. The stability of the method is analyzed with assumptions regarding the behavior of the solution. Through the exploration of Courant numbers and CFL conditions, the stability of the scheme for the advection equation is discussed, distinguishing between various spatial and temporal methods.
Download Presentation
Please find below an Image/Link to download the presentation.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author.If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.
You are allowed to download the files provided on this website for personal or commercial use, subject to the condition that they are used lawfully. All files are the property of their respective owners.
The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author.
E N D
Presentation Transcript
Finite difference IV ex) advection equation 2 Keywords leapfrog scheme physical and computational modes
2 Advection equation dh dt h t h x = + = 0 u It is assumed that u is positive constant for simplicity 1.2 t=0 t=30 t=70 1 0.8 move eastward h 0.6 0.4 0.2 0 10 20 30 40 50 x 60 70 80 90 100
3 2. Leapfrog method Here, we consider the center method in time and space + + ( , h x t ) 2 ( , h x t dt ) ( , ) ( , ) dt dt h x dx t h x dx t dx = u 2 To check the stability, it is again assumed that ( , ( , ) h x t + = ) ( )exp( h t ( , ) h x t ikx h x t dt = ) If the upwind and downwind methods in space are used, the center method in time is unstable (not shown)
4 2. Leapfrog method 2 ikdx ikdx 1 e e = u 2 2 dt dx If the Courant number C is defined as udt/dx, this equation becomes = i C 2 1 2 sin kdx = + 2 2 sin sin 1 iC kdx C kdx
5 2. Leapfrog method = + 2 2 sin sin 1 iC kdx C kdx It is necessary to separate the sign of the inside of the square root (i) sin C (ii) sin C 1 0 + 1 0 + 2 2 kdx kdx 2 2
6 2. Leapfrog method 1 0 + 2 2 (i) sin C kdx = + 2 2 sin sin real part 1 iC imaginary part kdx C kdx 2 = + + = 2 2 2 2 sin 1 sin 1 C kdx C kdx neutral
7 2. Leapfrog method What is 1 0? + 2 2 sin C kdx 2 0 sin 1 1 1 kdx C C CFL condition Because 2
8 2. Leapfrog method (ii) 1 0 + 2 2 sin C kdx = + 2 2 sin sin 1 iC kdx C kdx = 2 2 sin 1 sin iC kdx i C kdx = + 2 2 2 2 2 2 2 sin 1 sin 2 sin C 1 sin C kdx C kdx kdx C kdx 1 2 sin C = 2 2 1 sin kdx C kdx unstable Because C sin kdx > 1, | | > 1 Thus, if we satisfy the CFL condition, the center method in time and space is stable for the advection equation
9 Summary of stability dx dt u Upwind Downwind Center stable (numerical diffusion) Forward unstable unstable neutral Center unstable unstable
10 Note + + ( , h x t ) 2 ( , h x t dt ) ( , ) ( , ) dt dt h x dx t h x dx t dx = u 2 At the first timestep (t = 2dt), it is necessary to obtain the data at t = dt, as well as the initial condition (t = 0) In general, the forward method is used to calculate the data at t = dt
11 Example no numerical diffusion 1.4 1.2 1 0.8 h 0.6 0.4 0.2 0 -0.2 1. Why are the values negative? 2. Why is the shape changed? x -0.4 0 20 40 60 80 100
12 Why are the values negative? 2 = 1 1 Recalling that = =1 is a physical mode, and = 1 is a computational mode, which is unrealistic The forward and backward methods have only one mode, but the center method has two modes One cause of a computational mode is inconsistency between values at the current and previous timesteps
13 Example: physical mode Initial conditions: h(t = 0) = sin(2 x/16dx) h(t = dt) = sin(2 (x udt)/16dx) h x
14 Example: computational mode Initial conditions: h(t = 0) = sin(2 x/16dx) h(t = dt) = sin(2 (x udt)/16dx) h x
15 Example: computational mode Initial conditions: h(t = 0) = sin(2 x/16dx) h(t = dt) = sin(2 (x udt)/16dx) h 0 time
16 Asselin filter Several ways are proposed to reduce the effect of the computational mode *( , ) x t h = + + ) 2 ( , ) h x t + ( , ) h x t ( , h x t *( , ) dt h x t dt 2 h* is corrected values by the Asselin filter is from 0.1 to 0.3, and is around 0.1
17 Asselin filter + + ) 2 ( , ) h x t + ( , ) h x t ( , h x t ( , h x t ) dt dt 2 1 1 original correction sum 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 h 0 0 -0.2 -0.2 -0.4 -0.4 -0.6 -0.6 -0.8 -0.8 -1 -1 2 4 6 8 10 12 14 16 18 2 4 6 8 10 12 14 16 18 time time This example is = 0.4 If = 0.5, the summation becomes zero
18 Example: Asselin filter 1.2 a little bit better 1 0.8 0.6 h 0.4 0.2 0 -0.2 -0.4 0 20 40 60 80 100 x
19 Why is the shape changed? Wavenumber dependency of propagation speed is checked + ( , ) ( , ) h t h x dx t h x dx t dx ( )exp( h t = u 2 = ( , ) h x t ) ikx ( ) ikdx exp( ) exp( ) exp( dx ) ikx ikdx dh dt = exp( ) ikx uh 2 sin dh dt kdx dx = iu h
20 Why is the shape changed? dh dt sin kdx dx sin iu = iu h kdx dx ( ) = exp h t h t 0 where h0= h (0). Therefore, ( )exp( h t = ( , ) h x t ) ikx sin kdx kdx = exp h ik x u t 0
21 Why is the shape changed? sin kdx kdx = propagation speed u u The propagation speed is not equal to u, but depends on wavenumber k (i.e. wavelength ) When kdx << 1 (large wavelength), because sin kdx kdx, the propagation speed is equal to u However, when kdx = (i.e. = 2dx), the propagation speed is zero
22 Why is the shape changed? In this center difference method, = 2dx is the shortest wavelength 1000.01 1000.005 h 1000 999.995 999.99 0 2 4 6 8 x 10 12 14 16 This is because h(x+dx) h(x dx) = 0 Since the propagation speed depends on wavenumber, the shape of the waves is changed (dispersive)
23 Advection+damping equation h u t h x + = h advection Upwind Downwind Center neutral Leapfrog unstable unstable damping Forward Backward Center unstable stable stable
24 Advection+damping equation h x
25 Example %matplotlib inline import matplotlib as mpl import matplotlib.pyplot as plt plt.contourf(h.T) plt.colorbar() time x
