Cyclic-symmetry modal analysis#

A rotor built from N identical sectors spanning 360° admits a cyclic-symmetry modal reduction: every full-rotor eigenmode satisfies

\[u_{\text{sector } j+1}(\text{global}) = e^{i k \alpha}\, R(\alpha)\, u_{\text{sector } j}(\text{global})\]

with \(\alpha = 2\pi / N\) and \(R(\alpha)\) the spatial rotation about the symmetry axis. The integer \(k \in \{0, 1, \dots, N/2\}\) is the harmonic index (or nodal diameter). Solving one sector per k reproduces the full-rotor spectrum at a cost of \(1 / N\) per harmonic — for a 15-sector bladed rotor the eight needed sector solves already beat a single full-rotor eigensolve.

High-level API#

The cyclic-aware front door is femorph_solver.CyclicModel: it wraps a base-sector femorph_solver.Model plus the cyclic parameters (n_sectors and the rotation axis) and detects the low / high cyclic-face node pairing internally. Downstream code does not need to roll its own KD-tree or hand-build DOF index arrays.

import femorph_solver as fs

cm = fs.CyclicModel.from_pv(
    fs.examples.cyclic_bladed_rotor_sector_path(),
    n_sectors=15,
    axis="z",       # default; pass "x"/"y"/"z" or a 3-vector
)
results = cm.solve_modal(n_modes=4)
for r in results:
    print(f"k={r.harmonic_index}  f={r.frequency}")

Each entry of results is a CyclicModalResult for one harmonic index, holding omega_sq, frequency, and the complex base-sector mode_shapes. When you want the full-rotor frequency multiset (each intermediate-k mode counted twice), call cm.aggregated_frequencies(n_modes=4).

Automatic sector identification#

Leave n_sectors (and optionally axis) out of the constructor to let the wrapper infer them from the sector geometry alone:

cm = fs.CyclicModel.from_pv("rotor_sector.pv")          # auto N, axis="z"
cm = fs.CyclicModel.from_pv("rotor_sector.pv", axis="auto")  # auto N + axis

The algorithm is a port of the femorph.rotor.guess_sectors heuristic. For each candidate axis the wrapper

  1. builds a cylindrical (rho, phi) representation of the point cloud,

  2. uses same-radius nearest-neighbour pairs to bound the largest plausible sector angle (and therefore the smallest plausible N),

  3. probes increasing N values, rotating the cloud by \(360°/N\) and counting nearest-neighbour matches against the original points,

  4. accepts the first N whose match count exceeds both the 180° baseline and a 0.1 %-of-points noise floor.

When axis="auto" the candidate axes (z, y, x) are ranked by the strongest matched-count score, so an axis whose “matches” come from a single coincidence loses to the one with the real cyclic alignment. Tuning is via the guess_tol kwarg (default \(10^{-4}\)); pass n_sectors and axis explicitly to skip auto-identification entirely.

Auto-identification fails (with a clear ValueError) when no candidate N clears the noise floor — typically a model that isn’t actually cyclic, or one whose interface tolerance is too loose for guess_tol.

Loading directly from a MAPDL CDB sector#

CyclicModel.from_cdb() reads a single-sector MAPDL archive, builds the base Model via femorph_solver.interop.mapdl.from_cdb(), and (with the defaults) auto-identifies the sector count and cyclic-face pairing in one call. The bundled quarter_arc_sector.cdb deck (4 sectors, 90° each) is small enough for a runnable example:

import femorph_solver as fs

cm = fs.CyclicModel.from_cdb(fs.examples.quarter_arc_sector_cdb_path())
print(cm.n_sectors, cm.axis)   # 4, [0. 0. 1.]
results = cm.solve_modal(n_modes=4)

Pass n_sectors= and/or axis= to override the auto-detect, or low_face / high_face to force a specific cyclic-face pairing. **archive_kwargs are forwarded to mapdl_archive.Archive, so e.g. parse_vtk=True is honoured. Requires the optional mapdl extra (pip install femorph_solver[mapdl]).

Manual face override#

When auto-detection misclassifies the cyclic faces — for example, on a mesh whose interface nodes don’t sit within pair_tol of their rotated partners because of upstream meshing tolerance — pass the low / high faces explicitly as grid-point-index arrays:

cm = fs.CyclicModel.from_pv(
    fs.examples.cyclic_bladed_rotor_sector_path(),
    n_sectors=15,
    low_face=low_indices,
    high_face=high_indices,
)

The arrays must have the same length and pass an R(α)-image-within- pair_tol consistency check; misuses raise ValueError with a clear message naming the offending pair.

Loading a Model from a .pv file#

The bundled rotor sector ships as a single native .pv file: geometry plus element-type assignments, material properties, and unit-system stamp baked in. This is the same format save() writes, so any Model you serialise can come back through from_pv() ready-to-solve, with no APDL setup:

import femorph_solver as fs

m = fs.Model.from_pv("rotor_sector.pv")
r = m.solve_modal(n_modes=4)         # full-sector modal (no cyclic BC)

cm = fs.CyclicModel.from_pv("rotor_sector.pv", n_sectors=15)
r = cm.solve_modal(n_modes=4)        # cyclic sweep across harmonics

Underneath: the building blocks#

femorph_solver.solvers.cyclic.solve_cyclic() is the lower- level entry point if you have a model assembled outside this package or want to drive the cyclic reduction directly:

from femorph_solver.solvers.cyclic import solve_cyclic

results = solve_cyclic(
    K,                    # (N, N) sector stiffness, real SPD
    M,                    # (N, N) sector mass, real SPD
    low_node_dofs,        # (P, d) global DOFs at P paired low-face nodes
    high_node_dofs,       # (P, d) paired high-face DOFs (same labels, same order)
    n_sectors=N,
    n_modes=10,
    pair_rotation=R_alpha,
    harmonic_indices=None,
)

The constraint matrix is assembled directly in COO form: low-face and high-face DOFs are identified, the high-face DOFs are eliminated via the phase-rotation relation, and the reduced Hermitian generalised problem is solved with scipy.linalg.eigh() (dense, subset by index). For k = 0 and k = N/2 (even N) the phase is \(\pm 1\) and the reduced system is real; intermediate k give a complex Hermitian reduction.

CyclicModel does this assembly + DOF lookup for you: the inner Model exposes stiffness_matrix() and mass_matrix(), the wrapper translates its cached point-pair arrays into global DOF indices, and the spatial rotation is built from the user-supplied axis / angle. Call cm.solve_modal(...) and you get the same CyclicModalResult list back.

Harmonic-index counting#

The cyclic sweep is equivalent to the full-rotor eigenproblem — every full-rotor frequency falls in exactly one harmonic index — but modes are degenerate in pairs for intermediate k:

  • k = 0: single real eigenmodes (the “in-phase” family)

  • k = N/2 (even N only): single real eigenmodes (anti-phase)

  • 0 < k < N/2: each eigenvalue represents a travelling-wave pair in the full rotor; count it twice when aggregating to the full-rotor multiset.

Rigid-body modes are partitioned across k = 0 (axial translation + axial rotation) and k = 1 (the in-plane translation pair).

Why the spatial rotation matters#

Stiffness / mass are usually assembled in a global XYZ frame, so the phase relation between neighbouring sectors picks up not just the scalar exp(i k α) but also the rotation \(R(\alpha)\) that takes sector j to sector j+1. CyclicModel builds it from the axis / angle kwargs (default: α = 2π / n_sectors about +z) and threads it through the lower-level solver automatically; only when DOFs are already in a per-sector cylindrical frame would you bypass CyclicModel and pass pair_rotation=None to solve_cyclic() directly.

Elements with both translational and rotational DOFs (shells, beams) use a block-diagonal R of two copies of the spatial rotation — one for the translations, one for the rotations. The lower-level solver itself is element-agnostic: it reads only the (P, d) DOF-index arrays and the (d, d) rotation.

Performance#

The cyclic-modal solve is dominated by the per-harmonic Pardiso factor and the inner ARPACK shift-invert iterations. Roughly half the per-harmonic wall is the factor and half is the eigsh sweep, both of which scale super-linearly with reduced-DOF count.

Two install-time choices have outsized impact:

  • Install pypardiso (via the pardiso extra: pip install femorph_solver[pardiso]). Without it the complex cyclic path falls back to scipy’s bundled SuperLU, whose MMD column ordering produces 2-3× more fill than Pardiso’s METIS. On rotor-scale K_red (488 k DOFs, 75 M nnz) SuperLU complex factor exceeds ~50 GiB peak RSS and OOMs on a 64 GiB host; Pardiso fits at ~22 GiB.

  • Make sure MKL_NUM_THREADS (or OMP_NUM_THREADS) is set to the host’s physical core count. Pardiso scales well to 16 threads on rotor-scale problems.

The cyclic-symmetry analysis solver path has been progressively tuned across PRs #698, #702, #703, #705, #706, #708, #715:

  • Pardiso mtype=-4 (Hermitian indefinite) for the complex harmonics — Bunch-Kaufman pivoting on the upper-triangular CSR.

  • Symbolic-factor reuse across harmonics via refactor() (MKL phase 22).

  • Pre-computed triu of the cyclic plan’s three phase buckets so the per-harmonic sp.triu is amortised.

  • ARPACK ncv = 2k+10 for slack on clustered cyclic eigenvalues.

  • Eigenvector warm-start across harmonics — cyclic modes vary continuously with k, so the previous mode 0 is a good v0 seed for the next harmonic’s Lanczos.

Cumulative impact on the 488 k-DOF reference (pbs-sector-hd.cdb, N = 20, 4 modes per harmonic, full 0..N//2 sweep): wall ≈ 28 min vs the original ~30 min Pardiso baseline, peak RSS 21.5 GiB vs 34 GiB. See perf/cyclic-perf.md in the repository for the per-iteration breakdown.

To benchmark locally:

# Bundled quarter-arc sector, runs in ~2 s
python3 perf/bench_cyclic_modal.py --mode quick

# Larger workload (specify your own CDB)
python3 perf/bench_cyclic_modal.py --mode full --cdb /path/to/sector.cdb \
    > bench_run.json

The output JSON carries per-stage wall time, peak RSS, and the per-(k, mode) frequency table — easy to diff across runs.

Bundled example — bladed-rotor sector#

The femorph-solver wheel ships a 15-sector × 24° bladed-rotor base sector (101 HEX8 cells) at femorph_solver.examples.cyclic_bladed_rotor_sector_path(). The file is a complete save() snapshot: the HEX8 element-type and isotropic material assignments are baked in, so the demo above runs without any external download or APDL setup.

Cyclic faces on this mesh are not θ=const planes — they’re curved surfaces that map to each other under a 24° rotation about z. CyclicModel resolves the pairing uniformly on the rotated point-cloud match.

The deeper-history version of this example (driven directly off mapdl_archive.examples.sector_archive_file with hand-rolled cKDTree pair detection) is exercised in tests/elements/solid185/test_solid185_cyclic_sector_archive.py as a regression against the legacy code path.