From 1dd2c60535f8f73d33b9e7e92da543fc8a8ffeba Mon Sep 17 00:00:00 2001 From: Albert Kottke Date: Sun, 1 Mar 2015 13:00:36 -0800 Subject: [PATCH] Added function and tests for computing attenuation parameter. --- pyrvt/motions.py | 208 +++++++++++++++++++++++++++++------- pyrvt/peak_calculators.py | 1 + pyrvt/tests/test_motions.py | 166 ++++++++++++++++------------ pyrvt/tools.py | 4 +- 4 files changed, 269 insertions(+), 110 deletions(-) diff --git a/pyrvt/motions.py b/pyrvt/motions.py index cad31d4..ff9768a 100644 --- a/pyrvt/motions.py +++ b/pyrvt/motions.py @@ -20,10 +20,32 @@ """ Classes and functions used to define random vibration theory (RVT) based motions. + +References +---------- +.. [AH84] Anderson, J. G., & Hough, S. E. (1984). A model for the shape of the + Fourier amplitude spectrum of acceleration at high frequencies. Bulletin of + the Seismological Society of America, 74(5), 1969-1993. + +.. [AB11] Atkinson, G. M., & Boore, D. M. (2011). Modifications to existing + ground-motion prediction equations in light of new data. Bulletin of the + Seismological Society of America, 101(3), 1121-1135. + +.. [C03] Campbell, K. W. (2003). Prediction of strong ground motion using + the hybrid empirical method and its use in the development of ground-motion + (attenuation) relations in eastern North America. Bulletin of the + Seismological Society of America, 93(3), 1012-1033. + +.. [GV76] Gasparini, D. A., & Vanmarcke, E. H. (1976). Simulated earthquake + motions compatible with prescribed response spectra. Massachusetts + Institute of Technology, Department of Civil Engineering, Constructed + Facilities Division. + """ import numpy as np +from scipy.stats import linregress from scipy.interpolate import interp1d from . import peak_calculators @@ -31,6 +53,75 @@ DEFAULT_CALC = 'V75' +def increasing_x(x, y=None): + """Check if *x* is monotonically increasing. If not, reverse it. + + Parameters + ---------- + x : :class:`numpy.array` + X values, which are checked. + + y : :class:`numpy.array` or ``None``, default: ``None`` + Y values, which are reversed if *x* is reversed + + Returns + ------- + If *y* is not ``None``, then returns (*x*, *y*), otherwise only returns *x*. + x : :class:`numpy.array` + X values in monotonically increasing order. + + y : :class:`numpy.array` + Y values in same order as *x*. + + Raises + ------ + :class:`NotImplementedError` + If *x* is not monotonic. + """ + + diffs = np.diff(x) + if np.all(0 <= diffs): + # All increasing, do nothing + pass + elif np.all(diffs <= 0): + # All decreasing, reverse + x = x[::-1] + if y is not None: + y = y[::-1] + else: + raise NotImplementedError('Values are not regularly ordered.') + + if y is None: + return x + else: + return x, y + + +def log_spaced_values(min, max): + """Generate values with constant log-spacing. + + Values are generated with 512 points per decade. + + Parameters + ---------- + min : float + Minimum value of the range. + + max : float + Maximum value of the range. + + Returns + ------- + value : :class:`numpy.array` + Log-spaced values + + """ + lower = np.log10(min) + upper = np.log10(max) + count = np.ceil(512 * (upper - lower)) + return np.logspace(lower, upper, count) + + def compute_sdof_tf(freqs, osc_freq, osc_damping): """Compute the single-degree-of-freedom transfer function. When applied on the acceleration Fourier amplitude spectrum, it provides the @@ -61,12 +152,6 @@ def compute_sdof_tf(freqs, osc_freq, osc_damping): def compute_stress_drop(magnitude): """Compute the stress drop using [AB11]_ model. - References - ---------- - .. [AB11] Atkinson, G. M., & Boore, D. M. (2011). Modifications to existing - ground-motion prediction equations in light of new data. Bulletin of - the Seismological Society of America, 101(3), 1121-1135. - Parameters ---------- magnitude : float @@ -121,13 +206,13 @@ class RvtMotion(object): Parameters ---------- - freqs : :class:`numpy.array` + freqs : :class:`numpy.array` or ``None``, default: ``None`` Frequency array [Hz] - fourier_amps : :class:`numpy.array` + fourier_amps : :class:`numpy.array` or ``None``, default: ``None`` Absolute value of acceleration Fourier amplitudes. - duration : float + duration : float or ``None``, default: ``None`` Ground motion duration [dec]. peak_calculator : str or :class:`~.peak_calculators.Calculator`, default: ``None`` @@ -147,6 +232,10 @@ def __init__(self, freqs=None, fourier_amps=None, duration=None, self._fourier_amps = fourier_amps self._duration = duration + if self._freqs is not None: + self._freqs, self._fourier_amps = increasing_x( + self._freqs, self._fourier_amps) + if isinstance(peak_calculator, peak_calculators.Calculator): self.peak_calculator = peak_calculator else: @@ -161,7 +250,7 @@ def freqs(self): @property def fourier_amps(self): """Acceleration Fourier amplitude values.""" - return self._freqs + return self._fourier_amps @property def duration(self): @@ -222,17 +311,67 @@ def compute_peak(self, transfer_func=None, osc_freq=None, osc_freq=osc_freq, osc_damping=osc_damping, full_output=False) + def compute_attenuation(self, min_freq, max_freq=None, full=False): + """Compute the site attenuation (κ) based on a log-linear fit. + + This function computes the site attenuation defined by [AH84]_ as: + + .. math:: + a(f) = A_0 \exp(-\pi \kappa f) \\text( for ) f > f_E + + for a single Fourier amplitude spectrum + + + Parameters + ---------- + min_freq : float + minimum frequency of the fit + + max_freq : float or ``None``, default: ``None`` + maximum frequency of the fit. If ``None``, then the maximum + frequency range is used. + + full : bool, default: ``False`` + If the complete output should be returned + + Returns + ------- + If *full* is ``False``, then only *atten* is returned. Otherwise, + the tuple (*atten*, *r_value*, *freqs*, *fitted*) is return. + + atten : float + attenuation parameter + + r_sqrt : float + sqared correlation coefficient of the fit (R²). See + :function:`scipy.stats.linregress`. + + freqs : :class:`numpy.array` + selected frequencies + + fitted : :class:`numpy.array` + fitted values + + """ + max_freq = max_freq or self.freqs[-1] + mask = (min_freq <= self.freqs) & (self.freqs <= max_freq) + + slope, intercept, r_value, p_value, stderr = linregress( + self.freqs[mask], np.log(self.fourier_amps[mask])) + + atten = slope / -np.pi + + if full: + freqs = self.freqs[mask] + fitted = np.exp(intercept + slope * freqs) + return atten, r_value ** 2, freqs, fitted + else: + return atten + class SourceTheoryMotion(RvtMotion): """Single-corner source theory model with default parameters from [C03]_. - References - ---------- - .. [C03] Campbell, K. W. (2003). Prediction of strong ground motion using - the hybrid empirical method and its use in the development of - ground-motion (attenuation) relations in eastern North America. - Bulletin of the Seismological Society of America, 93(3), 1012-1033. - Parameters ---------- magnitude : float @@ -367,17 +506,21 @@ def compute_duration(self): return duration_source + duration_path - def compute_fourier_amps(self, freqs): + def compute_fourier_amps(self, freqs=None): """Compute the acceleration Fourier amplitudes for a frequency range. Parameters ---------- - freqs : numpy.array - Frequency range + freqs : :class:`numpy.array` or ``None``, default: ``None`` + Frequency range. If no frequency range is specified then + :function:`log_spaced_values(0.05, 200.)` is used. """ + if freqs is None: + self._freqs = log_spaced_values(0.05, 200.) + else: + self._freqs = increasing_x(np.asarray(freqs)) - self._freqs = np.asarray(freqs) self._duration = self.compute_duration() # Model component @@ -467,13 +610,8 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None, super(CompatibleRvtMotion, self).__init__( peak_calculator=peak_calculator) - osc_freqs = np.asarray(osc_freqs) - osc_accels_target = np.asarray(osc_accels_target) - - # Order by increasing frequency - ind = osc_freqs.argsort() - osc_freqs = osc_freqs[ind] - osc_accels_target = osc_accels_target[ind] + osc_freqs, osc_accels_target = increasing_x( + np.asarray(osc_freqs), np.asarray(osc_accels_target)) if duration: self._duration = duration @@ -485,12 +623,10 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None, osc_freqs, osc_accels_target, osc_damping) # The frequency needs to be extended to account for the fact that the - # osciallator transfer function has a width. The number of frequencies + # oscillator transfer function has a width. The number of frequencies # depends on the range of frequencies provided. - lower = np.log10(osc_freqs[0] / 2.) - upper = np.log10(2 * osc_freqs[-1]) - count = int(512 * (upper - lower)) - self._freqs = np.logspace(lower, upper, count) + self._freqs = log_spaced_values(osc_freqs[0] / 2., + 2. * osc_freqs[-1]) self._fourier_amps = np.empty_like(self._freqs) # Indices of the first and last point with the range of the provided @@ -501,7 +637,6 @@ def __init__(self, osc_freqs, osc_accels_target, duration=None, # last is extend one past the usable range to allow use of first:last # notation last = indices[-1, 0] + 1 - log_freqs = np.log(self._freqs) log_osc_freqs = np.log(osc_freqs) @@ -581,7 +716,7 @@ def _extrap(freq, freqs, fourier_amps, max_slope=None): self.iterations += 1 def _estimate_fourier_amps(self, osc_freqs, osc_accels, osc_damping): - """Compute an estimate of the FAS using the Vanmarcke methodology. + """Compute an estimate of the FAS using the [GV76]_ methodology. Parameters ---------- @@ -601,10 +736,9 @@ def _estimate_fourier_amps(self, osc_freqs, osc_accels, osc_damping): frequencies specifed by *osc_freqs*. """ - # TODO Add reference to the docstring # Compute initial value using Vanmarcke methodology. The response is - # first computed at the lowest frequency and then subsequently comptued + # first computed at the lowest frequency and then subsequently computed # at higher frequencies. peak_factor = 2.5 diff --git a/pyrvt/peak_calculators.py b/pyrvt/peak_calculators.py index 0693812..d75e4fa 100644 --- a/pyrvt/peak_calculators.py +++ b/pyrvt/peak_calculators.py @@ -852,6 +852,7 @@ def get_peak_calculator(method, calc_kwds): calculator : :class:`.Calculator` """ + calc_kwds = calc_kwds or dict() calculators = [ BooreJoyner1984, diff --git a/pyrvt/tests/test_motions.py b/pyrvt/tests/test_motions.py index ec4147d..b6a74f0 100644 --- a/pyrvt/tests/test_motions.py +++ b/pyrvt/tests/test_motions.py @@ -1,70 +1,96 @@ -#!/usr/bin/env python3 -# encoding: utf-8 - -# pyRVT: Seismological random vibration theory implemented with Python -# Copyright (C) 2013-2014 Albert R. Kottke albert.kottke@gmail.com -# -# This program is free software: you can redistribute it and/or modify -# it under the terms of the GNU 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 General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with this program. If not, see . - -import numpy as np -from numpy.testing import assert_allclose - -from .. import motions -from .. import peak_calculators - -import matplotlib.pyplot as plt - -def test_compatible_rvt_motion(): - # Compute the target from the point source model. - target = motions.SourceTheoryMotion( - 6., 20., 'wna', - peak_calculator=peak_calculators.DerKiureghian1985()) - target.compute_fourier_amps(np.logspace(-1.5, 2, 1024)) - - osc_freqs = np.logspace(-1, 2, num=50) - osc_accels_target = target.compute_osc_accels(osc_freqs, 0.05) - - compat = motions.CompatibleRvtMotion( - osc_freqs, osc_accels_target, - duration=target.duration, osc_damping=0.05, - peak_calculator=peak_calculators.DerKiureghian1985()) - - osc_accels_compat = compat.compute_osc_accels(osc_freqs, 0.05) - - # Might be off by a few percent because of difficulties with the inversion. - assert_allclose(osc_accels_target, osc_accels_compat, rtol=0.03, atol=0.05) - - # fig, axes = plt.subplots(2, 1) - # - # ax = axes.flat[0] - # - # ax.plot(target.freqs, target.fourier_amps, 'b-', label='Target') - # ax.plot(compat.freqs, compat.fourier_amps, 'r--', label='Compatible') - # - # ax.set_xlabel('Frequency (Hz)') - # ax.set_xscale('log') - # ax.set_ylabel('Fourier Ampl. (cm/s)') - # ax.set_yscale('log') - # - # ax = axes.flat[1] - # - # ax.plot(osc_freqs, osc_resp_target, 'b-', label='Target') - # ax.plot(osc_freqs, osc_resp_compat, 'r--', label='Compatible') - # - # ax.set_xlabel('Frequency (Hz)') - # ax.set_xscale('log') - # ax.set_ylabel('Spectral Accel. (cm/s²)') - # ax.set_yscale('log') - # - # fig.savefig('compatible_fas.png', dpi=300) +#!/usr/bin/env python3 +# encoding: utf-8 + +# pyRVT: Seismological random vibration theory implemented with Python +# Copyright (C) 2013-2014 Albert R. Kottke albert.kottke@gmail.com +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU 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 General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +import numpy as np +from numpy.testing import assert_allclose + +from .. import motions +from .. import peak_calculators + +import matplotlib.pyplot as plt + +def test_compute_attenuation(): + m = motions.SourceTheoryMotion(5.5, 0, 'cena', depth=1) + m.compute_fourier_amps() + + atten = m.compute_attenuation(50) + assert_allclose(0.006, atten, rtol=0.01) + +def test_compute_attenuation_full(): + m = motions.SourceTheoryMotion(5.5, 0, 'cena', depth=1) + m.compute_fourier_amps() + + atten, r_value, freqs, fitted = m.compute_attenuation(50, full=True) + + assert_allclose(0.006, atten, rtol=0.01) + assert_allclose(1.0, r_value, rtol=0.01) + + # fig = plt.figure() + # ax = fig.add_subplot(111) + # ax.plot(m.freqs, m.fourier_amps, 'b-') + # ax.set_xlabel('Frequency (Hz)') + # ax.set_xscale('log') + # ax.set_ylabel('Amplitude') + # ax.set_yscale('log') + # + # fig.savefig('test') + +def test_compatible_rvt_motion(): + # Compute the target from the point source model. + target = motions.SourceTheoryMotion( + 6., 20., 'wna', + peak_calculator=peak_calculators.DerKiureghian1985()) + target.compute_fourier_amps(np.logspace(-1.5, 2, 1024)) + + osc_freqs = np.logspace(-1, 2, num=50) + osc_accels_target = target.compute_osc_accels(osc_freqs, 0.05) + + compat = motions.CompatibleRvtMotion( + osc_freqs, osc_accels_target, + duration=target.duration, osc_damping=0.05, + peak_calculator=peak_calculators.DerKiureghian1985()) + + osc_accels_compat = compat.compute_osc_accels(osc_freqs, 0.05) + + # Might be off by a few percent because of difficulties with the inversion. + assert_allclose(osc_accels_target, osc_accels_compat, rtol=0.03, atol=0.05) + + # fig, axes = plt.subplots(2, 1) + # + # ax = axes.flat[0] + # + # ax.plot(target.freqs, target.fourier_amps, 'b-', label='Target') + # ax.plot(compat.freqs, compat.fourier_amps, 'r--', label='Compatible') + # + # ax.set_xlabel('Frequency (Hz)') + # ax.set_xscale('log') + # ax.set_ylabel('Fourier Ampl. (cm/s)') + # ax.set_yscale('log') + # + # ax = axes.flat[1] + # + # ax.plot(osc_freqs, osc_resp_target, 'b-', label='Target') + # ax.plot(osc_freqs, osc_resp_compat, 'r--', label='Compatible') + # + # ax.set_xlabel('Frequency (Hz)') + # ax.set_xscale('log') + # ax.set_ylabel('Spectral Accel. (cm/s²)') + # ax.set_yscale('log') + # + # fig.savefig('compatible_fas.png', dpi=300) diff --git a/pyrvt/tools.py b/pyrvt/tools.py index 893f730..a888d29 100644 --- a/pyrvt/tools.py +++ b/pyrvt/tools.py @@ -18,9 +18,7 @@ # along with this program. If not, see . """ -File: tools.py -Author: Albert Kottke -Description: Tools for reading/writing of files and performing operations. +Tools for reading/writing of files and performing operations. """ import csv