examples: Boundary reconstruction #1926

fffarias · 2022-05-23T21:00:11Z

Hello everybody.

This PR tries to suggest reconstructing the source wave field from the one saved in the border. This should be an option for the gradient computation, in which the forward field is reconstructed inside the backwards loop. It was necessary to put a boolean argument inside the SparseTimeFunction's inject method, since I needed the injected field not to be incremented. Given that I changed an important part of the code, I added a test to validate the option of not incrementing. As an example, I also add a demonstration tutorial and application of the method in FWI.

Thanks a lot

review-notebook-app · 2022-05-23T21:00:15Z

Check out this pull request on

See visual diffs & provide feedback on Jupyter Notebooks.

Powered by ReviewNB

FabioLuporini · 2022-05-24T07:50:20Z

devito/operations/interpolators.py

+                            implicit_dims=self.sfunction.dimensions)
+                        for b, vsub in zip(self._interpolation_coeffs, idx_subs)]
+            else: 
+                eqns = [Eq(field.xreplace(idx_subs[0]), _expr.xreplace(idx_subs[0]),


why are you picking idx_subs[0]? should you not traverse the whole list of idx_subs like in the increment case?

Could you do something like:

if increment: cls: Inc else: cls: Eq eqns = [cls(field.xreplace(....... .....) for b, vsub in zip(....)]

?

because I thought in this case very specifically to apply boundary conditions, not source injection, for example. The ` increment=False option as I suggested should not be used for off-grid points.

While this is true for your case, the generic case is still a standard potentially off the grid point injection without increment. So there is two option

Use the injection as is knowing the coefficient will still be zero for the interpolation if the point is on the grid

Define somehow an "on-the-grid" sparse function that ignores the interpolation.

The first case is definitely easier and should still work for your case and would boil down to what @FabioLuporini wrote earlier

I believe the first option is not possible for the case increment = False. Because even if the point is exactly on the grid, interpolation coefficients equal to zero will vanish the fields at points that should not be vanished. Which is not the case for injecting the source (when increment = True), for example, since accumulating zeros does not change the field at the other points involved in the interpolation.

The way I see it, we are left with the second option that @mloubout suggested, or use subdomains that are under development in another PR if I'm not mistaken.

Indeed good point I missed that gotta think about it

@mloubout prodding

Getting back here since been looking into sparse, I think that the "best" way to go for increment=Flase is to do what interpolation does i.e

Define sum and do the standard sum += inject(expr)

assign sum afterward field = sum

This would cover the off-the-grid case as well

FabioLuporini · 2022-05-24T07:50:47Z

tests/test_interpolation.py

+    ((10, 10, 10), 1., False),
+    ((10, 10, 10), 5., True)        
+])
+def test_inject_time_increment(shape, result, increment):


Thank you for attaching tests to this PR!

Would it be worth having another test to check that increment=False specified on SparseFunction with non-grid-aligned points is caught and handled appropriately?

FabioLuporini · 2022-05-24T07:53:02Z

Thank you, Fernanda! We will review this PR carefully. Looks nice!

Meanwhile, note that some of the commit messages would have to be edited, via rebasing, to use valid tags: https://github.com/devitocodes/devito/wiki/Tags-for-commit-messages-and-PR-titles

…ng interpolation weights.

…e increment option in the SparseTimeFunction inject method

mloubout

Thanks so much for the contribution @fffarias this is very much appreciated. I left a couple comments for the tutorial to make it a bit more pedagogic, happy to discuss what's doable without adding too much on your plate.

mloubout · 2022-05-24T19:21:30Z

devito/operations/interpolators.py

+                            implicit_dims=self.sfunction.dimensions)
+                        for b, vsub in zip(self._interpolation_coeffs, idx_subs)]
+            else: 
+                eqns = [Eq(field.xreplace(idx_subs[0]), _expr.xreplace(idx_subs[0]),


While this is true for your case, the generic case is still a standard potentially off the grid point injection without increment. So there is two option

Use the injection as is knowing the coefficient will still be zero for the interpolation if the point is on the grid

Define somehow an "on-the-grid" sparse function that ignores the interpolation.

The first case is definitely easier and should still work for your case and would boil down to what @FabioLuporini wrote earlier

mloubout · 2022-05-24T19:25:00Z

examples/seismic/tutorials/16_source_wavefield_reconstruction.ipynb

@@ -0,0 +1,1069 @@
+{


Nice tutorial. Some generic comments.

I think it would be beneficial to split and/or extend this tutorial for readability. I would

first focus purely on the wavefield itself and show snapshots/movie/... of the forward wavefield and its reconstruction (i.e check the snapshot tutorial for ploting)

Second detail the gradient computation

And on the into side, It would be beneficial to use slightly more explicit variable names, such as u and u_boundary or such. I would also completely skip the illumination part. While this is indeed very useful in practice, it adds to the complexity of the tutorial.

mloubout · 2022-05-24T19:26:17Z

tests/test_interpolation.py

+    a = SparseTimeFunction(name=name, grid=grid, npoint=npoints, nt=nt)
+    dims = tuple([np.linspace(0., 1., d) for d in grid.shape]) 
+    for i in range(len(grid.shape)):
+        a.coordinates.data[:,i] = np.meshgrid(*dims)[i].flatten()


can probably compute meshgrid only once before the for loop

EdCaunt · 2022-05-31T13:37:23Z

examples/seismic/tutorials/16_source_wavefield_reconstruction.ipynb

@@ -0,0 +1,1069 @@
+{


Might be good to have this presented as percentage of max wavefield value?

Reply via ReviewNB

mloubout · 2022-06-08T18:27:12Z

Generic question. Would it be easier to use 1D Function instead of Sparse functions for this?

fffarias · 2022-06-09T13:47:46Z

Generic question. Would it be easier to use 1D Function instead of Sparse functions for this?

possibly. I'm not sure how to get just the necessary points from and into the wavefield. Would I have to use subdomains you think? Also it may be more practical to use 2D function in space, to accommodate more layers according to the operator's spatial order.

FabioLuporini reviewed May 24, 2022

View reviewed changes

FabioLuporini changed the title ~~Boundary reconstruction~~ examples: Boundary reconstruction May 24, 2022

FabioLuporini requested review from mloubout and EdCaunt May 24, 2022 07:51

FabioLuporini added the examples examples label May 24, 2022

fffarias force-pushed the BoundaryReconstruction branch from ad7916e to 4952720 Compare May 24, 2022 14:21

fffarias added 4 commits May 24, 2022 14:38

compiler: Add increment boll to inject method in SparseTimeFunction

77aec6c

compiler: Add condition not to increment during injection, disregardi…

4feee17

…ng interpolation weights.

examples: Add source wavefield reconstruction strategy tutorial

e3d57c1

tests: Add tests referring to commits 3edf20 and 118ec4, which add th…

67d3f34

…e increment option in the SparseTimeFunction inject method

fffarias force-pushed the BoundaryReconstruction branch from 4952720 to 67d3f34 Compare May 24, 2022 17:42

mloubout reviewed May 24, 2022

View reviewed changes

EdCaunt reviewed May 31, 2022

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

examples: Boundary reconstruction #1926

examples: Boundary reconstruction #1926

fffarias commented May 23, 2022

review-notebook-app bot commented May 23, 2022

FabioLuporini May 24, 2022

fffarias May 24, 2022

mloubout May 24, 2022

fffarias May 27, 2022

mloubout May 27, 2022

FabioLuporini Jun 22, 2022

mloubout Jun 28, 2023 •

edited

FabioLuporini May 24, 2022

EdCaunt May 31, 2022

FabioLuporini commented May 24, 2022

mloubout left a comment

mloubout May 24, 2022

mloubout May 24, 2022

mloubout May 24, 2022

EdCaunt May 31, 2022

mloubout commented Jun 8, 2022

fffarias commented Jun 9, 2022

examples: Boundary reconstruction #1926

Are you sure you want to change the base?

examples: Boundary reconstruction #1926

Conversation

fffarias commented May 23, 2022

review-notebook-app bot commented May 23, 2022

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

mloubout Jun 28, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FabioLuporini commented May 24, 2022

mloubout 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

mloubout commented Jun 8, 2022

fffarias commented Jun 9, 2022

mloubout Jun 28, 2023 •

edited