Leapfrog Scheme for Advection Equation

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.

• Leapfrog Scheme
• Advection Equation
• Stability Analysis
• Courant Number
• CFL Condition

kri_whi Follow

Uploaded on Sep 19, 2024 | 0 Views

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. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.

Presentation Transcript

1. Finite difference IV ex) advection equation 2 Keywords leapfrog scheme physical and computational modes

2. 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. 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. 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. 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. 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. 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. 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. 9 Summary of stability dx dt u Upwind Downwind Center stable (numerical diffusion) Forward unstable unstable neutral Center unstable unstable

10. 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. 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. 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. 13 Example: physical mode Initial conditions: h(t = 0) = sin(2 x/16dx) h(t = dt) = sin(2 (x udt)/16dx) h x

14. 14 Example: computational mode Initial conditions: h(t = 0) = sin(2 x/16dx) h(t = dt) = sin(2 (x udt)/16dx) h x

15. 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. 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. 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. 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. 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. 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. 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. 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. 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. 24 Advection+damping equation h x

25. 25 Example %matplotlib inline import matplotlib as mpl import matplotlib.pyplot as plt plt.contourf(h.T) plt.colorbar() time x