diff --git a/CHANGELOG.md b/CHANGELOG.md index d357569..8dccd01 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/VALIDATION.md b/VALIDATION.md index e572b2f..99243cd 100644 --- a/VALIDATION.md +++ b/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 diff --git a/doc/openfast.rst b/doc/openfast.rst index 93828d6..5f579a3 100644 --- a/doc/openfast.rst +++ b/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. diff --git a/tests/test_supercycle_ab.py b/tests/test_supercycle_ab.py new file mode 100644 index 0000000..08a3dd6 --- /dev/null +++ b/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() diff --git a/validation/scripts/openfast_supercycle_ab.py b/validation/scripts/openfast_supercycle_ab.py new file mode 100644 index 0000000..fed2b80 --- /dev/null +++ b/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 .out): mean, std, max and the RMS of the time-series difference + relative to the reference std; + * fairlead and anchor tensions (CableDyn .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 (.out and .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())