-
Notifications
You must be signed in to change notification settings - Fork 228
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
examples: Boundary reconstruction #1926
base: master
Are you sure you want to change the base?
Conversation
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
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]), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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(....)]
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Indeed good point I missed that gotta think about it
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@mloubout prodding
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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 standardsum += inject(expr)
- assign sum afterward
field = sum
This would cover the off-the-grid case as well
((10, 10, 10), 1., False), | ||
((10, 10, 10), 5., True) | ||
]) | ||
def test_inject_time_increment(shape, result, increment): |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for attaching tests to this PR!
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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?
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 |
ad7916e
to
4952720
Compare
…ng interpolation weights.
…e increment option in the SparseTimeFunction inject method
4952720
to
67d3f34
Compare
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
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]), |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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
@@ -0,0 +1,1069 @@ | |||
{ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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.
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() |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
can probably compute meshgrid only once before the for loop
@@ -0,0 +1,1069 @@ | |||
{ |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Generic question. Would it be easier to use 1D |
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. |
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