pdr: coulomb failure stress modeling

src/gf/seismosizer.py

@@ -2612,7 +2612,7 @@ class PseudoDynamicRupture(SourceWithDerivedMagnitude):
 
     tractions = TractionField.T(
         optional=True,
-        help='Traction field the rupture plane is exposed to. See the'
+        help='Traction field the rupture plane is exposed to. See the '
              ':py:mod:`pyrocko.gf.tractions` module for more details. '
              'If ``tractions=None`` and ``rake`` is given'
              ' :py:class:`~pyrocko.gf.tractions.DirectedTractions` will'
@@ -4070,6 +4070,158 @@ class PseudoDynamicRupture(SourceWithDerivedMagnitude):
             magnitude=mt.magnitude,
             duration=t_max)
 
+    def get_coulomb_failure_stress(
+            self,
+            receiver_points,
+            friction,
+            pressure,
+            strike,
+            dip,
+            rake,
+            time=None,
+            *args,
+            **kwargs):
+        '''
+        Calculate Coulomb failure stress change CFS.
+
+        The function obtains the Coulomb failure stress change :math:`\\Delta
+        \\sigma_C` at arbitrary receiver points with a commonly oriented
+        receiver plane assuming:
+
+        .. math::
+
+            \\Delta \\sigma_C = \\sigma_S - \\mu (\\sigma_N - \\Delta p)
+
+        with the shear stress :math:`\\sigma_S`, the coefficient of friction
+        :math:`\\mu`, the normal stress :math:`\\sigma_N`, and the pore fluid
+        pressure change :math:`\\Delta p`. Each receiver point is characterized
+        by its geographical coordinates, and depth. The required receiver plane
+        orientation is defined by ``strike``, ``dip``, and ``rake``. The
+        Coulomb failure stress change is calculated for a given time after
+        rupture origin time.
+
+        :param receiver_points:
+            Location of the receiver points in Northing, Easting, and depth in
+            [m].
+        :type receiver_points:
+            :py:class:`~numpy.ndarray`: ``(n_receiver, 3)``
+
+        :param friction:
+            Coefficient of friction.
+        :type friction:
+            float
+
+        :param pressure:
+            Pore pressure change in [Pa].
+        :type pressure:
+            float
+
+        :param strike:
+            Strike of the receiver plane in [deg].
+        :type strike:
+            float
+
+        :param dip:
+            Dip of the receiver plane in [deg].
+        :type dip:
+            float
+
+        :param rake:
+            Rake of the receiver plane in [deg].
+        :type rake:
+            float
+
+        :param time:
+            Time after origin [s], for which the resulting :math:`\\Delta
+            \\Sigma_c` is computed. If not given, :math:`\\Delta \\Sigma_c` is
+            derived based on the final static slip.
+        :type time:
+            optional, float > 0.
+
+        :returns:
+            The Coulomb failure stress change :math:`\\Delta \\Sigma_c` at each
+            receiver point in [Pa].
+        :rtype:
+            :py:class:`~numpy.ndarray`: ``(n_receiver,)``
+        '''
+        # dislocation at given time
+        source_slip = self.get_slip(time=time, scale_slip=True)
+
+        # source planes
+        source_patches = num.array([
+            src.source_patch() for src in self.patches])
+
+        # earth model
+        lambda_mean = num.mean([src.lamb for src in self.patches])
+        mu_mean = num.mean([src.shearmod for src in self.patches])
+
+        # Dislocation and spatial derivatives from okada in NED
+        results = okada_ext.okada(
+            source_patches,
+            source_slip,
+            receiver_points,
+            lambda_mean,
+            mu_mean,
+            rotate_sdn=False,  # TODO Check
+            stack_sources=0,  # TODO Check
+            *args, **kwargs)
+
+        # resolve stress tensor (sum!)
+        diag_ind = [0, 4, 8]
+        kron = num.zeros(9)
+        kron[diag_ind] = 1.
+        kron = kron[num.newaxis, num.newaxis, :]
+
+        eps = 0.5 * (
+            results[:, :, 3:] +
+            results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
+
+        dilatation \
+            = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
+
+        stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps
+
+        # superposed stress of all sources at receiver locations
+        stress_sum = num.sum(stress, axis=0)
+
+        # get shear and normal stress from stress tensor
+        strike_rad = d2r * strike
+        dip_rad = d2r * dip
+        rake_rad = d2r * rake
+
+        n_rec = receiver_points.shape[0]
+        stress_normal = num.zeros(n_rec)
+        tau = num.zeros(n_rec)
+
+        # Get vectors in receiver fault normal (ns), strike (rst) and
+        # dip (rdi) direction
+        ns = num.zeros(3)
+        rst = num.zeros(3)
+        rdi = num.zeros(3)
+
+        ns[0] = num.sin(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
+        ns[1] = num.sin(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
+        ns[2] = -num.cos(dip_rad)
+
+        rst[0] = num.cos(strike_rad)
+        rst[1] = num.sin(strike_rad)
+        rst[2] = 0.0
+
+        rdi[0] = num.cos(dip_rad) * num.cos(strike_rad + 0.5 * num.pi)
+        rdi[1] = num.cos(dip_rad) * num.sin(strike_rad + 0.5 * num.pi)
+        rdi[2] = num.sin(dip_rad)
+
+        ts = rst * num.cos(rake_rad) - rdi * num.sin(rake_rad)
+
+        stress_normal = num.sum(
+            num.tile(ns, 3) * stress_sum * num.repeat(ns, 3), axis=1)
+
+        tau = num.sum(
+            num.tile(ts, 3) * stress_sum * num.repeat(ns, 3), axis=1)
+
+        # calculate cfs using formula above and return
+        return tau + friction * (stress_normal + pressure)
+
 
 class DoubleDCSource(SourceWithMagnitude):
     '''

src/modelling/ext/okada_ext.c

@@ -939,6 +939,12 @@ static okada_error_t dc3d_flexi(
     */
 
     rot_u(uokada, rotmat, u);
+
+    /*
+     * Optional rotation of displacement and strain vector/tensor into
+     * strike-dip-normal (sdn) coordinate system
+    */
+
     if (rot_sdn == 1) {
         euler_to_matrix((dip + 180.)*D2R, strike*D2R, 0.,  rotmat);
         rot_u(u, rotmat, u);

test/gf/test_gf_source_types.py

@@ -206,6 +206,126 @@ class GFSourceTypesTestCase(unittest.TestCase):
             plt.colorbar(im, label='Opening [m] after %.2f s' % deltat)
             plt.show()
 
+    @unittest.skipUnless(*have_store('iceland_reg_v2'))
+    def test_pseudo_dynamic_rupture_cfs(self):
+        from pyrocko.modelling import okada_ext
+        store_id = 'iceland_reg_v2'
+
+        engine = gf.get_engine()
+        store = engine.get_store(store_id)
+
+        source = gf.PseudoDynamicRupture(
+            length=20000., width=10000., depth=2000.,
+            anchor='top', gamma=0.8, dip=45., strike=60.,
+            slip=1., rake=0.,
+            nx=5, ny=3, smooth_rupture=False,
+            decimation_factor=1000)
+
+        source.discretize_patches(store)
+
+        target_strike = -92.
+        target_dip = 73.
+        target_rake = -8.
+
+        nnorths = 50
+        neasts = 50
+        norths = num.linspace(-50., 50., nnorths) * 1e3
+        easts = num.linspace(-50., 50., neasts) * 1e3
+        depth_target = 10e3
+
+        receiver_points = num.zeros((nnorths * neasts, 3))
+        receiver_points[:, 0] = num.repeat(norths, neasts)
+        receiver_points[:, 1] = num.tile(easts, nnorths)
+        receiver_points[:, 2] = num.ones(nnorths * neasts) * depth_target
+
+        friction = 0.6
+        pressure = 0.
+
+        cfs_init = source.get_coulomb_failure_stress(
+            receiver_points, friction=friction, pressure=pressure,
+            strike=target_strike, dip=target_dip, rake=target_rake,
+            nthreads=1)
+
+        # Calculate CFS from scratch
+        source_slip = source.get_slip(time=None, scale_slip=True)
+
+        # source planes
+        source_patches = num.array([
+            src.source_patch() for src in source.patches])
+
+        # earth model
+        lambda_mean = num.mean([src.lamb for src in source.patches])
+        mu_mean = num.mean([src.shearmod for src in source.patches])
+
+        # Dislocation and spatial derivatives from okada in NED
+        results = okada_ext.okada(
+            source_patches,
+            source_slip,
+            receiver_points,
+            lambda_mean,
+            mu_mean,
+            rotate_sdn=False,
+            stack_sources=0)
+
+        # resolve stress tensor (sum!)
+        diag_ind = [0, 4, 8]
+        kron = num.zeros(9)
+        kron[diag_ind] = 1.
+        kron = kron[num.newaxis, num.newaxis, :]
+
+        eps = 0.5 * (
+            results[:, :, 3:] +
+            results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
+
+        dilatation \
+            = eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
+
+        stress = kron*lambda_mean*dilatation + 2.*mu_mean*eps
+        # stress shape (n_sources, n_receivers, n_stress_components)
+
+        # superposed stress of all sources at receiver locations
+        stress_sum = num.sum(stress, axis=0)
+        # stress_sum shape: (n_receivers, n_stress_components)
+
+        # get shear and normal stress from stress tensor
+        st0 = d2r * target_strike
+        di0 = d2r * target_dip
+        ra0 = d2r * target_rake
+
+        n_rec = receiver_points.shape[0]
+        stress_normal = num.zeros(n_rec)
+        tau = num.zeros(n_rec)
+
+        for irec in range(n_rec):
+            ns = num.zeros(3)
+            rst = num.zeros(3)
+            rdi = num.zeros(3)
+
+            ns[0] = num.sin(di0) * num.cos(st0 + 0.5 * num.pi)
+            ns[1] = num.sin(di0) * num.sin(st0 + 0.5 * num.pi)
+            ns[2] = -num.cos(di0)
+
+            rst[0] = num.cos(st0)
+            rst[1] = num.sin(st0)
+            rst[2] = 0.0
+
+            rdi[0] = num.cos(di0) * num.cos(st0 + 0.5 * num.pi)
+            rdi[1] = num.cos(di0) * num.sin(st0 + 0.5 * num.pi)
+            rdi[2] = num.sin(di0)
+
+            ts = rst * num.cos(ra0) - rdi * num.sin(ra0)
+
+            for j in range(3):
+                for i in range(3):
+                    stress_normal[irec] += \
+                        ns[i] * stress_sum[irec, j*3 + i] * ns[j]
+                    tau[irec] += ts[i] * stress_sum[irec, j*3 + i] * ns[j]
+
+        # calculate cfs using formula above and return
+        cfs_comp = tau + friction * (stress_normal + pressure)
+
+        num.testing.assert_allclose(cfs_init, cfs_comp, atol=1.)
+
     def test_pseudo_dynamic_rupture_outline(self):
         length = 20000.
         width = 10000.