Programming 1D Linear Wave Equation with Upstream Scheme

Model Task #0A: Programming the 1D upstream scheme ATM 562 Fall 2021 Fovell 1

Overview Program the 1D linear wave equation using the upstream approximation in a periodic domain Task illustrates the basics of Model initialization Time stepping Application of boundary conditions Amplitude, phase, and dispersion errors Dependence on wavelength and time step Avoidance of roundoff errors Use programming language (e.g., Python, Fortran, C/C++) and visualization package (e.g., NCL, GrADS, GNUplot, Excel*) of your choice If you plan on using (modifying) NWP models, learning Fortran is very highly recommended 2 *Excel won t be useful for model tasks 1

Problem 1D linear wave equation When c > 0, wave translates to the right Upstream approximation is forward in time and upwind in space with 1st order accuracy 3

Explicit form This approximation can be written explicitly as The initial time = 0 seconds The time index n starts at 0 The space index runs from i = 1 to NX (Fortran) or i = 0 to NX-1 (Python, C/C++, other zero-based languages) The exact solution preserves its original shape when c is constant With periodic boundary conditions, the exact solution will return to its original position after a certain number of time steps (NT) 4

Program structure Declare needed arrays I will call present value u , forecast value up Initialize parameters and constants (NX, NT, n, dt, dx, trigonometric , etc.) Provide initial condition for u at n = 0 for all i Step ahead by time step dt, using upstream scheme to create up for all real points (see slide 8) n = n + 1 time = time + dt not the best choice (see next slide) Apply boundary conditions (BCs) on up Set for new time step (i.e., u = up) If time < timend, or n < NT, loop 5

Tips Define trigonometric to machine precision: Fortran: 4*arctan(1.0) Python: np.pi from numpy package To avoid confusion with nondimensional pressure (also called pi ), I call this trigpi To avoid roundoff error, advance time as time = (n-1)*dt instead of time = time + dt (see Appendix) 6

Periodic boundary conditions For NX points, NX-2 are real points, and the other two (fake) will facilitate application of the boundary conditions (BCs) For periodic BCs (Fortran) up(1) = up(nx-1) and up(nx) = up(2) So, we loop from i = 2, nx-1 (the real points) and then update the fake points For periodic BCs (zero-based) up(0) = up(nx-2) and up(nx-1) = up(1) - So, we loop from i = 1, nx-2 (the real points) and then update the fake points 7

Logic of fake points (Fortran) Part of time stepping loop ! solve for up over real points [rdx = 1/dx] do i=2,nx-1 up(i) = u(i)-c*dt*rdx*(u(i)-u(i-1)) ! needs u(1)!! enddo ! apply boundary conditions on up up(1) = up(nx-1) up(nx) = up(2) ! set for the next time step, over real AND FAKE points do i=1,nx u(i) = up(i) enddo ! creates u(1)!! Fake points help avoid IF statements that make code harder to read, more likely to contain errors When time loop repeats, don t need special handling for u(1), since it s been taken care of Facilitates application of all boundary condition types, not just periodic 8

Test problem Let c = c t/ x (by definition) Test: c = 1.0, dt = 1.0, dx = 1.0, so c = 1.0 NX = 52 (50 real points) NT = 50 (timend = 50 since dt = 1.0) Wavelength L = 50 x (one sine wave in domain), with amplitude 1.0 Execute for one revolution Plot on next slide shows initial condition (& exact solution) and upstream approximation after one revolution for test For c = 1.0, there should be no significant error (so this is your model test) We will manipulate c by changing t 9

After 1 revolution Upstream approximation to u_t = -cu_x 1.5 1.0 0.5 0.0 0 5 10 15 20 25 30 35 40 45 50 -0.5 -1.0 -1.5 exact upstream (cfl=1, 1 rev) (No significant error) 10 Abscissa: grid point number (Fortran)

Initial condition for u (Fortran) wavelength = 50 ! i.e., 50 x amp = 1.0 ! amplitude do i=2,nx-1 ! loop over real points xi=float(i-2) ! re-map the x-coordinate u(i) = amp*sin(2*xi*trigpi/wavelength) print *,' i ',i,' u ',u(i) enddo ! enforce periodic boundary conditions u(nx)= u(2) u(1) = u(nx-1) 11

Exact solution for u at time step n do i=2,nx-1 ! exact solution at any time n for the real points xei = float(i-2)-c*n*dt uexact(i) = amp*sin(2*xei*trigpi/wavelength) enddo ! then enforce the boundary conditions Wave movement to the right is realized by sliding the coordinate axis to the left as a function of time step 12

Errors Amplitude error occurs if the original wave amplitude is not preserved Phase error occurs if the simulated wave translates too quickly or slowly Plots on next slide show amplitude error (left) and phase error (right) as a function of wavelength (horizontal axis, longer waves at leftend) and c Shortest resolvable wavelength is 2 x (is that obvious?) Plots on next slide show: Amplitude decays if c < 1, grows if c > 1, and shorter waves handled worse than longer waves Waves propagate too slowly for all wavelengths; errors reasonable for long waves but large for short waves 13

Amplitude and phase errors (Fig. 4.3 from notes) c = c t/ x 0.9 > 1.0 0.96 1.0 > 1.0 0.9 0.6 Values < 1: damping Values < 1: too slow 14

Amplitude error For c = 1, upstream scheme has essentially zero amplitude error (apart from roundoff) For c 1, performance worst for shortest waves c > 1 unstable solution grows exponentially with time c < 1 stable but solution damps exponentially with time In both cases, growth or damping largest for shortest waves Sampling the wave propagation better(i.e., c < 1) makes the solution worse (counterintuitive!) {Frame rate} For c < 1, performance worst at c = 0.5 for all waves (even more counterintuitive and somewhat deceptive!) At wavelength of 4 x and c = 0.5, amplitude loss is 30% (0.7). That is per time step. That s why it is exponential. 15

Phase error For the exact solution, phase speed is independent of wavelength. All waves move at the same speed, c. The exact solution is nondispersive. Upstream scheme phase propagation is not absolutely perfect evenfor c = 1 For wavelengths > 4 x, propagation too slow by 4% (0.96) For c < 1, some wavelengths move too quickly, others too slowly Again, the shorter the wavelength, the poorer the performance 2 x waves are stationary. (Really?) Since phase error depends on wavelength, the numerical scheme is dispersive when there s more than one wave. 16

Experiments Experiments 1-5 check your code and understanding. Turn in result of Experiment 6 with a copy of your code 17

Experiment #1 Try values of t < 1.0 (so c < 1) and compare to exact solution for the 50 x wave Observe how amplitude and phase error vary with c after 1 or more revolutions Compare your results to Fig. 4.3 (left) If you plot maximum wave amplitude vs. time for c < 1, you should observe exponential decay. Plot on semi-log graph to confirm 18

Experiment #2 Try values of c > 1 (by increasing t). What happens? How does amplitude error after 1 revolution change with c ? You should observe exponential growth This is linear instability. 19

Experiment #3 Change the wavelength of the initial condition to, say, 5 x (wavelength = 5) and c = 0.5 ( t = 0.5). From Fig. 4.3 from the notes, we expect to lose 20% amplitude per time step. Plot maximum wave value vs. time. Do you? Warning: because of the periodic BCs, your initial condition has to have an integer number of waves in the domain, or the BC will mishandle the wave. To explore some wavelengths, it may be necessary to change NX. 20

Experiment #4 Compare solutions for the 50 x wave with c = 0.5 and 0.25. Run each simulation for 50 time units (e.g., seconds). According to Fig. 4.3 (left), the amplitude error is supposed to be worse for c = 0.5. Is it? If not, why not? Note 50 seconds requires 100 time steps with dt = 0.5 second, and 200 time steps for dt = 0.25 second. 21

Experiment #5 According to Fig. 4.3 (right), the 2 x wave isn t supposed to move with time for any value of c . Does it move? If it does, why does it? Note that because a 2 x wave is easily lost on the grid, you have to modify the code a little xi=float(i-2) if(wavelength.eq.2) xi=xi+0.5 u(i) = amp*sin(2*xi*trigpi/wavelength) 22

Experiment #6 An initial condition can consist of several waves, which combine to make a certain shape. The exact solution, being nondispersive, preserves that shape as the wave moves. The upstream scheme, however, is dispersive, as it has a different amplitude and phase error for different waves, so the combined shape will change. Demonstrate this with an initial condition with two waves. Specifically, additively combine a 50 and 10 x wave, each with initial amplitude of 1.0, and set c =0.5. Anticipate what the result would look like after 1 revolution before running the model and checking your guess. (If you wish to explore other wave combinations, keep in mind you need to have an integral number of waves in the initial condition.) For your report, it suffices to turn in a plot illustrating the result of this experiment, along with your model code. 23

Exp. 6 initial condition and exact solution Initial condition and exact solution for MT#0 Exp. 6 2.5 2.0 1.5 1.0 0.5 0.0 0 5 10 15 20 25 30 35 40 45 50 -0.5 -1.0 -1.5 -2.0 -2.5 exact 24

Appendix: roundoff error [You do not need to turn in suggested exercises] 25

Test #1 For Fortran, a real-valued (single precision) constant is 4 bytes long. A double precision value is 8 bytes. Write a program that assigns 10 entries of an single precision array to values of 0.1 to 0.9. Print these values to at least 10 decimal places. What do you observe? Now make the array elements double precision. What happens? 26

Test #2 Let dt, time1 and time2 be single precision, real variables Set dtto a value like 0.6, and time1and time2 to 0. Compare these after 20000 iterations (iter) time1 = time1 + dt time2 = float(iter)*dt (For a model with a 4-5 sec time step, 20000 iterations represents about 1 day of integration.) 27

Questions for thought Most NWP models are single precision. Why? You write separate codes with Fortran and Python, and don t obtain the same results. Why? 28