Fixed error on Reynolds ABS.

Working for both 2D and 3D.
Need to implement someday gaussian damping again.

cookbook/seismic_wavefd_scalar.py

@@ -21,19 +21,19 @@
 fc = 15.
 sources = [wavefd.GaussSource((125*ds, 75*ds), area, shape,  1., fc)]
 dt = wavefd.scalar_maxdt(area, shape, np.max(velocity))
-duration = 2.5
+duration = 1.9
 maxit = int(duration/dt)
 stations = [[75*ds, 125*ds]] # x, z coordinate of the seismometer
 snapshots = 3 # every 3 iterations plots one
 simulation = wavefd.scalar(velocity, area, dt, maxit, sources, stations, snapshots)
 
-# This part makes an animation using matplotlibs animation API
+# This part makes an animation using matplotlib animation API
 background = (velocity-4000)*10**-1
 fig = mpl.figure(figsize=(8, 6))
 mpl.subplots_adjust(right=0.98, left=0.11, hspace=0.5, top=0.93)
 mpl.subplot2grid((4, 3), (0,0), colspan=3,rowspan=3)
 wavefield = mpl.imshow(np.zeros_like(velocity), extent=area, cmap=mpl.cm.gray_r,
-                       vmin=-1000, vmax=1000)
+                       vmin=-300, vmax=300)
 mpl.points(stations, '^b', size=8)
 mpl.ylim(area[2:][::-1])
 mpl.xlabel('x (km)')

fatiando/seismic/_wavefd.pyx

@@ -14,13 +14,17 @@ ctypedef numpy.float_t double
 
 __all__ = [
     '_apply_damping',
+    '_apply_damping3',
     '_step_elastic_sh',
     '_step_elastic_psv',
     '_xz2ps',
     '_nonreflexive_sh_boundary_conditions',
     '_nonreflexive_psv_boundary_conditions',
-    '_reflexive_scalar_boundary_conditions',
-    '_step_scalar']
+    '_nonreflexive_scalar_boundary_conditions',
+    '_nonreflexive3_scalar_boundary_conditions',
+    '_step_scalar',
+    '_step_scalar3'
+    ]
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -89,6 +93,36 @@ def _apply_damping(double[:,::1] array not None,
         for j in range(nx):
             array[i,j] *= exp(-((decay*(i - nz + pad))**2))
 
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def _apply_damping3(double[:,:,::1] array not None,
+    unsigned int nx, unsigned int ny, unsigned int nz,
+    unsigned int pad, double decay):
+    """
+    Apply a decay factor to the values of the array in the padding region.
+    Damping 3D padding regions. x (north), y (east), z (down)
+    """
+    cdef:
+        unsigned int i, j, k
+    # Damping on the north/south
+    for k in range(nz-pad):
+        for j in range(pad,ny-pad):  # avoid double damping
+            for i in range(pad):
+                array[k, j, i] *= exp(-(decay*(pad - i))**2) # south
+                array[k, j, (nx-pad) + i] *= exp(-(decay*i)**2) # north
+    # Damping on the east/west
+    for k in range(nz-pad):
+        for i in range(nx):
+            for j in range(pad):
+                array[k, j, i] *= exp(-(decay*(pad - j))**2) # west
+                array[k, (ny-pad) + j, i] *= exp(-(decay*j)**2) # east
+    # Damping on the bottom
+    for i in range(nx):
+        for j in range(ny):
+            for k in range(pad):
+                array[(nz-pad) + k, i, j] *= exp(-(decay*k)**2)
+
+
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def _nonreflexive_psv_boundary_conditions(
@@ -271,26 +305,47 @@ def _step_elastic_psv(
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def _reflexive_scalar_boundary_conditions(
-    double[:,::1] u not None,
+def _nonreflexive_scalar_boundary_conditions(
+    double[:,::1] u_tp1 not None,
+    double[:,::1] u_t not None,
+    double[:,::1] u_tm1 not None,
+    double[:,::1] vel not None,
+    double dt, double ds,
     unsigned int nx, unsigned int nz):
     """
-    Apply the boundary conditions: free-surface at top, fixed on the others.
+    Apply the boundary conditions: free-surface at top, transparent in the borders
     4th order (+2-2) indexes
+
+    Uses Reynolds, A. C. - Boundary conditions for numerical solution of wave propagation problems
+    Geophysics p 1099-1110 - 1978
+    The finite difference approximation used by Reynolds for the transparent boundary condition is of first
+    order, though the scalar schema of propagation is of fourth order in space
+
     """
     cdef unsigned int i
-    # Top
+    # Top free surface: do I really need this? I might phase inversion
     for i in xrange(nx):
-        u[1, i] = u[2, i] #up
-        u[0, i] = u[1, i]
-        u[nz - 1, i] *= 0 #down
-        u[nz - 2, i] *= 0
-    # Sides
+        u_tp1[0, i] = 0.0 #up
+        u_tp1[1, i] = 0.0
+    # Transparent boundary condition applied exactly at
+    # The edge of the fourth order calculation. It's necessary to store t-1 pane for that
     for i in xrange(nz):
-        u[i, 0] *= 0 #left
-        u[i, 1] *= 0
-        u[i, nx - 1] *= 0 #right
-        u[i, nx - 2] *= 0
+        # left
+        for p in xrange(2):
+            u_tp1[i, p] = ( u_t[i, p] + u_t[i, p+1] - u_tm1[i,p+1] +
+                (vel[i, p]*dt/ds)*(u_t[i, p+1] - u_t[i, p] - u_tm1[i, p+2] + u_tm1[i, p+1])
+                )
+         #right
+        for p in xrange(2):
+            u_tp1[i, nx-2+p] = ( u_t[i, nx-2+p] + u_t[i, nx-3+p] - u_tm1[i, nx-3+p] -
+                (vel[i, nx-2+p]*dt/ds)*(u_t[i, nx-2+p] - u_t[i, nx-3+p] - u_tm1[i, nx-3+p] + u_tm1[i, nx-4+p])
+                )
+    # Down
+    for i in xrange(nx):
+        for p in xrange(2):
+            u_tp1[nz-2+p, i] = ( u_t[nz-2+p, i] + u_t[nz-3+p, i] - u_tm1[nz-3+p, i] -
+                    (vel[nz-2+p, i]*dt/ds)*(u_t[nz-2+p, i] - u_t[nz-3+p, i] - u_tm1[nz-3+p, i] + u_tm1[nz-4+p, i])
+                    )
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
@@ -314,3 +369,90 @@ def _step_scalar(
                      16.*u_t[i,j - 1] - u_t[i,j - 2])/12. +
                     (-u_t[i + 2,j] + 16.*u_t[i + 1,j] - 30.*u_t[i,j] +
                      16.*u_t[i - 1,j] - u_t[i - 2,j])/12.))
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def _nonreflexive3_scalar_boundary_conditions(
+    double[:,:,::1] u_tm1 not None,
+    double[:,:,::1] u_t not None,
+    double[:,:,::1] u_tp1 not None,
+    double dt, double dx, double dy, double dz,
+    double[:,:,::1] c not None,
+    unsigned int nx, unsigned int ny, unsigned int nz):
+    """
+    Apply the boundary conditions: fixed-surface at top, transparent in the borders
+    4th order (+2-2) indexes
+
+    Threat the field as it was just scalar and uses Reynolds ABS conditions.
+    Reynolds, A. C. - Boundary conditions for numerical solution of wave propagation problems
+    Geophysics p 1099-1110 - 1978
+    The finite difference approximation used by Reynolds for the transparent boundary condition is of first
+    order, though the scalar schema of propagation is of fourth order in space.
+
+    """
+    cdef unsigned int i, j, k
+
+    # Transparent boundary condition applied exactly at
+    # The edge of the fourth order calculation
+    # Sides
+    for j in xrange(ny):  #  faces west (left), east (right), down
+         for k in xrange(nz):
+             for p in xrange(2):
+                # left at x=0...1
+                u_tp1[k, j, 1-p] = ( u_t[k, j, 1-p] + u_t[k, j, 2-p] - u_tm1[k, j, 2-p] +
+                     (c[k, j, 1-p]*dt/dx)*(u_t[k, j, 2-p] - u_t[k, j, 1-p] - u_tm1[k, j, 3-p] + u_tm1[k, j, 2-p])
+                    )
+                #right at x=nx-1...nx-2
+                u_tp1[k, j, nx-2+p] = ( u_t[k, j, nx-2+p] + u_t[k, j, nx-3+p] - u_tm1[k, j, nx-3+p] -
+                     (c[k, j, nx-2+p]*dt/dx)*(u_t[k, j, nx-2+p] - u_t[k, j, nx-3+p] - u_tm1[k, j, nx-3+p] + u_tm1[k, j, nx-4+p])
+                    )
+         # Down at z=nz-2...nz-1
+         for i in xrange(nx):
+             for p in xrange(2):
+                u_tp1[nz-2+p, j, i] = ( u_t[nz-2+p,j, i] + u_t[nz-3+p,j, i] - u_tm1[nz-3+p,j, i] -
+                     (c[nz-2+p,j, i]*dt/dz)*(u_t[nz-2+p,j, i] - u_t[nz-3+p,j, i] - u_tm1[nz-3+p,j, i] + u_tm1[nz-4+p,j, i])
+                    )
+    # Faces front and back
+    for k in xrange(nz):
+         for i in xrange(nx):
+             for p in xrange(2):
+                #  back at y=1...0
+                u_tp1[k, 1-p, i] = ( u_t[k, 1-p, i] + u_t[k, 2-p, i] - u_tm1[k, 2-p, i] +
+                    (c[k, 1-p, i]*dt/dy)*(u_t[k, 2-p, i] - u_t[k, 1-p, i] - u_tm1[k, 3-p, i] + u_tm1[k, 2-p, i])
+                    )
+                #  front at y=ny-2..ny-1
+                u_tp1[k, ny-2+p, i] = ( u_t[k, ny-2+p, i] + u_t[k, ny-3+p, i] - u_tm1[k, ny-3+p, i] -
+                    (c[k, ny-2+p, i]*dt/dy)*(u_t[k, ny-2+p, i] - u_t[k, ny-3+p, i] - u_tm1[k, ny-3+p, i] + u_tm1[k, ny-4+p, i])
+                    )
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def _step_scalar3(
+    double[:,:,::1] u_tm1 not None,
+    double[:,:,::1] u_t not None,
+    double[:,:,::1] u_tp1 not None,
+    unsigned int x1, unsigned int x2,
+    unsigned int y1, unsigned int y2,
+    unsigned int z1, unsigned int z2,
+    double dt, double dx, double dy, double dz,
+    double[:,:, ::1] c not None):
+    """
+    Perform a single time step in the Finite Difference solution for scalar 3D
+    waves 4th order in space
+    """
+    cdef unsigned int k, i, j
+
+    for k in xrange(z1, z2):
+        for j in xrange(y1, y2):
+            for i in xrange(x1, x2):
+                u_tp1[k,j,i] = (2.*u_t[k,j,i] - u_tm1[k,j,i] +
+                    (((c[k,j,i]*dt/dx)**2)*(
+                        (-u_t[k,j,i+2] + 16.*u_t[k,j,i+1] - 30.*u_t[k,j,i] +
+                         16.*u_t[k,j,i-1] - u_t[k,j,i-2])/12.)
+                    + ((c[k,j,i]*dt/dz)**2)*(
+                        (-u_t[k+2,j,i] + 16.*u_t[k+1,j,i] - 30.*u_t[k,j,i] +
+                         16.*u_t[k-1,j,i] - u_t[k-2,j,i])/12.)
+                    + ((c[k,j,i]*dt/dy)**2)*(
+                        (-u_t[k,j+2,i] + 16.*u_t[k,j+1,i] - 30.*u_t[k,j,i] +
+                         16.*u_t[k,j-1,i] - u_t[k,j-2,i])/12.)))
+

fatiando/seismic/wavefd.py

@@ -635,42 +635,39 @@ def scalar(vel, area, dt, iterations, sources, stations=None,
     else:
         stations, seismograms = [], []
     # Add some padding to x and z. The padding region is where the wave is
-    # absorbed
+    # absorbed by gaussian dumping
     pad = int(padding)
     if pad == -1:   # default 5% percent nz
-        pad = int(0.05*nz) + 2
+        pad = int(0.05*nz) + 2  # plus 2 due 4th order
     nx += 2*pad
     nz += pad
     # Pad the velocity as well
     vel_pad = _add_pad(vel, pad, (nz, nx))
-    # Pack the particle position u at 2 different times in one 3d array
-    # u[0] = u(t-1)
-    # u[1] = u(t)
-    # The next time step overwrites the t-1 panel
-    u = numpy.zeros((2, nz, nx), dtype=numpy.float)
+    # Pack the particle position u at 3 different times in one 3d array
+    u = numpy.zeros((3, nz, nx), dtype=numpy.float)
     # Compute and yield the initial solutions
     for src in sources:
         i, j = src.indexes()
-        u[1, i, j + pad] += -((vel[i,j]*dt)**2)*src(0)
+        u[2, i, j + pad] += -((vel[i,j]*dt)**2)*src(0)
     # Update seismograms
     for station, seismogram in zip(stations, seismograms):
         i, j = station
-        seismogram[0] = u[1, i, j + pad]
+        seismogram[0] = u[2, i, j + pad]
     if snapshot is not None:
-        yield 0, u[1, :-pad, pad:-pad], seismograms
-
+        yield 0, u[2, :-pad, pad:-pad], seismograms
     for iteration in xrange(1, iterations):
-        t, tm1 = iteration%2, (iteration + 1)%2
-        tp1 = tm1
+        tm1, t, tp1 = iteration % 3, (iteration+1) % 3, (iteration+2) % 3  # to avoid copying between panels
+        # dumping not working with Reynolds need to fix apply dumping
+        # # Damp the regions in the padding to make waves go to infinity
+        # # apply dumping just outside the nonreflexive boundary conditions
+        # _apply_damping(u[tm1], nx-2, nz-2, pad-2, taper)
+        # _apply_damping(u[t], nx-2, nz-2, pad-2, taper)
         _step_scalar(u[tp1], u[t], u[tm1], 2, nx - 2, 2, nz - 2,
                      dt, ds, vel_pad)
+        # _apply_damping(u[tp1], nx-2, nz-2, pad-2, taper)
         # forth order +2-2 indexes needed
-        # Damp the regions in the padding to make waves go to infinity
-        #_apply_damping(u[t], nx, nz, pad, taper)
-        # not PML yet or anything similar
+        # apply Reynolds 1d plane wave absorbing condition
         _nonreflexive_scalar_boundary_conditions(u[tp1], u[t], u[tm1], vel_pad, dt, ds, nx, nz)
-        # Damp the regions in the padding to make waves go to infinity
-        #_apply_damping(u[tp1], nx, nz, pad, taper)
         for src in sources:
             i, j = src.indexes()
             u[tp1, i, j + pad] += -((vel[i,j]*dt)**2)*src(iteration*dt)
@@ -754,30 +751,24 @@ def scalar3(c, area, dt, iterations, sources, stations=None,
     ny += 2 * pad
     nz += pad
     c_pad = _add_pad(c, pad, (nz, ny, nx))
-    # Pack the particle position u at 2 different times in one 3d array
-    # u[0] = u(t-1)
-    # u[1] = u(t)
-    # The next time step overwrites the t-1 panel
-    u = numpy.zeros((2, nz, ny, nx), dtype=numpy.float)
+    # Pack the particle position u at 3 different times in one 3d array
+    u = numpy.zeros((3, nz, ny, nx), dtype=numpy.float)
     # Compute and yield the initial solutions
     for src in sources:
         k, j, i = src.indexes()
-        u[1, k, j + pad, i + pad] += - ((c[k, j, i]*dt)**2)*src(0)
+        u[2, k, j + pad, i + pad] += - ((c[k, j, i]*dt)**2)*src(0)
     # Update seismograms
     for station, seismogram in zip(stations, seismograms):
         k, j, i = station
-        seismogram[0] = u[1, k, j + pad, i + pad]
+        seismogram[0] = u[2, k, j + pad, i + pad]
     if snapshot is not None:
-        yield 0, u[1, :-pad, pad:-pad, pad:-pad], seismograms
+        yield 0, u[2, :-pad, pad:-pad, pad:-pad], seismograms
     for iteration in xrange(1, iterations):
-        t, tm1 = iteration % 2, (iteration + 1) % 2
-        tp1 = tm1
-        _step_scalar3(u[tp1], u[t], u[tm1], 2, nx - 2, 2, ny - 2,
+        tm1, t, tp1 = iteration % 3, (iteration + 1) % 3, (iteration + 2) % 3  # to avoid copying between panels
+        _step_scalar3(u[tm1], u[t], u[tp1], 2, nx - 2, 2, ny - 2,
                           2, nz - 2, dt, dx, dy, dz, c_pad)
-        _apply_damping3(u[t], nx, ny, nz, pad, taper)
-        # not PML yet or anything similar
+        # apply 1d Reynolds for plane waves
         _nonreflexive3_scalar_boundary_conditions(u[tm1], u[t], u[tp1], dt, dx, dy, dz, c_pad, nx, ny, nz)
-        _apply_damping3(u[tp1], nx, ny, nz, pad, taper)
         for src in sources:
             k, j, i = src.indexes()
             u[tp1, k, j + pad, i + pad] += -((c[k, j, i]*dt)**2)*src(iteration*dt)