Source code for openquake.hazardlib.mfd.truncated_gr

# The Hazard Library
# Copyright (C) 2012-2023 GEM Foundation
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as
# published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
Module :mod:`openquake.hazardlib.mfd.truncated_gr` defines a Truncated
Gutenberg-Richter MFD.
"""
import math
import numpy as np
from openquake.baselib.python3compat import round
from openquake.hazardlib.mfd.base import BaseMFD


[docs]class TruncatedGRMFD(BaseMFD): """ Truncated Gutenberg-Richter MFD is defined in a functional form. The annual occurrence rate for a specific bin (magnitude band) is defined as :: rate = 10 ** (a_val - b_val * mag_lo) - 10 ** (a_val - b_val * mag_hi) where * ``a_val`` is the cumulative ``a`` value (``10 ** a`` is the number of earthquakes per year with magnitude greater than or equal to 0), * ``b_val`` is Gutenberg-Richter ``b`` value -- the decay rate of exponential distribution. It describes the relative size distribution of earthquakes: a higher ``b`` value indicates a relatively larger proportion of small events and vice versa. * ``mag_lo`` and ``mag_hi`` are lower and upper magnitudes of a specific bin respectively. :param min_mag: The lowest possible magnitude for this MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` will be aligned to make its left border match this value. :param max_mag: The highest possible magnitude. The same as for ``min_mag``: the last bin in the histogram will correspond to the magnitude value equal to ``max_mag - bin_width / 2``. :param bin_width: A positive float value -- the width of a single histogram bin. Values for ``min_mag`` and ``max_mag`` don't have to be aligned with respect to ``bin_width``. They get rounded accordingly anyway so that both are divisible by ``bin_width`` just before converting a function to a histogram. See :meth:`_get_min_mag_and_num_bins`. """ MODIFICATIONS = {'increment_max_mag', 'set_max_mag', 'increment_b', 'set_ab', 'set_bGR', 'increment_max_mag_no_mo_balance'} def __init__(self, min_mag, max_mag, bin_width, a_val, b_val): self.min_mag = min_mag self.max_mag = max_mag self.bin_width = bin_width self.a_val = a_val self.b_val = b_val self.check_constraints()
[docs] def check_constraints(self): """ Checks the following constraints: * Bin width is greater than 0. * Minimum magnitude is positive. * Maximum magnitude is greater than minimum magnitude by at least one bin width (or equal to that value). * ``b`` value is positive. """ if not self.bin_width > 0: raise ValueError('bin width %g must be positive' % self.bin_width) if not self.min_mag >= 0: raise ValueError('minimum magnitude %g must be non-negative' % self.min_mag) if not self.max_mag >= self.min_mag + self.bin_width: raise ValueError('maximum magnitude %g must be higher than ' 'minimum magnitude %g by ' 'bin width %g at least' % (self.max_mag, self.min_mag, self.bin_width)) if self.b_val <= 0: raise ValueError('b-value %g must be non-negative' % self.b_val) if not np.isfinite(self.a_val): raise ValueError(self.a_val)
def _get_rate(self, mag): """ Calculate and return an annual occurrence rate for a specific bin. :param mag: Magnitude value corresponding to the center of the bin of interest. :returns: Float number, the annual occurrence rate calculated using formula described in :class:`TruncatedGRMFD`. """ mag_lo = mag - self.bin_width / 2.0 mag_hi = mag + self.bin_width / 2.0 return (10 ** (self.a_val - self.b_val * mag_lo) - 10 ** (self.a_val - self.b_val * mag_hi)) def _get_min_mag_and_num_bins(self): """ Estimate the number of bins in the histogram and return it along with the first bin center abscissa (magnitude) value. Rounds ``min_mag`` and ``max_mag`` with respect to ``bin_width`` to make the distance between them include integer number of bins. :returns: A tuple of two items: first bin center and total number of bins. """ min_mag = round(self.min_mag / self.bin_width) * self.bin_width max_mag = round(self.max_mag / self.bin_width) * self.bin_width if min_mag != max_mag: min_mag += self.bin_width / 2.0 max_mag -= self.bin_width / 2.0 # here we use math round on the result of division and not just # cast it to integer because for some magnitude values that can't # be represented as an IEEE 754 double precisely the result can # look like 7.999999999999 which would become 7 instead of 8 # being naively casted to int so we would lose the last bin. num_bins = int(round((max_mag - min_mag) / self.bin_width)) + 1 return min_mag, num_bins
[docs] def get_min_max_mag(self): """ Return the minum magnitude """ min_mag, num_bins = self._get_min_mag_and_num_bins() return min_mag, min_mag + self.bin_width * (num_bins - 1)
[docs] def get_annual_occurrence_rates(self): """ Calculate and return the annual occurrence rates histogram. The result histogram has only one bin if minimum and maximum magnitude values appear equal after rounding. :returns: See :meth: `openquake.hazardlib.mfd.base.BaseMFD.get_annual_occurrence_rates`. """ mag, num_bins = self._get_min_mag_and_num_bins() rates = [] for i in range(num_bins): rate = self._get_rate(mag) rates.append((mag, rate)) mag += self.bin_width return rates
def _get_total_moment_rate(self): """ Calculate total moment rate (total energy released per unit time) :: TMR = ((10**ai) / bi) * (10 ** (bi*max_mag) - 10 ** (bi*min_mag)) where ``ai = a + log10(b) + 9.05`` and ``bi = 1.5 - b``. In case of ``bi == 0`` the following formula is applied:: TMR = (10 ** ai) * (max_mag - min_mag) :returns: Float, calculated TMR value in ``N * m / year`` (Newton-meter per year). """ ai = 9.05 + self.a_val + math.log10(self.b_val) bi = 1.5 - self.b_val if bi == 0.0: return (10 ** ai) * (self.max_mag - self.min_mag) else: return (((10 ** ai) / bi) * (10 ** (bi * self.max_mag) - 10 ** (bi * self.min_mag))) def _set_a(self, tmr): """ Recalculate an ``a`` value preserving a total moment rate ``tmr`` :: a = (log10((tmr * bi) / (10 ** (bi*max_mag) - 10 ** (bi*min_mag))) - 9.05 - log10(b)) where ``bi = 1.5 - b``. If ``bi == 0`` the following formula is used: a = log10(tmr / (max_mag - min_mag)) - 9.05 - log10(b) """ bi = 1.5 - self.b_val if bi == 0.0: self.a_val = (math.log10(tmr / (self.max_mag - self.min_mag)) - 9.05 - math.log10(self.b_val)) else: self.a_val = (math.log10(tmr * bi / (10 ** (bi * self.max_mag) - 10 ** (bi * self.min_mag))) - 9.05 - math.log10(self.b_val))
[docs] def modify_increment_max_mag_no_mo_balance(self, value): """ Apply relative maximum magnitude modification. :param value: A float value to add to ``max_mag``. """ self.max_mag += value
[docs] def modify_increment_max_mag(self, value): """ Apply relative maximum magnitude modification. :param value: A float value to add to ``max_mag``. The Gutenberg-Richter ``a`` value is :meth:`recalculated <_set_a>` with respect to old :meth:`total moment rate <_get_total_moment_rate>`. """ tmr = self._get_total_moment_rate() self.max_mag += value # need to check constraints here because _set_a() would die # if new max_mag <= min_mag. self.check_constraints() self._set_a(tmr)
[docs] def modify_set_max_mag(self, value): """ Apply absolute maximum magnitude modification. :param value: A float value to assign to ``max_mag``. No specific recalculation of other Gutenberg-Richter parameters is done after assigning a new value to ``max_mag``. """ self.max_mag = value
[docs] def modify_set_bGR(self, b_val: float): """ Updates the b_value of the GR relationship. :param b_val: The value of the new maximum magnitude """ self.b_val = b_val
[docs] def modify_increment_b(self, value): """ Apply relative ``b``-value modification. :param value: A float value to add to ``b_val``. After changing ``b_val`` the ``a_val`` is recalculated the same way as for :meth:`modify_increment_max_mag` (with respect to TMR). """ tmr = self._get_total_moment_rate() self.b_val += value self.check_constraints() self._set_a(tmr)
[docs] def modify_set_ab(self, a_val, b_val): """ Apply absolute ``a`` and ``b`` values modification. :param a_val: A float value to use as a new ``a_val``. :param b_val: A float value to use as a new ``b_val``. No recalculation of other Gutenberg-Richter parameters is done. """ self.b_val = b_val self.a_val = a_val
[docs] @classmethod def from_moment(cls, min_mag, max_mag, bin_width, b_val, moment_rate, d_val=9.1): """ :param min_mag: The lowest possible magnitude for this MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` will be aligned to make its left border match this value. :param max_mag: The highest possible magnitude. The same as for ``min_mag``: the last bin in the histogram will correspond to the magnitude value equal to ``max_mag - bin_width / 2``. :param bin_width: A positive float value -- the width of a single histogram bin. :param b_val: The slope of the GR relationship :param moment_rate: The value of scalar seismic moment per year released by this MFD. Unit of measure is N ・m :param d_val: The constant term of the equation used to compute moment from magnitude. Set to 9.1 for backcompatibility. """ c_val = 1.5 tmp = 0 mou = 10**(c_val * max_mag + d_val) beta = b_val * np.log(10.0) term2 = np.exp(-beta*(max_mag - tmp)) rate = (moment_rate * (c_val - b_val) * (1 - term2) / (b_val * mou * term2)) a_val = np.log10(rate) self = cls(min_mag, max_mag, bin_width, a_val, b_val) return self
[docs] @classmethod def from_slip_rate(cls, min_mag, max_mag, bin_width, b_val, slip_rate, rigidity, area, constant_term=9.1): """ Calls .from_moment with moment = slip_rate * rigidity * area :param min_mag: The lowest possible magnitude for this MFD. The first bin in the :meth:`result histogram <get_annual_occurrence_rates>` will be aligned to make its left border match this value. :param max_mag: The highest possible magnitude. The same as for ``min_mag``: the last bin in the histogram will correspond to the magnitude value equal to ``max_mag - bin_width / 2``. :param bin_width: A positive float value -- the width of a single histogram bin. :param b_val: The slope of the GR relationship :param slip_rate: A float defining the slip rate [mm/yr] :param rigidity: A float defining the rigidity [GPa] :param area: A float defining the area of the fault surface [km^2] :param constant_term: A float defining the constant term of the equation used to compute the log M0 from magnitude. """ mm = 1E-3 # conversion millimiters -> meters moment_rate = (slip_rate * mm) * (rigidity * 1e9) * (area * 1e6) self = cls.from_moment(min_mag, max_mag, bin_width, b_val, moment_rate, constant_term) self.slip_rate = slip_rate self.rigidity = rigidity return self