examples: Born approximation for TTI media #1555

ofmla · 2021-01-15T14:21:02Z

At this point @FabioLuporini might be asking: where is @ofmla? xD well, I've been extra-busy with some other projects and preparing classes, but here is a small contribution. This PR adds new operators to the anisotropic solver. In the absence of Mathias @jkwashbourne could take a look.

ofmla · 2021-01-15T14:43:07Z

Note: the perturbation is only in the velocity in the Born approx.

codecov · 2021-01-15T15:30:32Z

Codecov Report

Merging #1555 (73b4a9f) into master (72023d2) will decrease coverage by 27.88%.
The diff coverage is 21.95%.

@@             Coverage Diff             @@
##           master    #1555       +/-   ##
===========================================
- Coverage   86.30%   58.42%   -27.89%     
===========================================
  Files         217      200       -17     
  Lines       31428    29921     -1507     
  Branches     4186     3983      -203     
===========================================
- Hits        27125    17482     -9643     
- Misses       3818    11849     +8031     
- Partials      485      590      +105

Impacted Files	Coverage Δ
tests/test_gradient.py	`21.83% <18.75%> (-76.77%)`	⬇️
tests/test_adjoint.py	`26.66% <20.00%> (-73.34%)`	⬇️
examples/seismic/tti/wavesolver.py	`38.31% <22.22%> (-58.51%)`	⬇️
examples/seismic/tti/operators.py	`34.76% <23.07%> (-61.08%)`	⬇️
examples/seismic/viscoacoustic/operators.py	`9.04% <0.00%> (-90.96%)`	⬇️
tests/test_dimension.py	`9.10% <0.00%> (-90.90%)`	⬇️
examples/seismic/self_adjoint/operators.py	`9.23% <0.00%> (-90.77%)`	⬇️
tests/test_caching.py	`13.51% <0.00%> (-86.25%)`	⬇️
tests/test_dse.py	`14.04% <0.00%> (-85.59%)`	⬇️
tests/test_lower_exprs.py	`11.34% <0.00%> (-85.57%)`	⬇️
... and 161 more

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 72023d2...2de2bfe. Read the comment docs.

jkwashbourne · 2021-01-15T16:14:45Z

examples/seismic/tti/operators.py

+    m = model.m
+    vp = model.vp
+
+    time_order = 1 if kernel == 'staggered' else 2


we don't handle time_order = 2 right? should maybe catch as an error?

aha, I thought this was a first order op, first time I really looked ...

jkwashbourne

Missing the "linearization test". For me the proof of correctness is passing the linearization and adjoint tests.

devito/tests/test_gradient.py

Line 147 in 4a52d5b

def test_gradientFWI(self, dtype, shape, kernel, space_order, checkpointing):

I guess for completeness could ensure the operator replicates acoustic modeling when isotropic. Not sure if people feel that is important or now, but it is like that acoustic op is nailed down via the accuracy notebook, so we can test that other ops reproduce acoustic isotropic wavefields.

jkwashbourne · 2021-01-15T16:16:46Z

examples/seismic/tti/operators.py

+    m = model.m
+    vp = model.vp
+
+    time_order = 1 if kernel == 'staggered' else 2


aha, I thought this was a first order op, first time I really looked ...

mloubout

Reparametrize in term of squared slowness m not velocity like for the rest of the operators and to match the examples and test setups.

examples/seismic/tti/operators.py

ofmla · 2021-01-26T16:43:22Z

Missing the "linearization test". For me the proof of correctness is passing the linearization and adjoint tests.

devito/tests/test_gradient.py

Line 147 in 4a52d5b

def test_gradientFWI(self, dtype, shape, kernel, space_order, checkpointing):

I guess for completeness could ensure the operator replicates acoustic modeling when isotropic. Not sure if people feel that is important or now, but it is like that acoustic op is nailed down via the accuracy notebook, so we can test that other ops reproduce acoustic isotropic wavefields.

I added that test as suggested

georgebisbas

Am I wrong, or it needs a rebase on latest master for easier review ?
Thanks Oscar! Hope you are doing well!

ofmla · 2021-01-27T12:24:10Z

Am I wrong, or it needs a rebase on latest master for easier review ?
Thanks Oscar! Hope you are doing well!

Hi George, all good, thx. You are correct, I made a mess. I hope that everything is in order now.

mloubout · 2021-01-29T13:27:54Z

examples/seismic/tti/wavesolver.py

+            The time-constant Azimuth angle (radians).
+        """
+        if kernel != 'centered':
+            raise RuntimeError('Only centered kernel is supported for the jacobian')


argument error not runtime

ValueError right?

it also was resolved.

examples/seismic/tti/wavesolver.py

mloubout · 2021-01-29T13:30:25Z

tests/test_gradient.py

+        (np.float64, 4, 'OT2', (70, 80), True, setup),
+        (np.float64, 4, 'OT2', (70, 80), False, setup),
+        (np.float64, 4, 'centered', (70, 80), True, tti_setup),
+        (np.float64, 4, 'centered', (70, 80), False, tti_setup),


no float32?

I forgot to add np.float32. I just tried it and the test failed (AssertionError) for both Checkpointing=True and False even for the first order error

Operator `initdamp` run in 0.01 s Operator `pad_vp` run in 0.01 s Operator `pad_epsilon` run in 0.01 s Operator `pad_delta` run in 0.01 s Operator `pad_theta` run in 0.01 s Operator `smoother` run in 0.01 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.03 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.03 s 1st order error, Phi(m0+dm)-Phi(m0): [ 0.87556505 -5.82210223] 2nd order error, Phi(m0+dm)-Phi(m0) - <J(m0)^T \delta d, dm>: [ 0.87554379 -5.82207074] FOperator `initdamp` run in 0.01 s Operator `pad_vp` run in 0.01 s Operator `pad_epsilon` run in 0.01 s Operator `pad_delta` run in 0.01 s Operator `pad_theta` run in 0.01 s Operator `smoother` run in 0.01 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `GradientTTI` run in 0.03 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s Operator `ForwardTTI` run in 0.02 s 1st order error, Phi(m0+dm)-Phi(m0): [ 0.87556505 -5.82210223] 2nd order error, Phi(m0+dm)-Phi(m0) - <J(m0)^T \delta d, dm>: [ 0.87554379 -5.82207074]

I don't really know what happened as with np.float64 the test runs without problems

Operator `initdamp` run in 0.01 s Operator `pad_vp` run in 0.01 s Operator `pad_epsilon` run in 0.01 s Operator `pad_delta` run in 0.01 s Operator `pad_theta` run in 0.01 s Operator `smoother` run in 0.01 s Operator `ForwardTTI` run in 0.03 s Operator `ForwardTTI` run in 0.04 s Operator `ForwardTTI` run in 0.12 s Operator `ForwardTTI` run in 0.07 s Operator `ForwardTTI` run in 0.07 s Operator `ForwardTTI` run in 0.07 s Operator `ForwardTTI` run in 0.06 s Operator `ForwardTTI` run in 0.06 s Operator `ForwardTTI` run in 0.07 s 1st order error, Phi(m0+dm)-Phi(m0): [ 0.92561436 -10.858962 ] 2nd order error, Phi(m0+dm)-Phi(m0) - <J(m0)^T \delta d, dm>: [ 1.97902976 -10.94126684] .Operator `initdamp` run in 0.01 s Operator `pad_vp` run in 0.01 s Operator `pad_epsilon` run in 0.01 s Operator `pad_delta` run in 0.01 s Operator `pad_theta` run in 0.01 s Operator `smoother` run in 0.01 s Operator `ForwardTTI` run in 0.03 s Operator `ForwardTTI` run in 0.05 s Operator `GradientTTI` run in 0.05 s Operator `ForwardTTI` run in 0.05 s Operator `ForwardTTI` run in 0.05 s Operator `ForwardTTI` run in 0.06 s Operator `ForwardTTI` run in 0.06 s Operator `ForwardTTI` run in 0.07 s Operator `ForwardTTI` run in 0.07 s Operator `ForwardTTI` run in 0.07 s 1st order error, Phi(m0+dm)-Phi(m0): [ 0.92561436 -10.858962 ] 2nd order error, Phi(m0+dm)-Phi(m0) - <J(m0)^T \delta d, dm>: [ 1.97902976 -10.94126684]

I'll have to look it up and that might take a while

I added a new commit with the test for np.float32 and several tests from CI failed. I review the code in operators.py and wavesolver.py and I couldn't find an error. Then, I tried without the optimization flags (-O3 and -ffast-math) but the test still fails. @mloubout if you have any suggestions or comments about this, please let me know.

There is definitely something off, I would understand if things were a bit off but that's way off. Can you plot the convergence curves and such to see how it behaves and maybe make sure model/data looks ok.

I discovered what the problem was. The tti_setup function uses a default value of tn=250.0, so only the direct wave is registered, doubling the time solves the problem, thx.

Ah, now I see why the bump to 500... OK. Instead of changing it from tti_setup and run, can you pass the value from the tests only? also, is there a smaller value such that the forward wave is recorded, such as 300 or 350? or a smaller grid? we need to bring down the turnaround time of these tests

I've reduced the grid dimensions as well as the recording time to 400.

georgebisbas · 2021-01-29T15:20:09Z

tests/test_gradient.py

+        (np.float64, 4, 'centered', (70, 80), True, tti_setup),
+        (np.float64, 4, 'centered', (70, 80), False, tti_setup),
+    ])
+    def test_gradientFWI(self, dtype, space_order, kernel, shape, ckp, setup_func):
        r"""


is this a leftover?

It seems so, and now I realize there are others in that file, thx.

FabioLuporini

There is a > 5 minutes increase in turnaround time in each build . Can we reduce this. Options incluse:

Bringing back that tn=500 to tn=250 should already have an impact, no?
Dropping either some of the so=4 or some of the so=8
Dropping some of the double precision?

examples/seismic/tti/operators.py

FabioLuporini · 2021-02-17T07:49:39Z

examples/seismic/tti/tti_example.py

@@ -21,7 +21,7 @@ def tti_setup(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
    return AnisotropicWaveSolver(model, geometry, space_order=space_order, **kwargs)


-def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
+def run(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=500.0,


I don't like this bump to 500. I don't see why it's necessary. I'm used to the turnaround time of python tti_example.py which I find not too short not too long, doubling it would make it slow for no evident reason

tn was set again to be 250.0. It only is changed in the specific test as you suggested.

FabioLuporini · 2021-02-17T07:49:52Z

examples/seismic/tti/tti_example.py

@@ -7,7 +7,7 @@
 from examples.seismic.tti import AnisotropicWaveSolver


-def tti_setup(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=250.0,
+def tti_setup(shape=(50, 50, 50), spacing=(20.0, 20.0, 20.0), tn=500.0,


Don't like this bump to 500. See comment below

tn was set again to be 250.0. It only is changed in the specific test as you suggested.

examples/seismic/tti/wavesolver.py

…ror if kernel != centered

…loat32

…nd minimize CI times

FabioLuporini

LGTM. Also approved offline by @mloubout

FabioLuporini · 2021-02-19T08:33:25Z

Thanks, merged.

jkwashbourne reviewed Jan 15, 2021

View reviewed changes

jkwashbourne requested changes Jan 15, 2021

View reviewed changes

mloubout added the examples examples label Jan 19, 2021

mloubout reviewed Jan 19, 2021

View reviewed changes

examples/seismic/tti/operators.py Outdated Show resolved Hide resolved

examples/seismic/tti/operators.py Outdated Show resolved Hide resolved

georgebisbas reviewed Jan 27, 2021

View reviewed changes

ofmla force-pushed the born_approx_tti branch from 721f4ad to 12ad366 Compare January 27, 2021 12:11

FabioLuporini assigned mloubout Jan 29, 2021

ofmla force-pushed the born_approx_tti branch from 12ad366 to 4f9543a Compare January 29, 2021 12:49

mloubout reviewed Jan 29, 2021

View reviewed changes

georgebisbas reviewed Jan 29, 2021

View reviewed changes

ofmla force-pushed the born_approx_tti branch 3 times, most recently from 8cefa00 to 29c0b43 Compare February 16, 2021 15:54

FabioLuporini requested changes Feb 17, 2021

View reviewed changes

ofmla force-pushed the born_approx_tti branch from 6165b83 to cc0c465 Compare February 18, 2021 15:18

ofmla added 10 commits February 18, 2021 21:05

examples: add born modeling and adjoint migration operators

e0e72ed

examples: add operators to solver

7b9c86c

examples: add Jacobian operator adjoint test for TTI

1494102

examples: fix flake8

11d79e5

examples: change parametrization & add test_gradientFWI

ea77186

examples: add born modeling and adjoint migration operators

84c0c05

examples: add np.float32 (tti_setup) test_gradientFWI & raise ValueEr…

450b74b

…ror if kernel != centered

examples: doubling tn to make sure that tes_gradientFWI work for np.f…

afa8f5e

…loat32

examples: small corrections/changes in order to reduce CI time

f520922

examples: size reduction in test_gradientFWI to maintain uniformity a…

2de2bfe

…nd minimize CI times

ofmla force-pushed the born_approx_tti branch from cc0c465 to 2de2bfe Compare February 19, 2021 00:23

FabioLuporini approved these changes Feb 19, 2021

View reviewed changes

FabioLuporini merged commit d8721eb into devitocodes:master Feb 19, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

examples: Born approximation for TTI media #1555

examples: Born approximation for TTI media #1555

ofmla commented Jan 15, 2021 •

edited

Loading

ofmla commented Jan 15, 2021

codecov bot commented Jan 15, 2021 •

edited

Loading

jkwashbourne Jan 15, 2021

jkwashbourne Jan 15, 2021

jkwashbourne left a comment •

edited

Loading

jkwashbourne Jan 15, 2021

mloubout left a comment

ofmla commented Jan 26, 2021

georgebisbas left a comment

ofmla commented Jan 27, 2021

mloubout Jan 29, 2021

ofmla Jan 29, 2021

ofmla Feb 17, 2021

mloubout Jan 29, 2021

ofmla Jan 29, 2021 •

edited

Loading

ofmla Feb 1, 2021

mloubout Feb 9, 2021

ofmla Feb 10, 2021

FabioLuporini Feb 17, 2021

ofmla Feb 17, 2021

georgebisbas Jan 29, 2021

ofmla Jan 29, 2021

FabioLuporini left a comment

FabioLuporini Feb 17, 2021

ofmla Feb 17, 2021

FabioLuporini Feb 17, 2021

ofmla Feb 17, 2021

FabioLuporini left a comment

FabioLuporini commented Feb 19, 2021

examples: Born approximation for TTI media #1555

examples: Born approximation for TTI media #1555

Conversation

ofmla commented Jan 15, 2021 • edited Loading

ofmla commented Jan 15, 2021

codecov bot commented Jan 15, 2021 • edited Loading

Codecov Report

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jkwashbourne left a comment • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mloubout left a comment

Choose a reason for hiding this comment

ofmla commented Jan 26, 2021

georgebisbas left a comment

Choose a reason for hiding this comment

ofmla commented Jan 27, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ofmla Jan 29, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FabioLuporini left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FabioLuporini left a comment

Choose a reason for hiding this comment

FabioLuporini commented Feb 19, 2021

ofmla commented Jan 15, 2021 •

edited

Loading

codecov bot commented Jan 15, 2021 •

edited

Loading

jkwashbourne left a comment •

edited

Loading

ofmla Jan 29, 2021 •

edited

Loading