Cyclic-symmetry modal analysis ============================== A rotor built from ``N`` identical sectors spanning 360° admits a *cyclic-symmetry* modal reduction: every full-rotor eigenmode satisfies .. math:: u_{\text{sector } j+1}(\text{global}) = e^{i k \alpha}\, R(\alpha)\, u_{\text{sector } j}(\text{global}) with :math:`\alpha = 2\pi / N` and :math:`R(\alpha)` the spatial rotation about the symmetry axis. The integer :math:`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 :math:`1 / N` per harmonic — for a 15-sector bladed rotor the eight needed sector solves already beat a single full-rotor eigensolve. .. contents:: :local: :depth: 2 High-level API -------------- The cyclic-aware front door is :class:`femorph_solver.CyclicModel`: it wraps a base-sector :class:`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. .. code-block:: python 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 :class:`~femorph_solver.solvers.cyclic.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: .. code-block:: python 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 :math:`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 :math:`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 ---------------------------------------- :meth:`CyclicModel.from_cdb` reads a single-sector MAPDL archive, builds the base :class:`Model` via :func:`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: .. code-block:: python 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 :class:`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: .. code-block:: python 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 :meth:`~femorph_solver.Model.save` writes, so any ``Model`` you serialise can come back through :meth:`~femorph_solver.Model.from_pv` ready-to-solve, with no APDL setup: .. code-block:: python 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 ------------------------------- :func:`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: .. code-block:: python 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 :func:`scipy.linalg.eigh` (dense, subset by index). For ``k = 0`` and ``k = N/2`` (even ``N``) the phase is :math:`\pm 1` and the reduced system is real; intermediate ``k`` give a complex Hermitian reduction. :class:`CyclicModel` does this assembly + DOF lookup for you: the inner :class:`~femorph_solver.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 :math:`R(\alpha)` that takes sector ``j`` to sector ``j+1``. :class:`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 :class:`CyclicModel` and pass ``pair_rotation=None`` to :func:`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 :doc:`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 :meth:`~femorph_solver.solvers.linear.LinearSolverBase.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: .. code-block:: bash # 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 :meth:`~femorph_solver.Model.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. :class:`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.