From bfdb39133d861e68bd345c8db5bee09f16e16a5e Mon Sep 17 00:00:00 2001 From: ofmla Date: Wed, 26 Jul 2023 14:38:50 -0300 Subject: [PATCH 1/9] examples: add fs arg for tti layered model --- examples/seismic/preset_models.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/seismic/preset_models.py b/examples/seismic/preset_models.py index 1b8ac3ea7d..0b79b102bb 100644 --- a/examples/seismic/preset_models.py +++ b/examples/seismic/preset_models.py @@ -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']: From f7dfe7e443ee291c9b469d1bd962d5d374bf13cf Mon Sep 17 00:00:00 2001 From: ofmla Date: Wed, 26 Jul 2023 14:40:45 -0300 Subject: [PATCH 2/9] examples: free surface boundary condition following the zero stress formulation for tti case --- examples/seismic/tti/operators.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py index 642593b38f..b30cf214c9 100644 --- a/examples/seismic/tti/operators.py +++ b/examples/seismic/tti/operators.py @@ -1,6 +1,8 @@ from devito import (Eq, Operator, Function, TimeFunction, NODE, Inc, solve, - cos, sin, sqrt) + cos, sin, sqrt, sign) from examples.seismic import PointSource, Receiver +from examples.seismic.acoustic.operators import freesurface +from devito.symbolics import retrieve_functions, INT def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True): @@ -20,10 +22,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))) + stencils.append(freesurface(model, Eq(vnext, stencilr))) + return stencils From 21990c7a97b037292a4d45113846b87825439aeb Mon Sep 17 00:00:00 2001 From: ofmla Date: Wed, 26 Jul 2023 14:43:41 -0300 Subject: [PATCH 3/9] examples: fs only implemented for centered kernel --- examples/seismic/tti/wavesolver.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py index 5aef427401..9b7b56451e 100644 --- a/examples/seismic/tti/wavesolver.py +++ b/examples/seismic/tti/wavesolver.py @@ -33,6 +33,12 @@ def __init__(self, model, geometry, space_order=4, kernel='centered', self.model._initialize_bcs(bcs="damp") self.geometry = geometry self.kernel = kernel + + if model.fs and kernel == 'staggered': + warning("Bad arguments in AnisotropicWaveSolver() - " + + "if model.fs is True, kernel should be centered. " + + "Switching to 'centered'") + self.kernel = 'centered' if space_order % 2 != 0: raise ValueError("space_order must be even but got %s" From a3808b51d7af1e52cf920ba4bc9f34f5da0c38df Mon Sep 17 00:00:00 2001 From: ofmla Date: Wed, 26 Jul 2023 15:10:59 -0300 Subject: [PATCH 4/9] examples: add some adjoint test for fs --- tests/test_adjoint.py | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/tests/test_adjoint.py b/tests/test_adjoint.py index 4f0e558697..9427f70e32 100644 --- a/tests/test_adjoint.py +++ b/tests/test_adjoint.py @@ -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}, } @@ -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', (200, 80), '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), @@ -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', (200, 80), '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), @@ -90,7 +96,7 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun location from data. This test uses the conventional dot test: < Fx, y> = """ - tn = 500. # Final time + tn = 1000. if mkey[-3:] == "-fs" else 500. # Final time (twice the time for fs) # Create solver from preset solver = setup_func(shape=shape, spacing=[15. for _ in shape], @@ -123,6 +129,8 @@ 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', (200, 80), '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), @@ -130,6 +138,8 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun # 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', (200, 80), '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), @@ -155,7 +165,7 @@ def test_adjoint_J(self, mkey, shape, kernel, space_order, time_order, setup_fun dot test: < Jx, y> = """ - tn = 500. # Final time + tn = 1000. if mkey[-3:] == "-fs" else 500. # Final time (twice the time for fs) nbl = 10 + space_order / 2 spacing = tuple([10.]*len(shape)) # Create solver from preset From 17397de550bb41160f802d6636aef76884388472 Mon Sep 17 00:00:00 2001 From: ofmla Date: Wed, 26 Jul 2023 15:22:49 -0300 Subject: [PATCH 5/9] examples: fix flake8 --- examples/seismic/tti/operators.py | 11 +++++------ examples/seismic/tti/wavesolver.py | 2 +- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/examples/seismic/tti/operators.py b/examples/seismic/tti/operators.py index b30cf214c9..b6e8687c04 100644 --- a/examples/seismic/tti/operators.py +++ b/examples/seismic/tti/operators.py @@ -1,8 +1,7 @@ from devito import (Eq, Operator, Function, TimeFunction, NODE, Inc, solve, - cos, sin, sqrt, sign) + cos, sin, sqrt) from examples.seismic import PointSource, Receiver from examples.seismic.acoustic.operators import freesurface -from devito.symbolics import retrieve_functions, INT def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True): @@ -26,12 +25,12 @@ def second_order_stencil(model, u, v, H0, Hz, qu, qv, forward=True): 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))) - stencils.append(freesurface(model, Eq(vnext, stencilr))) - + stencils.append(freesurface(model, Eq(unext, stencilp))) + stencils.append(freesurface(model, Eq(vnext, stencilr))) + return stencils diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py index 9b7b56451e..aa51665ac9 100644 --- a/examples/seismic/tti/wavesolver.py +++ b/examples/seismic/tti/wavesolver.py @@ -33,7 +33,7 @@ def __init__(self, model, geometry, space_order=4, kernel='centered', self.model._initialize_bcs(bcs="damp") self.geometry = geometry self.kernel = kernel - + if model.fs and kernel == 'staggered': warning("Bad arguments in AnisotropicWaveSolver() - " + "if model.fs is True, kernel should be centered. " + From 38b5048c6cecfae8264567b4d6d330e59ec6947e Mon Sep 17 00:00:00 2001 From: ofmla Date: Thu, 27 Jul 2023 17:10:02 -0300 Subject: [PATCH 6/9] examples: throw error instead of warning --- examples/seismic/tti/wavesolver.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py index aa51665ac9..e4aa30051c 100644 --- a/examples/seismic/tti/wavesolver.py +++ b/examples/seismic/tti/wavesolver.py @@ -35,10 +35,8 @@ def __init__(self, model, geometry, space_order=4, kernel='centered', self.kernel = kernel if model.fs and kernel == 'staggered': - warning("Bad arguments in AnisotropicWaveSolver() - " + - "if model.fs is True, kernel should be centered. " + - "Switching to 'centered'") - self.kernel = 'centered' + raise ValueError("Bad arguments in AnisotropicWaveSolver() - " + + "if model.fs is True, kernel must be centered.") if space_order % 2 != 0: raise ValueError("space_order must be even but got %s" From 5e1504abac9ef5ba2600e80abe3a3635a165a48a Mon Sep 17 00:00:00 2001 From: ofmla Date: Thu, 27 Jul 2023 17:13:02 -0300 Subject: [PATCH 7/9] examples: make sure the time is long enough to record some data --- tests/test_adjoint.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/test_adjoint.py b/tests/test_adjoint.py index 9427f70e32..3ec6d7ee32 100644 --- a/tests/test_adjoint.py +++ b/tests/test_adjoint.py @@ -96,7 +96,7 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun location from data. This test uses the conventional dot test: < Fx, y> = """ - tn = 1000. if mkey[-3:] == "-fs" else 500. # Final time (twice the time for fs) + tn = 500. # Final time # Create solver from preset solver = setup_func(shape=shape, spacing=[15. for _ in shape], @@ -165,7 +165,7 @@ def test_adjoint_J(self, mkey, shape, kernel, space_order, time_order, setup_fun dot test: < Jx, y> = """ - tn = 1000. if mkey[-3:] == "-fs" else 500. # Final time (twice the time for fs) + tn = 750. if mkey[-3:] == "-fs" else 500 # Final time nbl = 10 + space_order / 2 spacing = tuple([10.]*len(shape)) # Create solver from preset From e4b97c4e5487d90457e8fa2b8710fe7ce62ded7b Mon Sep 17 00:00:00 2001 From: ofmla Date: Fri, 4 Aug 2023 15:21:51 -0300 Subject: [PATCH 8/9] examples: be more specific in the message of the exception constructor --- examples/seismic/tti/wavesolver.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/seismic/tti/wavesolver.py b/examples/seismic/tti/wavesolver.py index e4aa30051c..982731f031 100644 --- a/examples/seismic/tti/wavesolver.py +++ b/examples/seismic/tti/wavesolver.py @@ -35,8 +35,7 @@ def __init__(self, model, geometry, space_order=4, kernel='centered', self.kernel = kernel if model.fs and kernel == 'staggered': - raise ValueError("Bad arguments in AnisotropicWaveSolver() - " + - "if model.fs is True, kernel must be centered.") + 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" From 4db0c98969c2f2194b5216f84e9833cbcc600f60 Mon Sep 17 00:00:00 2001 From: ofmla Date: Fri, 4 Aug 2023 15:25:27 -0300 Subject: [PATCH 9/9] examples: move to original confgs of tests --- tests/test_adjoint.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/test_adjoint.py b/tests/test_adjoint.py index 3ec6d7ee32..473e484c0e 100644 --- a/tests/test_adjoint.py +++ b/tests/test_adjoint.py @@ -30,7 +30,7 @@ class TestAdjoint(object): ('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', (200, 80), 'OT2', 4, 2, acoustic_setup), + ('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), @@ -47,7 +47,7 @@ class TestAdjoint(object): ('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', (200, 80), 'centered', 4, 2, tti_setup), + ('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), @@ -130,7 +130,7 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun ('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', (200, 80), 'OT2', 4, 2, acoustic_setup), + ('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), @@ -139,7 +139,7 @@ def test_adjoint_F(self, mkey, shape, kernel, space_order, time_order, setup_fun ('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', (200, 80), 'centered', 4, 2, tti_setup), + ('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), @@ -165,7 +165,7 @@ def test_adjoint_J(self, mkey, shape, kernel, space_order, time_order, setup_fun dot test: < Jx, y> = """ - tn = 750. if mkey[-3:] == "-fs" else 500 # Final time + tn = 500. # Final time nbl = 10 + space_order / 2 spacing = tuple([10.]*len(shape)) # Create solver from preset