|
# https://pyrocko.org - GPLv3 # # The Pyrocko Developers, 21st Century # ---|P------/S----------~Lg----------
''' Base class for analytical source models. '''
optional=True, default='')
default=0., help='Source origin time', optional=True)
default=0., help='Rupture velocity', optional=True)
def northing(self):
def easting(self):
''' Rectangular analytical source model.
Coordinates on the source plane are with respect to the origin point given by `(lat, lon, east_shift, north_shift, depth)`. '''
default=0.0, help='Strike direction in [deg], measured clockwise from north')
default=90.0, help='Dip angle in [deg], measured downward from horizontal')
default=0.0, help='Rake angle in [deg], measured counter-clockwise from ' 'right-horizontal in on-plane view')
default=0., help='Left edge source plane coordinate [m]')
default=0., help='Right edge source plane coordinate [m]')
default=0., help='Lower edge source plane coordinate [m]')
default=0., help='Upper edge source plane coordinate [m]')
default=0., help='Slip on the rectangular source area [m]', optional=True)
def length(self):
def width(self):
def area(self):
''' Rectangular Okada source model. '''
default=0., help='Opening of the plane in [m]', optional=True)
default=0.25, help='Poisson\'s ratio, default 0.25', optional=True)
help='First Lame\' s parameter [Pa]', optional=True)
default=32.0e9, help='Shear modulus along the plane [Pa]', optional=True)
def poisson(self): ''' Calculation of Poisson ratio (if not given).
The Poisson ratio :math:`\\nu` can be calculated from the Lame parameters :math`\\lambda` and :math:`\\mu` using :math:`\\nu = \\frac{\\lambda}{2(\\lambda + \\mu)}` (e.g. Mueller 2007). '''
raise ValueError('Shearmod and lambda are needed')
def poisson(self, poisson):
def lamb(self): ''' Calculation of first Lame parameter (if not given).
Poisson ratio and shear modulus must be available. '''
raise ValueError('Shearmod and poisson ratio are needed')
2. * self.poisson__ * self.shearmod__) / (1. - 2. * self.poisson__)
def lamb(self, lamb):
def shearmod(self): ''' Calculation of shear modulus (if not given).
Poisson ratio must be available.
.. important ::
We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`
Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`
'''
if self.poisson__ is None: raise ValueError('Poisson ratio is needed')
return (8. * (1. + self.poisson__)) / (1. - 2. * self.poisson__)
def shearmod(self, shearmod):
def seismic_moment(self): ''' Scalar Seismic moment.
Code copied from Kite Disregarding the opening (as for now) We assume :math:`M_0 = mu A D`
.. important ::
We assume a perfect elastic solid with :math:`K=\\frac{5}{3}\\mu`
Through :math:`\\mu = \\frac{3K(1-2\\nu)}{2(1+\\nu)}` this leads to :math:`\\mu = \\frac{8(1+\\nu)}{1-2\\nu}`
:return: Seismic moment release :rtype: float '''
def moment_magnitude(self): ''' Moment magnitude from Seismic moment.
We assume :math:`M_\\mathrm{w} = {\\frac{2}{3}}\\log_{10}(M_0) - 10.7`
:returns: Moment magnitude :rtype: float '''
''' Build source information array for okada_ext.okada input.
:return: array of the source data as input for okada_ext.okada :rtype: :py:class:`numpy.ndarray`, ``(1, 9)`` ''' self.northing, self.easting, self.depth, self.strike, self.dip, self.al1, self.al2, self.aw1, self.aw2])
''' Build source dislocation for okada_ext.okada input.
:return: array of the source dislocation data as input for okada_ext.okada :rtype: :py:class:`numpy.ndarray`, ``(1, 3)`` ''' num.cos(self.rake * d2r) * self.slip, num.sin(self.rake * d2r) * self.slip, self.opening])
''' Discretize fault into rectilinear grid of fault patches.
Fault orientation, slip and elastic parameters are passed to the sub-faults unchanged.
:param nlength: Number of patches in strike direction :type nlength: int :param nwidth: Number of patches in down-dip direction :type nwidth: int :return: Discrete fault patches :rtype: list of :py:class:`~pyrocko.modelling.okada.OkadaPatch` objects '''
mt.euler_to_matrix(self.dip*d2r, self.strike*d2r, 0.))
prop: getattr(self, prop) for prop in self.T.propnames if prop not in [ 'north_shift', 'east_shift', 'depth', 'al1', 'al2', 'aw1', 'aw2']}
[OkadaPatch( parent=self, ix=src_point[0], iy=src_point[1], north_shift=coord[0], east_shift=coord[1], depth=coord[2], al1=al1, al2=al2, aw1=aw1, aw2=aw2, **kwargs) for src_point, coord in zip(source_points, source_points_rot)], source_points)
''' Okada source with additional 2D indexes for bookkeeping. '''
source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1, variant='normal'):
''' Build coefficient matrix for given fault patches.
The BEM for a discretized fault and the determination of the slip distribution from stress drop is based on the relation :math:`stress = coefmat \\cdot displ`. Here the coefficient matrix is built, based on the displacements from Okada's solution and their partial derivatives.
:param source_patches_list: Source patches, to be used in BEM. :type source_patches_list: list of :py:class:`~pyrocko.modelling.okada.OkadaSource` objects :param pure_shear: If ``True``, only shear forces are taken into account. :type pure_shear: optional, bool :param rotate_sdn: If ``True``, rotate to strike, dip, normal. :type rotate_sdn: optional, bool :param nthreads: Number of threads. :type nthreads: optional, int
:return: Coefficient matrix for all source combinations. :rtype: :py:class:`numpy.ndarray`, ``(len(source_patches_list) * 3, len(source_patches_list) * 3)`` '''
return _make_okada_coefficient_matrix_slow( source_patches_list, pure_shear, rotate_sdn, nthreads)
src.source_patch() for src in source_patches_list])
else:
'strikeslip': { 'slip': unit_disl, 'opening': 0., 'rake': 0.}, 'dipslip': { 'slip': unit_disl, 'opening': 0., 'rake': 90.}, 'tensileslip': { 'slip': 0., 'opening': unit_disl, 'rake': 0.} }
else:
'strikeslip', 'dipslip', 'tensileslip'][:n_eq]): case['slip'] * num.cos(case['rake'] * d2r), case['slip'] * num.sin(case['rake'] * d2r), case['opening']])
source_patches, num.tile(source_disl, npoints).reshape(-1, 3), receiver_coords, lambda_mean, mu_mean, nthreads=nthreads, rotate_sdn=int(rotate_sdn), stack_sources=int(variant != 'normal'))
results[:, :, 3:] + results[:, :, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
= eps[:, :, diag_ind].sum(axis=-1)[:, :, num.newaxis]
.reshape(-1, npoints*3).T else: source[num.newaxis, :], source_disl[num.newaxis, :], receiver_coords, lambda_mean, mu_mean, nthreads=nthreads, rotate_sdn=int(rotate_sdn))
results[:, 3:] + results[:, (3, 6, 9, 4, 7, 10, 5, 8, 11)])
= num.sum(eps[:, diag_ind], axis=1)[:, num.newaxis] = kron * lambda_mean * dilatation+2. * mu_mean * eps
= stress_sdn[:, (2, 5, 8)].ravel()
source_patches_list, pure_shear=False, rotate_sdn=True, nthreads=1):
source_patches = num.array([ src.source_patch() for src in source_patches_list]) receiver_coords = source_patches[:, :3].copy()
npoints = len(source_patches_list)
if pure_shear: n_eq = 2 else: n_eq = 3
coefmat = num.zeros((npoints * 3, npoints * 3))
def ned2sdn_rotmat(strike, dip): rotmat = mt.euler_to_matrix( (dip + 180.) * d2r, strike * d2r, 0.).A return rotmat
lambda_mean = num.mean([src.lamb for src in source_patches_list]) shearmod_mean = num.mean([src.shearmod for src in source_patches_list])
unit_disl = 1. disl_cases = { 'strikeslip': { 'slip': unit_disl, 'opening': 0., 'rake': 0.}, 'dipslip': { 'slip': unit_disl, 'opening': 0., 'rake': 90.}, 'tensileslip': { 'slip': 0., 'opening': unit_disl, 'rake': 0.} } for idisl, case_type in enumerate([ 'strikeslip', 'dipslip', 'tensileslip'][:n_eq]): case = disl_cases[case_type] source_disl = num.array([ case['slip'] * num.cos(case['rake'] * d2r), case['slip'] * num.sin(case['rake'] * d2r), case['opening']])
for isource, source in enumerate(source_patches): results = okada_ext.okada( source[num.newaxis, :].copy(), source_disl[num.newaxis, :].copy(), receiver_coords, lambda_mean, shearmod_mean, nthreads=nthreads, rotate_sdn=int(rotate_sdn))
for irec in range(receiver_coords.shape[0]): eps = num.zeros((3, 3)) for m in range(3): for n in range(3): eps[m, n] = 0.5 * ( results[irec][m * 3 + n + 3] + results[irec][n * 3 + m + 3])
stress = num.zeros((3, 3)) dilatation = num.sum([eps[i, i] for i in range(3)])
for m, n in zip([0, 0, 0, 1, 1, 2], [0, 1, 2, 1, 2, 2]): if m == n: stress[m, n] = \ lambda_mean * \ dilatation + \ 2. * shearmod_mean * \ eps[m, n]
else: stress[m, n] = \ 2. * shearmod_mean * \ eps[m, n] stress[n, m] = stress[m, n]
normal = num.array([0., 0., -1.]) for isig in range(3): tension = num.sum(stress[isig, :] * normal) coefmat[irec * n_eq + isig, isource * n_eq + idisl] = \ tension / unit_disl
return coefmat
stress_field, coef_mat=None, source_list=None, pure_shear=False, epsilon=None, nthreads=1, **kwargs): ''' BEM least squares inversion to get fault dislocations given stress field.
Follows least squares inversion approach by Menke (1989) to calculate dislocations on a fault with several segments from a given stress field. The coefficient matrix connecting stresses and displacements of the fault patches can either be specified by the user (``coef_mat``) or it is calculated using the solution of Okada (1992) for a rectangular fault in a homogeneous half space (``source_list``).
:param stress_field: Stress change [Pa] for each source patch (as ``stress_field[isource, icomponent]`` where isource indexes the source patch and ``icomponent`` indexes component, ordered (strike, dip, tensile). :type stress_field: :py:class:`numpy.ndarray`, shape ``(nsources, 3)``
:param coef_mat: Coefficient matrix connecting source patch dislocations and the stress field. :type coef_mat: optional, :py:class:`numpy.ndarray`, shape ``(len(source_list) * 3, len(source_list) * 3)``
:param source_list: list of all source patches to be used for BEM. :type source_list: optional, list of :py:class:`~pyrocko.modelling.okada.OkadaSource` objects
:param epsilon: If given, values in ``coef_mat`` smaller than ``epsilon`` are set to zero. :type epsilon: optional, float
:param nthreads: Number of threads allowed. :type nthreads: int
:return: Inverted displacements as ``displacements[isource, icomponent]`` where isource indexes the source patch and ``icomponent`` indexes component, ordered (strike, dip, tensile). :rtype: :py:class:`numpy.ndarray`, ``(nsources, 3)`` '''
coef_mat = make_okada_coefficient_matrix( source_list, pure_shear=pure_shear, nthreads=nthreads, **kwargs)
coef_mat[coef_mat < epsilon] = 0.
stress_field = stress_field.ravel()
num.linalg.inv(num.dot(coef_mat_in.T, coef_mat_in)), coef_mat_in.T, stress_field[idx]]) except num.linalg.LinAlgError as e: logger.warning('Linear inversion failed!') logger.warning( 'coef_mat: %s\nstress_field: %s', coef_mat_in, stress_field[idx]) raise e
'AnalyticalSource', 'AnalyticalRectangularSource', 'OkadaSource', 'make_okada_coefficient_matrix', 'invert_fault_dislocations_bem'] |