Motor PSHA — revisión técnica

§ 1Resumen ejecutivo

El motor reproduce con precisión < 5 % los 12 casos del PEER PSHA Verification Project Set 1 (Thomas, Wong, Abrahamson, 2010) — el estándar de facto para calibrar códigos de PSHA. La validación cubre falla cortical vertical y buzante, cuatro MFDs (delta, trunc exp, trunc normal, Y&C 1985), variabilidad de área de ruptura, tres regímenes de truncación de sigma y sources areales/volumétricos.

§ 2Integral Cornell-McGuire

Para cada sitio el motor calcula la tasa anual media de excedencia

La integral se vectoriza sobre los IML en PSHACalculator.compute_hazard_curve(), dando ~30× speed-up frente al loop por nivel. Periodo de retorno vía λ = −ln(1 − P) / T.

§ 3Modelo de source

La clase central es SeismicSource (hazard_curve.py:216):

SourceDistance.probability siempre suma 1.0 en cada lista (base y por-magnitud), así que el motor ve una distribución propiamente normalizada.

Builder público para fallas planares

El camino recomendado para construir una source planar es build_planar_fault_source() (peer.py). Precomputa las distancias flotantes por magnitud y las conecta a mag_distances:

# Ejemplo: San Andreas segmento central
from pipeline.psha.peer import (
    PEERTruncatedExponentialMFD, build_planar_fault_source,
)

mfd = PEERTruncatedExponentialMFD(
    b_value=0.9, m_min=5.0, m_max=7.0,
    area_km2=fault.length_km * fault.width_km,
)
source = build_planar_fault_source(
    source_id="SAF-central", name="San Andreas central",
    mfd=mfd, mechanism="strike_slip", tectonic="crustal",
    site_lat=site.lat, site_lon=site.lon,
    trace_s=(fault.south_lat, fault.south_lon),
    trace_n=(fault.north_lat, fault.north_lon),
    fault_top_km=0.0, fault_bot_km=12.0,
    dip_deg=90.0, dip_direction="east",
    fault_length_km=fault.length_km,
)

Desde source_service.discretize_source() la selección es automática: si la traza tiene profundidad uniforme (falla cortical planar) se usa floating rupture; si varía con la profundidad (slab de subducción) cae al mesh 3-D existente.

§ 4Floating rupture

Para cada magnitud m el motor calcula las dimensiones de ruptura con las relaciones PEER:

La ruptura se mueve uniformemente sobre el plano con regla del punto medio — positions (i + 0.5)·Δ para i ∈ [0, n − 1] — tanto en strike como en down-dip. Cada posición recibe peso 1/n_total. Esto evita el sesgo de borde del trapezoidal endpoint-inclusive.

§ 5Distribuciones de magnitud

Todos los MFDs implementan la misma interfaz .nu (tasa anual total) y .sample_magnitudes() → (centers, probs), así que son drop-in entre sí.

Convenciones PEER bakeadas en los MFDs

• Balance de momento integrado desde M = 0, no desde M_min.
• log₁₀(M₀) = 1.5·M + 16.05 (no 16.1).
• μ (shear modulus) = 3 × 10¹¹ dyne/cm².
• Y&C 1985: el boxcar tiene densidad elevada por factor exp(β·ΔM₂) respecto a la extrapolación continua de la exp. Este fue otro bug corregido — la implementación previa usaba densidad continua en la frontera, dando ~3× menos tasa.

§ 6Ground Motion Prediction Equations

Las GMPEs implementan la signatura compute(magnitude, distance_km, vs30, period, mechanism, **kwargs) → (ln_mean, ln_sigma).

§ 7Truncación de sigma

Parámetro PSHACalculator(sigma_truncation_mode=...):

Convención R-CRISIS vs PEER

R-CRISIS usa one-sided (trunca arriba, no renormaliza). PEER 2010/106 Cases 1.8b/c usan symmetric con renormalización. Ambas dan respuestas diferentes cuando k es finito; el motor soporta las dos y deja la elección al caller vía parámetro. El default es one-sided para mantener compatibilidad con proyectos existentes calibrados contra R-CRISIS (p.ej. Manila 5BGC).

Para el caso σ = 0 (medianas sin variabilidad aleatoria) añadimos median_only porque ninguna de las dos fórmulas anteriores converge correctamente al límite σ → 0 con renormalización finita.

§ 8Validación PEER 2010/106

Los benchmarks vienen del Excel oficial del proyecto PEER (docs/Tests/Set 1 Results.xlsm), parseado por docs/Tests/parse_peer_tests.py. Estos valores superan a las tablas del PDF publicado — en particular, los Cases 1.10/1.11 fueron regenerados en el Excel con sigma untruncated aunque el texto del PDF aún dice "σ = 0", algo que confirmamos empíricamente comparando con todos los 13 códigos participantes.

Bugs detectados y corregidos durante la validación

1. Sadigh 1997 sin split M ≤ 6.5 / M > 6.5 — coeficientes unificados. Afecta aplicaciones con M_max > 6.5.
2. Rrup de falla buzante — el borde superior de la ruptura se pinaba en rx=0 sin incluir el offset horizontal cuando la ruptura flotaba down-dip. Error de hasta 170 % en footwall sites.
3. Youngs & Coppersmith 1985 sin elevación del boxcar — densidad continua en la frontera en vez de la elevada por exp(β·ΔM₂). Resultado: 3× menos tasa, 0 / 65 passing inicialmente.
4. Floating rupture endpoint-inclusive — sesgo en sitios de borde por usar trapezoidal en vez de midpoint rule. Corregido pasando a midpoint.
5. Truncación de sigma sin renormalización simétrica — faltaban los modos symmetric y median_only.
6. Piso de distancia 1 km forzado en compute_hazard_curve — mataba la cola de near-fault cuando la GMPE ya cap-ea a 0.1 km por sí sola. Ahora min_distance_km=0.1.
7. MFD truncated normal con balance de momento sobre [M_min, M_max] — PEER integra desde M = 0. Corregido en PEERTruncatedNormalMFD.
8. Grilla de area source en coord. esférica sin cos(lat) — ligera anisotropía en la densidad de puntos. Corregido en area_source_distances con spacing real 1 km.

Cómo correr la suite

# Suite PEER completa con tabla de salida
uv run --project packages/pipeline python docs/Tests/run_peer_tests.py

# Suite PEER via pytest (idéntica lógica, formato pytest)
pytest tests/pipeline/psha/test_peer_verification.py --run-slow -q

# Suite rápida (PEER oculto, smoke tests de producción incluidos)
pytest tests/pipeline -q

§ 9Regresión Manila 5BGC

Snapshot blindado contra drift silencioso en el caso con el que calibramos el motor: 5 BGC, Taguig, Metro Manila (14.54762°N, 121.04674°E, Vs30 = 520 m/s), modelo de sources PEM1 de Dynamis Associates × Appgile cargado directamente desde packages/pipeline/data/PSHA_PEM1.DAT (41 sources tras el filtro de catalog_row_to_source). Sin base de datos, sin GEM API, sin R-CRISIS — 100 % offline.

El motor calcula RotD50 (mediana direccional) y Dynamis reporta RotD100 (dirección de máxima respuesta). La ratio típica RotD100/RotD50 es 1.1–1.3, lo que explica la mayoría del 25-40 % de delta contra Dynamis. Más adelante se añadirán el término de directividad Bayless-Somerville (2013) y el floor del 80 % del código (TBI 2017) que Dynamis aplica.

Top 5 contribuciones (RP = 475 yr)

Test

3 aserciones en tests/pipeline/psha/test_manila_regression.py:

1. source_count — 41 sources (a drift aquí señalaría un cambio en el parser de PEM1 DAT).
2. design_pga — PGA a 43 / 475 / 975 / 2475 yr dentro de ±0.5 % del snapshot.
3. top_source_contributions — los 5 sources dominantes y su porcentaje dentro de ±1 punto porcentual.

Marcado como @pytest.mark.slow. Runtime end-to-end 22 s. Opt-in con pytest --run-slow.

§ 10Anchors de GMPE (snapshot regression)

47 puntos de anclaje que fijan el output exacto (ln_mean, ln_sigma) de cada GMPE en tuplas específicas de (M, R, Vs30, T, mechanism). Tolerancia 0.001 en unidades logarítmicas — aprox. 0.1 % en valor lineal. Detecta cualquier drift de coeficiente al instante.

Los tests en tests/pipeline/test_gmpes.py ya existían pero solo verifican rangos amplios ("PGA debería estar entre 0.02 y 0.5 g"). Un bug de coeficientes puede caber holgadamente dentro de esos rangos y pasar desapercibido. Los anchor tests complementan ese chequeo con valores exactos.

§ 11PEER Set 2 y Set 3 — estado

El PEER 2010/106 distribuye tres sets de casos de verificación. Set 1 está completo (§ 8). Set 2 y Set 3 están parcialmente cubiertos. Este es el estado actual y el roadmap.

Set 2 — sources múltiples, logic tree, hanging wall, intraslab

Las 13 hojas de Set 2 sí tienen benchmarks (5 core codes → mean como verdad aceptada, ±5 % por sheet — misma convención que Set 1). La Barrera principal no es el motor sino que la geometría de source está fuera del Excel y hay que reverse-engineer-la del PDF. Plan: tras la reunión con Carlos, abordar Test 2.1 (el más simple) y Test 2.2/2.3 juntos (misma source, distintas GMPEs).

Set 3 — cross-code comparison, sin benchmarks absolutos

Los 5 sheets (Test 3.1 Bending Fault, Test 3.2 Logic tree fractiles, Test 3.3 Intraslab zones, Test 3.4 Areal finite rupture, Test 3.1 GW) son ejercicios de comparación cruzada entre códigos — no hay una respuesta única aceptada. Útiles para detectar desacuerdos, no para certificar corrección absoluta. Por ese motivo Set 3 no forma parte de la suite de validación del motor.

Podrían usarse como chequeo de razonabilidad: si el motor cae dentro del cloud P5-P95 de los 13 códigos participantes en un modelo concreto, es señal de que el motor está en el rango. Esto queda registrado como mejora futura pero sin prioridad.

§ 12Mapa de archivos

packages/pipeline/pipeline/psha/
├── hazard_curve.py              Cornell-McGuire core + PSHACalculator
├── peer.py                      MFDs PEER + floating rupture + builder público
├── fault_mesh.py                3-D meshing para slabs de subducción
├── gmpes/
│   ├── sadigh_1997.py           Sadigh (referencia PEER)
│   ├── abrahamson_2014.py       ASK14
│   ├── boore_2014.py            BSSA14
│   ├── campbell_bozorgnia_2014.py  CB14
│   ├── chiou_youngs_2014.py     CY14
│   ├── bchydro_2016.py          BCHydro interface / intraslab
│   └── base.py                  Protocol GMPE
├── pga_calculator.py            helpers PGA site-level
└── gmpe_registry.py             lookup + pesos por ensemble

packages/pipeline/pipeline/services/
├── source_service.py            catalog → SeismicSource pipeline
│                                # discretize_source ahora auto-selecciona
│                                # floating rupture vs 3-D mesh
└── hazard_service.py            capa de orquestación

tests/pipeline/psha/
├── test_peer_verification.py   PEER Set 1 (12 casos, --run-slow)
├── test_peer_set2_case2.py     Set 2 Case 2 reach probe (rate ±2 %)
├── test_manila_regression.py   Manila 5BGC snapshot (--run-slow)
├── test_production_fault_source.py  smoke tests end-to-end
├── fixtures/
│   └── manila_5bgc_snapshot.json  baseline Apr-13 para regresión
└── conftest.py                  gate --run-slow

tests/pipeline/
└── test_gmpe_anchors.py         47 anchors (ln_mean, ln_sigma) por GMPE

docs/engine/
├── README.md                    overview técnico (este documento en MD)
└── psha-engine-review.html      versión HTML (este archivo)

docs/Tests/
├── run_peer_tests.py            CLI demo 910/910 bins
├── parsed_tests.json            benchmarks del Excel PEER
└── PEER_2010_106.pdf            reporte original

§ 13Limitaciones conocidas

• Fallas no-planares (multi-segmento con bends). floating_rupture_fault_discretize usa el chord entre el primer y último punto de la traza para definir el plano, conservando la longitud de arco real. Para fallas con bends fuertes la aproximación subestima Rrup en sitios cercanos al bend. Extensión futura: floating rupture piecewise-planar.
• Aspect ratio de ruptura fijo (L / W = 2). Las escalas PEER están hard-coded. Para fallas con aspect ratios atípicos (normal faults, fuertemente elongadas) se necesitaría parametrizar.
• Subducción sigue usando mesh 3-D. Para slabs con traza variable en profundidad el fault_mesh.mesh_fault_surface_with_ruptures (centroide + hack de reducción de profundidad) todavía no está PEER-verificado. Funciona razonablemente para eventos pequeños pero no tenemos un benchmark formal.
• Set 2 del PEER no validado. Non-planar faults, logic trees, fractiles, listric faults, zonas intraslab — los benchmarks existen en Set 2 Results.xlsm pero la definición de la fuente no está en los sheets. Reproducción pendiente.