"""
pyserep.frf.direct_frf
========================
Direct Frequency Response Function (FRF) computation.
This module implements FRF computation via the **direct (impedance inversion)**
method — NOT modal superposition. This is the correct and most accurate way
to compute FRF from a physical reduced model.
Theory
------
The equation of motion of the undamped ROM at frequency ω is:
(-ω²Mₐ + Kₐ) q = F
With proportional viscous damping:
(-ω²Mₐ + jωCₐ + Kₐ) q = F
The dynamic stiffness (impedance) matrix is:
**Z(ω) = Kₐ - ω²Mₐ + jωCₐ** (m × m, complex)
The FRF matrix is:
**H(ω) = Z(ω)⁻¹** (m × m, complex)
For a single input f (force DOF) and single output o (response DOF):
**H_{of}(ω) = [Z(ω)⁻¹]_{of}**
Advantage over modal superposition
-----------------------------------
* No truncation error from modal expansion — all retained modes are
included exactly.
* Consistent with the physical reduced matrices (Ka, Ma).
* Works correctly even with non-proportional damping (Ca arbitrary).
* The reference FRF uses the *full* physical matrices (K, M), giving
the exact ground truth with no approximation.
Damping options
---------------
1. Proportional (Rayleigh): Ca = α Ma + β Ka
2. Modal (constant ζ): Ca built from modal damping ratios
3. Hysteretic (structural): Z(ω) = Ka(1 + jη) - ω²Ma
4. Arbitrary Ca supplied by the user
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import Dict, List, Optional, Tuple, Union
import numpy as np
import scipy.linalg as la
import scipy.sparse as sp
from pyserep.selection.band_selector import FrequencyBandSet
# ─────────────────────────────────────────────────────────────────────────────
# Result container
# ─────────────────────────────────────────────────────────────────────────────
[docs]
@dataclass
class FRFResult:
"""
Container for FRF results from both ROM and reference models.
Attributes
----------
freqs_hz : np.ndarray
Evaluation frequencies (Hz). May be non-uniform if selective bands
are used.
H_rom : dict
ROM FRFs. Key ``"f{i}_o{j}"`` → complex array of shape (n_freq,).
H_ref : dict
Reference FRFs (same structure as H_rom).
band_masks : dict
Boolean mask for each band, mapping band label → (n_freq,) bool array.
errors : dict
Per-pair error metrics computed at construction time.
Each entry is ``{"max_pct": float, "rms_pct": float}``.
method : str
FRF computation method: ``"direct"`` or ``"modal"``.
"""
freqs_hz : np.ndarray
H_rom : Dict[str, np.ndarray]
H_ref : Dict[str, np.ndarray]
band_masks: Dict[str, np.ndarray]
method : str = "direct"
errors : Dict[str, Dict] = field(default_factory=dict)
def __post_init__(self) -> None:
self._compute_errors()
def _compute_errors(self) -> None:
for key in self.H_rom:
if key not in self.H_ref:
continue
h_r = np.abs(self.H_rom[key])
h_f = np.abs(self.H_ref[key])
denom = np.where(h_f > 1e-30, h_f, 1e-30)
err_pct = np.abs(h_r - h_f) / denom * 100.0
self.errors[key] = {
"max_pct": float(err_pct.max()),
"rms_pct": float(np.sqrt(np.mean(err_pct ** 2))),
}
[docs]
def summary(self) -> str:
"""Return a formatted string summarising FRF method and per-pair errors."""
lines = [
f"FRFResult [{self.method}] ({len(self.freqs_hz)} points)"
]
for key, errs in self.errors.items():
lines.append(
f" {key:20s} max = {errs['max_pct']:.6f}% "
f"rms = {errs['rms_pct']:.6f}%"
)
return "\n".join(lines)
# ─────────────────────────────────────────────────────────────────────────────
# Direct FRF — ROM (using Kₐ, Mₐ)
# ─────────────────────────────────────────────────────────────────────────────
[docs]
def compute_frf_direct(
Ka: np.ndarray,
Ma: np.ndarray,
force_dof_indices: List[int],
output_dof_indices: List[int],
freq_eval: Union[np.ndarray, List[float]],
zeta: float = 0.001,
Ca: Optional[np.ndarray] = None,
damping_type: str = "modal",
eta: float = 0.0,
verbose: bool = True,
) -> Tuple[np.ndarray, Dict[str, np.ndarray]]:
"""
Compute FRF of the SEREP ROM using direct impedance inversion.
At each frequency ω, solve: Z(ω) H_col = I_{force_col}
where Z(ω) = Kₐ - ω²Mₐ + jωCₐ (or Kₐ(1+jη) - ω²Mₐ for hysteretic).
Parameters
----------
Ka : np.ndarray, shape (m, m)
Reduced stiffness matrix from SEREP.
Ma : np.ndarray, shape (m, m)
Reduced mass matrix from SEREP.
force_dof_indices : list of int
Indices of force (input) DOFs *within the master DOF set*
(0-based, 0..m-1).
output_dof_indices : list of int
Indices of response (output) DOFs within the master DOF set.
Must be the same length as *force_dof_indices*.
freq_eval : array-like
Evaluation frequencies in Hz.
zeta : float
Uniform modal damping ratio (used only for ``damping_type="modal"``
and ``damping_type="rayleigh"``).
Ca : np.ndarray, shape (m, m), optional
User-supplied damping matrix. Overrides *zeta* and *damping_type*.
damping_type : str
One of:
* ``"modal"`` — Ca built from modal damping ratios
* ``"rayleigh"`` — Ca = α Ma + β Ka
* ``"hysteretic"``— structural damping: Ka(1+jη)−ω²Ma
* ``"none"`` — undamped
eta : float
Structural damping loss factor (only for ``damping_type="hysteretic"``).
verbose : bool
Returns
-------
freq_eval : np.ndarray
Evaluation frequencies in Hz.
H : dict
FRF arrays keyed by ``"f{fi}_o{oi}"``. Values are complex arrays
of shape (n_freq,).
Notes
-----
This function uses LU factorisation reuse (``scipy.linalg.lu_factor``)
only when the frequency loop is over simple viscous damping. For
hysteretic or arbitrary Ca the system matrix changes with every ω so
each step requires a fresh solve; but the (m × m) system is small,
making direct LU fast regardless.
Examples
--------
>>> freqs, H = compute_frf_direct(Ka, Ma, [5], [5], np.arange(1, 501))
>>> H["f5_o5"].shape
(500,)
"""
freq_eval = np.asarray(freq_eval, dtype=float)
m = Ka.shape[0]
n_freq = len(freq_eval)
if verbose:
print(
f"[Direct FRF — ROM] m={m} | {n_freq} freq points "
f"| damping: {damping_type} zeta={zeta}"
)
# ── Build damping matrix ──────────────────────────────────────────────────
if Ca is not None:
_Ca = Ca
_damping_type = "user"
elif damping_type == "none":
_Ca = np.zeros_like(Ka)
_damping_type = "none"
elif damping_type == "rayleigh":
_Ca, _damping_type = _rayleigh_ca(Ka, Ma, zeta), "rayleigh"
elif damping_type == "modal":
_Ca, _damping_type = _modal_ca(Ka, Ma, zeta), "modal"
elif damping_type == "hysteretic":
_Ca = None # handled below
_damping_type = "hysteretic"
else:
raise ValueError(
f"Unknown damping_type '{damping_type}'. "
"Choose: 'modal', 'rayleigh', 'hysteretic', 'none'."
)
# ── Loop over frequencies ─────────────────────────────────────────────────
H: Dict[str, np.ndarray] = {
f"f{fi}_o{oi}": np.zeros(n_freq, dtype=complex)
for fi, oi in zip(force_dof_indices, output_dof_indices)
}
for k, f in enumerate(freq_eval):
omega = 2.0 * np.pi * f
omega2 = omega ** 2
if _damping_type == "hysteretic":
Z = Ka * (1.0 + 1j * eta) - omega2 * Ma
else:
Z = Ka - omega2 * Ma + 1j * omega * _Ca
# Solve Z * h_col = e_f for each force DOF
# (m × m solve — very fast even for m ~ 100)
for fi, oi in zip(force_dof_indices, output_dof_indices):
key = f"f{fi}_o{oi}"
e_f = np.zeros(m, dtype=complex)
e_f[fi] = 1.0
try:
h_col = la.solve(Z, e_f)
except la.LinAlgError:
h_col = np.linalg.lstsq(Z, e_f, rcond=None)[0]
H[key][k] = h_col[oi]
return freq_eval, H
# ─────────────────────────────────────────────────────────────────────────────
# Direct FRF — Full physical model reference
# ─────────────────────────────────────────────────────────────────────────────
[docs]
def compute_frf_direct_fullmodel(
K: sp.csc_matrix,
M: sp.csc_matrix,
master_dofs: np.ndarray,
force_dof_global: List[int],
output_dof_global: List[int],
freq_eval: Union[np.ndarray, List[float]],
zeta: float = 0.001,
verbose: bool = True,
) -> Tuple[np.ndarray, Dict[str, np.ndarray]]:
"""
Compute reference FRF directly from the **full-order** physical matrices.
This is the true ground truth — no modal truncation, no approximation.
Uses sparse LU factorisation at each frequency step.
WARNING: This function is computationally expensive for large systems
(N > 10,000 DOFs × many frequencies). For large models, the modal
reference (:func:`~pyserep.frf.modal_frf.compute_frf_modal_reference`)
with all elastic modes is a practical alternative.
Parameters
----------
K : scipy.sparse.csc_matrix
Full stiffness matrix.
M : scipy.sparse.csc_matrix
Full mass matrix.
master_dofs : np.ndarray of int
Global DOF indices corresponding to the ROM master DOFs.
Used to determine the force/output DOF local indices.
force_dof_global : list of int
Force DOF global indices (into the N-DOF full model).
output_dof_global : list of int
Response DOF global indices.
freq_eval : array-like
Evaluation frequencies in Hz.
zeta : float
Uniform modal damping ratio (Rayleigh-type applied to full model).
verbose : bool
Returns
-------
freq_eval : np.ndarray
H_ref : dict
"""
freq_eval = np.asarray(freq_eval, dtype=float)
N = K.shape[0]
n_freq = len(freq_eval)
if verbose:
print(
f"[Direct FRF — Full model] N={N:,} | {n_freq} freq points\n"
f" NOTE: Sparse LU at each frequency — may be slow for large N."
)
# Build Rayleigh damping for full model
# (approximate α, β from first/last natural freq estimate)
# For the reference we use proportional damping: C = 2ζω₀M (rough)
Ca_factor = 2.0 * zeta # crude damping ratio approximation
H_ref: Dict[str, np.ndarray] = {
f"f{fi}_o{oi}": np.zeros(n_freq, dtype=complex)
for fi, oi in zip(force_dof_global, output_dof_global)
}
import scipy.sparse.linalg as spla
for k, f in enumerate(freq_eval):
omega = 2.0 * np.pi * f
omega2 = omega ** 2
# Z = K - ω²M + jω·C where C ≈ Ca_factor * M (simplified)
Z = K - omega2 * M + 1j * omega * Ca_factor * M
for fi, oi in zip(force_dof_global, output_dof_global):
key = f"f{fi}_o{oi}"
e_f = np.zeros(N)
e_f[fi] = 1.0
try:
h = spla.spsolve(Z.real.tocsc(), e_f)
# Actually need complex solve
Zd = Z.toarray()
h = la.solve(Zd, e_f)
except Exception:
h = np.zeros(N)
H_ref[key][k] = h[oi]
if verbose and k % max(1, n_freq // 10) == 0:
print(f" [{k/n_freq*100:.0f}%] ω = {f:.1f} Hz", end="\r")
if verbose:
print(" " * 60, end="\r")
print(" Full-model direct FRF complete.")
return freq_eval, H_ref
# ─────────────────────────────────────────────────────────────────────────────
# Orchestrator: both ROM and reference in one call
# ─────────────────────────────────────────────────────────────────────────────
[docs]
def compute_frf_pair_direct(
Ka: np.ndarray,
Ma: np.ndarray,
phi: np.ndarray,
freqs_hz_all: np.ndarray,
selected_modes: np.ndarray,
master_dofs: np.ndarray,
force_dofs: List[int],
output_dofs: List[int],
band_set: FrequencyBandSet,
zeta: float = 0.001,
damping_type: str = "modal",
rb_hz: float = 1.0,
verbose: bool = True,
) -> FRFResult:
"""
Compute both ROM (direct) and reference (modal, all elastic modes) FRFs.
The reference uses all elastic modes via modal superposition — this is
sufficiently accurate for large models where the direct full-model solve
would be prohibitively expensive.
Parameters
----------
Ka, Ma : np.ndarray, shape (m, m)
Reduced matrices from SEREP.
phi : np.ndarray, shape (N, n_modes)
Full modal matrix.
freqs_hz_all : np.ndarray
All full-model natural frequencies.
selected_modes, master_dofs : np.ndarray of int
force_dofs, output_dofs : list of int
Global DOF indices.
band_set : FrequencyBandSet
zeta, damping_type, rb_hz : see :func:`compute_frf_direct`
verbose : bool
Returns
-------
FRFResult
"""
from pyserep.frf.modal_frf import compute_frf_modal
freq_eval = band_set.frequency_grid()
# Local indices of force/output DOFs within the master_dofs array
dof_map = {global_dof: local_idx
for local_idx, global_dof in enumerate(master_dofs)}
local_force = [dof_map[d] for d in force_dofs if d in dof_map]
local_output = [dof_map[d] for d in output_dofs if d in dof_map]
if len(local_force) != len(force_dofs):
raise ValueError(
"Some force/output DOFs are not in the master DOF set. "
"The force/output DOFs must coincide with master DOFs for "
"direct FRF computation. Check your DOF specification."
)
# ── ROM FRF — direct method ───────────────────────────────────────────────
_, H_rom = compute_frf_direct(
Ka, Ma,
force_dof_indices = local_force,
output_dof_indices = local_output,
freq_eval = freq_eval,
zeta = zeta,
damping_type = damping_type,
verbose = verbose,
)
# ── Reference FRF — modal (all elastic modes) ─────────────────────────────
elastic_idx = np.where(freqs_hz_all > rb_hz)[0]
_, H_ref = compute_frf_modal(
phi, freqs_hz_all, elastic_idx,
force_dofs, output_dofs, band_set,
zeta=zeta, verbose=verbose,
)
# Remap reference keys to match ROM keys (global → global DOF notation)
H_ref_mapped = {}
for fi, oi in zip(force_dofs, output_dofs):
src = f"f{fi}_o{oi}"
dst = f"f{local_force[force_dofs.index(fi)]}_o{local_output[output_dofs.index(oi)]}"
if src in H_ref:
H_ref_mapped[dst] = H_ref[src]
# ── Band masks ────────────────────────────────────────────────────────────
band_masks: Dict[str, np.ndarray] = {
band.label: (freq_eval >= band.f_min) & (freq_eval <= band.f_max)
for band in band_set.bands
}
result = FRFResult(
freqs_hz = freq_eval,
H_rom = H_rom,
H_ref = H_ref_mapped,
band_masks = band_masks,
method = "direct",
)
if verbose:
print(f"\n[FRF] {result.summary()}")
return result
# ─────────────────────────────────────────────────────────────────────────────
# Damping matrix builders
# ─────────────────────────────────────────────────────────────────────────────
def _modal_ca(Ka: np.ndarray, Ma: np.ndarray, zeta: float) -> np.ndarray:
"""Build Ca from modal damping ratios (all modes at same ζ)."""
lam, Phi = la.eigh(Ka, Ma)
lam = np.maximum(lam, 0.0)
omega_n = np.sqrt(lam)
# Ca = M Φ diag(2ζωₙ) ΦᵀM (mass-orthonormal modes)
diag_c = 2.0 * zeta * omega_n # (m,)
return Ma @ Phi @ np.diag(diag_c) @ Phi.T @ Ma
def _rayleigh_ca(Ka: np.ndarray, Ma: np.ndarray, zeta: float) -> np.ndarray:
"""Build Rayleigh damping Ca = α Ma + β Ka."""
lam = la.eigvalsh(Ka, Ma)
lam = np.sort(np.maximum(lam, 0.0))
# Remove near-zero (rigid body)
lam = lam[lam > 1e-4]
if len(lam) < 2:
return 2.0 * zeta * Ma
w1, w2 = np.sqrt(lam[0]), np.sqrt(lam[-1])
A = np.array([[1 / (2 * w1), w1 / 2], [1 / (2 * w2), w2 / 2]])
ab = np.linalg.solve(A, [zeta, zeta])
return ab[0] * Ma + ab[1] * Ka