seismosizer.RectangularSource add opening_fraction argument

CHANGELOG.md

@@ -7,6 +7,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 ## [Unreleased]
 
 ### Added
+- RectangularSource: added opening_fraction to model tensile dislocations
 - New command line option for jackseis: `--record-length`
 - Timing definition offsets can now take `%` as suffix to scale phase
   traveltimes relatively.
@@ -24,6 +25,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/).
 - Fix typos in STA/LTA documentation.
 - Fomosto PSGRN/PSCMP backend: improved control of modelling parameters,
   fixes some accuracy issues, regarding the spacial sampling interval.
+- Fomosto PSGRN/PSCMP backend: fixed scaling of isotropic GF components
 - Improved handling of differing sampling rates and interpolation settings
   when modelling multiple targets in `gf.Engine.process`.
 - PyQt compat issues with MacOS Big Sur.

doc/source/library/examples/gf_forward.rst

@@ -24,7 +24,9 @@ Download :download:`gf_forward_example1.py </../../examples/gf_forward_example1.
 
 
 Calculate spatial surface displacement from a local GF store
--------------------------------------------------------------
+------------------------------------------------------------
+Shear dislocation - Earthquake
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 In this example we create a :class:`~pyrocko.gf.seismosizer.RectangularSource` and compute the spatial static displacement invoked by that rupture.
 
@@ -43,6 +45,28 @@ Download :download:`gf_forward_example2.py </../../examples/gf_forward_example2.
     :language: python
 
 
+Tensile dislocation - Sill/Dike
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+In this example we create a :class:`~pyrocko.gf.seismosizer.RectangularSource` and compute the spatial static displacement invoked by a
+magmatic contracting sill. The same model can be used to model a magmatic dike intrusion (changing the "dip" argument).
+
+We will utilize :class:`~pyrocko.gf.seismosizer.LocalEngine`, :class:`~pyrocko.gf.targets.StaticTarget` and :class:`~pyrocko.gf.targets.SatelliteTarget`.
+
+.. figure:: /static/gf_static_displacement_sill.png
+    :align: center
+    :width: 90%
+    :alt: Static displacement from a contracting sill calculated through pyrocko
+
+    Synthetic surface displacement from a contracting sill. The sill has a strike of 104° N. The surface displacements are shown in Line-of-sight (LOS), east, north and vertical directions. Envisat satellite has a look angle of 23° and heading -76°. The motion is positive towards the satellite LOS.
+
+Download :download:`gf_forward_example2_sill.py </../../examples/gf_forward_example2_sill.py>`
+
+.. literalinclude :: /../../examples/gf_forward_example2_sill.py
+    :language: python
+
+
+
 Calculate spatial surface displacement and export Kite scenes
 -------------------------------------------------------------
 

examples/gf_forward_example2_sill.py

@@ -0,0 +1,116 @@
+import os.path
+from pyrocko import gf
+import numpy as num
+
+# Download a Greens Functions store, programmatically.
+store_id = 'gf_abruzzo_nearfield_vmod_Ameri'
+if not os.path.exists(store_id):
+    gf.ws.download_gf_store(site='kinherd', store_id=store_id)
+
+# Setup the LocalEngine and point it to the fomosto store you just downloaded.
+# *store_superdirs* is a list of directories where to look for GF Stores.
+engine = gf.LocalEngine(store_superdirs=['.'])
+
+# We define an extended source, in this case a rectangular geometry
+# Centroid UTM position is defined relatively to geographical lat, lon position
+# Horizontal closing sill with an N104W azimuth.
+# Slip is split to shear and tensile slip where "opening_fraction" determines
+# the direction and amount of opening/closing defined from -1, 1
+# for a pure shear dislocation "opening_fraction" is 0.
+
+km = 1e3  # for convenience
+d2r = num.pi / 180.
+
+rect_source = gf.RectangularSource(
+    lat=0., lon=0.,
+    north_shift=0., east_shift=0., depth=2.5*km,
+    width=4*km, length=8*km,
+    dip=0., rake=0., strike=104.,
+    slip=3., opening_fraction=-1.)
+
+# We will define a grid of targets
+# number in east and north directions, and total
+ngrid = 80
+
+# extension from origin in all directions
+obs_size = 20.*km
+ntargets = ngrid**2
+
+# make regular line vector
+norths = num.linspace(-obs_size, obs_size, ngrid)
+easts = num.linspace(-obs_size, obs_size, ngrid)
+
+# make regular grid
+norths2d = num.repeat(norths, easts.size)
+easts2d = num.tile(easts, norths.size)
+
+# We initialize the satellite target and set the line of sight vectors
+# direction, example of the Envisat satellite
+look = 23.     # angle between the LOS and the vertical
+heading = -76  # angle between the azimuth and the east (anti-clock)
+theta = num.empty(ntargets)  # vertical LOS from horizontal
+theta.fill((90. - look) * d2r)
+phi = num.empty(ntargets)  # horizontal LOS from E in anti-clokwise rotation
+phi.fill((-90 - heading) * d2r)
+
+satellite_target = gf.SatelliteTarget(
+    north_shifts=norths2d,
+    east_shifts=easts2d,
+    tsnapshot=24. * 3600.,  # one day
+    interpolation='nearest_neighbor',
+    phi=phi,
+    theta=theta,
+    store_id=store_id)
+
+# The computation is performed by calling process on the engine
+result = engine.process(rect_source, [satellite_target])
+
+
+def plot_static_los_result(result, target=0):
+    '''Helper function for plotting the displacement'''
+
+    import matplotlib.pyplot as plt
+    import matplotlib.ticker as tick
+
+    # get target coordinates and displacements from results
+    N = result.request.targets[target].coords5[:, 2]
+    E = result.request.targets[target].coords5[:, 3]
+    synth_disp = result.results_list[0][target].result
+
+    # get the component names of displacements
+    components = synth_disp.keys()
+    fig, _ = plt.subplots(int(len(components)/2), int(len(components)/2))
+
+    vranges = [(synth_disp[k].max(),
+                synth_disp[k].min()) for k in components]
+
+    for comp, ax, vrange in zip(components, fig.axes, vranges):
+
+        lmax = num.abs([num.min(vrange), num.max(vrange)]).max()
+
+        # plot displacements at targets as colored points
+        cmap = ax.scatter(E, N, c=synth_disp[comp], s=10., marker='s',
+                          edgecolor='face', cmap='seismic',
+                          vmin=-1.5*lmax, vmax=1.5*lmax)
+
+        ax.set_title(comp+' [m]')
+        ax.set_aspect('equal')
+        ax.set_xlim(-obs_size, obs_size)
+        ax.set_ylim(-obs_size, obs_size)
+        # We plot the modeled fault
+        n, e = rect_source.outline(cs='xy').T
+        ax.fill(e, n, color=(0.5, 0.5, 0.5), alpha=0.5)
+
+        fig.colorbar(cmap, ax=ax, aspect=5)
+
+        # reduce number of ticks
+        yticker = tick.MaxNLocator(nbins=5)
+        yax = ax.get_yaxis()
+        xax = ax.get_xaxis()
+        yax.set_major_locator(yticker)
+        xax.set_major_locator(yticker)
+
+    plt.show()
+
+
+plot_static_los_result(result)

src/fomosto/psgrn_pscmp.py

@@ -648,6 +648,27 @@ MTDev = {
     }
 
 
+def get_nullification_factor(mu, lame_lambda):
+    '''
+    Factor that needs to be multiplied to 2 of the tensile sources to
+    eliminate two of the isotropic components.
+    '''
+    return -lame_lambda / (2. * mu + 2. * lame_lambda)
+
+
+def get_trace_normalization_factor(mu, lame_lambda, nullification):
+    '''
+    Factor to be multiplied to elementary GF trace to have unit displacement.
+    '''
+    return 1. / ((2. * mu) + lame_lambda + (2. * lame_lambda * nullification))
+
+
+def get_iso_scaling_factors(mu, lame_lambda):
+    nullf = get_nullification_factor(mu, lame_lambda)
+    scale = get_trace_normalization_factor(mu, lame_lambda, nullf)
+    return nullf, scale
+
+
 class PsCmpTensileSF(Location, gf.seismosizer.Cloneable):
     '''
     Compound dislocation of 3 perpendicular, rectangular sources to approximate
@@ -670,8 +691,7 @@ class PsCmpTensileSF(Location, gf.seismosizer.Cloneable):
         default='nn',
         help='Axis index for opening direction; "nn", "ee", or "dd"')
 
-    def to_rfs(self):
-        rf = -0.25
+    def to_rfs(self, nullification):
 
         cmpd = []
         for comp, mt in MTIso.items():
@@ -679,7 +699,7 @@ class PsCmpTensileSF(Location, gf.seismosizer.Cloneable):
 
             if comp != self.idx:
                 params = copy.deepcopy(mt)
-                params['opening'] *= rf
+                params['opening'] *= nullification
 
             kwargs = copy.deepcopy(self.__dict__)
             kwargs.update(params)
@@ -742,13 +762,13 @@ class PsCmpMomentTensor(Location, gf.seismosizer.Cloneable):
         help='Axis index for MT component; '
              '"nn", "ee", "dd", "ne", "nd", or "ed".')
 
-    def to_rfs(self):
+    def to_rfs(self, nullification=-0.25):
         kwargs = copy.deepcopy(self.__dict__)
         kwargs.pop('_latlon')
 
         if self.idx in MTIso:
             tsf = PsCmpTensileSF(**kwargs)
-            return tsf.to_rfs()
+            return tsf.to_rfs(nullification)
 
         elif self.idx in MTDev:
             ssf = PsCmpShearSF(**kwargs)
@@ -1409,7 +1429,13 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
             lat=dummy_lat,
             lon=dummy_lon,
             points=points,
-            interpolation='nearest_neighbor')
+            interpolation='multilinear')
+
+        self.lambda_moduli = storeconf.get_lambda_moduli(
+            lat=dummy_lat,
+            lon=dummy_lon,
+            points=points,
+            interpolation='multilinear')
 
         if step == 0:
             block_size = (1, storeconf.nsource_depths, storeconf.ndistances)
@@ -1472,6 +1498,7 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
             (sz, firstx), (sz, lastx), (ns, nx) = \
                 self.get_block_extents(iblock)
             mu = self.shear_moduli[iblock]
+            lame_lambda = self.lambda_moduli[iblock]
 
             rz = self.store.config.receiver_depth
         else:
@@ -1538,11 +1565,13 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
             mtsize = float(dx * cc.rectangular_fault_size_factor)
 
             Ai = 1. / num.power(mtsize, 2)
+            nullf, sf = get_iso_scaling_factors(mu=mu, lame_lambda=lame_lambda)
+
             mui = 1. / mu
 
             gfmapping = [
                 (('nn',),
-                 {'un': (0, 0.4 * Ai * mui), 'ud': (5, 0.4 * Ai * mui)}),
+                 {'un': (0, Ai * sf), 'ud': (5, Ai * sf)}),
                 (('ne',),
                  {'ue': (3, 1 * Ai * mui)}),
                 (('nd',),
@@ -1550,9 +1579,9 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
                 (('ed',),
                  {'ue': (4, 1 * Ai * mui)}),
                 (('dd',),
-                 {'un': (2, 0.4 * Ai * mui), 'ud': (7, 0.4 * Ai * mui)}),
+                 {'un': (2, Ai * sf), 'ud': (7, Ai * sf)}),
                 (('ee',),
-                 {'un': (8, 0.4 * Ai * mui), 'ud': (9, 0.4 * Ai * mui)}),
+                 {'un': (8, Ai * sf), 'ud': (9, Ai * sf)}),
                 ]
 
             for mt, gfmap in gfmapping:
@@ -1566,7 +1595,8 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
                         length=mtsize,
                         idx=idx)
 
-                    cc.rectangular_source_patches.extend(pmt.to_rfs())
+                    cc.rectangular_source_patches.extend(
+                        pmt.to_rfs(nullf))
 
                 runner.run(cc)
 

src/gf/meta.py

@@ -1734,6 +1734,52 @@ class Config(Object):
                                           parameter='shear_moduli',
                                           interpolation=interpolation)
 
+    def get_lambda_moduli(self, lat, lon, points,
+                          interpolation=None):
+        '''
+        Get lambda moduli at given points from contained velocity model.
+
+        :param lat: surface origin for coordinate system of ``points``
+        :param points: NumPy array of shape ``(N, 3)``, where each row is
+            a point ``(north, east, depth)``, relative to origin at
+            ``(lat, lon)``
+        :param interpolation: interpolation method. Choose from
+            ``('nearest_neighbor', 'multilinear')``
+        :returns: NumPy array of length N with extracted shear moduli at each
+            point
+
+        The default implementation retrieves and interpolates the lambda moduli
+        from the contained 1D velocity profile.
+        '''
+        return self.get_material_property(lat, lon, points,
+                                          parameter='lambda_moduli',
+                                          interpolation=interpolation)
+
+    def get_bulk_moduli(self, lat, lon, points,
+                        interpolation=None):
+        '''
+        Get bulk moduli at given points from contained velocity model.
+
+        :param lat: surface origin for coordinate system of ``points``
+        :param points: NumPy array of shape ``(N, 3)``, where each row is
+            a point ``(north, east, depth)``, relative to origin at
+            ``(lat, lon)``
+        :param interpolation: interpolation method. Choose from
+            ``('nearest_neighbor', 'multilinear')``
+        :returns: NumPy array of length N with extracted shear moduli at each
+            point
+
+        The default implementation retrieves and interpolates the lambda moduli
+        from the contained 1D velocity profile.
+        '''
+        lambda_moduli = self.get_material_property(
+            lat, lon, points, parameter='lambda_moduli',
+            interpolation=interpolation)
+        shear_moduli = self.get_material_property(
+            lat, lon, points, parameter='shear_moduli',
+            interpolation=interpolation)
+        return lambda_moduli + (2 / 3) * shear_moduli
+
     def get_vs(self, lat, lon, points, interpolation=None):
         '''
         Get Vs at given points from contained velocity model.
@@ -1831,6 +1877,22 @@ class Config(Object):
                 store_depth_profile, z_profile, rho_profile)
 
             profile = num.power(store_vs_profile, 2) * store_rho_profile
+
+        elif parameter == 'lambda_moduli':
+            vs_profile = earthmod.profile('vs')
+            vp_profile = earthmod.profile('vp')
+            rho_profile = earthmod.profile('rho')
+
+            store_vs_profile = num.interp(
+                store_depth_profile, z_profile, vs_profile)
+            store_vp_profile = num.interp(
+                store_depth_profile, z_profile, vp_profile)
+            store_rho_profile = num.interp(
+                store_depth_profile, z_profile, rho_profile)
+
+            profile = store_rho_profile * (
+                num.power(store_vp_profile, 2) -
+                num.power(store_vs_profile, 2) * 2)
         else:
             raise TypeError(
                 'parameter %s not available' % parameter)

src/gf/seismosizer.py

@@ -1968,6 +1968,13 @@ class RectangularSource(SourceWithDerivedMagnitude):
         optional=True,
         help='Slip on the rectangular source area [m]')
 
+    opening_fraction = Float.T(
+        default=0.,
+        help='Determines fraction of slip related to opening. '
+             '(``-1``: pure tensile closing, '
+             '``0``: pure shear, '
+             '``1``: pure tensile opening)')
+
     decimation_factor = Int.T(
         optional=True,
         default=1,
@@ -2016,7 +2023,12 @@ class RectangularSource(SourceWithDerivedMagnitude):
                     'interpolation method are available.')
 
             amplitudes = self._discretize(store, target)[2]
-            return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
+            if amplitudes.ndim == 2:
+                # CLVD component has no net moment, leave out
+                return float(pmt.moment_to_magnitude(
+                    num.sum(num.abs(amplitudes[0:2, :]).sum())))
+            else:
+                return float(pmt.moment_to_magnitude(num.sum(amplitudes)))
 
         else:
             return float(pmt.moment_to_magnitude(1.0))
@@ -2024,6 +2036,12 @@ class RectangularSource(SourceWithDerivedMagnitude):
     def get_factor(self):
         return 1.0
 
+    def get_slip_tensile(self):
+        return self.slip * self.opening_fraction
+
+    def get_slip_shear(self):
+        return self.slip - abs(self.get_slip_tensile)
+
     def _discretize(self, store, target):
         if self.nucleation_x is not None:
             nucx = self.nucleation_x * 0.5 * self.length
@@ -2059,13 +2077,38 @@ class RectangularSource(SourceWithDerivedMagnitude):
                 points=points,
                 interpolation=interpolation)
 
-            amplitudes *= dl * dw * shear_moduli * self.slip
+            tensile_slip = self.get_slip_tensile()
+            shear_slip = self.slip - abs(tensile_slip)
+
+            amplitudes_total = [shear_moduli * shear_slip]
+            if tensile_slip != 0:
+                bulk_moduli = store.config.get_bulk_moduli(
+                    self.lat, self.lon,
+                    points=points,
+                    interpolation=interpolation)
+
+                tensile_iso = bulk_moduli * tensile_slip
+                tensile_clvd = (2. / 3.) * shear_moduli * tensile_slip
+
+                amplitudes_total.extend([tensile_iso, tensile_clvd])
+
+            amplitudes_total = num.vstack(amplitudes_total).squeeze() * \
+                amplitudes * dl * dw
+
         else:
             # normalization to retain total moment
-            amplitudes /= num.sum(amplitudes)
-            amplitudes *= self.get_moment(store, target)
+            amplitudes_norm = amplitudes / num.sum(amplitudes)
+            moment = self.get_moment(store, target)
+
+            amplitudes_total = [
+                amplitudes_norm * moment * (1 - abs(self.opening_fraction))]
+            if self.opening_fraction != 0.:
+                amplitudes_total.append(
+                    amplitudes_norm * self.opening_fraction * moment)
+
+            amplitudes_total = num.vstack(amplitudes_total).squeeze()
 
-        return points, times, amplitudes, dl, dw
+        return points, times, num.atleast_1d(amplitudes_total), dl, dw
 
     def discretize_basesource(self, store, target=None):
 
@@ -2073,9 +2116,43 @@ class RectangularSource(SourceWithDerivedMagnitude):
 
         mot = pmt.MomentTensor(
             strike=self.strike, dip=self.dip, rake=self.rake)
-
         m6s = num.repeat(mot.m6()[num.newaxis, :], times.size, axis=0)
-        m6s[:, :] *= amplitudes[:, num.newaxis]
+
+        if amplitudes.ndim == 1:
+            m6s[:, :] *= amplitudes[:, num.newaxis]
+        elif amplitudes.ndim == 2:
+            # shear MT components
+            rotmat1 = pmt.euler_to_matrix(
+                d2r * self.dip, d2r * self.strike, d2r * -self.rake)
+            m6s[:, :] *= amplitudes[0, :][:, num.newaxis]
+
+            if amplitudes.shape[0] == 2:
+                # tensile MT components - moment/magnitude input
+                tensile = pmt.symmat6(1., 1., 3., 0., 0., 0.)
+                rot_tensile = pmt.to6(rotmat1.T * tensile * rotmat1)
+
+                m6s_tensile = rot_tensile[
+                    num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
+                m6s += m6s_tensile
+
+            elif amplitudes.shape[0] == 3:
+                # tensile MT components - slip input
+                iso = pmt.symmat6(1., 1., 1., 0., 0., 0.)
+                clvd = pmt.symmat6(-1., -1., 2., 0., 0., 0.)
+
+                rot_iso = pmt.to6(rotmat1.T * iso * rotmat1)
+                rot_clvd = pmt.to6(rotmat1.T * clvd * rotmat1)
+
+                m6s_iso = rot_iso[
+                    num.newaxis, :] * amplitudes[1, :][:, num.newaxis]
+                m6s_clvd = rot_clvd[
+                    num.newaxis, :] * amplitudes[2, :][:, num.newaxis]
+                m6s += m6s_iso + m6s_clvd
+            else:
+                raise ValueError('Unknwown amplitudes shape!')
+        else:
+            raise ValueError(
+                'Unexpected dimension of {}'.format(amplitudes.ndim))
 
         ds = meta.DiscretizedMTSource(
             lat=self.lat,

test/gf/test_gf_psgrn_pscmp.py

@@ -7,6 +7,8 @@ from tempfile import mkdtemp
 import numpy as num
 import os
 import shutil
+from copy import deepcopy
+from multiprocessing import cpu_count
 
 from pyrocko import orthodrome as ortd
 from pyrocko import util, gf, cake  # noqa
@@ -61,20 +63,41 @@ def statics(engine, source, starget):
 class GFPsgrnPscmpTestCase(unittest.TestCase):
 
     tempdirs = []
+    pscmp_store_dir = None
 
     def __init__(self, *args, **kwargs):
-        self.pscmp_store_dir = None
         unittest.TestCase.__init__(self, *args, **kwargs)
 
+    @classmethod
+    def setUpClass(cls):
+        cls.pscmp_store_dir, cls.psgrn_config = \
+            cls._create_psgrn_pscmp_store(cls)
+
     @classmethod
     def tearDownClass(cls):
         for d in cls.tempdirs:
             shutil.rmtree(d)
 
+    def test_isolated_isotropic_component(self):
+        shear = 31126160000.0
+        lame_lambda = 25211680000.0
+
+        cdm = num.array([
+            [(2 * shear + lame_lambda), lame_lambda, lame_lambda],
+            [lame_lambda, (2 * shear + lame_lambda), lame_lambda],
+            [lame_lambda, lame_lambda, (2 * shear + lame_lambda)]])
+
+        nullf = psgrn_pscmp.get_nullification_factor(
+            mu=shear, lame_lambda=lame_lambda)
+
+        cdm[0:2, :] *= nullf
+        null_cdm = cdm.sum(0)
+
+        num.testing.assert_allclose(null_cdm[0:2], num.zeros(2))
+
     def get_pscmp_store_info(self):
         if self.pscmp_store_dir is None:
-            self.pscmp_store_dir, self.psgrn_config = \
-                self._create_psgrn_pscmp_store()
+            raise ValueError('Store does not exist!')
 
         return self.pscmp_store_dir, self.psgrn_config
 
@@ -100,6 +123,7 @@ mantle
  660. 10.79 5.965 4.229 1349. 549.6 5e17 1e19 1'''.lstrip()))
 
         store_dir = mkdtemp(prefix='gfstore')
+        # store_dir = '/tmp/pscmp_gfstore'
         self.tempdirs.append(store_dir)
         store_id = 'psgrn_pscmp_test'
 
@@ -124,60 +148,76 @@ mantle
             receiver_depth=0. * km,
             source_depth_min=0. * km,
             source_depth_max=15. * km,
-            source_depth_delta=.5 * km,
+            source_depth_delta=.25 * km,
             distance_min=0. * km,
             distance_max=50. * km,
-            distance_delta=.5 * km)
+            distance_delta=.25 * km)
 
         config.validate()
-        gf.store.Store.create_editables(
-            store_dir, config=config, extra={'psgrn_pscmp': c})
-
-        store = gf.store.Store(store_dir, 'r')
-        store.close()
-
-        # build store
-        try:
-            psgrn_pscmp.build(store_dir, nworkers=1)
-        except psgrn_pscmp.PsCmpError as e:
-            if str(e).find('could not start psgrn/pscmp') != -1:
-                logger.warn('psgrn/pscmp not installed; '
-                            'skipping test_pyrocko_gf_vs_pscmp')
-                return
-            else:
-                raise e
+        if not os.path.exists(os.path.join(store_dir, 'config')):
+            gf.store.Store.create_editables(
+                store_dir, config=config, extra={'psgrn_pscmp': c})
+
+            store = gf.store.Store(store_dir, 'r')
+            store.close()
+
+            # build store
+            try:
+                psgrn_pscmp.build(store_dir, nworkers=cpu_count())
+            except psgrn_pscmp.PsCmpError as e:
+                if str(e).find('could not start psgrn/pscmp') != -1:
+                    logger.warn('psgrn/pscmp not installed; '
+                                'skipping test_pyrocko_gf_vs_pscmp')
+                    return
+                else:
+                    raise e
 
         return store_dir, c
 
-    def test_fomosto_vs_psgrn_pscmp(self):
+    def fomosto_vs_psgrn_pscmp(self, pscmp_sources, gf_sources, atol=2*mm):
+
+        def plot_components_compare(fomosto_comps, psgrn_comps):
+            import matplotlib.pyplot as plt
+
+            fig, axes = plt.subplots(4, 3)
+
+            for i, (fcomp, pscomp, cname) in enumerate(
+                    zip(fomosto_comps, psgrn_comps, ['N', 'E', 'D'])):
+                fdispl = fcomp.reshape(nnorth, neast)
+                pdispl = pscomp.reshape(nnorth, neast)
+                pcbound = num.max([num.abs(pdispl.min()), pdispl.max()])
+                # fcbound = num.max([num.abs(fdispl.min()), fdispl.max()])
+
+                axes[0, i].imshow(
+                    pdispl, cmap='seismic', vmin=-pcbound, vmax=pcbound)
+                axes[1, i].imshow(
+                    fdispl, cmap='seismic', vmin=-pcbound, vmax=pcbound)
+                diff = pdispl - fdispl
+                rdiff = pdispl / fdispl
+                axes[2, i].imshow(diff, cmap='seismic')
+                axes[3, i].imshow(rdiff, cmap='seismic')
+
+                axes[0, i].set_title('PSCMP %s' % cname)
+                axes[1, i].set_title('Fomosto %s' % cname)
+                axes[2, i].set_title(
+                    'abs diff min max %f, %f' % (diff.min(), diff.max()))
+                axes[3, i].set_title(
+                    'rel diff min max %f, %f' % (rdiff.min(), rdiff.max()))
+
+            plt.show()
+
         store_dir, c = self.get_pscmp_store_info()
 
         origin = gf.Location(
             lat=10.,
             lon=-15.)
 
-        # test GF store
-        TestRF = dict(
-            lat=origin.lat,
-            lon=origin.lon,
-            depth=2. * km,
-            width=0.2 * km,
-            length=0.5 * km,
-            rake=uniform(-90., 90.),
-            dip=uniform(0., 90.),
-            strike=uniform(0., 360.),
-            slip=uniform(1., 5.))
-
-        source_plain = gf.RectangularSource(**TestRF)
-        source_with_time = gf.RectangularSource(time=123.5, **TestRF)
-
         N, E = num.meshgrid(num.linspace(-20. * km, 20. * km, nnorth),
                             num.linspace(-20. * km, 20. * km, neast))
 
         # direct pscmp output
         lats, lons = ortd.ne_to_latlon(
             origin.lat, origin.lon, N.flatten(), E.flatten())
-        pscmp_sources = [psgrn_pscmp.PsCmpRectangularSource(**TestRF)]
 
         cc = c.pscmp_config
         cc.observation = psgrn_pscmp.PsCmpScatter(lats=lats, lons=lons)
@@ -218,18 +258,120 @@ mantle
 
         engine = gf.LocalEngine(store_dirs=[store_dir])
 
-        for source in [source_plain, source_with_time]:
+        for source in gf_sources:
             t0 = time()
             r = engine.process(source, [starget_nn, starget_ml])
             logger.info('pyrocko stacking time %f' % (time() - t0))
-            for static_result in r.static_results():
+            for i, static_result in enumerate(r.static_results()):
                 un_fomosto = static_result.result['displacement.n']
                 ue_fomosto = static_result.result['displacement.e']
                 ud_fomosto = static_result.result['displacement.d']
 
-                num.testing.assert_allclose(un_fomosto, un_pscmp, atol=1*mm)
-                num.testing.assert_allclose(ue_fomosto, ue_pscmp, atol=1*mm)
-                num.testing.assert_allclose(ud_fomosto, ud_pscmp, atol=1*mm)
+                if show_plot:
+                    fomosto_comps = [un_fomosto, ue_fomosto, ud_fomosto]
+                    psgrn_comps = [un_pscmp, ue_pscmp, ud_pscmp]
+                    plot_components_compare(fomosto_comps, psgrn_comps)
+
+                num.testing.assert_allclose(un_fomosto, un_pscmp, atol=atol)
+                num.testing.assert_allclose(ue_fomosto, ue_pscmp, atol=atol)
+                num.testing.assert_allclose(ud_fomosto, ud_pscmp, atol=atol)
+
+    def test_fomosto_vs_psgrn_pscmp_shear(self):
+
+        origin = gf.Location(
+            lat=10.,
+            lon=-15.)
+
+        # test GF store
+        TestRF = dict(
+            lat=origin.lat,
+            lon=origin.lon,
+            depth=2. * km,
+            width=0.2 * km,
+            length=7. * km,
+            rake=uniform(-90., 90.),
+            dip=uniform(0., 90.),
+            strike=uniform(0., 360.),
+            slip=uniform(1., 5.))
+
+        source_plain = gf.RectangularSource(**TestRF)
+        source_with_time = gf.RectangularSource(time=123.5, **TestRF)
+
+        gf_sources = [source_plain, source_with_time]
+        pscmp_sources = [psgrn_pscmp.PsCmpRectangularSource(**TestRF)]
+
+        self.fomosto_vs_psgrn_pscmp(
+            pscmp_sources=pscmp_sources, gf_sources=gf_sources, atol=5*mm)
+
+    def test_fomosto_vs_psgrn_pscmp_tensile(self):
+
+        origin = gf.Location(
+            lat=10.,
+            lon=-15.)
+
+        # test GF store
+        TestRF = dict(
+            lat=origin.lat,
+            lon=origin.lon,
+            depth=2. * km,
+            width=2. * km,
+            length=5. * km,
+            rake=uniform(-90., 90.),
+            dip=uniform(0., 90.),
+            strike=uniform(0., 360.))
+
+        slip = uniform(1., 2.)
+
+        for open_mode in [-1., 1.]:   # closing, opening
+
+            source_plain = gf.RectangularSource(**TestRF)
+            source_plain.update(slip=slip, opening_fraction=open_mode)
+
+            source_with_time = deepcopy(source_plain)
+            source_with_time.update(time=123.5)
+
+            gf_sources = [source_plain, source_with_time]
+
+            pscmp_sources = [psgrn_pscmp.PsCmpRectangularSource(
+                opening=slip * open_mode, slip=0., **TestRF)]
+
+            self.fomosto_vs_psgrn_pscmp(
+                pscmp_sources=pscmp_sources, gf_sources=gf_sources, atol=7*mm)
+
+    def test_fomosto_vs_psgrn_pscmp_tensile_shear(self):
+
+        origin = gf.Location(
+            lat=10.,
+            lon=-15.)
+
+        # test GF store
+        TestRF = dict(
+            lat=origin.lat,
+            lon=origin.lon,
+            depth=3. * km,
+            width=2. * km,
+            length=5. * km,
+            rake=uniform(-90., 90.),
+            dip=uniform(0., 90.),
+            strike=uniform(0., 360.))
+
+        opening_fraction = 0.4
+        slip = uniform(1., 5.)
+        opening = slip * opening_fraction
+        pscmp_slip = slip - opening
+
+        source_plain = gf.RectangularSource(**TestRF)
+        source_plain.update(slip=slip, opening_fraction=opening_fraction)
+
+        source_with_time = deepcopy(source_plain)
+        source_with_time.update(time=123.5)
+
+        gf_sources = [source_plain, source_with_time]
+        pscmp_sources = [psgrn_pscmp.PsCmpRectangularSource(
+            opening=opening, slip=pscmp_slip, **TestRF)]
+
+        self.fomosto_vs_psgrn_pscmp(
+            pscmp_sources=pscmp_sources, gf_sources=gf_sources, atol=7*mm)
 
     def plot_gf_distance_sampling(self):
         origin = gf.Location(
@@ -250,8 +392,8 @@ mantle
 
         source_plain = gf.RectangularSource(**TestRF)
 
-        N, E = num.meshgrid(num.linspace(-40. * km, 40. * km, nnorth),
-                            num.linspace(-40. * km, 40. * km, neast))
+        N, E = num.meshgrid(num.linspace(-20. * km, 20. * km, nnorth),
+                            num.linspace(-20. * km, 20. * km, neast))
 
         starget_ml = gf.StaticTarget(
             lats=num.full(N.size, origin.lat),
@@ -364,5 +506,5 @@ mantle
 
 
 if __name__ == '__main__':
-    util.setup_logging('test_gf_psgrn_pscmp', 'warning')
+    util.setup_logging('test_gf_psgrn_pscmp', 'info')
     unittest.main()

test/gf/test_gf_sources.py

@@ -403,6 +403,8 @@ class GFSourcesTestCase(unittest.TestCase):
 
         store = self.dummy_homogeneous_store()
 
+        depth = 10 * km
+        # shear
         rect1 = gf.RectangularSource(
             depth=10*km,
             magnitude=5.0,
@@ -410,10 +412,10 @@ class GFSourcesTestCase(unittest.TestCase):
             length=5*km)
 
         rect2 = gf.RectangularSource(
-            depth=10*km,
+            depth=depth,
             slip=pmt.magnitude_to_moment(5.0) / (
                 5*km * 5*km * store.config.earthmodel_1d.material(
-                    10*km).shear_modulus()),
+                    depth).shear_modulus()),
             width=5*km,
             length=5*km)
 
@@ -422,6 +424,54 @@ class GFSourcesTestCase(unittest.TestCase):
             rect2.get_magnitude(
                 store, gf.Target(interpolation='nearest_neighbor')))
 
+        # tensile
+        rect3 = gf.RectangularSource(
+            depth=depth,
+            magnitude=5.0,
+            width=5*km,
+            length=5*km,
+            opening_fraction=1.)
+
+        rect4 = gf.RectangularSource(
+            depth=depth,
+            slip=pmt.magnitude_to_moment(5.0) / (
+                5*km * 5*km * store.config.earthmodel_1d.material(
+                    depth).bulk()),
+            width=5*km,
+            length=5*km,
+            opening_fraction=1.)
+
+        self.assertAlmostEqual(
+            rect3.get_magnitude(),
+            rect4.get_magnitude(
+                store, gf.Target(interpolation='nearest_neighbor')))
+
+        # mixed
+        of = -0.4
+        rect5 = gf.RectangularSource(
+            depth=depth,
+            magnitude=5.0,
+            width=5*km,
+            length=5*km,
+            opening_fraction=of)
+
+        rect6 = gf.RectangularSource(
+            depth=depth,
+            slip=pmt.magnitude_to_moment(5.0) / (
+                5*km * 5*km * (
+                    store.config.earthmodel_1d.material(
+                        depth).bulk() * abs(of) +
+                    store.config.earthmodel_1d.material(
+                        depth).shear_modulus() * (1 - abs(of)))),
+            width=5*km,
+            length=5*km,
+            opening_fraction=of)
+
+        self.assertAlmostEqual(
+            rect5.get_magnitude(),
+            rect6.get_magnitude(
+                store, gf.Target(interpolation='nearest_neighbor')))
+
     def test_discretize_rect_source(self):
 
         store = self.dummy_homogeneous_store()