Skip to content
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: free surface for coupled tti eqs #2173

Merged
merged 9 commits into from
Aug 9, 2023
3 changes: 2 additions & 1 deletion examples/seismic/preset_models.py
Original file line number Diff line number Diff line change
Expand Up @@ -216,7 +216,8 @@ def demo_model(preset, **kwargs):

model = SeismicModel(space_order=space_order, vp=v, origin=origin, shape=shape,
dtype=dtype, spacing=spacing, nbl=nbl, epsilon=epsilon,
delta=delta, theta=theta, phi=phi, bcs="damp", **kwargs)
delta=delta, theta=theta, phi=phi, bcs="damp",
fs=fs, **kwargs)

if kwargs.get('smooth', False):
if len(shape) > 2 and preset.lower() not in ['layers-tti-noazimuth']:
Expand Down
11 changes: 9 additions & 2 deletions examples/seismic/tti/operators.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
from devito import (Eq, Operator, Function, TimeFunction, NODE, Inc, solve,
cos, sin, sqrt)
from examples.seismic import PointSource, Receiver
from examples.seismic.acoustic.operators import freesurface


def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True):
Expand All @@ -20,10 +21,16 @@ def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True):
stencilp = solve(m * u.dt2 - H0 - qu + damp * udt, unext)
stencilr = solve(m * v.dt2 - Hz - qv + damp * vdt, vnext)

first_stencil = Eq(unext, stencilp)
second_stencil = Eq(vnext, stencilr)
first_stencil = Eq(unext, stencilp, subdomain=model.grid.subdomains['physdomain'])
second_stencil = Eq(vnext, stencilr, subdomain=model.grid.subdomains['physdomain'])

stencils = [first_stencil, second_stencil]

# Add free surface
if model.fs:
stencils.append(freesurface(model, Eq(unext, stencilp)))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No need to recreate the eq, just freesurface(first_stencil)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

firs_stencil is defined with the "physdomain" subdomain, for that reason I used Eq again

stencils.append(freesurface(model, Eq(vnext, stencilr)))

return stencils


Expand Down
3 changes: 3 additions & 0 deletions examples/seismic/tti/wavesolver.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,9 @@ def __init__(self, model, geometry, space_order=4, kernel='centered',
self.geometry = geometry
self.kernel = kernel

if model.fs and kernel == 'staggered':
raise ValueError("Free surface only supported for centered TTI kernel")

if space_order % 2 != 0:
raise ValueError("space_order must be even but got %s"
% space_order)
Expand Down
10 changes: 10 additions & 0 deletions tests/test_adjoint.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,9 @@
presets = {
'constant': {'preset': 'constant-isotropic'},
'layers': {'preset': 'layers-isotropic', 'nlayers': 2},
'layers-fs': {'preset': 'layers-isotropic', 'nlayers': 2, 'fs': True},
'layers-tti': {'preset': 'layers-tti', 'nlayers': 2},
'layers-tti-fs': {'preset': 'layers-tti', 'nlayers': 2, 'fs': True},
'layers-viscoacoustic': {'preset': 'layers-viscoacoustic', 'nlayers': 2},
}

Expand All @@ -27,6 +29,8 @@ class TestAdjoint(object):
('layers', (60, 70), 'OT2', 8, 2, acoustic_setup),
('layers', (60, 70), 'OT2', 4, 2, acoustic_setup),
('layers', (60, 70), 'OT4', 2, 2, acoustic_setup),
# 2D test with 2 layers and freesurface
('layers-fs', (60, 70), 'OT2', 4, 2, acoustic_setup),
# 3D tests with varying time and space orders
('layers', (60, 70, 80), 'OT2', 8, 2, acoustic_setup),
('layers', (60, 70, 80), 'OT2', 6, 2, acoustic_setup),
Expand All @@ -42,6 +46,8 @@ class TestAdjoint(object):
('layers-tti', (30, 35), 'centered', 4, 2, tti_setup),
('layers-tti', (30, 35), 'staggered', 8, 1, tti_setup),
('layers-tti', (30, 35), 'staggered', 4, 1, tti_setup),
# 2D TTI test with 2 layers and freesurface
('layers-tti-fs', (30, 35), 'centered', 4, 2, tti_setup),
# 3D TTI tests with varying space orders
('layers-tti', (30, 35, 40), 'centered', 8, 2, tti_setup),
('layers-tti', (30, 35, 40), 'centered', 4, 2, tti_setup),
Expand Down Expand Up @@ -123,13 +129,17 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun
('layers', (60, 70), 'OT2', 12, 2, acoustic_setup),
('layers', (60, 70), 'OT2', 8, 2, acoustic_setup),
('layers', (60, 70), 'OT2', 4, 2, acoustic_setup),
# 2D test with 2 layers and freesurface
('layers-fs', (60, 70), 'OT2', 4, 2, acoustic_setup),
# 3D tests with varying time and space orders
('layers', (40, 50, 30), 'OT2', 12, 2, acoustic_setup),
('layers', (40, 50, 30), 'OT2', 8, 2, acoustic_setup),
('layers', (40, 50, 30), 'OT2', 4, 2, acoustic_setup),
# 2D TTI tests with varying space orders
('layers-tti', (20, 25), 'centered', 8, 2, tti_setup),
('layers-tti', (20, 25), 'centered', 4, 2, tti_setup),
# 2D TTI test with 2 layers and freesurface
('layers-tti-fs', (20, 25), 'centered', 4, 2, tti_setup),
# 3D TTI tests with varying space orders
('layers-tti', (20, 25, 30), 'centered', 8, 2, tti_setup),
('layers-tti', (20, 25, 30), 'centered', 4, 2, tti_setup),
Expand Down