docs(validation): record the supercycling verification

CHANGELOG.md

@@ -31,6 +31,10 @@ All notable changes to CableDyn are recorded in this file. The format follows
   `dtM = DT` (the default) do not change; the previous solve's loads are checkpointed, so a restart
   continues the extrapolation.
 
+- `validation/scripts/openfast_supercycle_ab.py`: compares coupled runs at a coarser `dtM` with
+  the same model at `dtM = DT` (platform time series, tension statistics, rainflow
+  damage-equivalent loads); `VALIDATION.md` records the supercycling verification it scores.
+
 ### Changed
 
 - OpenFAST: `CalcOutput` no longer evaluates the reactions of the finite-EI cables when it saves

VALIDATION.md

@@ -321,6 +321,23 @@ scripts, results, and source provenance are in
     tension statistics match CableDyn's own per-glue-step solution to 0.01% mean /
     ~3% std. Tension agreement with MoorDyn is scored in the settled case above and
     in the DLC matrix below.
+  - *Supercycling* (`dtM` = 4 `DT` = 0.1 s): the chain mooring in a DLC-1.6
+    environment (NTM 10.6 m/s turbulent wind, JONSWAP Hs 8 m / Tp 12 s, 3800 s), scored
+    over 200–3800 s on the common 0.1 s grid against the same model at `dtM = DT`
+    (`validation/scripts/openfast_supercycle_ab.py`). With the coupled loads
+    extrapolated between CableDyn solves (the default), every platform DOF's time series
+    differs by at most 0.14 % of its standard deviation (0.41 % with held loads); the
+    fairlead and anchor tensions keep their mean within 0.01 %, std within 0.12 % and
+    maximum within 0.16 %; their damage-equivalent loads (rainflow, slope 3) change by
+    at most 0.21 % at the fairleads and 1.5 % at the anchors (held: 0.27 % / 1.1 %), the
+    anchor changes being of the order of the reference's own sampling effect (its DEL
+    from 0.025 s against 0.1 s samples, up to 0.94 %). Against stock MoorDyn (same
+    seeds, 200–1200 s) the in-plane platform standard deviations agree within 0.03 %
+    (0.04 % at `dtM = DT`, 0.14 % held). The tension time series differ from those at
+    `dtM = DT` by 5–7 % of their standard deviation on the leeward lines, because the
+    lines step on the coarser `dtM`; statistics and damage are unaffected. A finite-EI
+    cable needs a finer `dtM` (the initialisation advisory; 0.05 s for the mixed
+    example's power cable).
   - *Mixed deck*: the shipped example `iea15mw_umaine_mixed_cabledyn.dat` (three
     chain moorings and a 248-element finite-EI power cable that touches down on the
     200 m seabed) runs to completion under full coupled physics with finite, bounded

doc/openfast.rst

@@ -107,7 +107,9 @@ peaks, and damage. If CableDyn needs a smaller step than the glue step, reduce `
    the ``dtM = DT`` run by at most 0.12 % of their standard deviation with the extrapolation, on
    the chain mooring as on the chain and power-cable mooring, where the held loads give up to
    0.42 %. The chain-only run took a fifth less time (116 s against 141 s for ``dtM = DT``
-   on the same machine). In the three-seed DLC-1.1 comparison
+   on the same machine). Over a one-hour run in the same environment the fairlead and
+   anchor tension statistics stay within 0.16 % and their fatigue damage-equivalent loads
+   within 1.5 % (see :doc:`validation`). In the three-seed DLC-1.1 comparison
    of :doc:`validation`, made with held loads, the platform motions with ``dtM`` above ``DT``
    agree with MoorDyn within 0.6 % in mean and standard deviation, but the fairlead tension's
    standard deviation differs by about 3 % from CableDyn's own solution at every glue step.

tests/test_supercycle_ab.py

@@ -0,0 +1,85 @@
+# SPDX-License-Identifier: Apache-2.0
+"""The rainflow count and damage-equivalent load of validation/scripts/openfast_supercycle_ab.py."""
+
+from __future__ import annotations
+
+import importlib.util
+import math
+import unittest
+from collections import Counter
+from pathlib import Path
+
+_SCRIPT = Path(__file__).resolve().parents[1] / "validation" / "scripts" / "openfast_supercycle_ab.py"
+_SPEC = importlib.util.spec_from_file_location("openfast_supercycle_ab", _SCRIPT)
+ab = importlib.util.module_from_spec(_SPEC)
+_SPEC.loader.exec_module(ab)
+
+# ASTM E1049-85 (2017), Fig. 6: the rainflow counts of this history
+_ASTM_HISTORY = [-2.0, 1.0, -3.0, 5.0, -1.0, 3.0, -4.0, 4.0, -2.0]
+_ASTM_COUNTS = {3.0: 0.5, 4.0: 1.5, 6.0: 0.5, 8.0: 1.0, 9.0: 0.5}
+
+
+def _counts(history):
+    c = Counter()
+    for r, n in ab.rainflow(history):
+        c[r] += n
+    return dict(c)
+
+
+class RainflowTest(unittest.TestCase):
+    def test_astm_example(self):
+        self.assertEqual(_counts(_ASTM_HISTORY), _ASTM_COUNTS)
+
+    def test_plateaus_count_once(self):
+        # equal neighbours (a tension written to seven digits) neither add nor drop a turning point
+        plateaued = [-2.0, 1.0, 1.0, -3.0, -3.0, 5.0, 5.0, 5.0, -1.0, 3.0, -4.0, 4.0, -2.0]
+        self.assertEqual(_counts(plateaued), _ASTM_COUNTS)
+        self.assertEqual(ab.turning_points([1.0, 2.0, 2.0, 1.0]), [1.0, 2.0, 1.0])
+
+    def test_damage_equivalent_load(self):
+        # a sine of amplitude A over N whole periods in T seconds: N full cycles of range 2A
+        a, n, t = 3.0, 40, 200.0
+        history = [a * ((-1.0) ** k) for k in range(2 * n + 1)]
+        expected = (n * (2.0 * a) ** 3 / t) ** (1.0 / 3.0)
+        self.assertAlmostEqual(ab.damage_equivalent_load(history, t, 3.0), expected, places=12)
+
+
+class CompareChannelTest(unittest.TestCase):
+    @staticmethod
+    def _grid(values):
+        return {k: [0.1 * k, v] for k, v in enumerate(values)}
+
+    def test_truncated_twin_is_compared_over_its_overlap(self):
+        # a run identical to the reference but ending at a tenth of it: no difference at all,
+        # and the partial coverage is reported (the DEL is not set against the whole reference)
+        ref = [3.0 * ((-1.0) ** k) * (1.0 + 0.01 * k) for k in range(1000)]
+        c = ab.compare_channel(self._grid(ref[:100]), 1, self._grid(ref), 1, 0.1, 3.0)
+        self.assertEqual(c["rel"], 0.0)
+        self.assertEqual(c["ddel"], 0.0)
+        self.assertAlmostEqual(c["covered"], 0.1)
+
+    def test_a_scaled_run_changes_its_del_by_the_scale(self):
+        ref = [3.0 * ((-1.0) ** k) for k in range(400)]
+        c = ab.compare_channel(self._grid([1.02 * v for v in ref]), 1, self._grid(ref), 1, 0.1, 3.0)
+        self.assertAlmostEqual(c["ddel"], 2.0, places=9)
+        self.assertEqual(c["covered"], 1.0)
+
+    def test_zero_damage_reference(self):
+        # a constant tension has no cycles: identical constant histories change by 0 %, and a run
+        # with cycles against a constant reference has an undefined change (not a crash)
+        flat = [1000.0] * 200
+        same = ab.compare_channel(self._grid(flat), 1, self._grid(flat), 1, 0.1, 3.0)
+        self.assertEqual(same["del"], 0.0)
+        self.assertEqual(same["ddel"], 0.0)
+        wavy = ab.compare_channel(self._grid([1000.0 + (-1.0) ** k for k in range(200)]), 1,
+                                  self._grid(flat), 1, 0.1, 3.0)
+        self.assertGreater(wavy["del"], 0.0)
+        self.assertTrue(math.isnan(wavy["ddel"]))
+
+    def test_no_common_samples(self):
+        self.assertIsNone(ab.compare_channel({0: [0.0, 1.0]}, 1, self._grid([1.0, 2.0, 3.0]), 1, 0.1, 3.0))
+        self.assertEqual(ab.damage_equivalent_load([1.0, 2.0], 0.0, 3.0), 0.0)
+
+
+if __name__ == "__main__":
+    unittest.main()

validation/scripts/openfast_supercycle_ab.py

@@ -0,0 +1,205 @@
+#!/usr/bin/env python3
+# SPDX-License-Identifier: Apache-2.0
+"""Supercycling A/B: coupled CableDyn runs at a coarser dtM against the same model at dtM = DT.
+
+Twin OpenFAST runs of the SAME model and environment differ only in the CableDyn time step (and,
+with dtM larger than DT, in how the coupled loads are served between CableDyn solves: deck
+OPTIONS predictLoads True extrapolates them, False holds them). The reference is the run at
+dtM = DT. Over the window t >= t_start, on the common grid of the coarsest committed CableDyn
+samples (--grid, default 0.1 s), the script reports per channel:
+
+  * platform DOFs (glue <root>.out): mean, std, max and the RMS of the time-series difference
+    relative to the reference std;
+  * fairlead and anchor tensions (CableDyn <root>.CD.out, committed samples): the same, plus the
+    damage-equivalent load DEL = (sum n_i S_i^m / T)^(1/m) from a rainflow count (ASTM E1049,
+    four-point method, residue counted as half cycles; a plateau of equal samples at the
+    output precision counts as one turning point), slope m (--slope, default 3 for steel chain),
+    1 Hz equivalent count, and its change against the reference DEL over the same samples.
+
+Every statistic of a run is taken over the samples it shares with the reference; a run that
+covers fewer of them (a shorter run, missing rows) is flagged in the report.
+
+It also reports the reference DEL at its native sampling against the common grid, the part of
+a DEL difference that is sampling alone. It reports; it does not gate.
+
+Usage:
+    python validation/scripts/openfast_supercycle_ab.py REF_ROOT LABEL=ROOT [LABEL=ROOT ...]
+        [--t-start 200] [--grid 0.1] [--slope 3]
+"""
+
+from __future__ import annotations
+
+import argparse
+import math
+from pathlib import Path
+
+PLATFORM = ["PtfmSurge", "PtfmSway", "PtfmHeave", "PtfmRoll", "PtfmPitch", "PtfmYaw"]
+TENSION_PREFIXES = ("FairTen", "AnchTen")
+
+
+def read_table(path: Path) -> tuple[list[str], list[list[float]]]:
+    """An OpenFAST or CableDyn text table: the header row starting with Time, the units row,
+    then numeric rows (rows that do not parse to the header width are skipped)."""
+    lines = path.read_text(encoding="utf-8", errors="replace").splitlines()
+    hi = next((i for i, ln in enumerate(lines) if ln.split()[:1] == ["Time"]), -1)
+    if hi < 0:
+        raise ValueError(f"{path}: no Time header row")
+    names = lines[hi].split()
+    rows = []
+    for ln in lines[hi + 2:]:
+        parts = ln.split()
+        if len(parts) != len(names):
+            continue
+        try:
+            rows.append([float(x) for x in parts])
+        except ValueError:
+            continue
+    return names, rows
+
+
+def on_grid(rows: list[list[float]], step: float, t_start: float) -> dict[int, list[float]]:
+    """The rows at whole multiples of step from t_start on, keyed by the multiple."""
+    out = {}
+    for row in rows:
+        k = round(row[0] / step)
+        if abs(row[0] - step * k) < 1e-6 and row[0] >= t_start - 1e-9:
+            out[k] = row
+    return out
+
+
+def turning_points(x: list[float]) -> list[float]:
+    """Peaks and valleys of x; a plateau of equal neighbours counts once."""
+    y = [v for i, v in enumerate(x) if i == 0 or v != x[i - 1]]
+    if len(y) < 3:
+        return y
+    out = [y[0]]
+    for i in range(1, len(y) - 1):
+        if (y[i] - y[i - 1]) * (y[i + 1] - y[i]) < 0.0:
+            out.append(y[i])
+    out.append(y[-1])
+    return out
+
+
+def rainflow(x: list[float]) -> list[tuple[float, float]]:
+    """(range, count) pairs: full cycles by the four-point method, the residue as half cycles."""
+    stack: list[float] = []
+    out: list[tuple[float, float]] = []
+    for v in turning_points(x):
+        stack.append(v)
+        while len(stack) >= 4:
+            s0, s1, s2, s3 = stack[-4:]
+            inner = abs(s2 - s1)
+            if inner <= abs(s1 - s0) and inner <= abs(s3 - s2):
+                out.append((inner, 1.0))
+                del stack[-3:-1]
+            else:
+                break
+    out.extend((abs(b - a), 0.5) for a, b in zip(stack[:-1], stack[1:]))
+    return out
+
+
+def damage_equivalent_load(x: list[float], duration: float, slope: float) -> float:
+    """DEL at a 1 Hz equivalent count over duration seconds; zero for a history without cycles
+    or without duration."""
+    total = sum(n * r ** slope for r, n in rainflow(x))
+    if total <= 0.0 or duration <= 0.0:
+        return 0.0
+    return (total / duration) ** (1.0 / slope)
+
+
+def relative_change(value: float, reference: float) -> float:
+    """100 (value/reference - 1); 0 when both are zero, NaN (undefined) when only the reference is."""
+    if reference == 0.0:
+        return 0.0 if value == 0.0 else math.nan
+    return 100.0 * (value / reference - 1.0)
+
+
+def mean_std(x: list[float]) -> tuple[float, float]:
+    m = sum(x) / len(x)
+    return m, math.sqrt(sum((v - m) ** 2 for v in x) / len(x))
+
+
+def compare_channel(grid: dict[int, list[float]], j: int, ref_grid: dict[int, list[float]], jr: int,
+                    step: float, slope: float | None) -> dict[str, float] | None:
+    """One channel of a run against the reference over the samples BOTH have: mean, std and max of
+    the run, the RMS of the difference over the reference std, and (slope given) the run's DEL and
+    its change against the reference DEL over the same samples and duration (relative_change).
+    covered is the fraction of the reference samples the run has; None when the two share fewer
+    than two samples."""
+    keys = sorted(set(grid) & set(ref_grid))
+    if len(keys) < 2:
+        return None
+    x = [grid[k][j] for k in keys]
+    xr = [ref_grid[k][jr] for k in keys]
+    m, s = mean_std(x)
+    _, sr = mean_std(xr)
+    out = {"mean": m, "std": s, "max": max(x), "covered": len(keys) / len(ref_grid),
+           "rel": math.sqrt(sum((p - q) ** 2 for p, q in zip(x, xr)) / len(x)) / max(sr, 1e-30)}
+    if slope is not None:
+        duration = (keys[-1] - keys[0]) * step
+        out["del"] = damage_equivalent_load(x, duration, slope)
+        out["ddel"] = relative_change(out["del"], damage_equivalent_load(xr, duration, slope))
+    return out
+
+
+def main() -> int:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("ref", help="root of the dtM = DT run (<root>.out and <root>.CD.out)")
+    ap.add_argument("runs", nargs="+", help="LABEL=ROOT of a run to compare (split at the last =)")
+    ap.add_argument("--t-start", type=float, default=200.0)
+    ap.add_argument("--grid", type=float, default=0.1)
+    ap.add_argument("--slope", type=float, default=3.0)
+    a = ap.parse_args()
+
+    runs = [("dtM = DT (reference)", a.ref)] + [tuple(r.rsplit("=", 1)) for r in a.runs]
+    tables = {}
+    for label, root in runs:
+        of_names, of_rows = read_table(Path(root + ".out"))
+        cd_names, cd_rows = read_table(Path(root + ".CD.out"))
+        tables[label] = (of_names, on_grid(of_rows, a.grid, a.t_start), cd_names,
+                         on_grid(cd_rows, a.grid, a.t_start), cd_rows)
+    ref_label = runs[0][0]
+    ref = tables[ref_label]
+    tension = [n for n in ref[2] if n.startswith(TENSION_PREFIXES)]
+    print(f"window t >= {a.t_start:g} s on the common {a.grid:g} s grid; DEL slope m = {a.slope:g}, "
+          "1 Hz equivalent count")
+    print(f"{'channel':<11} {'run':<28} {'mean':>12} {'std':>11} {'max':>12} {'RMSdiff/std':>11} "
+          f"{'DEL':>11} {'dDEL':>8}")
+    for ch, kind in [(c, 0) for c in PLATFORM] + [(c, 2) for c in tension]:
+        rnames, rgrid = ref[kind], ref[kind + 1]
+        if ch not in rnames:
+            continue
+        jr = rnames.index(ch)
+        for label, tab in tables.items():
+            names, grid = tab[kind], tab[kind + 1]
+            if ch not in names:
+                continue
+            c = compare_channel(grid, names.index(ch), rgrid, jr, a.grid, a.slope if kind == 2 else None)
+            if c is None:
+                print(f"{ch:<11} {label:<28} no common samples with the reference")
+                continue
+            line = f"{ch:<11} {label:<28} {c['mean']:12.5e} {c['std']:11.4e} {c['max']:12.5e} {c['rel']:11.4f}"
+            if kind == 2:
+                change = "    n/a" if math.isnan(c["ddel"]) else f"{c['ddel']:7.2f}%"
+                line += f" {c['del']:11.4e} {change}"
+            if c["covered"] < 1.0:
+                line += f"  (covers {100.0 * c['covered']:.1f} % of the reference samples; compared over those)"
+            print(line)
+        print()
+    print("reference DEL at its native sampling against the common grid (sampling alone):")
+    names, grid, rows = ref[2], ref[3], ref[4]
+    window = [r for r in rows if r[0] >= a.t_start - 1e-9]
+    if len(window) < 2 or len(grid) < 2:
+        print("  (fewer than two reference samples in the window)")
+        return 0
+    for ch in tension:
+        j = names.index(ch)
+        native = damage_equivalent_load([r[j] for r in window], window[-1][0] - window[0][0], a.slope)
+        keys = sorted(grid)
+        coarse = damage_equivalent_load([grid[k][j] for k in keys], (keys[-1] - keys[0]) * a.grid, a.slope)
+        print(f"  {ch:<9} native {native:.4e}  grid {coarse:.4e}")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())