Examples: Added metric for geometry-compensated energy concentration …

…to the infrasound example

examples/infrasound/2d_infrasound_example.py

@@ -10,11 +10,15 @@
 from propagator import InfrasoundPropagator
 
 
-src_coords = np.array([500., 500.])[np.newaxis, :]
-rec_coords = np.array([[250., 250.], [750., 750.]])
-t0, tn, dt = 0., 2.25, 0.007  # Courant number ~0.5
-src_f = 2.
-model = InfrasoundModel(dims=2, shape=(201, 201), extent=(1000., 1000.),
+src_coords = np.array([4800., 2550.])[np.newaxis, :]
+rec_coords = np.array([[2000., 1900.], [7505., 3025.],
+                       [7200., 1275.], [2150., 3300.],
+                       [3000., 1000.], [6505, 4500.]])
+# rec_coords = np.array([[2400., 1275.], [7200., 3825.],
+#                        [7200., 1275.], [2400., 3825.]])
+t0, tn, dt = 0., 13., 0.021  # Courant number ~0.25
+src_f = 1.
+model = InfrasoundModel(dims=2, shape=(321, 171), extent=(9600., 5100.),
                         src_coords=src_coords, rec_coords=rec_coords,
                         t0=t0, tn=tn, dt=dt, src_f=src_f)
 
@@ -23,6 +27,30 @@
 
 print(np.linalg.norm(model.p.data[-1]))
 
-plt.imshow(model.p.data[-1])
+plt.imshow(model.p.data[-1].T)
+plt.colorbar()
+plt.show()
+
+print(np.linalg.norm(model.rec.data[:, 0]))
+plt.plot(model.rec.data[:, 0])
+plt.show()
+
+# Reset the fields
+model.p.data[:] = 0
+model.p_aux[0].data[:] = 0
+model.p_aux[1].data[:] = 0
+model.A[0].data[:] = 0
+model.A[1].data[:] = 0
+
+bpropagator = InfrasoundPropagator(model=model, mode='adjoint',
+                                   track_zsc=True)
+bpropagator.run()
+
+print(np.linalg.norm(model.p.data[-1]))
+plt.imshow(model.p.data[-1].T)
+plt.colorbar()
+plt.show()
+
+plt.imshow(model.zsc.data.T)
 plt.colorbar()
 plt.show()

examples/infrasound/model.py

@@ -134,14 +134,19 @@ def _setup_fields(self):
         # Signed distance function for boundary
         self._sdf = dv.Function(name='sdf', grid=self.grid)
 
+        # Z score (measure of energy concentration, Kim and Lees 2014)
+        self._zsc = dv.TimeFunction(name='zsc', grid=self.grid,
+                                    time_order=0,
+                                    staggered=dv.NODE)
+
     def _setup_damping(self):
         """Initialise the damping mask"""
         d = self.damp
         grid = self.grid
         # Hardcoded for 10 layers of PMLs
         eqs = []
         default_names = ('interior', 'domain')
-        d_par = 10.  # Damping parameter
+        d_par = 2.  # Damping parameter
         for name, subdomain in self.grid.subdomains.items():
             for i in range(self._ndims):
                 if name[i] == 'l' and name not in default_names:
@@ -240,6 +245,11 @@ def sdf(self):
         """Signed distance function for boundary"""
         return self._sdf
 
+    @property
+    def zsc(self):
+        """The z-score as defined by Kim and Lees 2014"""
+        return self._zsc
+
     @property
     def t0(self):
         """Start time"""

examples/infrasound/propagator.py

@@ -9,6 +9,7 @@
 
 import devito as dv
 import sympy as sp
+import numpy as np
 
 from cached_property import cached_property
 
@@ -25,6 +26,9 @@ class InfrasoundPropagator:
     mode : str
         Mode of propagation. Can be 'forward' or 'adjoint'. Default is
         'forward'.
+    track_zsc : bool
+        Track the z-score. This is a measure of energy concentration devised
+        by Kim and Lees 2014
     """
     def __init__(self, *args, **kwargs):
         self._model = kwargs.get('model')
@@ -33,6 +37,8 @@ def __init__(self, *args, **kwargs):
         self._mode = kwargs.get('mode', 'forward')
         self._setup_dense()
         self._setup_sparse()
+        self._track_zsc = kwargs.get('track_zsc', False)
+        self._setup_heuristics()
 
     def _setup_dense(self):
         """Set up the update equations (for dense kernels)"""
@@ -105,12 +111,55 @@ def _setup_sparse(self):
         if self.mode == 'forward':
             # Timestep to update
             pf = p.forward
+            if self.model.src is not None:
+                self._sparse.append(self.model.src.inject(field=pf,
+                                                          expr=self.model.src))
+
+            if self.model.rec is not None:
+                self._sparse.append(self.model.rec.interpolate(expr=pf))
         else:
             pf = p.backward
+            if self.model.rec is not None:
+                self._sparse.append(self.model.rec.inject(field=pf,
+                                                          expr=self.model.rec))
+
+    def _setup_heuristics(self):
+        """Set up tracking terms"""
+        if self.mode == 'forward':
+            def fwd(f):
+                return f.forward
+        else:
+            def fwd(f):
+                return f.backward
+
+        p = self.model.p
+        zsc = self.model.zsc
+        c = self.model.c
+        dt = self.model.dt
+        self._heuristics = []
+
+        main_name = 'm'*self.model._ndims
 
-        if self.model.src is not None:
-            self._sparse.append(self.model.src.inject(field=pf,
-                                                      expr=self.model.src))
+        if self._track_zsc:
+            grid_size = np.prod(self.model.grid.shape)
+            # Mean of field
+            mu = dv.TimeFunction(name='mu', dimensions=(p.time_dim,),
+                                 time_order=2, shape=(self.model.src.nt,))
+            self._heuristics.append(dv.Inc(mu, fwd(p)/grid_size))
+            # Squared standard deviation
+            sigma2 = dv.TimeFunction(name='sigma2', dimensions=(p.time_dim,),
+                                     time_order=2, shape=(self.model.src.nt,))
+            self._heuristics.append(dv.Inc(sigma2,
+                                           (fwd(p)-mu)**2/grid_size))
+            sigma = dv.sqrt(sigma2)
+            # Scaling for geometric attenuation (approximate)
+            scale = dv.TimeFunction(name='scale', dimensions=(p.time_dim,),
+                                    time_order=2, shape=(self.model.src.nt,))
+            self._heuristics.append(dv.Eq(fwd(scale), scale + c*dt))
+            zsc_eq = dv.Eq(fwd(zsc),
+                           dv.Max(zsc, scale**2*(fwd(p)-mu)/sigma),
+                           subdomain=self.model.grid.subdomains[main_name])
+            self._heuristics.append(zsc_eq)
 
     def run(self):
         """Run the propagator"""
@@ -131,5 +180,6 @@ def mode(self):
     def operator(self):
         """The finite-difference operator"""
         return dv.Operator(self._A_eqs + self._p_aux_eqs
-                           + self._p_eqs + self._sparse,
+                           + self._p_eqs + self._sparse
+                           + self._heuristics,
                            name=self.mode)