variable order density acoustic wave equation

[krober10nd]

• we should implement the same spatial orders (i.e., up to P=16) for the variable density acoustic wave equation.
• Stopping at order 2 is not sufficient for the article in my opinion.

[jaimesouza]

To do that, we need discretization for higher space orders. Apparently the adaptation is not straightforward as in the case of constant density.

[krober10nd]

Yea, the arbitrary order discretization is the old Overleaf document.

[jaimesouza]

There was an arbitrary order discretization (for variable density), but it had some issues. We have tried to implement it, but the discretization was wrong.

[krober10nd]

I don't recall that. The discretization is straight from a textbook.

[joao-bapdm]

There was an arbitrary order discretization (for variable density), but it had some issues. We have tried to implement it, but the discretization was wrong.

where did you take it from? Was it the one I passed you? Is it on Overleaf?

[jaimesouza]

@joao-bapdm did it, right @joao-bapdm ? It was on https://www.overleaf.com/project/5ef60c808491e40001ee1749

[joao-bapdm]

Is it currently in the code?

[jaimesouza]

No, It is not.

[krober10nd]

The actual equation that needs to be implemented is this:

which is essentially just one more term than the standard constant density wave propagator.

[krober10nd]

keep in mind that the gradient term is a dot product so you have to expand the x and y derivatives and compute the dot product between them.

[jaimesouza]

Can you write a discretized equation? Like in https://www.researchgate.net/publication/313765589_3D_Isotropic_Acoustic_Finite-Difference_Wave_Equation_Code_A_Many-Core_Processor_Implementation_and_Analysis

[krober10nd]

Sure thing.

[jaimesouza]

Great. Thanks!

[krober10nd]

Here you go:

It's also now on our Overleaf document "StencilCode work".

I tried to be faithful to the document you linked above in notation.

[jaimesouza]

Awesome. I am gonna implement it.

[jaimesouza]

Just one detail:

For second-order derivative Viis mirrored:

>>> findiff.coefficients(deriv=2, acc=2)['center']['coefficients']
array([ 1., -2.,  1.])

>>> findiff.coefficients(deriv=2, acc=4)['center']['coefficients']
array([-0.08333333,  1.33333333, -2.5       ,  1.33333333, -0.08333333])

Every neighbor (withidistance from central) has the same coefficient Vi.

For first-order derivative, Wi is is inverted (+ or -):

>>> findiff.coefficients(deriv=1, acc=2)['center']['coefficients']
array([-0.5,  0. ,  0.5])

>>> findiff.coefficients(deriv=1, acc=4)['center']['coefficients']
array([ 0.08333333, -0.66666667,  0.        ,  0.66666667, -0.08333333])

Every neighbor (withi distance from central) has the it's own coefficient Wi.

[krober10nd]

yes, nice catch so the sign is negative then.

[jaimesouza]

I think the positive sign is correct, but each neighbor must have it's own W. For example:

(Wi1 * A(x+i,y)) + (Wi2 * A(x-i,y))

Actually, It (mirrored or inverted) is not general for higher orders.

>>> findiff.coefficients(deriv=2, acc=40)['center']['coefficients']
array([-1.13053274e-18,  5.10002009e-17, -1.12795295e-15,  1.62971722e-14,
       -1.72903714e-13,  1.43487153e-12, -9.68135853e-12,  5.44309419e-11,
       -2.58655786e-10,  1.04260583e-09, -3.51796416e-09,  9.41282242e-09,
       -1.59424003e-08, -1.20980925e-08,  2.03509442e-07, -4.46634191e-09,
       -2.23897479e-05,  6.07536917e-04, -2.44297594e-02,  6.21361254e+00,
       -1.23777807e+01,  6.21196796e+00, -2.41143293e-02,  7.01280168e-04,
       -4.01799412e-05,  3.09901734e-06, -2.62774410e-07,  2.08958937e-08,
       -1.21531619e-09, -3.10406330e-12,  1.29815547e-11, -1.87570750e-12,
        1.23813981e-13,  8.36333182e-16, -4.45754518e-16, -1.15077553e-16,
        3.09318071e-17, -2.59401813e-18, -1.08898448e-20,  1.64577467e-20,
       -8.54423500e-22])

>>> findiff.coefficients(deriv=1, acc=40)['center']['coefficients']
array([ 2.35587290e-15, -1.05892978e-13,  2.33062032e-12, -3.34533738e-11,
        3.51757155e-10, -2.88306272e-09,  1.91107809e-08, -1.04653275e-07,
        4.77069794e-07, -1.78908021e-06,  5.20332490e-06, -8.84837384e-06,
       -1.73282766e-05,  2.48622091e-04, -1.35895605e-03,  3.57929131e-03,
        3.84132748e-02, -1.40421277e+00,  6.02842212e+01, -9.43893980e+03,
        1.87563245e+04, -9.43495928e+03,  5.98624413e+01, -1.72732638e+00,
        9.81282851e-02, -7.46391379e-03,  6.16078210e-04, -4.59687635e-05,
        2.08588568e-06,  1.42593476e-07, -4.92723098e-08,  6.42163915e-09,
       -4.43856326e-10,  4.55246526e-12,  9.74690928e-13,  2.79858242e-13,
       -7.05515600e-14,  4.16147002e-15,  4.19654506e-16, -7.55345208e-17,
        3.41672095e-18])

I am gonna have to fix that on constant density code as well. Each neighbor must have it's own coefficient. even when they are the same for lower orders.

[krober10nd]

I wasn't aware that the schemes weren't symmetric past a certain point.

[jaimesouza]

Me neither. I have to check that.

[krober10nd]

What about up to order 16? Do they remain symmetric?

[jaimesouza]

Yes, they do.

[jaimesouza]

>>> findiff.coefficients(deriv=2, acc=16)['center']['coefficients']
array([-2.42812743e-06,  5.07429079e-05, -5.18000518e-04,  3.48096348e-03,
       -1.76767677e-02,  7.54208754e-02, -3.11111111e-01,  1.77777778e+00,
       -3.05484410e+00,  1.77777778e+00, -3.11111111e-01,  7.54208754e-02,
       -1.76767677e-02,  3.48096348e-03, -5.18000518e-04,  5.07429079e-05,
       -2.42812743e-06])

>>> findiff.coefficients(deriv=1, acc=16)['center']['coefficients']
array([ 9.71250971e-06, -1.77600178e-04,  1.55400155e-03, -8.70240870e-03,
        3.53535353e-02, -1.13131313e-01,  3.11111111e-01, -8.88888889e-01,
       -3.35961124e-10,  8.88888889e-01, -3.11111111e-01,  1.13131313e-01,
       -3.53535354e-02,  8.70240870e-03, -1.55400155e-03,  1.77600178e-04,
       -9.71250971e-06])

[krober10nd]

I think findiff has trouble finding the coefficients for very high orders. They should remain symmetric.

[krober10nd]

>>> findiff.coefficients(deriv=1, acc=40)['center']['coefficients']
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/Keith/Desktop/PostDoctoralWork/codes/firedrake/lib/python3.7/site-packages/findiff/coefs.py", line 64, in coefficients
    ret["center"] = calc_coefs(deriv, offsets, symbolic)
  File "/Users/Keith/Desktop/PostDoctoralWork/codes/firedrake/lib/python3.7/site-packages/findiff/coefs.py", line 99, in calc_coefs
    coefs = np.linalg.solve(matrix, rhs)
  File "<__array_function__ internals>", line 6, in solve
  File "/Users/Keith/Desktop/PostDoctoralWork/codes/firedrake/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 393, in solve
    r = gufunc(a, b, signature=signature, extobj=extobj)
  File "/Users/Keith/Desktop/PostDoctoralWork/codes/firedrake/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 88, in _raise_linalgerror_singular
    raise LinAlgError("Singular matrix")
numpy.linalg.LinAlgError: Singular matrix

[jaimesouza]

Yes, thats it. findiff has this problem.

[krober10nd]

Realistically, up to order 20 is as high as most people would go. Could you put a warning in alerting users?

[jaimesouza]

Sure! I am gonna limit it up to 20 and add an assert.

[jaimesouza]

https://github.com/HPCSys-Lab/simwave/pull/35