psgrn_pscmp: upd timing and gf_spacing

doc/source/apps/fomosto/backends.rst

@@ -97,6 +97,55 @@ The stresses, tilts and geoidal changes can still be computed but are (not yet)
 supported by a Fomosto Green's function store and the respective stacking
 functions.
 
+Viscoelastic media are defined by Burger's parameters, these are columns 
+7, 8 and 9 in the stores' ``earthmodel_1d``. See this example with a viscoelastic
+mantle rheology.
+
+.. code-block :: yaml
+
+    earthmodel_1d: |2
+          0.             2.5            1.2            2.1           50.            50.
+          1.             2.5            1.2            2.1           50.            50.
+          1.             6.2            3.6            2.8          600.           400.
+         17.             6.2            3.6            2.8          600.           400.
+         17.             6.6            3.7            2.9         1432.           600.
+         32.             6.6            3.7            2.9         1432.           600.
+         32.             7.3            4.             3.1         1499.           600.
+         41.             7.3            4.             3.1         1499.           600.
+      mantle
+         41.             8.2            4.7            3.4         1370.           600.            5.000E+17      1.000E+19      1.
+         91.             8.2            4.7            3.4         1370.           600.            5.000E+17      1.000E+19      1.
+
+
+The time span which is covered in the static displacement series is defined in
+``extra/psgrn_pscmp`` file (default is only coseismic displacement).
+Change `tmin_days` and `tmax_days` for the time of interest.
+The sampling rate is taken from the ``config`` file's definition of `sampling_rate`.
+
+.. note ::
+
+    Static stores define the sampling rate in Hz.
+    ``sampling_rate: 1.157e-05 Hz`` is a sampling rate of one day!
+
+.. code-block :: yaml
+
+    --- !pf.PsGrnPsCmpConfig
+    tmin_days: 0.0
+    tmax_days: 1000.0
+    gf_outdir: psgrn_functions
+    psgrn_config: !pf.PsGrnConfig
+      version: 2008a
+      sampling_interval: 1.0
+      gf_depth_spacing: -1.0
+      gf_distance_spacing: -1.0
+      observation_depth: 0.0
+    pscmp_config: !pf.PsCmpConfig
+      version: 2008a
+      observation: !pf.PsCmpScatter {}
+      rectangular_fault_size_factor: 1.0
+      rectangular_source_patches: []
+
+
 The current version of PSGRN/PSCMP is ``2008a`` (at the time of writing,
 2017-02-14), and can be downloaded from
 https://git.pyrocko.org/pyrocko/fomosto-psgrn-pscmp .

doc/source/library/examples/gf_forward.rst

@@ -104,7 +104,12 @@ In this advanced example we leverage the viscoelastic forward modelling capabili
 
 Viscoelastic static GF store forward-modeling the transient effects of a deep dislocation source, mimicking a transform plate boundary. Together with a shallow seismic source. The cross denotes the tracked pixel location. (Top) Displacement of the tracked pixel in time.
 
-The static store has to be setup with Burger material describing the viscoelastic properties of the medium, see this `config` for the fomosto store:
+The static store has to be setup with Burger material describing the viscoelastic properties of the medium, see this ``config`` for the fomosto store:
+
+.. note ::
+
+    Static stores define the sampling rate in Hz.
+    ``sampling_rate: 1.157e-06 Hz`` is a sampling rate of 10 days!
 
 .. code-block:: yaml
 
@@ -138,29 +143,25 @@ The static store has to be setup with Burger material describing the viscoelasti
     distance_delta: 1000.0
 
 
-In the `extra/psgrn_pscmp` configruation file we have to define the number and spacing of timesteps. Here we use 64 timesteps within 600 days after time of the dislocation for `psgrn` and 600 days with a deltat of 10 days for `pscmp`.
+In the ``extra/psgrn_pscmp`` configruation file we have to define the timespan from `tmin_days` to `tmax_days`, covered by the `sampling_rate` (see above)
 
 .. code-block:: yaml
 
     --- !pf.PsGrnPsCmpConfig
+    tmin_days: 0.0
+    tmax_days: 600.0
+    gf_outdir: psgrn_functions
     psgrn_config: !pf.PsGrnConfig
       version: 2008a
       sampling_interval: 1.0
-      n_time_samples: 64
-      max_time: 600
-      gf_depth_spacing: 1.0
-      gf_distance_spacing: 10.0
+      gf_depth_spacing: -1.0
+      gf_distance_spacing: -1.0
       observation_depth: 0.0
     pscmp_config: !pf.PsCmpConfig
       version: 2008a
       observation: !pf.PsCmpScatter {}
       rectangular_fault_size_factor: 1.0
-      snapshots: !pf.PsCmpSnapshots
-        tmin: 0.0
-        tmax: 600
-        deltatdays: 10
       rectangular_source_patches: []
-    gf_outdir: psgrn_functions
 
 
 Download :download:`gf_forward_viscoelastic.py </../../examples/gf_forward_viscoelastic.py>`

src/fomosto/psgrn_pscmp.py

@@ -5,6 +5,7 @@
 from __future__ import absolute_import, division
 
 import logging
+import warnings
 import os
 from os.path import join as pjoin
 import signal
@@ -35,6 +36,7 @@ psgrn_gravity_names = ('gd', 'gr')
 psgrn_components = 'ep ss ds cl'.split()
 
 km = 1000.
+day = 3600. * 24
 
 guts_prefix = 'pf'
 logger = logging.getLogger('pyrocko.fomosto.psgrn_pscmp')
@@ -98,20 +100,23 @@ class PsGrnConfig(Object):
         help='exponential sampling 1.- equidistant, > 1. decreasing sampling'
              ' with distance')
     n_time_samples = Int.T(
-        default=1,
+        optional=True,
         help='Number of temporal GF samples up to max_time. Has to be equal'
              ' to a power of 2! If not, next power of 2 is taken.')
     max_time = Float.T(
-        default=1.,
+        optional=True,
         help='Maximum time [days] after seismic event.')
 
     gf_depth_spacing = Float.T(
-        default=1.,
-        help='spatial depth spacing [km] for the medium properties.')
+        default=-1.,
+        help='depth spacing [km] for the observation points. '
+             '-1 is using the distance spacing from the store\'s config')
     gf_distance_spacing = Float.T(
-        default=10.,
-        help='spatial distance spacing [km] for the medium properties.')
-    observation_depth = Float.T(default=0.)
+        default=-1.,
+        help='spatial spacing [km] for the observation points. '
+             '-1 is using the distance spacing from the store\'s config')
+    observation_depth = Float.T(
+        default=0.)
 
     def items(self):
         return dict(self.T.inamevals(self))
@@ -760,7 +765,8 @@ class PsCmpConfig(Object):
              ' spacing of the GF database. This factor is applied to the'
              ' relationship in i.e. fault length & width = f * dx.')
 
-    snapshots = PsCmpSnapshots.T(PsCmpSnapshots.D())
+    snapshots = PsCmpSnapshots.T(
+        optional=True)
 
     rectangular_source_patches = List.T(PsCmpRectangularSource.T())
 
@@ -1113,11 +1119,17 @@ class PsGrnPsCmpConfig(Object):
     '''
     Combined config Object of PsGrn and PsCmp.
     '''
+    tmin_days = Float.T(
+        default=0.,
+        help='Min. time in days')
+    tmax_days = Float.T(
+        default=1.,
+        help='Max. time in days')
+    gf_outdir = String.T(default='psgrn_functions')
+
     psgrn_config = PsGrnConfig.T(PsGrnConfig.D())
     pscmp_config = PsCmpConfig.T(PsCmpConfig.D())
 
-    gf_outdir = String.T(default='psgrn_functions')
-
 
 class PsCmpError(gf.store.StoreError):
     pass
@@ -1354,7 +1366,30 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
         baseconf = self.store.get_extra('psgrn_pscmp')
 
         cg = PsGrnConfigFull(**baseconf.psgrn_config.items())
+        if cg.n_time_samples is None or cg.max_time is None:
+            deltat_days = 1. / storeconf.sample_rate / day
+
+            cg.n_time_samples = int(baseconf.tmax_days // deltat_days)
+            cg.max_time = baseconf.tmax_days
+        else:
+            warnings.warn(
+                'PsGrnConfig is defining n_times_samples and max_time,'
+                ' this will be replaced by PsGrnPsCmpConfig tmin and tmax.',
+                FutureWarning)
+
         cc = PsCmpConfigFull(**baseconf.pscmp_config.items())
+        if cc.snapshots is None:
+            deltat_days = 1. / storeconf.sample_rate / day
+
+            cc.snapshots = PsCmpSnapshots(
+                tmin=baseconf.tmin_days,
+                tmax=baseconf.tmax_days,
+                deltatdays=deltat_days)
+        else:
+            warnings.warn(
+                'PsCmpConfig is defining snapshots,'
+                ' this will be replaced by PsGrnPsCmpConfig tmin and tmax.',
+                FutureWarning)
 
         cg.earthmodel_1d = storeconf.earthmodel_1d
 
@@ -1394,10 +1429,17 @@ class PsGrnCmpGFBuilder(gf.builder.Builder):
             (self.step + 1, self.nsteps, iblock + 1, self.nblocks))
 
         if self.step == 0:
+            gf_depth_spacing = cg.gf_depth_spacing * km
+            if cg.gf_depth_spacing == -1.:
+                gf_depth_spacing = fc.source_depth_delta
             n_steps_depth = int((fc.source_depth_max - fc.source_depth_min) /
-                                (cg.gf_depth_spacing * km)) + 1
+                                gf_depth_spacing) + 1
+
+            gf_distance_spacing = cg.gf_distance_spacing * km
+            if cg.gf_distance_spacing == -1.:
+                gf_distance_spacing = fc.distance_delta
             n_steps_distance = int((fc.distance_max - fc.distance_min) /
-                                   (cg.gf_distance_spacing * km)) + 1
+                                   gf_distance_spacing) + 1
 
             cg.depth_grid = PsGrnSpatialSampling(
                 n_steps=n_steps_depth,

src/gf/meta.py

@@ -1851,8 +1851,9 @@ class Config(Object):
             raise OutOfBounds()
 
     def is_static(self):
-        if self.modelling_code_id in ('psgrn_pscmp', 'poel'):
-            return True
+        for code in ('psgrn_pscmp', 'poel'):
+            if self.modelling_code_id.startswith(code):
+                return True
         return False
 
     def is_dynamic(self):

test/gf/test_gf_psgrn_pscmp.py

@@ -15,6 +15,7 @@ from ..common import Benchmark
 
 logger = logging.getLogger('pyrocko.test.test_gf_psgrn_pscmp')
 benchmark = Benchmark()
+uniform = num.random.uniform
 
 km = 1000.
 
@@ -23,6 +24,9 @@ d2r = 1.0 / r2d
 km = 1000.
 mm = 1e-3
 
+neast = 40
+nnorth = 40
+
 
 def statics(engine, source, starget):
     store = engine.get_store(starget.store_id)
@@ -74,7 +78,7 @@ class GFPsgrnPscmpTestCase(unittest.TestCase):
 
         return self.pscmp_store_dir, self.psgrn_config
 
-    def _create_psgrn_pscmp_store(self):
+    def _create_psgrn_pscmp_store(self, extra_config=None):
 
         mod = cake.LayeredModel.from_scanlines(cake.read_nd_model_str('''
  0. 5.8 3.46 2.6 1264. 600.
@@ -99,10 +103,13 @@ mantle
         self.tempdirs.append(store_dir)
         store_id = 'psgrn_pscmp_test'
 
-        version = '2008a'
+        if not extra_config:
+            c = psgrn_pscmp.PsGrnPsCmpConfig()
+            c.psgrn_config.sampling_interval = 1.
+        else:
+            c = extra_config
 
-        c = psgrn_pscmp.PsGrnPsCmpConfig()
-        c.psgrn_config.sampling_interval = 1.
+        version = '2008a'
         c.psgrn_config.version = version
         c.pscmp_config.version = version
 
@@ -112,10 +119,10 @@ mantle
             sample_rate=1. / (3600. * 24.),
             receiver_depth=0. * km,
             source_depth_min=0. * km,
-            source_depth_max=5. * km,
+            source_depth_max=15. * km,
             source_depth_delta=0.1 * km,
             distance_min=0. * km,
-            distance_max=40. * km,
+            distance_max=60. * km,
             distance_delta=0.1 * km,
             modelling_code_id='psgrn_pscmp.%s' % version,
             earthmodel_1d=mod,
@@ -142,7 +149,6 @@ mantle
         return store_dir, c
 
     def test_fomosto_vs_psgrn_pscmp(self):
-
         store_dir, c = self.get_pscmp_store_info()
 
         origin = gf.Location(
@@ -156,15 +162,14 @@ mantle
             depth=2. * km,
             width=0.2 * km,
             length=0.5 * km,
-            rake=90., dip=45., strike=45.,
-            slip=num.random.uniform(1., 5.))
+            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)
 
-        neast = 40
-        nnorth = 40
-
         N, E = num.meshgrid(num.linspace(-20. * km, 20. * km, nnorth),
                             num.linspace(-20. * km, 20. * km, neast))
 
@@ -185,8 +190,7 @@ mantle
         runner = psgrn_pscmp.PsCmpRunner(keep_tmp=False)
         runner.run(ccf)
         ps2du = runner.get_results(component='displ')[0]
-        t3 = time()
-        logger.info('pscmp stacking time %f' % (t3 - t2))
+        logger.info('pscmp stacking time %f s' % (time() - t2))
 
         un_pscmp = ps2du[:, 0]
         ue_pscmp = ps2du[:, 1]
@@ -212,8 +216,7 @@ mantle
         for source in [source_plain, source_with_time]:
             t0 = time()
             r = engine.process(source, [starget_nn, starget_ml])
-            t1 = time()
-            logger.info('pyrocko stacking time %f' % (t1 - t0))
+            logger.info('pyrocko stacking time %f' % (time() - t0))
             for static_result in r.static_results():
                 un_fomosto = static_result.result['displacement.n']
                 ue_fomosto = static_result.result['displacement.e']
@@ -223,6 +226,142 @@ mantle
                 num.testing.assert_allclose(ue_fomosto, ue_pscmp, atol=1*mm)
                 num.testing.assert_allclose(ud_fomosto, ud_pscmp, atol=1*mm)
 
+    def test_gf_distance_sampling(self):
+        origin = gf.Location(
+            lat=10.,
+            lon=-15.)
+
+        # test GF store
+        TestRF = dict(
+            lat=origin.lat,
+            lon=origin.lon,
+            depth=2. * km,
+            width=1. * km,
+            length=3. * km,
+            rake=uniform(-90., 90.),
+            dip=uniform(0., 90.),
+            strike=uniform(0., 360.),
+            slip=uniform(1., 5.))
+
+        source_plain = gf.RectangularSource(**TestRF)
+
+        N, E = num.meshgrid(num.linspace(-40. * km, 40. * km, nnorth),
+                            num.linspace(-40. * km, 40. * km, neast))
+
+        starget_ml = gf.StaticTarget(
+            lats=num.full(N.size, origin.lat),
+            lons=num.full(N.size, origin.lon),
+            north_shifts=N.flatten(),
+            east_shifts=E.flatten(),
+            interpolation='multilinear')
+
+        # Getting reference gf_distance_sampling = 10.
+
+        def get_displacements(source):
+            store_dir, c = self.get_pscmp_store_info()
+            engine = gf.LocalEngine(store_dirs=[store_dir])
+            r = engine.process(source, starget_ml)
+            ref_result = r.static_results()[0]
+
+            compare_results = {}
+            for gf_dist_spacing in (0.25, .5, 1., 2., 4., 8., 10.,):
+                extra_config = psgrn_pscmp.PsGrnPsCmpConfig()
+                extra_config.psgrn_config.gf_distance_spacing = gf_dist_spacing
+                extra_config.psgrn_config.gf_depth_spacing = .5
+
+                store_dir, c = self._create_psgrn_pscmp_store(extra_config)
+                engine = gf.LocalEngine(store_dirs=[store_dir])
+                t0 = time()
+                r = engine.process(source, starget_ml)
+                logger.info('pyrocko stacking time %f s' % (time() - t0))
+
+                static_result = r.static_results()[0]
+                compare_results[gf_dist_spacing] = static_result
+
+                num.testing.assert_allclose(
+                    ref_result.result['displacement.n'],
+                    static_result.result['displacement.n'],
+                    atol=1*mm)
+                num.testing.assert_allclose(
+                    ref_result.result['displacement.e'],
+                    static_result.result['displacement.e'],
+                    atol=1*mm)
+                num.testing.assert_allclose(
+                    ref_result.result['displacement.d'],
+                    static_result.result['displacement.d'],
+                    atol=1*mm)
+
+            return ref_result, compare_results
+
+        # ref_result, compare_results = get_displacements(source_plain)
+        # self.plot_displacement_differences(ref_result, compare_results)
+
+        if True:
+            import matplotlib
+            matplotlib.use('Qt5Agg')
+            import matplotlib.pyplot as plt
+
+            fig = plt.figure()
+            ax = fig.gca()
+
+            for depth in (2., 4., 8., 12.):
+                source_plain.depth = depth * km
+                ref_result, compare_results = get_displacements(source_plain)
+                self.plot_differences(
+                    ax, ref_result, compare_results,
+                    label='Source Depth %.2f km' % depth)
+
+            ax.legend()
+            plt.show()
+
+    def plot_displacement_differences(self, reference, compare):
+        import matplotlib
+        matplotlib.use('Qt5Agg')
+        import matplotlib.pyplot as plt
+
+        ncompare = len(compare)
+        fig, axes = plt.subplots(ncompare + 1, 1)
+
+        ax = axes[0]
+        ref_displ = reference.result['displacement.d'].reshape(nnorth, neast)
+        ax.imshow(ref_displ, cmap='seismic')
+        ax.set_title('Reference displacement')
+
+        for (gf_distance, comp), ax in zip(compare.items(), axes[1:]):
+            displ = comp.result['displacement.d'].reshape(nnorth, neast)
+            displ -= ref_displ
+            ax.imshow(displ, cmap='seismic')
+            ax.set_title('GF Distance %.2f' % gf_distance)
+
+        ax.legend()
+        plt.show()
+
+    def plot_differences(self, ax, reference, compare, **kwargs):
+        from matplotlib.ticker import FuncFormatter
+
+        for key, comp in compare.items():
+            print(key, comp.result['displacement.n'].shape)
+
+        ref_displ = num.linalg.norm(tuple(reference.result.values()), axis=0)
+        print(ref_displ.shape)
+
+        differences = num.array([
+            num.linalg.norm(tuple(comp.result.values()), axis=0)
+            for comp in compare.values()])
+
+        differences = num.abs(differences)
+        differences /= num.abs(ref_displ)
+        differences -= 1.
+
+        ax.plot(tuple(compare.keys()), differences.max(axis=1)*100., **kwargs)
+        ax.set_ylabel('% Difference')
+        ax.set_xlabel('gf_distance_spacing [km]')
+        # ax.xaxis.set_major_formatter(
+        #     FuncFormatter(lambda x, v: '%.1f%' % v))
+        # ax.set_xticks(tuple(compare.keys()))
+        ax.grid(alpha=.3)
+        return ax
+
 
 if __name__ == '__main__':
     util.setup_logging('test_gf_psgrn_pscmp', 'warning')

test/gf/test_gf_static.py

@@ -225,6 +225,46 @@ mantle
         # self.plot_static_los_result(ml)
         # print benchmark
 
+    def test_pseudo_dyn_performance(self):
+        engine = gf.LocalEngine(store_dirs=[self.get_store_dir('pscmp')])
+        store = engine.get_store('psgrn_pscmp_test')
+        ntargets = 250
+        interpolation = 'nearest_neighbor'
+
+        def calc_dyn_rupt(nx=4, ny=4):
+            from pyrocko.plot.dynamic_rupture import RuptureMap
+            dyn_rupture = gf.PseudoDynamicRupture(
+                nx=nx, ny=ny,
+                tractions=(1.e4, 0.4e4, 0.1),
+                north_shift=2*km,
+                east_shift=2*km,
+                depth=6.5*km,
+                width=10.*km,
+                length=40*km,
+                dip=random.uniform(0., 90.),
+                strike=random.uniform(-180., 180.),
+                magnitude=7.,
+                anchor='top',
+                decimation_factor=4)
+
+            static_target = gf.StaticTarget(
+                north_shifts=(random.rand(ntargets)-.5) * 25. * km,
+                east_shifts=(random.rand(ntargets)-.5) * 25. * km,
+                tsnapshot=20,
+                interpolation=interpolation)
+
+            t = time.time()
+            # dyn_rupture.discretize_patches(store)
+            engine.process(dyn_rupture, static_target)
+            map = RuptureMap(source=dyn_rupture, lat=0., lon=0., radius=40*km)
+            map.draw_patch_parameter('traction')
+            map.save('/tmp/test.pdf')
+            return dyn_rupture.nx*dyn_rupture.ny, time.time() - t
+
+        for n in (10, 20, 30):
+            npatches, t = calc_dyn_rupt(n, n)
+            
+
     @staticmethod
     def plot_static_los_result(result):
         import matplotlib.pyplot as plt