From 9ee654fbc85d92b2de4df30a2b4f43a826d7df95 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 22 Oct 2026 11:24:23 +0900 Subject: [PATCH 1/3] docs(validation): record the supercycling verification The coupled time-step work (#38 to #40) lets CableDyn solve on a dtM larger than the OpenFAST DT. This records the verification of that change on the IEA-15 MW VolturnUS-S chain mooring in a DLC-1.6 environment (NTM 10.6 m/s turbulent wind, JONSWAP Hs 8 m / Tp 12 s, 3800 s), at dtM = 4 DT = 0.1 s against the same model at dtM = DT and against stock MoorDyn. With the coupled loads extrapolated between CableDyn solves, every platform DOF's time series differs from the dtM = DT run by at most 0.14 % of its standard deviation (0.41 % held), the tension statistics stay within 0.16 %, and the rainflow damage-equivalent loads change by at most 0.21 % at the fairleads and 1.5 % at the anchors, of the order of the reference's own sampling effect (0.94 %). Against stock MoorDyn the in-plane platform standard deviations agree within 0.03 %. validation/scripts/openfast_supercycle_ab.py computes these figures from the twin runs; its rainflow count (ASTM E1049, four-point method, half-cycle residue, a plateau of equal samples counted once) is tested against the standard's worked example. Signed-off-by: Jae Hoon Seo --- CHANGELOG.md | 4 + VALIDATION.md | 17 ++ doc/openfast.rst | 4 +- tests/test_supercycle_ab.py | 47 ++++++ validation/scripts/openfast_supercycle_ab.py | 169 +++++++++++++++++++ 5 files changed, 240 insertions(+), 1 deletion(-) create mode 100644 tests/test_supercycle_ab.py create mode 100644 validation/scripts/openfast_supercycle_ab.py 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..562fb75 --- /dev/null +++ b/tests/test_supercycle_ab.py @@ -0,0 +1,47 @@ +# 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 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) + + +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..6114c34 --- /dev/null +++ b/validation/scripts/openfast_supercycle_ab.py @@ -0,0 +1,169 @@ +#!/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. + +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.""" + total = sum(n * r ** slope for r, n in rainflow(x)) + return (total / duration) ** (1.0 / slope) + + +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 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) + del_ref = None + for label, tab in tables.items(): + names, grid = tab[kind], tab[kind + 1] + if ch not in names: + continue + j = names.index(ch) + keys = sorted(set(grid) & set(rgrid)) + x = [grid[k][j] for k in keys] + xr = [rgrid[k][jr] for k in keys] + m, s = mean_std(x) + _, sr = mean_std(xr) + rel = math.sqrt(sum((p - q) ** 2 for p, q in zip(x, xr)) / len(x)) / max(sr, 1e-30) + line = f"{ch:<11} {label:<28} {m:12.5e} {s:11.4e} {max(x):12.5e} {rel:11.4f}" + if kind == 2: + d = damage_equivalent_load(x, (keys[-1] - keys[0]) * a.grid, a.slope) + del_ref = d if del_ref is None else del_ref + line += f" {d:11.4e} {100.0 * (d / del_ref - 1.0):7.2f}%" + 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] + 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()) From 3e1909be314723adb02918fe3e4ab2bba52b6acf Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 22 Oct 2026 11:31:02 +0900 Subject: [PATCH 2/3] fix(validation): compare each run's DEL over its overlap with the reference Signed-off-by: Jae Hoon Seo --- tests/test_supercycle_ab.py | 21 ++++++++++ validation/scripts/openfast_supercycle_ab.py | 42 ++++++++++++++------ 2 files changed, 50 insertions(+), 13 deletions(-) diff --git a/tests/test_supercycle_ab.py b/tests/test_supercycle_ab.py index 562fb75..8407e6d 100644 --- a/tests/test_supercycle_ab.py +++ b/tests/test_supercycle_ab.py @@ -43,5 +43,26 @@ def test_damage_equivalent_load(self): 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) + + if __name__ == "__main__": unittest.main() diff --git a/validation/scripts/openfast_supercycle_ab.py b/validation/scripts/openfast_supercycle_ab.py index 6114c34..973a805 100644 --- a/validation/scripts/openfast_supercycle_ab.py +++ b/validation/scripts/openfast_supercycle_ab.py @@ -14,7 +14,10 @@ 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. + 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. @@ -106,6 +109,26 @@ def mean_std(x: list[float]) -> tuple[float, float]: 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]: + """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. covered is the + fraction of the reference samples the run has.""" + keys = sorted(set(grid) & set(ref_grid)) + 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"] = 100.0 * (out["del"] / damage_equivalent_load(xr, duration, slope) - 1.0) + 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)") @@ -134,23 +157,16 @@ def main() -> int: if ch not in rnames: continue jr = rnames.index(ch) - del_ref = None for label, tab in tables.items(): names, grid = tab[kind], tab[kind + 1] if ch not in names: continue - j = names.index(ch) - keys = sorted(set(grid) & set(rgrid)) - x = [grid[k][j] for k in keys] - xr = [rgrid[k][jr] for k in keys] - m, s = mean_std(x) - _, sr = mean_std(xr) - rel = math.sqrt(sum((p - q) ** 2 for p, q in zip(x, xr)) / len(x)) / max(sr, 1e-30) - line = f"{ch:<11} {label:<28} {m:12.5e} {s:11.4e} {max(x):12.5e} {rel:11.4f}" + c = compare_channel(grid, names.index(ch), rgrid, jr, a.grid, a.slope if kind == 2 else None) + 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: - d = damage_equivalent_load(x, (keys[-1] - keys[0]) * a.grid, a.slope) - del_ref = d if del_ref is None else del_ref - line += f" {d:11.4e} {100.0 * (d / del_ref - 1.0):7.2f}%" + line += f" {c['del']:11.4e} {c['ddel']:7.2f}%" + 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):") From 9b56a67de326bd2db80949de3a5adfbcfea40c21 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Thu, 22 Oct 2026 11:42:16 +0900 Subject: [PATCH 3/3] fix(validation): report a zero-damage reference instead of dividing by it Signed-off-by: Jae Hoon Seo --- tests/test_supercycle_ab.py | 17 +++++++++++ validation/scripts/openfast_supercycle_ab.py | 32 ++++++++++++++++---- 2 files changed, 43 insertions(+), 6 deletions(-) diff --git a/tests/test_supercycle_ab.py b/tests/test_supercycle_ab.py index 8407e6d..08a3dd6 100644 --- a/tests/test_supercycle_ab.py +++ b/tests/test_supercycle_ab.py @@ -4,6 +4,7 @@ from __future__ import annotations import importlib.util +import math import unittest from collections import Counter from pathlib import Path @@ -63,6 +64,22 @@ def test_a_scaled_run_changes_its_del_by_the_scale(self): 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 index 973a805..fed2b80 100644 --- a/validation/scripts/openfast_supercycle_ab.py +++ b/validation/scripts/openfast_supercycle_ab.py @@ -99,23 +99,36 @@ def rainflow(x: list[float]) -> list[tuple[float, float]]: def damage_equivalent_load(x: list[float], duration: float, slope: float) -> float: - """DEL at a 1 Hz equivalent count over duration seconds.""" + """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]: + 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. covered is the - fraction of the reference samples the run has.""" + 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) @@ -125,7 +138,7 @@ def compare_channel(grid: dict[int, list[float]], j: int, ref_grid: dict[int, li if slope is not None: duration = (keys[-1] - keys[0]) * step out["del"] = damage_equivalent_load(x, duration, slope) - out["ddel"] = 100.0 * (out["del"] / damage_equivalent_load(xr, duration, slope) - 1.0) + out["ddel"] = relative_change(out["del"], damage_equivalent_load(xr, duration, slope)) return out @@ -162,9 +175,13 @@ def main() -> int: 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: - line += f" {c['del']:11.4e} {c['ddel']:7.2f}%" + 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) @@ -172,6 +189,9 @@ def main() -> int: 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)