4
$\begingroup$

I am trying to compare my finite difference's solution of the scalar (or simple acoustic) wave equation with an analytic solution.

For that purpose I am using the following analytic solution presented in the old paper Accuracy of the finite-difference modeling of the acoustic wave equation - Geophysics 1974 - R. M. Alford et. al.

The solution, equation (9) in that paper, for a constant velocity medium is given in cylindrical coordinates by:

\begin{eqnarray} u_s(\rho, \phi, t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \{ -i \pi H_0^{(2)}\left(k | \sigma - \sigma_s| \right) F(\omega) e^{i\omega t} d\omega\} \end{eqnarray} (A)

$H_0^{(2)}$ is the Hankel function of the second kind with order zero.

$| \sigma - \sigma_s| $ is the distance from the source to the observation point.

Also $ k = \frac{\omega}{C_0}$ come from $C_0$ the medium velocity converted to the Fourier domain, defined in the original scalar wave equation as:

\begin{eqnarray} \left[ \nabla^2 - \frac{1}{C_0^2}\frac{\partial ^2}{\partial t^2}\right] u(\rho, \phi, t) = -4 \pi \frac{\delta(\rho-\rho_s)\delta(\phi-\phi_s)f(t)}{\rho} \end{eqnarray}

$ (\rho_s, \phi_s) $ defines the source location in the $\phi$ x $\rho$ plane.

$ F(\omega)$ comes from the source function $f(t)$

The problem:

  • When implementing numerically A I'm not getting meaningful results. The velocity obtained using A is totally different from the $C_0$ used as input

My code example is bellow where I am evaluating A, using python. I calculate the solution at zero solution0 and at a distance D=8000 solutionD. Then I use the global maximum to calculate an approximation to velocity $C_0$

Probably it's a quite simple error that I am doing. If someone can point it out for me. The velocity I am getting shows that maybe the equation (9) is wrong or I am implementing it wrongly (what I believe to be the case).

from scipy.special import hankel2
import numpy as np

# a simple source function f(t) gauss derivative normalized
def gaussource(time, wlength, delay=None):
        if delay == None: # enough delay time
            delay = 3*wlength
        t = time - delay
        return ((2*np.sqrt(np.e)/(wlength))
               *t*np.exp(-2*(t**2)/(wlength**2)))

# analytic solution equation (A) paper
def waves(source, distance, c):
    n = len(source)
    sourcew = np.fft.fft(source) # source in the frequency domain
    hnk_factor = distance/c # distance to source divided by velocity
    if not distance == 0.0:
        hankelshift = np.complex(0,-np.pi)*np.array([hankel2(0,hnk_factor*omega) for omega in xrange(n)])
        hankelshift[0] = 0. # i infinity in limit
        return np.real(np.fft.ifft(hankelshift*sourcew))
    return source

c = 4000.
fc = 5.0
dt = 0.01
t = np.arange(0., 2., dt) # 2 seconds to evaluate the solution
source = gaussource(t, 1./fc) # source
solution0 = waves(source, 0., c) # solution at zero source (source itself)
solutionD = waves(source, 8000., c) #solution at certain distance D=8000 metres e.g.   
t0 = solution0.argmax(); # get max index time at zero
tD = solutionD.argmax(); # get max index time at distance
print "phase velocity aproximation ", 8000./((tD-t0)*dt) 
$\endgroup$
0
$\begingroup$

After some weeks and no one found the answer to my problem. I found myself the errors and the solution. To help the community I'am sharing it (meta advice).

It was indeed a very simple error.

Really basic sampling theorem, maximum frequency based on sample rate is:

$$ f_{max} = 1/\Delta t $$

or using $ \omega $

$$ \omega_{max} = 2\pi/\Delta t $$

I changed the def waves(source, distance, c): function definition. I used $\Delta t$ to calculate $ \Delta \omega $ to sample the solution in frequency. Also removed the calculation at zero distance since the source is the solution itself.

# analytic solution equation (A) paper
def waves(source, distance, c, dt):
    n = len(source)
    sourcew = np.fft.fft(source) # source in the frequency domain
    hnk_factor = distance/c # distance to source divided by velocity
    dw = np.pi*2.0*(1./dt)*1/n
    hankelshift = np.complex(0,-np.pi)*np.array([hankel2(0.0,hnk_factor*k*dw) for k in xrange(n)])
    hankelshift[0] = 0. # i infinity in limit
    return np.real(np.fft.ifft(hankelshift*sourcew))

t = np.arange(0., 3.5, dt) # 3.5 seconds to evaluate the solution
source = gaussource(t, 1./fc) # source and solution at 0
solutionD = waves(source, 8000., c, dt) # solution at certain distance     D=8000 metres e.g.   
t0 = source.argmax(); # get max index time at zero
tD = solutionD.argmax(); # get max index time at distance

I also increased the solution time spam to 3.5 seconds.

And bellow is how it did work, it even can give closer results changing the sample rate.

phase velocity approximation 3883.49514563

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.