Skip to content

Commit

Permalink
Merge pull request #260 from rmodrak/master
Browse files Browse the repository at this point in the history
Miscellaneous improvements
  • Loading branch information
rmodrak authored Mar 27, 2024
2 parents 578f0ce + dbb2fc1 commit 5ebdd79
Show file tree
Hide file tree
Showing 16 changed files with 271 additions and 166 deletions.
2 changes: 1 addition & 1 deletion docs/user_guide/03.rst
Original file line number Diff line number Diff line change
Expand Up @@ -94,7 +94,7 @@ To generate a full set of Green's functions for a given hypocenter and depth, si
.. code ::
{event_id}/
{depth_in_m}/
{depth_in_km}/
{net}.{sta}.{loc}.Z.Mrr.sac
{net}.{sta}.{loc}.Z.Mtt.sac
{net}.{sta}.{loc}.Z.Mpp.sac
Expand Down
61 changes: 37 additions & 24 deletions docs/user_guide/03/source_side.rst
Original file line number Diff line number Diff line change
Expand Up @@ -4,35 +4,40 @@
This page is still under construction. To help improve the
documentation, feel free to submit a pull request.


Source-side 3D Green's functions
================================

Generating source-side 3D Green's functions using SPECFEM3D/3D_GLOBE
--------------------------------------------------------------------

In principle, any 3D solver can be used to generate Green's functions, as long as the `requirements <https://uafgeotools.github.io/mtuq/user_guide/03/source_side.html#requirements-for-mtuq-source-side-green-s-functions>`_ below are satisfied.
In principle, any 3D solver can be used to generate source-side Green's functions, as long as the `requirements <https://uafgeotools.github.io/mtuq/user_guide/03/source_side.html#requirements-for-mtuq-source-side-green-s-functions>`_ below are satisfied.

So far, however, the source-side Green's function machinery has been tested using only SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary.
So far, however, the machinery has only been tested using SPECFEM3D/3D_GLOBE. To convert SPECFEM3D/3D_GLOBE output to MTUQ-compliant Green's functions, the following steps are necessary.


**Generate SAC binary files**
Generate SAC binary files
.........................

SAC binary output format is natively supported by SPECFEM3D_GLOBE (in the parameter file, set `OUTPUT_SEISMOS_SAC_BINARY = .true.`).

Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it becomes necessary to manually convert SPECFEM3D output to SAC binary format.
Unfortunately, SAC binary output format is not natively supported by SPECFEM3D, so it is necessary to manually convert SPECFEM3D output.


**Convert to SI units**
Convert to SI units
...................

MTUQ uses the fully SI convention described below. In contrast, SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions.
SPECFEM3D/3D_GLOBE use a mixed SI and CGS convention, in which moment tensor elements are input in terms of dynes and centimeters, and seismograms are output in meters. In contrast, MTUQ uses a fully SI :ref:`convention<Units convention>`. As a result, it is necessary to scale SPECFEM3D/3D_GLOBE seismograms by 10^7 prior to using them as MTUQ Green's functions.


**Additional amplitude scaling**
Other amplitude scaling
.......................

In addition to converting to SI units, it is also necessary to account for any scaling factors in the SPECFEM3D/3D_GLOBE input files. Such scaling factors can enter, for example, through the `M_rr,M_tt,M_pp,M_rt,M_rp,M_tp` values in the moment tensor input file or through the `scaling factor <https://github.com/SPECFEM/specfem3d/blob/bf45798f3af9d792326a829de920fd944cf7c7dd/EXAMPLES/applications/homogeneous_halfspace_HEX27_elastic_no_absorbing/DATA/FORCESOLUTION#L8>`_ in the force input file.


**Rotate to vertical, radial, transverse components**
Rotate traces
.............

Conveniently, SPECFEM3D_GLOBE can be made to automatically rotate output seismograms into vertical, radial and transverse components (set `ROTATE_SEISMOGRAMS_RT = .true.` in the parameter file).

Expand All @@ -44,15 +49,16 @@ No modifications are necessary on account of moment tensor basis convention, sin
Requirements for MTUQ source-side 3D Green's functions
------------------------------------------------------

**File format**
File format
...........

Individual Green's functions must be written to SAC binary files.

A total 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion).

A total of 18 SAC binary files are required to represent the response between a given hypocenter and station (corresponding to 6 moment tensor elements times 3 directions of motion).


**Units convention**
Units convention
................

For a moment tensor inversion, each SAC binary file must give the response in meters to a 1 Newton-meter force couple.

Expand All @@ -62,9 +68,10 @@ In both cases, MTUQ uses a fully SI units convention (compare with SPECFEM3D/3D_



**Basis convention**
Basis convention
................

MTUQ uses an `up-south-east` basis convention in which `r` denotes up, `t` denotes south, and `p` denotes east.
MTUQ uses an `up-south-east` basis convention in which `r` denotes vertically upward, `t` denotes south, and `p` denotes east.

Green's functions must be rotated into into vertical `Z`, radial `R` and transverse `T` components relative to the source-receiver backazimuth.

Expand Down Expand Up @@ -96,23 +103,29 @@ Place all seismograms for the same hypocenter in a single directory as follows:
The corresponding convention for force responses is:

.. warning::

Not yet tested; subject to change without notice.


.. code ::
{event_id}/
{depth_in_km}/
{net}.{sta}.{loc}.Z.Fr.sac
{net}.{sta}.{loc}.Z.Ft.sac
{net}.{sta}.{loc}.Z.Fp.sac
{net}.{sta}.{loc}.R.Fr.sac
{net}.{sta}.{loc}.R.Ft.sac
{net}.{sta}.{loc}.R.Fp.sac
{net}.{sta}.{loc}.T.Fr.sac
{net}.{sta}.{loc}.T.Fp.sac
{net}.{sta}.{loc}.T.Fp.sac
{net}.{sta}.{loc}.Z.Fz.sac
{net}.{sta}.{loc}.Z.Fs.sac
{net}.{sta}.{loc}.Z.Fe.sac
{net}.{sta}.{loc}.R.Fz.sac
{net}.{sta}.{loc}.R.Fs.sac
{net}.{sta}.{loc}.R.Fe.sac
{net}.{sta}.{loc}.T.Fz.sac
{net}.{sta}.{loc}.T.Fs.sac
{net}.{sta}.{loc}.T.Fe.sac
...
**Origin time convention**
Origin time convention
......................

For origin time, MTUQ uses a centroid convention (`more details <https://github.com/uafgeotools/mtuq/issues/140>`_), so that `t=0` in the `GreensTensor` time discretization corresponds to mean source excitation time.

Expand Down
28 changes: 8 additions & 20 deletions examples/GridSearch.Force.py
Original file line number Diff line number Diff line change
Expand Up @@ -18,19 +18,7 @@

if __name__=='__main__':
#
# Carries out grid search over 64,000 double couple moment tensors
#
# USAGE
# mpirun -n <NPROC> python GridSearch.DoubleCouple.py
#
# For a simpler example, see SerialGridSearch.DoubleCouple.py,
# which runs the same inversion in serial
#


#
# We will investigate the source process of an Mw~4 earthquake using data
# from a regional seismic array
# Carries out grid search over force orientation and magnitude
#

path_data= fullpath('data/examples/20210809074550/*[ZRT].sac')
Expand All @@ -40,7 +28,7 @@


#
# We are only using surface waves in this example. Check out the DC or FMT examples for multi-mode inversions.
# We are only using surface waves in this example. Check out the DC or FMT examples for multi-mode inversions
#

process_sw = ProcessData(
Expand Down Expand Up @@ -77,7 +65,7 @@


#
# Next, we specify the moment tensor grid and source-time function
# Next, we specify the force grid and source-time function
#

grid = ForceGridRegular(magnitudes_in_N=10**np.arange(9,12.1,0.1), npts_per_axis=90)
Expand Down Expand Up @@ -189,17 +177,17 @@

print('Generating figures...\n')

plot_data_greens1(event_id+'FMT_waveforms1.png',
plot_data_greens1(event_id+'_force_waveforms.png',
data_sw, greens_sw, process_sw, misfit_sw, stations, origin, best_force, direction_dict)

plot_misfit_force(event_id+'DC_misfit_force.png', results)
plot_misfit_force(event_id+'_force_misfit.png', results)

print('Saving results...\n')

# save best-fitting source
save_json(event_id+'Force_solution.json', merged_dict)
save_json(event_id+'_force_solution.json', merged_dict)

# save misfit surface
results.save(event_id+'DC_misfit.nc')
results.save(event_id+'_force_misfit.nc')

print('\nFinished\n')
print('\nFinished\n')
Loading

0 comments on commit 5ebdd79

Please sign in to comment.