Example 3: Wave propagation calculators

Use RVT input motion with:

1. Linear elastic

2. Equivalent linear (e.g., SHAKE)

3. Frequency-dependent equivalent linear

import matplotlib.pyplot as plt
import numpy as np

import pystrata

%matplotlib inline

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

Create a point source theory RVT motion

motion = pystrata.motion.SourceTheoryRvtMotion(7.0, 30, "wna")
motion.calc_fourier_amps()

Create site profile

Create a simple soil profile with a single soil layer with nonlinear properties defined by the Darendeli nonlinear model.

profile = pystrata.site.Profile(
    [
        pystrata.site.Layer(
            pystrata.site.DarendeliSoilType(
                18.0, plas_index=30, ocr=1, stress_mean=200
            ),
            30,
            400,
        ),
        pystrata.site.Layer(pystrata.site.SoilType("Rock", 24.0, None, 0.01), 0, 1200),
    ]
).auto_discretize()

Create the site response calculator

calcs = [
    ("LE", pystrata.propagation.LinearElasticCalculator()),
    ("EQL", pystrata.propagation.EquivalentLinearCalculator(strain_ratio=0.65)),
    (
        "FDM",
        pystrata.propagation.FrequencyDependentEqlCalculator(method="ka02"),
    ),
]

Specify the output

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

outputs = pystrata.output.OutputCollection(
    [
        pystrata.output.AccelTransferFunctionOutput(
            # Frequency
            freqs,
            # Location in (denominator),
            pystrata.output.OutputLocation("outcrop", index=-1),
            # Location out (numerator)
            pystrata.output.OutputLocation("outcrop", index=0),
        ),
        pystrata.output.ResponseSpectrumOutput(
            # Frequency
            freqs,
            # Location of the output
            pystrata.output.OutputLocation("outcrop", index=0),
            # Damping
            0.05,
        ),
        pystrata.output.ResponseSpectrumRatioOutput(
            # Frequency
            freqs,
            # Location in (denominator),
            pystrata.output.OutputLocation("outcrop", index=-1),
            # Location out (numerator)
            pystrata.output.OutputLocation("outcrop", index=0),
            # Damping
            0.05,
        ),
        pystrata.output.MaxStrainProfile(),
        pystrata.output.InitialVelProfile(),
        pystrata.output.CompatVelProfile(),
    ]
)

Perform the calculation

Compute the response of the site, and store the state within the calculation object. Use the calculator, to compute the outputs.

calcs

[('LE', <pystrata.propagation.LinearElasticCalculator at 0x7f698e8e70b0>),
 ('EQL', <pystrata.propagation.EquivalentLinearCalculator at 0x7f698e8e6b10>),
 ('FDM',
  <pystrata.propagation.FrequencyDependentEqlCalculator at 0x7f698e8e6a80>)]

len(profile)

20

properties = {}

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

Plot the properties of the EQL and EQL-FDM methods

for key in properties["EQL"].keys():
    fig, ax = plt.subplots()

    ax.axhline(properties["EQL"][key], label="EQL", color="C0")
    ax.plot(motion.freqs, properties["FDM"][key], label="FDM", color="C1")

    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.

for o in outputs:
    o.plot(style="indiv")

Show how to manually create a plot

for o in outputs[-1:]:
    fig, ax = plt.subplots()

    for name, refs, values in o.iter_results():
        if name == "LE":
            # No strain results for LE analyses
            continue
        ax.plot(values, refs, label=name)

    ax.set(xlabel=o.xlabel, xscale="log", ylabel=o.ylabel)
    ax.invert_yaxis()
    ax.legend()
    fig.tight_layout()