Geol 455/655 - Lab 2, Courant Condition

For the second exercise, the objective is to begin exploring the parameters of the numerical simulation. The first one is the Courant Condition, which affects the dt parameter.

Aliasing and the Courant Condition

aliasing graph
The blue curve above is a sine wave with wavelength lambda. On the left, the green lines represent a discrete approximation to the continuous curve, with two sample points (white circles) per wavelength. The green lines are clearly an adequate approximation of the sine wave, with insignificant differences between the two.

On the right the red lines represent a sampling of the sine wave at only 1.5 points per wavelength. This is clearly not enough samples per wavelength, as the red lines do not approximate the blue curve well at all, and there are large errors between them.

In E3D's computations, the time-step dt must be short enough to properly represent the waves being modeled. The equation v = f (lambda) above shows that the velocity at any point in the model is related to the wavelength.

The Nyquist Criterion states that a discrete representation of a wave in time must have a sampling frequency fs that provides two points per period, or fs = 2 fmax , where fmax is the highest frequency found in the waves, fs = 1/dt , and dt is the time sampling interval.

The Courant Condition states that the time sampling dt must be small enough that the longest wavelengths, propagating at the highest velocity Vmax , do not outrun the spatial grid sampling dh. There is an additional adjustment to the Courant Condition that we will not worry about here, since it comes from the nature of the finite-difference numerical approximation used by E3D.

Combining the Nyquist Criterion with v = f (lambda) will show you how small you have to make dt to keep two samples per spatial wavelength.

Tasks:

  1. Re-run the computation from Lab 1, to make sure you still have all the input and output files in your folder. As you watch it run, note what E3D says about the Courant Condition, ct/dt, and the MFLOPS achieved.
  2. Re-name the sac files from the Lab 1 computation, to indicate the value of the dt parameter used on the ``time'' line of the ``test.in'' and ``test1.in'' files: 
    	cp LSM.000.x LSM.12dt.x
	cp UNLV.001.x UNLV.12dt.x
  3. Edit the ``test1.in'' file to change the ``dt'' and ``t'' parameters, making dt half of what it was, and doubling t so the same amount of model time is run. Don't bother with the ``test.in'' file, because for this lab we will not need to run ModelAssembler. The time line will look like: 
    	# time       dt=0.12  t=1000
	 time       dt=0.06  t=2000
    So the previous values are just commented out. This should provide finer time sampling, and lead to more accuracy in the computation.
  4. Run the model with the command: 
    	/quake/rg/e3d/bin/e3d test1.in &
    We don't need to run "csh run.sh" because we already have all the medium volume files, and don't need to run ModelAssembler now.
  5. Observe the wave propagation in the pop-up image window. What is different about this run compared to the Lab 1 run? The Courant Condition? How much longer (wall-clock time) does it take? Why doesn't the MFLOPS change?
  6. Re-name the sac files from this new computation, to indicate the new value of the dt parameter: 
    	cp LSM.000.x LSM.06dt.x
	cp UNLV.001.x UNLV.06dt.x
  7. Use ``sac'' to plot and compare the four seismograms: the two LSM and UNLV x-components from the Lab 1 run with dt=0.12 s; and the two with dt=0.06 s. In sac you can get them with the command: 
    	read LSM.12dt.x LSM.06dt.x UNLV.12dt.x UNLV.06dt.x
    Note that the seismogram files made with dt=0.06 s take up twice as much space as those made with dt=0.12 s, and they took longer to compute. Since they all cover the same 120-sec time period, sac plots all four seismograms on the same horizontal axis. Are the more expensive seismograms using dt=0.06 s any different or better?
  8. Now, just as in preparing for the dt=0.06 s run as above, edit test1.in and make another E3D run using dt=0.24 and t=500. What happens? What is the Courant Condition for this run?

What to turn in:

  1. Answer: Given constant frequency, how do wavelengths change as velocity changes from place to place in the grid? Using this concept, explain how wavelengths on the surface of the computation change across the visual display while E3D is running.
  2. Answer: Given the usual relation between P and S velocities, will P waves or S waves have longer wavelengths?
  3. Answer: Solve v = f (lambda) for fmax given the maximum velocity reported in the grids by E3D, and using dh for lambda. According to the Nyquist criterion, what dt is required? Is this similar to E3D's reported dt as required by the Courant Condition?
  4. Answer the questions in italics above.
  5. Turn in the plot of the 4 seismograms from the 2 stations and the 2 different dt values, properly and completely labeled on all axes (with a pencil if necessary).
  6. Answer: Given velocity grids and a fixed dh such as we are using here, could we save any computational effort by removing structural features such as the basins? Can you propose a feature of the models that we could remove to enable us to make dt larger, and thus make a synthetic seismogram covering the same time with less computational effort? (Are there any structural features in the grids other than the basins?)

Resources: