Example 11 : Time series SRA using FDM

Time series analysis to acceleration transfer functions and spectral ratios.



In [1]:

    
import matplotlib.pyplot as plt
import numpy as np

import pysra

%matplotlib inline



In [2]:

    
# Increased figure sizes
plt.rcParams['figure.dpi'] = 120

Load time series data



In [3]:

    
fname = 'data/NIS090.AT2'
with open(fname) as fp:
    next(fp)
    description = next(fp).strip()
    next(fp)
    parts = next(fp).split()
    time_step = float(parts[1])

    accels = [float(p) for l in fp for p in l.split()]

    ts = pysra.motion.TimeSeriesMotion(
        fname, description, time_step, accels)



In [4]:

    
ts.accels









    Out[4]:





array([2.33833e-07, 2.99033e-07, 5.15835e-07, ..., 4.90601e-05,
       4.94028e-05, 4.96963e-05])

There are a few supported file formats. AT2 files can be loaded as follows:



In [5]:

    
ts = pysra.motion.TimeSeriesMotion.load_at2_file(fname)
ts.accels









    Out[5]:





array([2.33833e-07, 2.99033e-07, 5.15835e-07, ..., 4.90601e-05,
       4.94028e-05, 4.96963e-05])



In [6]:

    
fig, ax = plt.subplots()
ax.plot(ts.times, ts.accels)
ax.set(xlabel='Time (sec)', ylabel='Accel (g)')
fig.tight_layout();

Create site profile

This is about the simplest profile that we can create. Linear-elastic soil and rock.



In [7]:

    
profile = pysra.site.Profile([
    pysra.site.Layer(
        pysra.site.DarendeliSoilType(
            18., plas_index=0, ocr=1, stress_mean=200),
        30, 400
    ),
    pysra.site.Layer(
        pysra.site.SoilType(
            'Rock', 24., None, 0.01
        ),
        0, 1200
    ),
])

Create the site response calculator



In [8]:

    
calcs = [
    ('EQL', pysra.propagation.EquivalentLinearCalculator()),
    ('FDM (KA)', pysra.propagation.FrequencyDependentEqlCalculator(use_smooth_spectrum=True)),
    ('FDM (ZR)', pysra.propagation.FrequencyDependentEqlCalculator(use_smooth_spectrum=False))
]

Specify the output



In [9]:

    
freqs = np.logspace(-1, np.log10(50.), num=500)

outputs = pysra.output.OutputCollection([
    pysra.output.ResponseSpectrumOutput(
        # Frequency
        freqs,
        # Location of the output
        pysra.output.OutputLocation('outcrop', index=0),
        # Damping
        0.05),
    pysra.output.ResponseSpectrumRatioOutput(
        # Frequency
        freqs,
        # Location in (denominator),
        pysra.output.OutputLocation('outcrop', index=-1),
        # Location out (numerator)
        pysra.output.OutputLocation('outcrop', index=0),
        # Damping
        0.05),
    pysra.output.AccelTransferFunctionOutput(
        # Frequency
        freqs,
        # Location in (denominator),
        pysra.output.OutputLocation('outcrop', index=-1),
        # Location out (numerator)
        pysra.output.OutputLocation('outcrop', index=0),
    ),
    pysra.output.AccelTransferFunctionOutput(
        # Frequency
        freqs,
        # Location in (denominator),
        pysra.output.OutputLocation('outcrop', index=-1),
        # Location out (numerator)
        pysra.output.OutputLocation('outcrop', index=0),
        ko_bandwidth=30,
    )
])

Perform the calculation

Compute the response of the site, and store the state within the calculation object, which is then used along with the output collection to compute the desired outputs. Also, extract the computed properties for comparison.



In [10]:

    
properties = {}
for name, calc in calcs:
    calc(ts, profile, profile.location('outcrop', index=-1))
    outputs(calc, name)
    
    properties[name] = {
        key: getattr(profile[0], key)
        for key in ['shear_mod_reduc', 'damping'] 
    }









    



/home/albert/Documents/programs/pysra/pysra/propagation.py:669: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  a, b = np.linalg.lstsq(A, np.log(strain_fas[~mask]))[0]
/home/albert/Documents/programs/pysra/pysra/propagation.py:672: RuntimeWarning: divide by zero encountered in true_divide
  shape = np.minimum(1, np.exp(-a * freqs) / np.power(freqs, b))

Plot the final properties



In [11]:

    
for key in properties['EQL'].keys():
    fig, ax = plt.subplots()
    
    for i, (k, p) in enumerate(properties.items()):
        if k == 'EQL':
            ax.axhline(p[key], label=k, color=f'C{i}')
        else:
            ax.plot(ts.freqs, p[key], label=k, color=f'C{i}')
    
    ax.set(
        ylabel={'damping': 'Damping (dec)', 
                'shear_mod_reduc': r'$G/G_{max}$'}[key],
        xlabel='Frequency (Hz)', xscale='log'
    )
    ax.legend()

Plot the outputs

Create a few plots of the output.



In [12]:

    
for output in outputs:
    fig, ax = plt.subplots()
    for name, refs, values in output.iter_results():
        ax.plot(refs, values, label=name)
        
    ax.set(xlabel=output.xlabel, xscale='log', ylabel=output.ylabel)
    ax.legend()
    fig.tight_layout();



In [ ]: