seg / tutorials-2018

Geophysical Tutorials column for 2018
118 stars 71 forks source link

Acoustic wave equation is not self-adjoint #3

Closed ediaz closed 6 years ago

ediaz commented 6 years ago

I've been following these very exciting tutorials. The aim for simplicity is superb as it translates the equations to "human-readable'' implementations.

I propose an amendment for one statement of your last tutorial and it's corresponding implementation.

The acoustic wave equation is self-adjoint only if the background model is constant (which is not a very exciting medium).

If we think about the wave equation without absorbing boundaries we need to implement the following steps at each time:

  1. apply Laplacian operator to u[x,z]
  2. scale output 1. by v^2[x,z]
  3. Add source term
  4. Do time-stepping

In adjoint mode, we need to perform this process in backward order:

  1. Do time-stepping (self-adjoint)
  2. Add source term
  3. Scale 2. by v^2[x,z]
  4. apply laplacian to output from 3.

Only if v^2[x,z] is constant, you can interchange 4. and 3.

In terms of matrices, steps 3. and 4. of the forward operator can be represented by: C=VL, with V being a diagonal matrix with the squared velocity in its diagonal, and L being the Laplacian operator. In adjoint mode, we need to apply C^t= LV (L and V are self-adjoint).

So, the adjoint equation is;

To test if your forward and adjoint implementation is correct you can do the dot-product test

I hope this point helps,

Cheers, Esteban

kwinkunks commented 6 years ago

Thank you for this remark, Esteban. I agree, these tutorials are really exciting. Nice to have some discussion!

I am pinging the authors @philippwitte and @mloubout so they can respond.

mloubout commented 6 years ago

Hi @ediaz, thank you for the input.

TO answer your question, if you look at the paper/notebook, you'll see that the squared slowness is the coefficient in front of the time derivative, not the laplacian. This is the conventional way to write it (e.g https://en.wikipedia.org/wiki/Acoustic_wave_equation)

From this, as the squared slowness is not time dependent, the equation is then self-adjoint. It is actually tested via the dottest at https://github.com/opesci/devito/blob/master/tests/test_adjointA.py.

In the case you would write it with the squared slowness in front of the laplacian, it would indeed not be self-adjoint anymore, however the source term would be scaled by the squared slowness too and would required to be taken into account for the derivation of gradient of the FWI objective (see part2 and coming part 3 of the tutorial).

Thank you again for your interest in the tutorial and feel free to ask any question you would have.

ediaz commented 6 years ago

Thanks for the follow-up, guys.

I need to go into the details of the implementation, but regardless of how you write the equation when you solve for u[t+dt] you would need to scale the laplacian (and the source) by v2.

I looked at your tests, and I see that you use several velocity models (non-homogeneous), so it seems that the forward-adjoint pair is correctly implemented. It doesn't prove the self-adjointness claim.

Here is a test if you're up for it:

  1. use a constant velocity model with some random noise
  2. Create an acoustic wave equation without absorbing boundary conditions.
  3. forward model a random source
  4. time-reverse another random source, forward model with it (this is your self-adjoint claim)
  5. try the dot product test

If the acoustic equation is self-adjoint as you claim, then if you turn off the boundary conditions, your adjoint should be equal to the time-reversed forward operator.

Cheers, Esteban