Geol 455/655 - Lab 3, Grid Dispersion

The second equation defines the Laplacian, which is just a second spatial derivative of the scalar field P(x,y,z) (partial derivatives over the spatial axes only, ignoring the time t axis). The Laplacian appears in our differential equation for heat flow, and it also appears in differential equations for wave propagation, called wave equations.

Now, one way to solve such differential equations is to find a closed-form analytic expression that solves a simplified case of the differential equation. We looked at an analytic solution for the heat flow equation in one dimension. The way to solve these equations without simplifying them is to finite difference them. This works well for fields that are discretized and exist as discrete samples in computer memory. We start with the definition of the analytic derivative, the third equation at right.

The fourth equation simply removes the limit and keeps dx finite, or greater than zero. The bold D_x is an operator, one that approximates the analytic partial-derivative operator. It says ``apply the single finite difference in x to the field P(x)'' that may vary in the x direction.

The fifth equation shows how to derive the second finite-difference operator D_x², which is just applying D_x twice to take the difference of the difference and approximate the second spatial derivative.

The sixth line is the result, defining the second-difference operator D_x² that we can just plug into differential equations containing the Laplacian, along with the first difference, to solve them explicitly. Finite-difference computer programs just loop along time t, in steps dt, to build up a computed wavefield (or temperature field) for the time range we desire. As you saw in lab 2, if we make dt too large in an attempt to save computation time, the faster velocities make the discretized wavefield aliased in time, and the program blows up.

At right is a blue sinusoid of amplitude one and wavelength lambda equal to 4 in x. On the left it is sampled with the green line at dx=1 for 4 samples per cycle (or wavelength). On the right it is sampled with the red line at dx=2 for 2 samples per cycle. The coarser sampling and the red line are not aliased, since the Nyquist Criterion says 2 samples per cycle perfectly represent the wave.

The equations on the lower right of the figure here show that the second derivative of this sinusoid is a scaled sinusoid itself. At a value of x equal to three-quarters of the wavelength lambda, which is x=3, the analytic 2^nd derivative has a value of 2.47.

Now, let's apply the formula for the 2^nd-difference approximation at the same point, x=3. The estimate of the 2^nd derivative withD_x² = [P(x+dx)-2P(x)+P(x-dx)]/dx² is centered at a value of x, which will be 3.

On the left, the 3^rd sample on the green line is at a minimum in the sine function. The 2^nd derivative is supposed to be about 2.47, and D_x² produces an estimate of 2.0. For geophysical work, this is sometimes accurate enough.

The second point on the red line, sampling the sine at 2 points per cycle, is at an identical minimum in the sine wave, and ought to still give a 2^nd derivative close to 2.47. However, this coarser sampling does not produce nearly as accurate a D_x² estimate, coming back with a value of 1. The 2^nd difference estimate is highly inaccurate, despite the sampling being sufficient for the Nyquist Criterion.

Note that computing D_x² at x=2, at the sine wave's zero crossing, will always yield the correct value, zero, for the 2^nd derivative there. This is true for both the green and red samplings. In fact, D_x² is near-perfect where the 2^nd derivative is small, and underestimates the 2^nd derivative more and more at it gets larger.

This factor of two error at the maxima will have drastic effects on the wave propagation in E3D. Waves with high frequencies, having the shortest wavelengths, will seem to slow down and even stop propagating once they get below six spatial samples per wavelength. This slowing of high-frequency parts of the wavefield is reminiscent of the dispersion (spread) of the velocities of surface waves according to their frequency. The computational artifact is thus known as grid dispersion, since it depends on the spatial sampling dh of the model grid.

For reliable computations, Shawn Larsen and other authors of finite-difference wave modeling codes suggest that the shortest wavelengths in the calculation cover at least 10 grid points-- in other words, that 10dh < lambda_min . At left is an example E3D computation visual (a map with m=0) that has at least 10 spatial samples per wavelength. At right is a visual from the same model, where the source was set to pump a higher frequency into the grid that results in fewer than 6 points per cycle. The grid dispersion generates false long durations of shaking. You can also notice a ``moire'' pattern of waves that seem to be propagating parallel to the axes rather than radiating circularly from the source. As the computation continues, you will see waves propagating so slowly they appear to be stationary.

Assignment:

1. Given the criterion for avoiding grid dispersion, 10dh < lambda_min , derive an expression that gives you the dh needed to meet this criterion for a highest desired frequency f_max. This expression will involve a velocity-- is this the minimum or the maximum velocity found in the model?

2. If one had a code that was computing only P-wave propagation in a crustal model, with a dh that avoided grid dispersion, and switched to computing only shear-wave propagation across the same realistic model, would dh have to increase or decrease? Could you guess by how much?

3. Briefly describe a few examples of lithospheric geologic features that affect the Courant Condition on dt, and then describe the types of features that affect the need to set dh small enough to avoid grid dispersion.

4. To compute a seismogram representing a given duration, the cost of computation goes up as dt gets smaller, since more time steps are needed. Each time step in general requires computation effort proportional to the number of grid nodes in the model, which is nx times ny times nz. Thus the computational effort required for a model of a certain size goes up as dh decreases, as well. Can you think of a set or sets of lithospheric features, which, if included in a model, cause the computational cost to skyrocket by demanding both small dt and small dh?

5. Download the test3.in E3D input file. Examine it in a text editor. Note that we are not using ModelAssembler for this lab, so all velocities are set in test3.in . Given the frequency specified on the source line of the input file, and the other settings in the file, what is the minimum number of grid samples per wavelength for the computation set up here?

   Run this model with the command:

   	/quake/rg/e3d/bin/e3d test3.in &
   

   Watch the visual to verify that no obvious grid dispersion occurs. As you did in Lab 2, rename the LSM and UNLV x-component output files to save them.

6. Edit the source line of the test3.in file to increase the source frequency by a factor of five. How many points per wavelength will this computation have as a minimum? Run this model and observe the grid dispersion artifacts, which become more obvious as the computation proceeds. Again rename the LSM and UNLV output files to save them.

7. Plot the four seismogram files for the two runs on the LSM and UNLV stations in SAC, label them completely, including which are from which station and which run, and point out some possible features of grid dispersion on the plot.

8. The original source line in test3.in is not a efficient as it might be. How high can you raise the frequency and still avoid grid dispersion? Set this value and make another run. Print out a comparison of the seismograms from these two runs (both without grid dispersion).

Resources:

• E3D modeling code manual, defining inputs and outputs
• test3.in E3D input file, defining how E3D will run.