Geometrical Optics¶

Setup and coding instructions:¶

This session involves little coding from you part. In order to comply with pysit and the conventions of seismic imaging, the z-axis is now pointing downwards.

In order for the Helmholtz computations to run in a short enough time, it is highly recommended to not increase the resolution to more than 200 x 200. (if the 150 x 150 default already take too long, don't think about it !)

If your computer has trouble with the Helmholtz solve, download the precomputed Helmholtz solutions and set the variable MY_COMPUTER_CAN_TAKE_THE_PAIN to False. Unzip the file to your directory. This file contains precomputed solutions for n = [150, 200], Cchoice = [0,1], and nu = [5.0, 10.0, 2.0] Hz.

The parts of the code that you may play with are

• the DASHBOARD part, containing the parameters you can play with
• the YOUR CODE STARTS HERE part, where you may code the remaining part of the exercises.

IMPORTANT: in your dashboard, update the problem_number variable with the number of the problem you are working on. If this variable is set on another problem, it will simply not run.

Part 2: Approximation of the Helmholtz equation¶

The present session aims at testing the validity of the approximation

\[ \widehat{G}(x,x_s,\omega) = a(x,x_s) e^{-i\omega\tau(x,x_s)} + \text{smoother}, \]

where \(\widehat{G}(x,x_s,\omega)\) solves the Helmholtz equation in free space

\[ -\left(\Delta_x + \frac{\omega^2}{c(x)^2}\right)\widehat{G} = \delta(x-x_s), \]

\(\tau\) solves the eikonal equation

\[ c(x) | \nabla_x\tau | = 1, \\\\ \tau(x_s,x_s) = 0. \]

and the amplitude \(a\) solves the transport equation

\[ \nabla a\cdot\nabla\tau = -\frac{1}{2} a\Delta \tau \]

For both equations, the solvers are already coded. The computational domain is \([0,1]\times [0,1]\).

• The Helmholtz solver comes from the pysit package and solves the equivalent of a problem in free space by implementing PMLs at the boundary. The frequency \(\omega\) is fixed and is such that \(\nu = 20 Hz\). By default, the point source position is \((0.1,0.1)\)
• The eikonal solver uses the fast sweeping algorithm introduced in cp7.
• The geodesic solver implemented in cp7 is also provided.

Problem 1 (approximation of the phase): The first thing to compare is the phase of the approximation with that of the exact solution to Helmholtz equation. Pick the sound speed of 'gradient' type and on two subplots of the same figure, visualize

• on the left, \(\arg{\widehat{G(x,x_s,\omega)}}\) (requires running a Helmholtz solve)
• on the right, \(\arg(\exp (-i\omega\tau))\) (requires running an eikonal solve)

Problem 2 (breakdown of the phase approximation): repeat the same experiment with the speed of 'lens' type. Are both phases in good agreement everywhere ?

Problem 3 (ripple effect on Helholtz equation): For the speed of 'lens' type,

• visualize the amplitude and the phase of the Helmoltz solution
• superimpose geodesics to both plots
• on the amplitude plot, what is the meaning of the ripple effects at/around the caustic ?

Problem 4 (breakdown of the amplitude approximation): For the speed of 'lens' type,

• Complete the laplacian function so as to compute the laplacian of a 2D array using a 5-point stencil (see class), or any method of your choice.
• Compute the Laplacian of \(\tau\) (obtained after running an eikonal solve) for various resolutions, say n=200,300,400,500. Based on the values displayed by the plt.colorbar() command for these various resolutions, does the function \(\Delta \tau\) seem to converge to a function in the classical sense ? Where does the trouble come from ?
• Based on the fact that the amplitude \(a(x,x_s)\) of the GO approximation satisfies, along a characteristic curve,

\[ \frac{d}{dt}a(X(t), x_s) = - \frac{c^2(X(t))a(X(t),x_s)}{2} \Delta\tau(X(t)), \]

explain how the model for the amplitude of the GO approximation also breaks down after the caustics.