Merge pull request #525 from opesci/reresample

Added support for resampling time series data in PointSource.

.travis.yml

@@ -109,8 +109,8 @@ script:
   - if [[ $DEVITO_BACKEND != 'yask' ]]; then python examples/checkpointing/checkpointing_example.py; fi
 
   # Test tutorial notebooks for the website using nbval
-  - if [[ $DEVITO_BACKEND != 'yask' ]]; then py.test -vs --nbval examples/seismic/tutorials; fi
   - if [[ $DEVITO_BACKEND != 'yask' ]]; then py.test -vs --nbval examples/cfd; fi
+  - if [[ $DEVITO_BACKEND != 'yask' ]]; then py.test -vs --nbval examples/seismic/tutorials; fi
 
   # Code coverage
   - codecov

examples/checkpointing/checkpointing_example.py

@@ -38,7 +38,7 @@ def gradient(self, m0, maxmem=None):
                                       v=self.adjoint_field, m=m0, rec=self.rec_g,
                                       grad=self.grad, dt=self.dt)
         # Subtracting 2 in the following line assuming time_order=2
-        wrp = Revolver(cp, wrap_fw, wrap_rev, n_checkpoints, self.nt-2)
+        wrp = Revolver(cp, wrap_fw, wrap_rev, n_checkpoints, self.src._time_range.num-2)
 
         wrp.apply_forward()
 

examples/seismic/acoustic/acoustic_example.py

@@ -4,7 +4,7 @@
 from devito.logger import info
 from devito import Constant
 from examples.seismic.acoustic import AcousticWaveSolver
-from examples.seismic import demo_model, RickerSource, Receiver
+from examples.seismic import demo_model, TimeAxis, RickerSource, Receiver
 
 
 # Velocity models
@@ -34,16 +34,15 @@ def acoustic_setup(shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0),
     # Derive timestepping from model spacing
     dt = model.critical_dt * (1.73 if kernel == 'OT4' else 1.0)
     t0 = 0.0
-    nt = int(1 + (tn-t0) / dt)  # Number of timesteps
-    time = np.linspace(t0, tn, nt)  # Discretized time axis
+    time_range = TimeAxis(start=t0, stop=tn, step=dt)
 
     # Define source geometry (center of domain, just below surface)
-    src = RickerSource(name='src', grid=model.grid, f0=0.01, time=time)
+    src = RickerSource(name='src', grid=model.grid, f0=0.01, time_range=time_range)
     src.coordinates.data[0, :] = np.array(model.domain_size) * .5
     src.coordinates.data[0, -1] = model.origin[-1] + 2 * spacing[-1]
 
     # Define receiver geometry (spread across x, just below surface)
-    rec = Receiver(name='rec', grid=model.grid, ntime=nt, npoint=nrec)
+    rec = Receiver(name='rec', grid=model.grid, time_range=time_range, npoint=nrec)
     rec.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
     rec.coordinates.data[:, 1:] = src.coordinates.data[0, 1:]
 

examples/seismic/acoustic/gradient_example.py

@@ -3,7 +3,7 @@
 from cached_property import cached_property
 
 from devito import TimeFunction, Function
-from examples.seismic import Receiver, RickerSource, demo_model
+from examples.seismic import TimeAxis, Receiver, RickerSource, demo_model
 from examples.seismic.acoustic import ForwardOperator, GradientOperator, smooth10
 
 
@@ -12,7 +12,8 @@ def __init__(self, shape=(50, 50, 50), spacing=(15.0, 15.0, 15.0), tn=500.,
                  kernel='OT2', space_order=4, nbpml=10):
         self.kernel = kernel
         self.space_order = space_order
-        self._setup_model_and_acquisition(space_order, shape, spacing, nbpml, tn)
+        self._setup_model_and_acquisition(space_order, shape, spacing,
+                                          nbpml+int(space_order/2), tn)
         self._true_data()
 
     @cached_property
@@ -26,31 +27,33 @@ def _setup_model_and_acquisition(self, space_order, shape, spacing, nbpml, tn):
                            shape=shape, spacing=spacing, nbpml=nbpml)
         self.model = model
         t0 = 0.0
-        self.nt = int(1 + (tn-t0) / self.dt)  # Number of timesteps
-        time = np.linspace(t0, tn, self.nt)  # Discretized time axis
+        time_range = TimeAxis(start=t0, stop=tn, step=self.dt)
+        self.nt = time_range.num
 
         # Define source geometry (center of domain, just below surface)
-        src = RickerSource(name='src', grid=model.grid, f0=0.01, time=time)
+        src = RickerSource(name='src', grid=model.grid, f0=0.01, time_range=time_range)
         src.coordinates.data[0, :] = np.array(model.domain_size) * .5
         src.coordinates.data[0, -1] = model.origin[-1] + 2 * spacing[-1]
 
         self.src = src
 
         # Define receiver geometry (spread across x, just below surface)
         # We need two receiver fields - one for the true (verification) run
-        rec_t = Receiver(name='rec', grid=model.grid, ntime=self.nt, npoint=nrec)
-        rec_t.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0], num=nrec)
+        rec_t = Receiver(name='rec', grid=model.grid, time_range=time_range,
+                         npoint=nrec)
+        rec_t.coordinates.data[:, 0] = np.linspace(0., model.domain_size[0],
+                                                   num=nrec)
         rec_t.coordinates.data[:, 1:] = src.coordinates.data[0, 1:]
 
         self.rec_t = rec_t
 
         # and the other for the smoothed run
-        self.rec = Receiver(name='rec', grid=model.grid, ntime=self.nt, npoint=nrec,
-                            coordinates=rec_t.coordinates.data)
+        self.rec = Receiver(name='rec', grid=model.grid, time_range=time_range,
+                            npoint=nrec, coordinates=rec_t.coordinates.data)
 
         # Receiver for Gradient
         self.rec_g = Receiver(name="rec", coordinates=self.rec.coordinates.data,
-                              grid=model.grid, dt=self.dt, ntime=self.nt)
+                              grid=model.grid, time_range=time_range)
 
         # Gradient symbol
         self.grad = Function(name="grad", grid=model.grid)
@@ -90,7 +93,7 @@ def temp_field(self):
 
     @cached_property
     def forward_field(self):
-        return TimeFunction(name="u", grid=self.model.grid, save=self.nt,
+        return TimeFunction(name="u", grid=self.model.grid, save=self.src._time_range.num,
                             time_order=2, space_order=self.space_order)
 
     @cached_property
@@ -123,10 +126,11 @@ def verify(self, m0, gradient, rec, dm):
         G = np.dot(gradient.reshape(-1), dm.reshape(-1))
         # FWI Gradient test
         H = [0.5, 0.25, .125, 0.0625, 0.0312, 0.015625, 0.0078125]
-        error1 = np.zeros(7)
-        error2 = np.zeros(7)
 
-        for i in range(0, 7):
+        error1 = np.zeros(len(H))
+        error2 = np.zeros(len(H))
+
+        for i in range(0, len(H)):
             # Add the perturbation to the model
             def initializer(data):
                 data[:] = m0.data + H[i] * dm
@@ -148,6 +152,7 @@ def initializer(data):
         # Test slope of the  tests
         p1 = np.polyfit(np.log10(H), np.log10(error1), 1)
         p2 = np.polyfit(np.log10(H), np.log10(error2), 1)
+
         assert np.isclose(p1[0], 1.0, rtol=0.1)
         assert np.isclose(p2[0], 2.0, rtol=0.1)
 

examples/seismic/acoustic/operators.py

@@ -66,9 +66,9 @@ def ForwardOperator(model, source, receiver, space_order=4,
     u = TimeFunction(name='u', grid=model.grid,
                      save=source.nt if save else None,
                      time_order=2, space_order=space_order)
-    src = PointSource(name='src', grid=model.grid, ntime=source.nt,
+    src = PointSource(name='src', grid=model.grid, time_range=source.time_range,
                       npoint=source.npoint)
-    rec = Receiver(name='rec', grid=model.grid, ntime=receiver.nt,
+    rec = Receiver(name='rec', grid=model.grid, time_range=receiver.time_range,
                    npoint=receiver.npoint)
 
     s = model.grid.stepping_dim.spacing
@@ -101,9 +101,9 @@ def AdjointOperator(model, source, receiver, space_order=4,
 
     v = TimeFunction(name='v', grid=model.grid, save=None,
                      time_order=2, space_order=space_order)
-    srca = PointSource(name='srca', grid=model.grid, ntime=source.nt,
+    srca = PointSource(name='srca', grid=model.grid, time_range=source.time_range,
                        npoint=source.npoint)
-    rec = Receiver(name='rec', grid=model.grid, ntime=receiver.nt,
+    rec = Receiver(name='rec', grid=model.grid, time_range=receiver.time_range,
                    npoint=receiver.npoint)
 
     s = model.grid.stepping_dim.spacing
@@ -140,8 +140,8 @@ def GradientOperator(model, source, receiver, space_order=4, save=True,
                      else None, time_order=2, space_order=space_order)
     v = TimeFunction(name='v', grid=model.grid, save=None,
                      time_order=2, space_order=space_order)
-    rec = Receiver(name='rec', grid=model.grid, ntime=receiver.nt,
-                   npoint=receiver.npoint)
+    rec = Receiver(name='rec', grid=model.grid,
+                   time_range=receiver.time_range, npoint=receiver.npoint)
 
     s = model.grid.stepping_dim.spacing
     eqn = iso_stencil(v, m, s, damp, kernel, forward=False)
@@ -176,9 +176,9 @@ def BornOperator(model, source, receiver, space_order=4,
     m, damp = model.m, model.damp
 
     # Create source and receiver symbols
-    src = PointSource(name='src', grid=model.grid, ntime=source.nt,
+    src = PointSource(name='src', grid=model.grid, time_range=source.time_range,
                       npoint=source.npoint)
-    rec = Receiver(name='rec', grid=model.grid, ntime=receiver.nt,
+    rec = Receiver(name='rec', grid=model.grid, time_range=receiver.time_range,
                    npoint=receiver.npoint)
 
     # Create wavefields and a dm field

examples/seismic/acoustic/wavesolver.py

@@ -88,7 +88,7 @@ def forward(self, src=None, rec=None, u=None, m=None, save=None, **kwargs):
         # Create a new receiver object to store the result
         if rec is None:
             rec = Receiver(name='rec', grid=self.model.grid,
-                           ntime=self.receiver.nt,
+                           time_range=self.receiver.time_range,
                            coordinates=self.receiver.coordinates.data)
 
         # Create the forward wavefield if not provided
@@ -123,7 +123,7 @@ def adjoint(self, rec, srca=None, v=None, m=None, **kwargs):
         # Create a new adjoint source and receiver symbol
         if srca is None:
             srca = PointSource(name='srca', grid=self.model.grid,
-                               ntime=self.source.nt,
+                               time_range=self.source.time_range,
                                coordinates=self.source.coordinates.data)
 
         # Create the adjoint wavefield if not provided
@@ -188,7 +188,7 @@ def born(self, dmin, src=None, rec=None, u=None, U=None, m=None, **kwargs):
         # Create a new receiver object to store the result
         if rec is None:
             rec = rec or Receiver(name='rec', grid=self.model.grid,
-                                  ntime=self.receiver.nt,
+                                  time_range=self.receiver.time_range,
                                   coordinates=self.receiver.coordinates.data)
 
         # Create the forward wavefields u and U if not provided