Source code for pyrotd

#!/usr/bin/python
import sys

import numpy as np
from pkg_resources import get_distribution
from numpy.lib import NumpyVersion

if sys.version_info >= (3, 6):
    import functools
    import multiprocessing

    processes = max(multiprocessing.cpu_count() - 1, 1)
else:
    processes = 1

__author__ = "Albert Kottke"
__copyright__ = "Copyright 2016-23 Albert Kottke"
__license__ = "MIT"
__title__ = "pyrotd"
__version__ = get_distribution("pyrotd").version

G_TO_CMPS = 980.665


[docs]def calc_oscillator_resp( freq, fourier_amp, osc_damping, osc_freq, max_freq_ratio=5.0, peak_resp_only=False, osc_type="psa", ): """Compute the time series response of an oscillator. Parameters ---------- freq : array_like frequency of the Fourier acceleration spectrum [Hz] fourier_amp : array_like Fourier acceleration spectrum [g-sec] osc_damping : float damping of the oscillator [decimal] osc_freq : float frequency of the oscillator [Hz] max_freq_ratio : float, default=5 minimum required ratio between the oscillator frequency and then maximum frequency of the time series. It is recommended that this value be 5. peak_resp_only : bool, default=False If only the peak response is returned. osc_type : str, default='psa' type of response. Options are: 'sd': spectral displacement 'sv': spectral velocity 'sa': spectral acceleration 'psv': psuedo-spectral velocity 'psa': psuedo-spectral acceleration Returns ------- response : :class:`numpy.ndarray` or float time series response of the oscillator """ ang_freq = 2 * np.pi * freq osc_ang_freq = 2 * np.pi * osc_freq # Single-degree of freedom transfer function h = 1 / ( ang_freq**2.0 - osc_ang_freq**2 - 2.0j * osc_damping * osc_ang_freq * ang_freq ) if osc_type == "sd": pass elif osc_type == "sv": h *= 1.0j * ang_freq elif osc_type == "sa": h *= 1 + (1.0j * ang_freq) ** 2 elif osc_type == "psa": h *= -(osc_ang_freq**2) elif osc_type == "psv": h *= -osc_ang_freq else: raise RuntimeError # Adjust the maximum frequency considered. The maximum frequency is 5 # times the oscillator frequency. This provides that at the oscillator # frequency there are at least tenth samples per wavelength. n = len(fourier_amp) m = max(n, int(max_freq_ratio * osc_freq / freq[1])) scale = float(m) / float(n) # Scale factor is applied to correct the amplitude of the motion for the # change in number of points resp = scale * np.fft.irfft(fourier_amp * h, 2 * (m - 1)) if peak_resp_only: resp = np.abs(resp).max() return resp
[docs]def calc_rotated_percentiles(accels, angles, percentiles=None, method="optimized"): """Compute the response spectrum for a time series. This function computes the period-dependent rotated percentiles. If `optimized` is *True*, then the percentiles are only computed at times greater than a threshold. This optimized procedure is described in Equations 2.4 and 2.5 of Stewart et al. (2017) PEER report [1]_. .. _[1]: https://peer.berkeley.edu/sites/default/files/2017_09_stewart_9.10.18.pdf Parameters ---------- accels : list of array_like pair of acceleration time series angles : array_like angles to which to compute the rotated time series percentiles : array_like or None percentiles to return method : str, default of `optimized` If `optimized`, then the rotation is done on a subset of values. If `rigorous`, then all values are considered. Returns ------- rotated_resp : :class:`np.recarray` Percentiles of the rotated response. Records have keys: 'percentile', 'spec_accel', and 'angle'. """ accels = np.asarray(accels) if method == "optimized": # Use a reduced set of accelerations level = 0.7 * np.min(np.max(np.abs(accels), axis=1)) vsum = np.sqrt(np.sum(accels**2, axis=0)) indices = vsum >= level accels = accels[:, indices] percentiles = ( np.array([0, 50, 100]) if percentiles is None else np.asarray(percentiles) ) angles = np.arange(0, 180, step=1) if angles is None else np.asarray(angles) # Compute rotated time series radians = np.radians(angles) coeffs = np.c_[np.cos(radians), np.sin(radians)] rotated_time_series = np.dot(coeffs, accels) # Sort this array based on the response peak_responses = np.abs(rotated_time_series).max(axis=1) rotated = np.rec.fromarrays([angles, peak_responses], names="angle,peak_resp") rotated.sort(order="peak_resp") # Get the peak response at the requested percentiles if NumpyVersion(np.__version__) < "1.22.0": p_peak_resps = np.percentile( rotated.peak_resp, percentiles, interpolation="linear" ) else: p_peak_resps = np.percentile(rotated.peak_resp, percentiles, method="linear") # Can only return the orientations for the minimum and maximum value as the # orientation is not unique (i.e., two values correspond to the 50% # percentile). p_angles = np.select( [np.isclose(percentiles, 0), np.isclose(percentiles, 100), True], [rotated.angle[0], rotated.angle[-1], np.nan], ) return np.rec.fromarrays( [percentiles, p_peak_resps, p_angles], names="percentile,spec_accel,angle" )
[docs]def calc_rotated_oscillator_resp( angles, percentiles, freqs, fourier_amps, osc_damping, osc_freq, max_freq_ratio=5.0, osc_type="psa", method="optimized", ): """Compute the percentiles of response of a rotated oscillator. Parameters ---------- percentiles : array_like percentiles to return. angles : array_like angles to which to compute the rotated time series. freq : array_like frequency of the Fourier acceleration spectrum [Hz] fourier_amps : [array_like, array_like] pair of Fourier acceleration spectrum [g-sec] osc_damping : float damping of the oscillator [decimal] osc_freq : float frequency of the oscillator [Hz] max_freq_ratio : float, default=5 minimum required ratio between the oscillator frequency and then maximum frequency of the time series. It is recommended that this value be 5. peak_resp_only : bool, default=False If only the peak response is returned. osc_type : str, default='psa' type of response. Options are: 'sd': spectral displacement 'sv': spectral velocity 'sa': spectral acceleration 'psv': psuedo-spectral velocity 'psa': psuedo-spectral acceleration method : str, default of `optimized` If `optimized`, then the rotation is done on a subset of values. If `rigorous`, then all values are considered. Returns ------- response : :class:`numpy.ndarray` or float time series response of the oscillator """ # Compute the oscillator responses osc_ts = np.vstack( [ calc_oscillator_resp( freqs, fa, osc_damping, osc_freq, max_freq_ratio=max_freq_ratio, peak_resp_only=False, osc_type=osc_type, ) for fa in fourier_amps ] ) # Compute the rotated values of the oscillator response rotated_percentiles = calc_rotated_percentiles( osc_ts, angles, percentiles, method=method ) # Stack all of the results return [(osc_freq,) + rp.tolist() for rp in rotated_percentiles]
[docs]def calc_spec_accels( time_step, accel_ts, osc_freqs, osc_damping=0.05, max_freq_ratio=5, osc_type="psa" ): """Compute the psuedo-spectral accelerations. Parameters ---------- time_step : float time step of the time series [s] accel_ts : array_like acceleration time series [g] osc_freqs : array_like natural frequency of the oscillators [Hz] osc_damping : float damping of the oscillator [decimal]. Default of 0.05 (i.e., 5%) max_freq_ratio : float, default=5 minimum required ratio between the oscillator frequency and then maximum frequency of the time series. It is recommended that this value be 5. osc_type : str, default='psa' type of response. Options are: 'sd': spectral displacement 'sv': spectral velocity 'sa': spectral acceleration 'psv': psuedo-spectral velocity 'psa': psuedo-spectral acceleration Returns ------- resp_spec : :class:`np.recarray` computed pseudo-spectral acceleration [g]. Records have keys: 'osc_freq', and 'spec_accel' """ fourier_amp = np.fft.rfft(accel_ts) freq = np.linspace(0, 1.0 / (2 * time_step), num=fourier_amp.size) if processes > 1: with multiprocessing.Pool(processes=processes) as pool: spec_accels = pool.map( functools.partial( calc_oscillator_resp, freq, fourier_amp, osc_damping, max_freq_ratio=max_freq_ratio, peak_resp_only=True, osc_type=osc_type, ), osc_freqs, ) else: # Single process spec_accels = [ calc_oscillator_resp( freq, fourier_amp, osc_damping, of, max_freq_ratio=max_freq_ratio, peak_resp_only=True, osc_type=osc_type, ) for of in osc_freqs ] return np.rec.fromarrays([osc_freqs, spec_accels], names="osc_freq,spec_accel")
[docs]def calc_rotated_spec_accels( time_step, accel_a, accel_b, osc_freqs, osc_damping=0.05, percentiles=None, angles=None, max_freq_ratio=5, osc_type="psa", method="optimized", ): """Compute the rotated psuedo-spectral accelerations. Parameters ---------- time_step : float time step of the time series [s] accel_a : array_like acceleration time series of the first motion [g] accel_b : array_like acceleration time series of the second motion that is perpendicular to the first motion [g] osc_freqs : array_like natural frequency of the oscillators [Hz] osc_damping : float damping of the oscillator [decimal]. Default of 0.05 (i.e., 5%) percentiles : array_like or None percentiles to return. Default of [0, 50, 100], angles : array_like or None angles to which to compute the rotated time series. Default of np.arange(0, 180, step=1) (i.e., 0, 1, 2, .., 179). max_freq_ratio : float, default=5 minimum required ratio between the oscillator frequency and then maximum frequency of the time series. It is recommended that this value be 5. method : str, default of `optimized` If `optimized`, then the rotation is done on a subset of values. If `rigorous`, then all values are considered. Returns ------- rotated_resp : :class:`np.recarray` computed pseudo-spectral acceleration [g] at each of the percentiles. Records have keys: 'osc_freq', 'percentile', 'spec_accel', and 'angle' """ percentiles = [0, 50, 100] if percentiles is None else np.asarray(percentiles) angles = np.arange(0, 180, step=1) if angles is None else np.asarray(angles) assert len(accel_a) == len(accel_b), "Time series not equal lengths!" # Compute the Fourier amplitude spectra fourier_amps = [np.fft.rfft(accel_a), np.fft.rfft(accel_b)] freqs = np.linspace(0, 1.0 / (2 * time_step), num=fourier_amps[0].size) if processes > 1: with multiprocessing.Pool(processes=processes) as pool: groups = pool.map( functools.partial( calc_rotated_oscillator_resp, angles, percentiles, freqs, fourier_amps, osc_damping, max_freq_ratio=max_freq_ratio, osc_type=osc_type, method=method, ), osc_freqs, ) else: # Single process groups = [ calc_rotated_oscillator_resp( angles, percentiles, freqs, fourier_amps, osc_damping, of, max_freq_ratio=max_freq_ratio, osc_type=osc_type, method=method, ) for of in osc_freqs ] records = [g for group in groups for g in group] # Reorganize the arrays grouping by the percentile rotated_resp = np.rec.fromrecords( records, names="osc_freq,percentile,spec_accel,angle" ) return rotated_resp
[docs]def calc_rotated_peaks( time_step, accel_a, accel_b, peak_type="pga", percentiles=None, angles=None, method="optimized", ): """Compute the rotated peak values of input accelerations. Parameters ---------- time_step : float time step of the time series [s] accel_a : array_like acceleration time series of the first motion [g] accel_b : array_like acceleration time series of the second motion that is perpendicular to the first motion [g] peak_type : str Either: pga, pgv, or pgd percentiles : array_like or None percentiles to return. Default of [0, 50, 100], angles : array_like or None angles to which to compute the rotated time series. Default of np.arange(0, 180, step=1) (i.e., 0, 1, 2, .., 179). method : str, default of `optimized` If `optimized`, then the rotation is done on a subset of values. If `rigorous`, then all values are considered. Returns ------- rotated_peak : float If `pga`, then the units are in *g*. Otherwise, they are reported in *cm/s* or *cm* for `pgv` or `pgd`, respectively. """ percentiles = [0, 50, 100] if percentiles is None else np.asarray(percentiles) angles = np.arange(0, 180, step=1) if angles is None else np.asarray(angles) assert len(accel_a) == len(accel_b), "Time series not equal lengths!" if peak_type in ["pgv", "pgd"]: # Compute the Fourier amplitude spectra fourier_amps = np.vstack([np.fft.rfft(accel_a), np.fft.rfft(accel_b)]) fourier_amps *= G_TO_CMPS # Convert to velocity or displacement using Fourier domain integration power = 2 if peak_type == "pgd" else 1 freqs = np.linspace(0, 1.0 / (2 * time_step), num=fourier_amps[0].size) fourier_amps[:, 1:] *= (2j * np.pi * freqs[1:]) ** (-power) values = [np.fft.irfft(fourier_amps[0, :]), np.fft.irfft(fourier_amps[1, :])] elif peak_type == "pga": values = [accel_a, accel_b] else: raise NotImplementedError ret = calc_rotated_percentiles(values, angles, percentiles, method=method) ret.dtype.names = "percentile", peak_type, "angle" return ret