Example 15: Catalog of published nonlinear curves¶

[1]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import pystrata

[2]:

pystrata.site._load_published_curves()
models = pystrata.site.PUBLISHED_CURVES

[3]:

df = pd.DataFrame(models).T.reset_index(drop=True)
df[["source", "suffix"]] = df.name.str.split(",", expand=True)

Published nonlinear curve library¶

[4]:

print("\n".join(df["source"].unique()))

Vucetic & Dobry (91)
EPRI (93)
GEI (83)
GeoMatrix (1990)
Idriss (1990)
Imperial Valley Soils
Iwasaki
Peninsular Range
Seed & Idriss

[5]:

def plot_group(group, title="", cycle_colors=False, label=False):
    fig, axes = plt.subplots(nrows=2, sharex=True, subplot_kw=dict(xscale="log"))

    for i, row in group.reset_index(drop=True).iterrows():
        for ax, prop in zip(axes, ["mod_reduc", "damping"]):
            if pd.isna(row[prop]):
                continue

            ax.plot(
                row[prop]["strains"],
                row[prop]["values"],
                color=f"C{i}" if cycle_colors else "C0",
                label=row["suffix"].strip() if label else None,
            )

    axes[0].set(ylabel="G/G$_{max}$")
    if title:
        axes[0].set_title(title)

    axes[1].set(ylabel="Damping (dec)", xlabel="Strain (dec)")

    if label:

        def get_labels(ax):
            return sorted([line.get_label() for line in ax.get_lines()])

        labels = [get_labels(ax) for ax in axes]
        if len(labels[1]):
            # Prefer to label damping
            axes[1].legend(
                loc="upper left",
                fontsize="x-small",
                ncols=2 if len(labels[1]) > 7 else 1,
            )
        # Add to shear-mod if different
        if labels[0] != labels[1] and len(labels[0]):
            axes[0].legend(
                loc="lower left",
                fontsize="x-small",
                ncols=2 if len(labels[0]) > 7 else 1,
            )

    fig.tight_layout()
    return fig, axes

[6]:

fig, axes = plot_group(df, cycle_colors=False, label=False)
plt.show(fig)

[7]:

for name, group in df.groupby("source"):
    fig, axes = plot_group(group, name, cycle_colors=True, label=True)
    plt.show(fig)

Empirical model comparisons¶

Fine-grained soils: Darendeli (2001) vs Wang & Stokoe (2022)¶

[8]:

plas_index = np.r_[0, 15, 30, 60]
strains = np.geomspace(0.0001, 1)

[9]:

darendelis = [
    pystrata.site.DarendeliSoilType(unit_wt=16, plas_index=pi, ocr=1, stress_mean=101.3)
    for pi in plas_index
]

wangs = [
    pystrata.site.WangSoilType(
        "clayey_soil",
        unit_wt=16,
        plas_index=pi,
        ocr=1,
        stress_mean=101.3,
        void_ratio=0.7,
        fines_cont=70,
        water_cont=25,
    )
    for pi in plas_index
]

[10]:

fig, axes = plt.subplots(nrows=2, sharex=True, subplot_kw={"xscale": "log"})

for color, models in [("C0", darendelis), ("C1", wangs)]:
    for pi, m in zip(plas_index, models):
        name = "D01" if isinstance(m, pystrata.site.DarendeliSoilType) else "WS22"

        for ax, attr in zip(axes, ["mod_reduc", "damping"]):
            scale = 100 if attr == "damping" else 1
            ax.plot(
                100 * getattr(m, attr).strains,
                scale * getattr(m, attr).values,
                label=f"{name}, PI={pi}",
                color=color,
                ls="-",
                alpha=0.8,
            )

axes[0].set(ylabel="$G/G_{\\max}$")
axes[1].legend(fontsize="x-small", loc="upper left")
axes[1].set(xlabel="Shear Strain (%)", ylabel="Damping (%)")

[10]:

[Text(0.5, 0, 'Shear Strain (%)'), Text(0, 0.5, 'Damping (%)')]

Transitional silts: Darendeli (2001) vs Alemu et al. (2025)¶

The AlemuEtAlSoilType model (Alemu et al. 2025) is valid for Pacific Northwest transitional silts with 0 ≤ PI ≤ 32, 10 ≤ p′ ≤ 125 kPa, and 1 ≤ OCR ≤ 9.1.

[12]:

plas_index_silt = np.r_[0, 10, 20, 30]
strains = np.geomspace(1e-6, 0.03)

darendelis_silt = [
    pystrata.site.DarendeliSoilType(unit_wt=18, plas_index=pi, ocr=1, stress_mean=100)
    for pi in plas_index_silt
]

alemu_silts = [
    pystrata.site.AlemuEtAlSoilType(
        unit_wt=18,
        plas_index=pi,
        ocr=1,
        stress_mean=100,
        fines_cont=0.9,
        strains=strains,
    )
    for pi in plas_index_silt
]

[13]:

fig, axes = plt.subplots(nrows=2, sharex=True, subplot_kw={"xscale": "log"})

for color, models in [("C0", darendelis_silt), ("C1", alemu_silts)]:
    for pi, m in zip(plas_index_silt, models):
        label = m.name.split(" - ")[0] if "Alemu" in m.name else "Darendeli (2001)"
        label = f"{label}, PI={pi:.0f}"

        for ax, attr in zip(axes, ["mod_reduc", "damping"]):
            scale = 100 if attr == "damping" else 1
            ax.plot(
                getattr(m, attr).strains,
                getattr(m, attr).values * scale,
                color=color,
                label=label if ax is axes[0] else None,
            )

axes[0].set_ylabel("G/Gmax")
axes[1].set_ylabel("Damping (%)")
axes[1].set_xlabel("Shear Strain (dec)")
axes[0].legend(fontsize=8, ncols=2)
fig.suptitle(
    "Darendeli (2001) [blue] vs Alemu et al. 2025 [orange] - Transitional Silts"
)
plt.tight_layout()
plt.show(fig)

Gravels: Rollins et al. (2020) vs Menq (2003)¶

Both models target gravelly soils and require the uniformity coefficient Cu. RollinsEtAlSoilType (Rollins et al. 2020, JGGE) uses large-scale field testing; MenqSoilType (Menq 2003) is derived from laboratory resonant-column/torsional-shear tests. Curves are shown for Cu = 7 and Cu = 20 at σ₀′ = 50 and 200 kPa.

[14]:

cu_vals = [7, 20]
stress_vals_gravel = [50, 200]
strains_gravel = np.geomspace(1e-6, 0.05)

rollins_models = [
    pystrata.site.RollinsEtAlSoilType(
        unit_wt=21, stress_mean=s, coef_unif=cu, strains=strains_gravel
    )
    for cu in cu_vals
    for s in stress_vals_gravel
]

menq_models = [
    pystrata.site.MenqSoilType(
        unit_wt=21,
        coef_unif=cu,
        diam_mean=10,
        stress_mean=s,
        num_cycles=10,
        strains=strains_gravel,
    )
    for cu in cu_vals
    for s in stress_vals_gravel
]

[15]:

fig, axes = plt.subplots(nrows=2, sharex=True, subplot_kw={"xscale": "log"})

ls_map = {7: "-", 20: "--"}
for color, model_list, tag in [
    ("C0", rollins_models, "Rollins (2020)"),
    ("C1", menq_models, "Menq (2003)"),
]:
    for (cu, s), m in zip(
        [(cu, s) for cu in cu_vals for s in stress_vals_gravel],
        model_list,
    ):
        label = f"{tag}, Cu={cu}, \u03c3\u2080\u2032={s} kPa"
        ls = ls_map[cu]
        for ax, attr in zip(axes, ["mod_reduc", "damping"]):
            scale = 100 if attr == "damping" else 1
            ax.plot(
                getattr(m, attr).strains,
                getattr(m, attr).values * scale,
                color=color,
                ls=ls,
                label=label if ax is axes[0] else None,
            )

axes[0].set_ylabel("G/G$_{max}$")
axes[1].set_ylabel("Damping (%)")
axes[1].set_xlabel("Shear Strain (dec)")
axes[0].legend(fontsize=7, ncols=2, loc="lower left")
fig.suptitle("Gravels: Rollins (2020) [blue] vs Menq (2003) [orange]")
plt.tight_layout()
plt.show(fig)

Example usage of the published nonlinear curve library¶

To create a soil profile, use the pystrata.site.SoilType.from_published() method. The protential curve names are here:

[16]:

print("\n".join(pystrata.site.PUBLISHED_CURVES.keys()))

Vucetic & Dobry (91), PI=0
Vucetic & Dobry (91), PI=15
Vucetic & Dobry (91), PI=30
Vucetic & Dobry (91), PI=50
Vucetic & Dobry (91), PI=100
Vucetic & Dobry (91), PI=200
EPRI (93), PI=10
EPRI (93), PI=30
EPRI (93), PI=50
EPRI (93), PI=70
EPRI (93), 0-20 ft
EPRI (93), 20-50 ft
EPRI (93), 50-120 ft
EPRI (93), 120-250 ft
EPRI (93), 250-500 ft
EPRI (93), 500-1000 ft
GEI (83), 0-50 ft
GEI (83), 50-100 ft
GEI (83), 100-250 ft
GEI (83), 250-500 ft
GEI (83), >500 ft
GeoMatrix (1990), 0-50 ft
GeoMatrix (1990), 50-150 ft
GeoMatrix (1990), >150 ft
GeoMatrix (1990), 150-300 ft
GeoMatrix (1990), >300 ft
Idriss (1990), Clay
Idriss (1990), Sand
Imperial Valley Soils, 0-300 ft
Imperial Valley Soils, >300 ft
Iwasaki, 0.25 atm
Iwasaki, 1.0 atm
Peninsular Range, Cohesionless 0-50 ft
Peninsular Range, Cohesionless 50-500 ft
Seed & Idriss, Sand Mean
Seed & Idriss, Sand Upper
Seed & Idriss, Sand Lower

[17]:

# Here we can create it with one model for both
st = pystrata.site.SoilType.from_published(unit_wt=18, model="Seed & Idriss, Sand Mean")
st

[17]:

<pystrata.site.SoilType at 0x7fe0e55e6420>

[18]:

# Here we can create it with two different models
st = pystrata.site.SoilType.from_published(
    unit_wt=18,
    model="Seed & Idriss, Sand Mean",
    model_damping="Seed & Idriss, Sand Lower",
)
st

[18]:

<pystrata.site.SoilType at 0x7fe0a4f799a0>

Site response using published nonlinear curves¶

A simple three-layer sand profile (Seed & Idriss, Sand Mean) is used to demonstrate equivalent-linear site-response analysis with pystrata.

[19]:

profile = pystrata.site.Profile(
    [
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            10,
            400,
        ),
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            10,
            450,
        ),
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            30,
            600,
        ),
        pystrata.site.Layer(pystrata.site.SoilType("Rock", 24.0, None, 0.01), 0, 1200),
    ]
).auto_discretize()

[20]:

profile = pystrata.site.Profile(
    [
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            10,
            400,
        ),
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            10,
            450,
        ),
        pystrata.site.Layer(
            pystrata.site.SoilType.from_published(
                unit_wt=18.0, model="Seed & Idriss, Sand Mean"
            ),
            30,
            600,
        ),
        pystrata.site.Layer(pystrata.site.SoilType("Rock", 24.0, None, 0.01), 0, 1200),
    ]
).auto_discretize()

[21]:

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

[22]:

calc = pystrata.propagation.EquivalentLinearCalculator()

[23]:

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

outputs = pystrata.output.OutputCollection(
    [
        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.InitialVelProfile(),
    ]
)

[24]:

calc(m, profile, profile.location("outcrop", index=-1))
outputs(calc)

[25]:

for o in outputs[:-1]:
    ax = o.plot()