Freesurface with subdomain (iso-acoustic) #1344

mloubout · 2020-06-09T13:59:36Z

POC for freesurface with a subdomain. Only implemented for isotropic acoustic as elastic is another story. Can be added to TTI as well.

FYI generated code with sign wrap around (((i2z - 2) > 0) - ((i2z - 2) < 0)):

          #pragma omp simd aligned(damp,u,vp:32)
          for (int i1z = i1z_ltkn + z_m; i1z <= -i1z_rtkn + z_M; i1z += 1)
          {
            float r13 = 1.0/(vp[x + 4][y + 4][i1z + 4]*vp[x + 4][y + 4][i1z + 4]);
            u[t1][x + 4][y + 4][i1z + 4] = (r13*(2*u[t0][x + 4][y + 4][i1z + 4] - u[t2][x + 4][y + 4][i1z + 4]) + dt*damp[x + 1][y + 1][i1z + 1]*u[t0][x + 4][y + 4][i1z + 4] + (dt*dt)*(3.70370379e-4F*(-u[t0][x + 2][y + 4][i1z + 4] - u[t0][x + 4][y + 2][i1z + 4] - u[t0][x + 4][y + 4][i1z + 2] - u[t0][x + 4][y + 4][i1z + 6] - u[t0][x + 4][y + 6][i1z + 4] - u[t0][x + 6][y + 4][i1z + 4]) + 5.92592607e-3F*(u[t0][x + 3][y + 4][i1z + 4] + u[t0][x + 4][y + 3][i1z + 4] + u[t0][x + 4][y + 4][i1z + 3] + u[t0][x + 4][y + 4][i1z + 5] + u[t0][x + 4][y + 5][i1z + 4] + u[t0][x + 5][y + 4][i1z + 4]) - 3.33333341e-2F*u[t0][x + 4][y + 4][i1z + 4]))/(r13 + dt*damp[x + 1][y + 1][i1z + 1]);
          }
          #pragma omp simd aligned(damp,u,vp:32)
          for (int i2z = z_m; i2z <= i2z_ltkn + z_m - 1; i2z += 1)
          {
            u[t1][x + 4][y + 4][i2z + 4] = (dt*damp[x + 1][y + 1][i2z + 1]*u[t0][x + 4][y + 4][i2z + 4] + (dt*dt)*(-3.70370379e-4F*(((i2z - 2) > 0) - ((i2z - 2) < 0))*u[t0][x + 4][y + 4][(int)(fabs(i2z - 2)) + 4] + 5.92592607e-3F*(((i2z - 1) > 0) - ((i2z - 1) < 0))*u[t0][x + 4][y + 4][(int)(fabs(i2z - 1)) + 4] - 3.70370379e-4F*u[t0][x + 2][y + 4][i2z + 4] + 5.92592607e-3F*u[t0][x + 3][y + 4][i2z + 4] - 3.70370379e-4F*u[t0][x + 4][y + 2][i2z + 4] + 5.92592607e-3F*u[t0][x + 4][y + 3][i2z + 4] - 3.33333341e-2F*u[t0][x + 4][y + 4][i2z + 4] + 5.92592607e-3F*u[t0][x + 4][y + 4][i2z + 5] - 3.70370379e-4F*u[t0][x + 4][y + 4][i2z + 6] + 5.92592607e-3F*u[t0][x + 4][y + 5][i2z + 4] - 3.70370379e-4F*u[t0][x + 4][y + 6][i2z + 4] + 5.92592607e-3F*u[t0][x + 5][y + 4][i2z + 4] - 3.70370379e-4F*u[t0][x + 6][y + 4][i2z + 4]) + (2*u[t0][x + 4][y + 4][i2z + 4] - u[t2][x + 4][y + 4][i2z + 4])/((vp[x + 4][y + 4][i2z + 4]*vp[x + 4][y + 4][i2z + 4])))/(dt*damp[x + 1][y + 1][i2z + 1] + 1.0/(vp[x + 4][y + 4][i2z + 4]*vp[x + 4][y + 4][i2z + 4]));
          }

codecov · 2020-06-09T14:02:54Z

Codecov Report

Merging #1344 into master will decrease coverage by 0.04%.
The diff coverage is 54.34%.

@@            Coverage Diff             @@
##           master    #1344      +/-   ##
==========================================
- Coverage   86.69%   86.65%   -0.05%     
==========================================
  Files         183      183              
  Lines       26254    26300      +46     
  Branches     3605     3613       +8     
==========================================
+ Hits        22762    22790      +28     
- Misses       3067     3083      +16     
- Partials      425      427       +2

Impacted Files	Coverage Δ
examples/seismic/acoustic/acoustic_example.py	`38.09% <0.00%> (-1.91%)`	⬇️
examples/seismic/elastic/elastic_example.py	`71.42% <0.00%> (ø)`
examples/seismic/skew_self_adjoint/example_iso.py	`47.72% <0.00%> (ø)`
examples/seismic/tti/tti_example.py	`45.45% <0.00%> (-4.55%)`	⬇️
examples/seismic/utils.py	`80.00% <0.00%> (ø)`
...les/seismic/viscoacoustic/viscoacoustic_example.py	`70.37% <0.00%> (ø)`
...mples/seismic/viscoelastic/viscoelastic_example.py	`71.42% <0.00%> (ø)`
examples/seismic/acoustic/operators.py	`79.74% <60.60%> (-15.34%)`	⬇️
examples/seismic/model.py	`84.91% <72.50%> (-4.83%)`	⬇️
examples/seismic/preset_models.py	`57.55% <100.00%> (+0.24%)`	⬆️
... and 3 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 85ea263...20fdea1. Read the comment docs.

jkwashbourne · 2020-06-09T18:00:40Z

should we make a unit test with a mirror condition?

this is friggin great! once we get to actually using devito the first thing we will need is this.

examples/seismic/acoustic/operators.py

examples/seismic/utils.py

FabioLuporini · 2020-06-10T12:46:46Z

examples/seismic/acoustic/operators.py

+    z = model.grid.dimensions[-1]
+    zfs = model.grid.subdomains['fsdomain'].dimensions[-1]
+
+    funcs = retrieve_functions(rhs)


this whole function is OKish for this PR, but I feel the user will need something easier than this

Yeah I know, I have a TODO in it, this just for "it works" but could be done better at some point.

FabioLuporini · 2020-06-10T12:49:05Z

examples/seismic/model.py

-        val = dampcoeff * (pos - sin(2*np.pi*pos)/(2*np.pi))
-        val = -val if abc_type == "mask" else val
-        eqs += [Inc(damp.subs({d: dim_l}), val/d.spacing)]
+        if not fs or d is not damp.dimensions[-1]:


I don't like this. This is a symptom that this code needs some refactoring to distinguish the fs vs no-fs cases

Not sure I follow. This pads damp where it's define, which is every side, except the top for the free surface case. That's all that's needed and refactoring would probably not make it easier to follow IMO.

ok, I think I understand now...it's a bit contrieved the way it is ATM, IMHO

so what you mean is :

pad_dims = list(damp.dimensions) if fs: pad_dims.pop(-1) for d in pad_dims: <padding code, just as in master>

?

But you still have to pad half of the last dim, this would remove ABC at the bottom.

just fission the loop into two loops ? it would be much much clearer I feel

examples/seismic/acoustic/acoustic_example.py

examples/seismic/utils.py

examples/seismic/model.py

examples/seismic/acoustic/acoustic_example.py

examples/seismic/acoustic/operators.py

FabioLuporini · 2020-06-10T15:08:08Z

examples/seismic/model.py

-        val = dampcoeff * (pos - sin(2*np.pi*pos)/(2*np.pi))
-        val = -val if abc_type == "mask" else val
-        eqs += [Inc(damp.subs({d: dim_l}), val/d.spacing)]
+        if not fs or d is not damp.dimensions[-1]:


ok, I think I understand now...it's a bit contrieved the way it is ATM, IMHO

so what you mean is :

pad_dims = list(damp.dimensions) if fs: pad_dims.pop(-1) for d in pad_dims: <padding code, just as in master>

?

examples/seismic/model.py

FabioLuporini · 2020-06-12T07:26:55Z

examples/seismic/tti/tti_example.py

+    parser = seismic_args(description)
+    parser.add_argument('--noazimuth', dest='azi', default=False, action='store_true',
+                        help="Whether or not to use an azimuth angle")
+    args = parser.parse_args()


blank line below

FabioLuporini · 2020-06-12T07:29:26Z

examples/seismic/viscoacoustic/viscoacoustic_example.py


    # Preset parameters
    ndim = args.ndim
    shape = args.shape[:args.ndim]
    spacing = tuple(ndim * [10.0])
-    tn = 750. if ndim < 3 else 1250.
+    tn = args.tn if args.tn > 0 else (750. if ndim < 3 else 1250.)


I'm writing this here, but it applies to all our examples. It feels like that much of the preprocessing could be moved inside seismic_args once and for all, instead of being repeated each time in the individual examples?

not necessarily for this PR (would be out of scope admittedly), but you know, just saying..

I looked into it, and it's not possible if we return the parser from seismic_args. We either return the args and put all the possible ones in seismic_args ir we return the parser and this can only be done afer parse_args() which is done here if we want these different cases.

I imagined... perhaps that's one reason to switch to click -- in which case, we could provide callbacks to each argument! as I said, just food for thought...

rhodrin · 2020-06-12T16:23:57Z

Just starting to look through this now, but first comment is 'no tutorial notebook? 😲' 😛

examples/seismic/acoustic/operators.py

rhodrin · 2020-06-12T17:14:08Z

examples/seismic/acoustic/operators.py

+
+    # Add free surface
+    if model.fs:
+        eqns.append(freesurface(model, Eq(unext, eq_time)))


Passing Eq(unext, eq_time) and then splitting it up again in freesurface?

Yes, this apply the free surface to the equation directly

examples/seismic/model.py

rhodrin · 2020-06-12T17:38:01Z

examples/seismic/model.py

        self.shape = shape
        self.space_order = space_order
        self.nbl = int(nbl)
        self.origin = tuple([dtype(o) for o in origin])
+        self.fs = fs
+        if self.fs:
+            warning("Freesurface is only supported for isotropic acoustic solver")


Should we not:

Not have a warning if using iso-acoustic

Throw a not implemented error if not using isco-acoustic
?

I can try, the model doesn't really know which physics it is at this point so would need to move it to each example. I can just remove it as would just be False for all other examples

rhodrin · 2020-06-12T17:48:12Z

examples/seismic/acoustic/operators.py

+    mapper = {}
+    # Mirror negative indices
+    # TODO: Make a proper "mirror_indices" tool function
+    for f in funcs:


~~I see the TODO - but if I'm reading this correctly it's messing with values in the halo to implement this atm? Is that something we really want to do? (Please correct me if I'm misunderstanding).~~

Oh, or is the sponge still there and values are being manipulated in that?

This is not touching the halo at all no. This is just turning any indices smaller than symbolic min into its mirror i.e u[x-1] => u[ abs(x-1)]. The halo is actually unnecessary then since nothing is ever accessed in it.

Oh, ok, I was initially thrown by the if ... < 0. So you're creating a 'boundary layer' of length so to mirror the values in?

Yes, the top so gridpoints are the layer in which indices are mirrored if negative. Ultimately (after some work) this would be basically autowrapping at the edge of any domain which would remove the need of a halo.

And the boundary layer could also be completely removed by modifying the fd stencil appropriately -- easily doable for iso-acoustic, but probably too much pain to do nicely once this is generalized.

And the boundary layer could also be completely removed by modifying the fd stencil appropriately

That'd be really ideal yes but not really easy on the symbolic side as now you would have a different stencil (coeffs and radius) per distance to the boundary so wouldn't really be doable with a subdomain that expects the same equation for all points in the domain. Even for iso-acoustic, it's fairly trivial to do it by hand but not as a generic symbolic rule (unless I missed something).

Even for iso-acoustic, it's fairly trivial to do it by hand but not as a generic symbolic rule (unless I missed something).

For flat boundaries we can make use of custom FD coefficients to do it pretty easily.

rhodrin · 2020-06-12T18:24:10Z

Another comment is why is this being restricted to the 'left' of the last dimension? I suppose the argument would be that the Model class if just a set of 'helper' wrappers anyway and if you want something different you're able to implement it yourself?

rhodrin · 2020-06-12T18:45:33Z

Also, shouldn't we have a norm test for a run with an FS?

mloubout · 2020-06-15T13:36:56Z

Also, shouldn't we have a norm test for a run with an FS?

I can add one to #1348 that adds lots of tests to the examples once this one is in.

mloubout requested review from FabioLuporini, rhodrin and jkwashbourne June 9, 2020 13:59

mloubout added enhancement examples examples labels Jun 9, 2020

mloubout force-pushed the fsdomain branch from 37fd587 to 1f78591 Compare June 9, 2020 14:33

georgebisbas reviewed Jun 10, 2020

View reviewed changes

FabioLuporini reviewed Jun 10, 2020

View reviewed changes

mloubout force-pushed the fsdomain branch 5 times, most recently from db669b1 to 2bf7992 Compare June 11, 2020 14:31

georgebisbas reviewed Jun 12, 2020

View reviewed changes

examples/seismic/model.py Outdated Show resolved Hide resolved

FabioLuporini approved these changes Jun 12, 2020

View reviewed changes

mloubout added the POC label Jun 12, 2020

mloubout force-pushed the fsdomain branch from 2bf7992 to 5f0a1ab Compare June 12, 2020 12:25

rhodrin reviewed Jun 12, 2020

View reviewed changes

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

rhodrin reviewed Jun 12, 2020

View reviewed changes

examples/seismic/model.py Outdated Show resolved Hide resolved

rhodrin reviewed Jun 12, 2020

View reviewed changes

mloubout force-pushed the fsdomain branch from 5f0a1ab to 00b2921 Compare June 15, 2020 13:40

rhodrin approved these changes Jun 15, 2020

View reviewed changes

mloubout added 2 commits June 15, 2020 13:09

isoacoustic: add freesurface with subdomain

7190bdb

reviewer comments first pass

20fdea1

mloubout force-pushed the fsdomain branch from 00b2921 to 20fdea1 Compare June 15, 2020 17:09

mloubout merged commit 377a2fb into master Jun 15, 2020

FabioLuporini deleted the fsdomain branch June 15, 2020 20:30

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Freesurface with subdomain (iso-acoustic) #1344

Freesurface with subdomain (iso-acoustic) #1344

mloubout commented Jun 9, 2020 •

edited

Loading

codecov bot commented Jun 9, 2020 •

edited

Loading

jkwashbourne commented Jun 9, 2020

FabioLuporini Jun 10, 2020

mloubout Jun 11, 2020

FabioLuporini Jun 10, 2020

mloubout Jun 10, 2020

FabioLuporini Jun 10, 2020

mloubout Jun 10, 2020

FabioLuporini Jun 12, 2020

FabioLuporini Jun 10, 2020

FabioLuporini Jun 12, 2020

FabioLuporini Jun 12, 2020

FabioLuporini Jun 12, 2020

mloubout Jun 12, 2020

FabioLuporini Jun 12, 2020

rhodrin commented Jun 12, 2020 •

edited

Loading

rhodrin Jun 12, 2020 •

edited

Loading

mloubout Jun 15, 2020

rhodrin Jun 12, 2020

mloubout Jun 15, 2020

rhodrin Jun 12, 2020 •

edited

Loading

rhodrin Jun 12, 2020

mloubout Jun 12, 2020

rhodrin Jun 12, 2020

mloubout Jun 12, 2020

rhodrin Jun 12, 2020

mloubout Jun 15, 2020 •

edited

Loading

rhodrin Jun 15, 2020

rhodrin commented Jun 12, 2020

rhodrin commented Jun 12, 2020

mloubout commented Jun 15, 2020

Freesurface with subdomain (iso-acoustic) #1344

Freesurface with subdomain (iso-acoustic) #1344

Conversation

mloubout commented Jun 9, 2020 • edited Loading

codecov bot commented Jun 9, 2020 • edited Loading

Codecov Report

jkwashbourne commented Jun 9, 2020

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

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

rhodrin commented Jun 12, 2020 • edited Loading

rhodrin Jun 12, 2020 • 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

rhodrin Jun 12, 2020 • 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

mloubout Jun 15, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rhodrin commented Jun 12, 2020

rhodrin commented Jun 12, 2020

mloubout commented Jun 15, 2020

mloubout commented Jun 9, 2020 •

edited

Loading

codecov bot commented Jun 9, 2020 •

edited

Loading

rhodrin commented Jun 12, 2020 •

edited

Loading

rhodrin Jun 12, 2020 •

edited

Loading

rhodrin Jun 12, 2020 •

edited

Loading

mloubout Jun 15, 2020 •

edited

Loading