From cf672fdc8683fa227d40e6ba57e18ab6dc19efad Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sat, 7 Nov 2026 21:34:32 +0900 Subject: [PATCH 1/5] feat(waves): per-row directions and a reference point for wave records WaveKin 7 refused rows of different directions, though MoorDyn-C reads a beta per row of wave_frequencies.txt and gives each its own direction. Rows that differ now make a short-crested table in ascending frequency: CD_Component_Wave_Kinematics takes the direction cosines of each component (comp_cos, comp_sin) and sums x and y apart, in a branch of its own, so a long-crested table runs as before. Vessel RAOs take each row's heading. The beta of the omega = 0 row has no effect. The new options WaveRefX and WaveRefY give the point that an elevation record or a component file (WaterKin WaveKinMod 1, WaveKin 3 and 7) refers to: each phase uses the position relative to it, and a vessel RAO the vessel's. The default is the origin. They are an extension of CableDyn, documented as such, and a deck with any other sea, or none, is refused by name; the Python deck validator mirrors the rule. Existing decks give the same bits as before. Two directions equal the two superposed seas to 1e-15, at the depth stretched by the total elevation, across the quadrants; a shifted reference point equals the deck translated. Signed-off-by: Jae Hoon Seo --- CMakeLists.txt | 2 +- doc/driver_format.md | 14 +- doc/file_formats.rst | 3 +- doc/options.rst | 10 ++ doc/theory.rst | 13 ++ examples/cabledyn_options_reference.dat | 2 + python/cabledyn/deck_file.py | 29 ++- python/tests/test_deck_native_parity.py | 10 ++ src/CableDyn_DeckDriver.f90 | 73 ++++++-- src/CableDyn_Hydro.f90 | 143 ++++++++++++--- tests/test_vessel_motion.f90 | 65 +++++++ tests/test_wave_component_sum.f90 | 86 ++++++++- tests/test_wavekin_modes.f90 | 223 +++++++++++++++++++++++- 13 files changed, 613 insertions(+), 60 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index d337db75..29d48c7f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1166,7 +1166,7 @@ set_tests_properties(cable_statics_modes PROPERTIES LABELS "fortran;driver;hermi # clamped hang-off moment under vessel pitch. add_executable(test_vessel_motion tests/test_vessel_motion.f90) target_link_libraries(test_vessel_motion PRIVATE cabledyn_core) -file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion) +file(MAKE_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion/vm_wk7) add_test(NAME vessel_motion COMMAND test_vessel_motion WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/vessel_motion) set_tests_properties(vessel_motion PROPERTIES LABELS "fortran;unit;dynamic;validation") diff --git a/doc/driver_format.md b/doc/driver_format.md index 922859b0..5dcc1a2e 100644 --- a/doc/driver_format.md +++ b/doc/driver_format.md @@ -1329,7 +1329,7 @@ every node's position and time as MoorDyn-C interpolates them. Standalone decks | Option | File | Meaning | |---|---|---| | `WaveKin 3` | `wave_elevation.txt` | `time elevation` rows from `t = 0` (no header rows), resampled linearly at `dtWave` to `floor(t_last/dtWave) + 1` samples (the last dropped when the count is odd) and reduced to its Fourier components; components above `WvHiCOff` (default 0.5 Hz, MoorDyn-C's fixed cut) are dropped and the mean is a mean-level offset, as MoorDyn-C does. Waves travel along +x | -| `WaveKin 7` | `wave_frequencies.txt` | `omega Re Im [beta]` rows [rad/s, m, m, rad], the first at `omega = 0`: each row is the component `\|c\| sin(omega t - k d + arg c)`, `d = x cos(beta) + y sin(beta)`; the `omega = 0` row is a mean-level offset. All rows must share one `beta` | +| `WaveKin 7` | `wave_frequencies.txt` | `omega Re Im [beta]` rows [rad/s, m, m, rad], the first at `omega = 0`: each row is the component `\|c\| sin(omega t - k d + arg c)`, `d = x cos(beta) + y sin(beta)` with the row's own `beta` (0 when the column is left out); the `omega = 0` row is a mean-level offset and its `beta` has no effect. Rows of different directions describe a short-crested sea, summed in ascending frequency | | `Currents 1` | `current_profile.txt` | `z ux uy uz` rows after its header rows, strictly increasing `z`; linear in `z` and held beyond the end rows | | `Currents 2` | `current_profile_dynamic.txt` | a depth profile in time: lines 1-4 and 6 are free text, line 5 lists the depths `z1 ... zn` (strictly increasing, read in single precision as MoorDyn-C does), and each later row is `t ux1 ... uxn [uy1 ... uyn [uz1 ... uzn]]` (missing components zero), at least two rows with strictly increasing `t` | | `Currents 5` | `current_profile_4d.txt` | a grid in x, y, z and time: line 1 is `nx ny nz nt`, lines 2-5 the x, y, z and t coordinates (each strictly increasing, at least two times), then one `x y z t ux uy uz` row per grid point, in any order, matched to the grid by value; a text line 6 is a header | @@ -1368,16 +1368,20 @@ where it is sampled on `water_grid.txt`; CableDyn adds the current at the nodes, Rejected by name: `WaveKin 1` (kinematics set node by node through the MoorDyn-C API; a deck has no source for them), `WaveKin 2` (the FFT grid: `wave_frequencies.txt` resampled to an even spectrum; `WaveKin 7` evaluates the same file's components exactly), `WaveKin 4`-`6` and -`Currents 3`-`4` (no kinematics source in MoorDyn-C v2.7.1), and a `wave_frequencies.txt` whose -rows have different directions. A deck that also declares a `waves` option, a WaterKin file, or a -second current source is rejected as double counting. On shared inputs the `wavekin_modes` test -holds `WaveKin 7` and `Currents 1` identical to MoorDyn-C v2.7.1 (`5e-12` relative) and `WaveKin 3` +`Currents 3`-`4` (no kinematics source in MoorDyn-C v2.7.1). A deck that also declares a `waves` +option, a WaterKin file, or a second current source is rejected as double counting. On shared +inputs the `wavekin_modes` test holds `WaveKin 7` and `Currents 1` identical to MoorDyn-C v2.7.1 (`5e-12` relative) and `WaveKin 3` within `3e-7` m in elevation and `1e-4` in velocity and acceleration at MoorDyn-C's grid points, the remainder being MoorDyn-C's grid interpolation and its approximate wave number. The `mdc_currents` test holds `Currents 2` and `5` identical to MoorDyn-C v2.7.1 at 739 points and times (including past the last sample), and `mdc_currents_runs` matches the tension change of a time-varying and of a sheared current to MoorDyn-C's within 0.9 % on a 70 m rope. +`wave_frequencies.txt` is the file MoorDyn-C v2.7.1 reads, a `beta` per row included: its +`SUM_COMPONENTS` mode gives each row its own direction too. `WaveRefX` and `WaveRefY` +([options](options.rst)) are an extension of CableDyn and no part of MoorDyn-C or MoorDyn-F: a +deck with them is not a plain MoorDyn deck, and a MoorDyn solver would not place the phases there. + (waterkin-file-modes)= ### MoorDyn-F `WaterKin` file modes diff --git a/doc/file_formats.rst b/doc/file_formats.rst index eeb19ded..65232a19 100644 --- a/doc/file_formats.rst +++ b/doc/file_formats.rst @@ -549,7 +549,8 @@ two columns: * - 2 - ``elevation`` - m - - Free-surface elevation η at the reference location. + - Free-surface elevation η at the reference location: the origin, or the point of the + ``WaveRefX`` and ``WaveRefY`` options, an extension of CableDyn. * Blank lines are skipped. Lines before the first data row whose first token is not a number, such as a ``time elevation`` header, are skipped. Once data has diff --git a/doc/options.rst b/doc/options.rst index 10750981..84cd525c 100644 --- a/doc/options.rst +++ b/doc/options.rst @@ -761,6 +761,16 @@ Ambient fluid and prescribed motion component at or below the cut-off is refused - standalone - ``6.2832 WvHiCOff`` + * - ``WaveRefX``, ``WaveRefY`` + - number [m], at most 10\ :sup:`6` in magnitude + - ``0`` + - the horizontal point that the phases of an elevation record or a component file + (WaterKin ``WaveKinMod 1``, ``WaveKin 3`` and ``7``) refer to: the record was measured + there, and each phase uses the position relative to it (:ref:`theory-waves`). An + extension of CableDyn, which neither MoorDyn-C nor MoorDyn-F reads. A deck with any other + sea, or none, is refused by name + - standalone + - ``250.0 WaveRefX`` * - ``Currents`` - ``0``, ``1``, ``2`` or ``5`` - ``0`` diff --git a/doc/theory.rst b/doc/theory.rst index 2d006363..8b9f9ed0 100644 --- a/doc/theory.rst +++ b/doc/theory.rst @@ -1112,6 +1112,19 @@ commensurate, the record does not repeat (an equally spaced comb would repeat af its component table (amplitude, frequency, phase), for example a WaterKin ``WaveKinMod`` 1 elevation history, which is reduced once to its Fourier components on the run time base. +The phase of a component of an elevation history or a ``WaveKin 7`` table is + +.. math:: + + \theta_i = k_i\big[(x - x_r)\cos\beta_i + (y - y_r)\sin\beta_i\big] - \omega_i t + \phi_i, + +with :math:`(x_r, y_r)` the ``WaveRefX``, ``WaveRefY`` point where the record was measured, the +origin by default. The sea of a record referred to a point is therefore the same sea with the +deck translated, and a vessel RAO takes the vessel's position relative to that point. A history +has one direction :math:`\beta_i = \beta`. The rows of a ``WaveKin 7`` table may each have +their own: a short-crested sea, whose elevation and kinematics are the sums of those of its +directions, stretched by the one total elevation. + Every point sums every component of the elevation. The sums run in interleaved partial sums added up in a fixed order, so a result does not depend on the vector width or the thread count, and each point is a function of its position and time alone: nothing carries from one diff --git a/examples/cabledyn_options_reference.dat b/examples/cabledyn_options_reference.dat index 85a57e0a..b9fc3128 100644 --- a/examples/cabledyn_options_reference.dat +++ b/examples/cabledyn_options_reference.dat @@ -99,6 +99,8 @@ none waves - Wave model and parameters (m, s, deg) {none; airy; -- 0 Currents - MoorDyn-C current mode {0=none; 1=current_profile.txt} [default: 0] -- 0.25 dtWave - MoorDyn-C WaveKin 3 resampling step (s) [default: 0.25] -- 3.14159 WvHiCOff - High cut-off of an elevation record, WaveKinMod 1 or WaveKin 3 (rad/s) {value; none} [default: 3.14159] +-- 0.0 WaveRefX - x of the point an elevation record or component file refers to (m) [default: 0] +-- 0.0 WaveRefY - y of the point an elevation record or component file refers to (m) [default: 0] -- Keyword aliases (equivalent spellings, case-insensitive): -- g = gravity; rhoW = WtrDnsty = water_density; WtrDpth = water_depth; kBot = kb; cBot = cb -- bathymetryFile = bathymetry_file = seafloorFile = seafloor_file; frictionMu = mu = frictionCoefficient diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index c2a54e18..107a242e 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -277,6 +277,8 @@ "bodysubstep", "dtwave", "wvhicoff", + "waverefx", + "waverefy", "bodywetting", "bodyhydro", "rodhydro", @@ -1229,6 +1231,8 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None: "dtout", "dtwave", "wvhicoff", + "waverefx", + "waverefy", "inertialf", "inertialf_rampt", }: @@ -2428,8 +2432,9 @@ def validate(self) -> None: and the MoorDyn-C kinematics files) are not read here, so the checks that depend on their contents stay with the solver: for example the anchor height on a bathymetry surface, whether a WaterKin file carries - waves (``WaveKinMod``) for ``vesselRAO``, and whether its current - doubles a ``current`` row or ``Currents 1``, ``2`` or ``5``. + waves (``WaveKinMod``) for ``vesselRAO``, ``WaveRefX`` or ``WaveRefY``, + and whether its current doubles a ``current`` row or ``Currents 1``, + ``2`` or ``5``. Raises ------ @@ -2675,6 +2680,26 @@ def effective(key: str, default: float) -> float: raise DeckFormatError( f"{self._label}: WaveKin {mdc_wave_kin} waves require dtM, TMax and WtrDpth" ) + # native validation: WaveRefX and WaveRefY place the phases of component waves alone. A + # WaterKin file may carry them (its WaveKinMod is read at run time). + wave_ref_rows = [o for o in options if o.keyword.lower() in {"waverefx", "waverefy"}] + if wave_ref_rows: + if any(abs(_native_float(o.values[0])) > _NATIVE_MAX_COORDINATE for o in wave_ref_rows): + raise DeckFormatError( + f"{self._label}: WaveRefX and WaveRefY must be at most 1e6 m in magnitude" + ) + waterkin_ref = option_by_key.get("waterkin") + waterkin_text = "" if waterkin_ref is None else waterkin_ref.values[0] + has_record = bool(mdc_wave_kin) or ( + waterkin_text.lower() not in {"", "none", "seastate"} + and any(character in _WATERKIN_FILENAME_LETTERS for character in waterkin_text) + ) + if self.caller_driven or has_waves or trains or not has_record: + raise DeckFormatError( + f"{self._label}: WaveRefX and WaveRefY place the reference point of an " + "elevation record or a component file (WaterKin WaveKinMod 1, WaveKin 3 or 7); " + "this deck has none" + ) if mdc_currents: if current is not None and current.values[0].lower() != "none": raise DeckFormatError( diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index 4f2ad44e..66c529ff 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -105,6 +105,16 @@ def test_dtwave_must_be_positive(value): (("uniform 0.1 0 0 current", "1 Currents"), "both a current OPTION and Currents 1"), (("uniform 0.1 0 0 current", "2 Currents"), "both a current OPTION and Currents 2"), (("uniform 0.1 0 0 current", "5 Currents"), "both a current OPTION and Currents 5"), + (("7 WaterKin", "12.5 WaveRefX", "-3 WaveRefY"), None), + (("3 WaveKin", "0.1 dtWave", "0 WaveRefY"), None), + (("wk.dat WaterKin", "12.5 WaveRefX"), None), + (("12.5 WaveRefX",), "reference point of an elevation record"), + (("0 WaveRefY",), "reference point of an elevation record"), + (("airy 2 8 0 waves", "1 WaveRefX"), "reference point of an elevation record"), + (("airy 2 8 0 wavetrain", "1 WaveRefY"), "reference point of an elevation record"), + (("7 WaterKin", "2e6 WaveRefX"), "at most 1e6 m"), + (("7 WaterKin", "nan WaveRefY"), "must be"), + (("7 WaterKin", "north WaveRefY"), "must be numeric"), ], ) def test_moordyn_c_kinematics_follow_the_native_rules(rows, message): diff --git a/src/CableDyn_DeckDriver.f90 b/src/CableDyn_DeckDriver.f90 index d4c02c2c..9c43c028 100644 --- a/src/CableDyn_DeckDriver.f90 +++ b/src/CableDyn_DeckDriver.f90 @@ -664,7 +664,14 @@ MODULE CableDyn_DeckDriver REAL(wp), ALLOCATABLE :: wave_component_amplitude(:), wave_component_phase(:) ! exp(-2 k (wtrdpth + wave_eta_mean)) and the tail_sums of the table (CD_Component_Wave_Kinematics) REAL(wp), ALLOCATABLE :: wave_component_e2kd(:), wave_component_tail(:, :) + ! a short-crested table (WaveKin 7 rows of different directions): each component's direction + ! cosine and sine, in ascending frequency; unallocated when all share wave_direction + REAL(wp), ALLOCATABLE :: wave_component_cos(:), wave_component_sin(:) REAL(wp) :: wave_eta_mean = CD_ZERO + ! OPTION WaveRefX, WaveRefY: the point the phases of an elevation record or a component file + ! refer to (the record was measured there); the origin by default + REAL(wp) :: wave_ref_x = CD_ZERO, wave_ref_y = CD_ZERO + LOGICAL :: has_wave_ref = .FALSE. ! N-level vertical current profile (ascending z): from the two-point deck ! OPTION row or a MoorDyn WaterKin file's CurrentMod-1 depth table. REAL(wp), ALLOCATABLE :: current_profile_z(:) @@ -13802,8 +13809,10 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg) !! the mean-level offset, as MoorDyn-C does. Waves travel along +x. !! - WaveKin 7, wave_frequencies.txt (`omega Re Im [beta]` rows, the first at omega = 0): !! each row is the component |c| sin(omega t - k d + arg c) of MoorDyn-C's SUM_COMPONENTS - !! mode, d = x cos(beta) + y sin(beta); the omega = 0 row is the mean-level offset. The rows - !! must share one direction beta. + !! mode, d = x cos(beta) + y sin(beta) with the row's own direction beta; the omega = 0 row + !! is the mean-level offset and its beta has no effect. Rows that share one direction are + !! a long-crested table in file order; rows that differ a short-crested table in ascending + !! frequency, the order the tail bound of the sums needs. !! - Currents 1, current_profile.txt (`z ux uy uz` rows after its header rows): the steady !! profile, linear in z and held beyond its ends, as MoorDyn-C's current grid. !! - Currents 2, current_profile_dynamic.txt (constructDynamicCurrentGrid of MoorDyn-C's @@ -13823,9 +13832,11 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg) TYPE(DeckOptions), INTENT(INOUT) :: opts INTEGER, INTENT(OUT) :: ErrStat CHARACTER(*), INTENT(OUT) :: ErrMsg - REAL(wp), ALLOCATABLE :: tab(:, :), eta(:), comp_w(:), comp_a(:), comp_p(:) + REAL(wp), ALLOCATABLE :: tab(:, :), eta(:), comp_w(:), comp_a(:), comp_p(:), comp_b(:) REAL(wp) :: dt, tl, frac, ca, sa, angle, scale, beta, mean_level, arg INTEGER :: n, nt, j, k, nc, bracket + INTEGER, ALLOCATABLE :: order(:) + LOGICAL :: one_direction INTEGER(KIND=8) :: kj CHARACTER(32) :: fname ! the longest row of a MoorDyn-C current table (a time and 3 nz velocities per row) @@ -13933,15 +13944,15 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg) IF (ABS(tab(1, 1)) > CD_ZERO) THEN CALL fail(ErrStat, ErrMsg, 'WaveKin 7: the first wave_frequencies.txt row must be at omega = 0'); RETURN END IF + ! the omega = 0 row has no direction of its own: the first component's stands for the table beta = tab(4, 1) + IF (ABS(tab(4, 2) - beta) > 1.0e-12_wp) beta = tab(4, 2) + one_direction = .TRUE. DO j = 1, n IF (ABS(tab(4, j)) > TWO_PI) THEN CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt directions are radians within +-2 pi'); RETURN END IF - IF (ABS(tab(4, j) - beta) > 1.0e-12_wp) THEN - CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt rows with different directions '// & - '(short-crested seas) are not supported; give every row one direction'); RETURN - END IF + IF (j > 1 .AND. ABS(tab(4, j) - beta) > 1.0e-12_wp) one_direction = .FALSE. IF (j > 1) THEN IF (.NOT. tab(1, j) > CD_ZERO) THEN CALL fail(ErrStat, ErrMsg, 'WaveKin 7: wave_frequencies.txt frequencies after the first row '// & @@ -13954,6 +13965,14 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg) comp_w = tab(1, 2:n) comp_a = SQRT(tab(2, 2:n)**2 + tab(3, 2:n)**2) comp_p = 0.25_wp*TWO_PI - ATAN2(tab(3, 2:n), tab(2, 2:n)) + IF (.NOT. one_direction) THEN + ALLOCATE (order(n - 1)) + CALL sort_index(comp_w, order) + comp_w = comp_w(order) + comp_a = comp_a(order) + comp_p = comp_p(order) + comp_b = tab(4, order + 1) + END IF CALL install_components() IF (ErrStat /= CD_DECKDRV_OK) RETURN END SELECT @@ -14378,7 +14397,7 @@ END FUNCTION real_to_str SUBROUTINE install_components() !! The component waves (omega, amplitude, phase in CableDyn's cos(k d - omega t + phase) - !! convention), their direction beta [rad] and the mean level. + !! convention), their direction beta [rad], or each their own comp_b, and the mean level. IF (.NOT. (ALL(IEEE_IS_FINITE(comp_w)) .AND. ALL(IEEE_IS_FINITE(comp_a)) .AND. & ALL(IEEE_IS_FINITE(comp_p)) .AND. IEEE_IS_FINITE(mean_level))) THEN CALL fail(ErrStat, ErrMsg, 'WaveKin '//TRIM(int_to_str(opts%mdc_wave_kin))//': '//TRIM(fname)// & @@ -14400,6 +14419,11 @@ SUBROUTINE install_components() opts%wave_component_phase = comp_p opts%wave_eta_mean = mean_level opts%wave_direction = beta*45.0_wp/ATAN(CD_ONE) + IF (ALLOCATED(opts%wave_component_cos)) DEALLOCATE (opts%wave_component_cos, opts%wave_component_sin) + IF (ALLOCATED(comp_b)) THEN + opts%wave_component_cos = COS(comp_b) + opts%wave_component_sin = SIN(comp_b) + END IF opts%wave_kind = WAVE_COMPONENTS opts%has_wave = .TRUE. END SUBROUTINE install_components @@ -16620,6 +16644,12 @@ SUBROUTINE apply_option(line, opts, ErrStat, ErrMsg) ELSE CALL rd(val, opts%wv_hi_cutoff, ErrStat, ErrMsg) END IF + CASE ('waverefx') + CALL rd(val, opts%wave_ref_x, ErrStat, ErrMsg) + opts%has_wave_ref = .TRUE. + CASE ('waverefy') + CALL rd(val, opts%wave_ref_y, ErrStat, ErrMsg) + opts%has_wave_ref = .TRUE. CASE ('ramptime', 'ramp_time', 'tramp') CALL rd(val, opts%ramp_time, ErrStat, ErrMsg) CASE ('maxstrain', 'max_strain') @@ -17124,7 +17154,7 @@ SUBROUTINE check_option_row_value(kw, val, ErrStat, ErrMsg) rule = 'non-negative' CASE ('aslaidtension', 'as_laid_tension') rule = 'non-negative' - CASE ('layazimuth', 'lay_azimuth') + CASE ('layazimuth', 'lay_azimuth', 'waverefx', 'waverefy') rule = 'finite' CASE ('rhoinf', 'rho_inf') rule = 'unit' @@ -17758,6 +17788,15 @@ SUBROUTINE finalize_checks(types, bodies, rod_types, rods, points, lines, sectio END IF END IF END IF + IF (opts%has_wave_ref) THEN + IF (ABS(opts%wave_ref_x) > MAX_COORDINATE .OR. ABS(opts%wave_ref_y) > MAX_COORDINATE) THEN + CALL fail(ErrStat, ErrMsg, 'WaveRefX and WaveRefY must be at most 1e6 m in magnitude'); RETURN + END IF + IF (.NOT. (opts%has_wave .AND. opts%wave_kind == WAVE_COMPONENTS)) THEN + CALL fail(ErrStat, ErrMsg, 'WaveRefX and WaveRefY place the reference point of an elevation record or a '// & + 'component file (WaterKin WaveKinMod 1, WaveKin 3 or 7); this deck has none'); RETURN + END IF + END IF IF (opts%has_wave) THEN IF (.NOT. opts%has_dtm) THEN CALL fail(ErrStat, ErrMsg, 'deck waves require dtM and TMax'); RETURN @@ -32232,21 +32271,23 @@ SUBROUTINE deck_wave_kinematics(opts, x, y, z, time, eta, velocity, acceleration pdyn) CASE (WAVE_COMPONENTS) IF (ramp < CD_ONE) THEN - CALL CD_Component_Wave_Kinematics(x, y, z - opts%wave_eta_mean, time, & + CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - opts%wave_eta_mean, time, & opts%wtrdpth + opts%wave_eta_mean, opts%wave_direction, .TRUE., & opts%wave_component_omega, opts%wave_component_k, & opts%wave_component_amplitude, opts%wave_component_phase, & eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated=.TRUE., & amplitude_scale=ramp, pdyn=pdyn, exp_m2kd=opts%wave_component_e2kd, & - tail_sums=opts%wave_component_tail) + tail_sums=opts%wave_component_tail, comp_cos=opts%wave_component_cos, & + comp_sin=opts%wave_component_sin) ELSE - CALL CD_Component_Wave_Kinematics(x, y, z - opts%wave_eta_mean, time, & + CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - opts%wave_eta_mean, time, & opts%wtrdpth + opts%wave_eta_mean, opts%wave_direction, .TRUE., & opts%wave_component_omega, opts%wave_component_k, & opts%wave_component_amplitude, opts%wave_component_phase, & eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated=.TRUE., & pdyn=pdyn, exp_m2kd=opts%wave_component_e2kd, & - tail_sums=opts%wave_component_tail) + tail_sums=opts%wave_component_tail, comp_cos=opts%wave_component_cos, & + comp_sin=opts%wave_component_sin) END IF eta = eta + opts%wave_eta_mean CASE DEFAULT @@ -35097,6 +35138,10 @@ SUBROUTINE synthesize_vessel_rao(path, opts, motion, rref, vref, aref, ErrStat, CASE (WAVE_COMPONENTS) ALLOCATE (om, SOURCE=opts%wave_component_omega); ALLOCATE (kk, SOURCE=opts%wave_component_k) ALLOCATE (amp, SOURCE=opts%wave_component_amplitude); ALLOCATE (phs, SOURCE=opts%wave_component_phase) + IF (ALLOCATED(opts%wave_component_cos)) THEN + ALLOCATE (dir(SIZE(om))) + dir = ATAN2(opts%wave_component_sin, opts%wave_component_cos)*R2D + END IF CASE (WAVE_SEA) nc = opts%sea%n_comp ALLOCATE (om, SOURCE=opts%sea%omega(1:nc)); ALLOCATE (kk, SOURCE=opts%sea%k(1:nc)) @@ -35122,7 +35167,7 @@ SUBROUTINE synthesize_vessel_rao(path, opts, motion, rref, vref, aref, ErrStat, DO i = 1, nc cb = COS(dir(i)*QUARTER_PI/45.0_wp) sb = SIN(dir(i)*QUARTER_PI/45.0_wp) - eps(i) = kk(i)*(opts%vessel_ref(1)*cb + opts%vessel_ref(2)*sb) + phs(i) + eps(i) = kk(i)*((opts%vessel_ref(1) - opts%wave_ref_x)*cb + (opts%vessel_ref(2) - opts%wave_ref_y)*sb) + phs(i) IF (.NOT. (amp(i) > CD_ZERO)) CYCLE CALL CD_Vessel_RAO_Eval(rao, TWO_PI/om(i), dir(i), cre(:, i), cim(:, i), clamped, es, em) IF (es /= CD_VESSEL_OK) THEN diff --git a/src/CableDyn_Hydro.f90 b/src/CableDyn_Hydro.f90 index 487e5c5b..2a6789ad 100644 --- a/src/CableDyn_Hydro.f90 +++ b/src/CableDyn_Hydro.f90 @@ -2239,7 +2239,7 @@ END SUBROUTINE CD_JONSWAP_Wave_Kinematics_Precomputed SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretch, & omega_comp, k_comp, amplitude, phase, & eta, velocity, acceleration, ErrStat, ErrMsg, inputs_validated, & - amplitude_scale, pdyn, exp_m2kd, tail_sums) + amplitude_scale, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin) !! Evaluate linear irregular-wave kinematics from a CALLER-SUPPLIED component table: !! the general (arbitrary length, explicit per-component phase) form of the fixed-size !! deterministic-phase JONSWAP evaluator above, sharing its conventions exactly -- @@ -2275,6 +2275,9 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc !! once that bound is below half an ulp of the magnitude already summed: a deep point drops !! the short components that have decayed to nothing there, and the result changes at !! round-off only. The elevation always sums every component. + !! comp_cos and comp_sin, given together, hold the cosine and sine of each component's own + !! direction, in place of the one direction_deg of a long-crested table: a short-crested + !! table. The bound of tail_sums holds for any direction, since it bounds the magnitude. REAL(wp), INTENT(IN) :: x, y, z, t, depth, direction_deg LOGICAL, INTENT(IN) :: stretch REAL(wp), INTENT(IN) :: omega_comp(:), k_comp(:), amplitude(:), phase(:) @@ -2286,7 +2289,7 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc ! the linear dynamic pressure per unit density (sum over the components of ! (omega^2/k) a cosh(k(z+h))/sinh(kh) cos(theta)), as in CD_Airy_Wave_Kinematics_Precomputed REAL(wp), INTENT(OUT), OPTIONAL :: pdyn - REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(:), tail_sums(:, :) + REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(:), tail_sums(:, :), comp_cos(:), comp_sin(:) INTEGER :: n REAL(wp) :: cb, sb, scl LOGICAL :: trust_inputs @@ -2316,6 +2319,22 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc RETURN END IF END IF + IF (PRESENT(comp_cos) .NEQV. PRESENT(comp_sin)) THEN + CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must be given together') + RETURN + END IF + IF (PRESENT(comp_cos)) THEN + IF (SIZE(comp_cos) /= n .OR. SIZE(comp_sin) /= n) THEN + CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must match the component-wave arrays') + RETURN + END IF + IF (.NOT. trust_inputs) THEN + IF (.NOT. (CD_All_Finite(comp_cos) .AND. CD_All_Finite(comp_sin))) THEN + CALL fail(ErrStat, ErrMsg, 'comp_cos and comp_sin must be finite') + RETURN + END IF + END IF + END IF ! Fortran does not require short-circuit evaluation: keep the O(n) scans inside a ! separate branch so the setup-validated hot path portably skips them at -O0 too. IF (.NOT. trust_inputs) THEN @@ -2345,7 +2364,7 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc ! a unit scale is exact, so the unscaled table is reproduced bit for bit scl = CD_ONE IF (PRESENT(amplitude_scale)) scl = amplitude_scale - IF (n == 1 .AND. ABS(phase(1)) <= CD_ZERO .AND. scl*amplitude(1) > CD_ZERO) THEN + IF (n == 1 .AND. ABS(phase(1)) <= CD_ZERO .AND. scl*amplitude(1) > CD_ZERO .AND. .NOT. PRESENT(comp_cos)) THEN CALL CD_Airy_Wave_Kinematics_Precomputed(x, y, z, t, 2.0_wp*(scl*amplitude(1)), omega_comp(1), k_comp(1), & depth, direction_deg, stretch, eta, velocity, acceleration, & ErrStat, ErrMsg, pdyn) @@ -2353,25 +2372,28 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc END IF cb = COS(direction_deg*DEG2RAD) sb = SIN(direction_deg*DEG2RAD) - CALL component_sums(n, x*cb + y*sb, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, & - velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums) + CALL component_sums(n, x, y, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, & + velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin) END SUBROUTINE CD_Component_Wave_Kinematics - SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, & - velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums) + SUBROUTINE component_sums(n, x, y, z, t, depth, stretch, scl, cb, sb, omega_comp, k_comp, amplitude, phase, eta, & + velocity, acceleration, ErrStat, ErrMsg, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin) !! The sums of CD_Component_Wave_Kinematics, on explicit-shape tables: their blocks pass to - !! the block kernels without a copy, and the loops run at unit stride. + !! the block kernels without a copy, and the loops run at unit stride. A long-crested table + !! forms each phase from the one distance along its direction and sums the horizontal + !! velocity and acceleration once, for cb and sb to multiply; a short-crested one (comp_cos, + !! comp_sin) forms each phase from the component's own direction and sums x and y apart. INTEGER, INTENT(IN) :: n - REAL(wp), INTENT(IN) :: x_along, z, t, depth, scl, cb, sb + REAL(wp), INTENT(IN) :: x, y, z, t, depth, scl, cb, sb LOGICAL, INTENT(IN) :: stretch REAL(wp), INTENT(IN) :: omega_comp(n), k_comp(n), amplitude(n), phase(n) REAL(wp), INTENT(INOUT) :: eta, velocity(3), acceleration(3) INTEGER, INTENT(INOUT) :: ErrStat CHARACTER(*), INTENT(INOUT) :: ErrMsg REAL(wp), INTENT(INOUT), OPTIONAL :: pdyn - REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(n), tail_sums(3, n) + REAL(wp), INTENT(IN), OPTIONAL :: exp_m2kd(n), tail_sums(3, n), comp_cos(n), comp_sin(n) INTEGER :: nb, n_small, i0, i1, j, nk, m, lanes - REAL(wp) :: z_eval + REAL(wp) :: z_eval, x_along ! the velocity sums run in blocks of NLANE and the elevation in chunks of NCHUNK (NBUF a ! multiple of NCHUNK and NCHUNK of NLANE, so a block lies wholly inside or past the buffer) INTEGER, PARAMETER :: NBUF = 4096, NCHUNK = 512 @@ -2379,10 +2401,12 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c REAL(wp) :: denom, e2, first_ratio REAL(wp) :: cth(NBUF), sth(NBUF), th(NCHUNK), bc(NCHUNK), bs(NCHUNK), ch(NLANE), sh(NLANE), e_bed(NLANE) REAL(wp) :: lane_eta(NLANE), lane_vc(NLANE), lane_vs(NLANE), lane_as(NLANE), lane_ac(NLANE) - REAL(wp) :: lane_p(NLANE), lane_m(NLANE, 3), sum_vc, sum_vs, sum_as, sum_ac + REAL(wp) :: lane_p(NLANE), lane_m(NLANE, 3), lane_vy(NLANE), lane_ay(NLANE), sum_vc, sum_vs, sum_as, sum_ac REAL(wp) :: mag(3), rest(3) - LOGICAL :: at_bed, use_tail, want_p + LOGICAL :: at_bed, use_tail, want_p, spread + spread = PRESENT(comp_cos) + x_along = x*cb + y*sb nb = MIN(n, NBUF) lanes = MIN(n, NLANE) ! the elevation sums every component (Wheeler stretching needs the total eta); the phase @@ -2391,9 +2415,15 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c DO i0 = 1, n, NCHUNK nk = MIN(n - i0 + 1, NCHUNK) i1 = i0 + nk - 1 - DO j = 1, nk - th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1) - END DO + IF (spread) THEN + DO j = i0, i1 + th(j - i0 + 1) = k_comp(j)*(x*comp_cos(j) + y*comp_sin(j)) - omega_comp(j)*t + phase(j) + END DO + ELSE + DO j = 1, nk + th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1) + END DO + END IF IF (i0 <= nb) THEN CALL CD_Wave_SinCos(nk, th, sth(i0), cth(i0)) CALL add_elevation_lanes(nk, amplitude(i0:i1), cth(i0), lane_eta) @@ -2429,6 +2459,10 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c lane_ac(1:lanes) = CD_ZERO lane_p(1:lanes) = CD_ZERO lane_m(1:lanes, :) = CD_ZERO + IF (spread) THEN + lane_vy(1:lanes) = CD_ZERO + lane_ay(1:lanes) = CD_ZERO + END IF first_ratio = CD_ZERO DO i0 = 1, n, NLANE nk = MIN(n - i0 + 1, NLANE) @@ -2479,20 +2513,42 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c END DO END IF END IF - IF (i0 <= nb) THEN + IF (i0 > nb) THEN + IF (spread) THEN + DO j = i0, i1 + th(j - i0 + 1) = k_comp(j)*(x*comp_cos(j) + y*comp_sin(j)) - omega_comp(j)*t + phase(j) + END DO + ELSE + DO j = 1, nk + th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1) + END DO + END IF + CALL CD_Wave_SinCos(nk, th, bs, bc) + END IF + IF (spread) THEN + IF (i0 <= nb) THEN + CALL add_spread_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), comp_cos(i0:i1), comp_sin(i0:i1), ch, sh, & + cth(i0), sth(i0), lane_vc, lane_vy, lane_vs, lane_as, lane_ay, lane_ac, lane_m(1, 1), & + lane_m(1, 2)) + ELSE + CALL add_spread_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), comp_cos(i0:i1), comp_sin(i0:i1), ch, sh, & + bc, bs, lane_vc, lane_vy, lane_vs, lane_as, lane_ay, lane_ac, lane_m(1, 1), & + lane_m(1, 2)) + END IF + ELSE IF (i0 <= nb) THEN CALL add_velocity_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), ch, sh, cth(i0), sth(i0), lane_vc, & lane_vs, lane_as, lane_ac, lane_m(1, 1), lane_m(1, 2)) - IF (want_p) CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, cth(i0), & - lane_p, lane_m(1, 3)) ELSE - DO j = 1, nk - th(j) = k_comp(i0 + j - 1)*x_along - omega_comp(i0 + j - 1)*t + phase(i0 + j - 1) - END DO - CALL CD_Wave_SinCos(nk, th, bs, bc) CALL add_velocity_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), ch, sh, bc, bs, lane_vc, lane_vs, lane_as, & lane_ac, lane_m(1, 1), lane_m(1, 2)) - IF (want_p) CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, bc, lane_p, & - lane_m(1, 3)) + END IF + IF (want_p) THEN + IF (i0 <= nb) THEN + CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, cth(i0), lane_p, & + lane_m(1, 3)) + ELSE + CALL add_pressure_lanes(nk, amplitude(i0:i1), omega_comp(i0:i1), k_comp(i0:i1), ch, bc, lane_p, lane_m(1, 3)) + END IF END IF IF (i0 == 1) first_ratio = ch(1) IF (use_tail) THEN @@ -2515,8 +2571,13 @@ SUBROUTINE component_sums(n, x_along, z, t, depth, stretch, scl, cb, sb, omega_c sum_vs = scl*CD_Wave_Sum_Lanes(lane_vs, lanes) sum_as = scl*CD_Wave_Sum_Lanes(lane_as, lanes) sum_ac = scl*CD_Wave_Sum_Lanes(lane_ac, lanes) - velocity = [sum_vc*cb, sum_vc*sb, sum_vs] - acceleration = [sum_as*cb, sum_as*sb, -sum_ac] + IF (spread) THEN + velocity = [sum_vc, scl*CD_Wave_Sum_Lanes(lane_vy, lanes), sum_vs] + acceleration = [sum_as, scl*CD_Wave_Sum_Lanes(lane_ay, lanes), -sum_ac] + ELSE + velocity = [sum_vc*cb, sum_vc*sb, sum_vs] + acceleration = [sum_as*cb, sum_as*sb, -sum_ac] + END IF IF (want_p) pdyn = scl*CD_Wave_Sum_Lanes(lane_p, lanes) END SUBROUTINE component_sums @@ -2557,6 +2618,34 @@ PURE SUBROUTINE add_velocity_lanes(n, amplitude, omega, ch, sh, c, s, lane_vc, l END DO END SUBROUTINE add_velocity_lanes + PURE SUBROUTINE add_spread_lanes(n, amplitude, omega, cosb, sinb, ch, sh, c, s, lane_vx, lane_vy, lane_vs, lane_ax, & + lane_ay, lane_ac, lane_mv, lane_ma) + !! add_velocity_lanes for components each with its own direction: the horizontal sums go to + !! x and y apart. + INTEGER, INTENT(IN) :: n + REAL(wp), INTENT(IN) :: amplitude(n), omega(n), cosb(n), sinb(n), ch(n), sh(n), c(n), s(n) + REAL(wp), INTENT(INOUT) :: lane_vx(n), lane_vy(n), lane_vs(n), lane_ax(n), lane_ay(n), lane_ac(n), lane_mv(n), & + lane_ma(n) + INTEGER :: j + REAL(wp) :: wv, wa, gv, ga, uh + DO j = 1, n + wv = amplitude(j)*omega(j) + wa = wv*omega(j) + gv = wv*ch(j) + ga = wa*ch(j) + uh = gv*c(j) + lane_vx(j) = lane_vx(j) + uh*cosb(j) + lane_vy(j) = lane_vy(j) + uh*sinb(j) + lane_vs(j) = lane_vs(j) + (wv*sh(j))*s(j) + uh = ga*s(j) + lane_ax(j) = lane_ax(j) + uh*cosb(j) + lane_ay(j) = lane_ay(j) + uh*sinb(j) + lane_ac(j) = lane_ac(j) + (wa*sh(j))*c(j) + lane_mv(j) = lane_mv(j) + gv + lane_ma(j) = lane_ma(j) + ga + END DO + END SUBROUTINE add_spread_lanes + PURE SUBROUTINE add_pressure_lanes(n, amplitude, omega, k, ch, c, lane_p, lane_mp) !! One block (n <= NLANE) of the pressure sum of CD_Component_Wave_Kinematics. INTEGER, INTENT(IN) :: n diff --git a/tests/test_vessel_motion.f90 b/tests/test_vessel_motion.f90 index 2dd92d10..1f2616f6 100644 --- a/tests/test_vessel_motion.f90 +++ b/tests/test_vessel_motion.f90 @@ -25,6 +25,9 @@ PROGRAM test_vessel_motion !! 7. Fail-closed option and file grammar. !! 8. Deck, vesselRAO in a spread sea: each spread component drives the vessel with its own !! heading and phase; a nonlinear stream wave is rejected by name. + !! 9. Deck, vesselRAO in a WaveKin 7 sea of rows with different directions referred to a point + !! (WaveRefX, WaveRefY): each row drives the vessel with its own heading, and its phase + !! from the position of the vessel relative to that point. USE CableDyn_Precision, ONLY: wp USE CableDyn_Vessel, ONLY: CD_VesselRAOType, CD_Vessel_Euler_DCM, CD_Vessel_Euler_Angular, & CD_Vessel_Point_Kinematics, CD_Vessel_RAO_Set, CD_Vessel_RAO_Eval, & @@ -49,6 +52,7 @@ PROGRAM test_vessel_motion CALL case_deck_rao_regular() CALL case_deck_rao_irregular() CALL case_deck_rao_spread() + CALL case_deck_rao_wavekin7() CALL case_deck_rao_stream_rejected() CALL case_deck_rejections() @@ -701,6 +705,67 @@ SUBROUTINE case_deck_rao_spread() CALL require(err < 1.0e-5_wp, 'rao spread: every spread component drives the vessel with its own phase') END SUBROUTINE case_deck_rao_spread + SUBROUTINE case_deck_rao_wavekin7() + !! A heave-only RAO over the full circle of headings: the fairlead below the reference point + !! heaves with the elevation of the file's rows there. + INTEGER, PARAMETER :: NC = 3 + REAL(wp), PARAMETER :: W(NC) = [0.6_wp, 0.75_wp, 0.9_wp], RE(NC) = [0.5_wp, -0.2_wp, 0.3_wp] + REAL(wp), PARAMETER :: IM(NC) = [0.1_wp, 0.4_wp, -0.25_wp], BETA(NC) = [0.65_wp, 2.9_wp, -2.2_wp] + REAL(wp), PARAMETER :: X0 = 31.5_wp, Y0 = -12.0_wp + CHARACTER(40) :: opts(9) + CHARACTER(512) :: msg + LOGICAL :: ok + REAL(wp), ALLOCATABLE :: d(:, :) + REAL(wp) :: kk(NC), t, eta, err, dist, r, rd, rdd + INTEGER :: nr, n, i, es, u + CHARACTER(200) :: em + OPEN (NEWUNIT=u, FILE='vm_wk7/vm_rao.txt', STATUS='REPLACE', ACTION='WRITE') + DO i = 1, 2 + WRITE (u, '(A,F7.1)') 'HEADING ', MERGE(-180.0_wp, 179.0_wp, i == 1) + WRITE (u, '(A)') '0.5 0 0 0 0 1 0 0 0 0 0 0 0' + WRITE (u, '(A)') '100 0 0 0 0 1 0 0 0 0 0 0 0' + END DO + CLOSE (u) + OPEN (NEWUNIT=u, FILE='vm_wk7/wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE') + WRITE (u, '(A)') '0.0 0.0 0.0 0.0' + DO i = 1, NC + WRITE (u, '(4ES24.16)') W(i), RE(i), IM(i), BETA(i) + END DO + CLOSE (u) + opts(1) = '0.05 dtM' + opts(2) = '20.0 TMax' + opts(3) = '60.0 WtrDpth' + opts(4) = '7 WaveKin' + opts(5) = 'vm_rao.txt vesselRAO' + opts(6) = '90.0|0.0|0.0 vesselRef' + opts(7) = '31.5 WaveRefX' + opts(8) = '-12.0 WaveRefY' + opts(9) = '6.0 rampTime' + CALL write_deck('vm_wk7/vm_wk7.dat', opts) + CALL run('vm_wk7/vm_wk7.dat', 'vm_wk7/out', ok, msg) + CALL require(ok, 'rao WaveKin 7: deck runs: '//TRIM(msg)) + CALL read_out('vm_wk7/out.out', d, nr, 10) + CALL require(nr == 401, 'rao WaveKin 7: full record written') + IF (nr /= 401) RETURN + DO i = 1, NC + CALL CD_Solve_Dispersion_Wavenumber(W(i), 60.0_wp, 9.80665_wp, kk(i), es, em) + CALL require(es == 0, 'rao WaveKin 7: dispersion solve: '//TRIM(em)) + END DO + err = 0.0_wp + DO n = 1, nr + t = d(1, n) + CALL CD_Vessel_Ramp(6.0_wp, t, r, rd, rdd) + eta = 0.0_wp + DO i = 1, NC + dist = (90.0_wp - X0)*COS(BETA(i)) + (0.0_wp - Y0)*SIN(BETA(i)) + eta = eta + SQRT(RE(i)**2 + IM(i)**2)*SIN(W(i)*t - kk(i)*dist + ATAN2(IM(i), RE(i))) + END DO + err = MAX(err, ABS(d(7, n) - (-14.0_wp + r*eta))) + END DO + WRITE (*, '(A,ES12.4)') ' rao WaveKin 7: max fairlead heave error [m] = ', err + CALL require(err < 1.0e-5_wp, 'rao WaveKin 7: each row drives the vessel with its own heading and phase') + END SUBROUTINE case_deck_rao_wavekin7 + SUBROUTINE case_deck_rao_stream_rejected() CHARACTER(40) :: opts(6) CHARACTER(512) :: msg diff --git a/tests/test_wave_component_sum.f90 b/tests/test_wave_component_sum.f90 index 1992e966..eddb3d3c 100644 --- a/tests/test_wave_component_sum.f90 +++ b/tests/test_wave_component_sum.f90 @@ -13,7 +13,11 @@ PROGRAM test_wave_component_sum !! sum of squares: they are independent) and, for the stretched depth, of the elevation; !! 2. the sum is a function of its arguments alone: the same call after others gives the same !! bits, and so do the tables in any thread count (the sums keep a fixed order); - !! 3. amplitude_scale equals a scaled table to round-off. + !! 3. amplitude_scale equals a scaled table to round-off; + !! 4. a short-crested table (a direction per component) equals the two long-crested tables it + !! is made of to round-off, for directions in every quadrant, one direction everywhere + !! equals the long-crested sum, and the stretched sum equals the unstretched one at the + !! stretched depth; half a table of cosines fails closed. USE, INTRINSIC :: ISO_FORTRAN_ENV, ONLY: INT64 USE CableDyn_Precision, ONLY: wp USE CableDyn_Hydro, ONLY: CD_Component_Wave_Kinematics, CD_Solve_Dispersion_Wavenumber, CD_HYDRO_OK @@ -45,6 +49,7 @@ PROGRAM test_wave_component_sum END DO WRITE (*, '(A,F6.3)') ' largest error over its round-off bound: ', worst CALL check_pure() + CALL check_short_crested() IF (nfail > 0) THEN WRITE (*, '(A,I0,A)') 'FAIL: ', nfail, ' assertion(s) failed' @@ -263,6 +268,85 @@ SUBROUTINE check_pure() 'amplitude_scale equals a scaled table to round-off') END SUBROUTINE check_pure + SUBROUTINE check_short_crested() + !! A table past the phase buffer in which every other component goes in direction a or b. + INTEGER, PARAMETER :: N = 4200, NH = N/2, NPAIR = 5 + REAL(wp), PARAMETER :: PAIRS(2, NPAIR) = RESHAPE([0.0_wp, 90.0_wp, 37.5_wp, 127.0_wp, 200.0_wp, 311.3_wp, & + -45.0_wp, 233.3_wp, 180.0_wp, 271.0_wp], [2, NPAIR]) + REAL(wp), PARAMETER :: X = 41.0_wp, Y = -17.5_wp, T = 1234.5_wp, D = 60.0_wp + REAL(wp) :: om(N), k(N), amp(N), ph(N), e2kd(N), tail(3, N), cb(N), sb(N), half(4, NH, 2) + REAL(wp) :: e0, u0(3), a0(3), p0, e1, u1(3), a1(3), p1, e2, u2(3), a2(3), p2, tol, ze + REAL(wp), PARAMETER :: ZS(3) = [-0.8_wp, -7.0_wp, -33.0_wp] + INTEGER :: ip, iz, es + CHARACTER(200) :: em + + CALL make_table(N, D, 4.0_wp, om, k, amp, ph, e2kd, tail) + tol = 32.0_wp*EPSILON(1.0_wp)*SUM(amp*MAX(1.0_wp, om*om, om*om/k)) + DO iz = 1, 2 + half(1, :, iz) = om(iz:N:2) + half(2, :, iz) = k(iz:N:2) + half(3, :, iz) = amp(iz:N:2) + half(4, :, iz) = ph(iz:N:2) + END DO + DO ip = 1, NPAIR + cb(1:N:2) = COS(PAIRS(1, ip)*PI/180.0_wp) + sb(1:N:2) = SIN(PAIRS(1, ip)*PI/180.0_wp) + cb(2:N:2) = COS(PAIRS(2, ip)*PI/180.0_wp) + sb(2:N:2) = SIN(PAIRS(2, ip)*PI/180.0_wp) + DO iz = 1, SIZE(ZS) + CALL CD_Component_Wave_Kinematics(X, Y, ZS(iz), T, D, 0.0_wp, .FALSE., om, k, amp, ph, e0, u0, a0, es, em, & + pdyn=p0, exp_m2kd=e2kd, tail_sums=tail, comp_cos=cb, comp_sin=sb) + CALL require(es == CD_HYDRO_OK, 'short-crested sum: '//TRIM(em)) + CALL long_crested(X, Y, ZS(iz), T, D, PAIRS(1, ip), half(:, :, 1), e1, u1, a1, p1, es, em) + CALL require(es == CD_HYDRO_OK, 'first long-crested half: '//TRIM(em)) + CALL long_crested(X, Y, ZS(iz), T, D, PAIRS(2, ip), half(:, :, 2), e2, u2, a2, p2, es, em) + CALL require(es == CD_HYDRO_OK, 'second long-crested half: '//TRIM(em)) + CALL require(ABS(e0 - e1 - e2) <= tol .AND. MAXVAL(ABS(u0 - u1 - u2)) <= tol .AND. & + MAXVAL(ABS(a0 - a1 - a2)) <= tol .AND. ABS(p0 - p1 - p2) <= tol, & + 'a short-crested table is the sum of its long-crested halves') + ! stretched: the sum at the stretched depth of the one total elevation + CALL CD_Component_Wave_Kinematics(X, Y, ZS(iz), T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, & + exp_m2kd=e2kd, tail_sums=tail, comp_cos=cb, comp_sin=sb) + CALL require(es == CD_HYDRO_OK, 'stretched short-crested sum: '//TRIM(em)) + IF (ZS(iz) > e1) CYCLE ! above this trough the point is dry + ze = (ZS(iz) - e1)*D/(D + e1) + CALL CD_Component_Wave_Kinematics(X, Y, ze, T, D, 0.0_wp, .FALSE., om, k, amp, ph, e2, u2, a2, es, em, & + exp_m2kd=e2kd, comp_cos=cb, comp_sin=sb) + CALL require(es == CD_HYDRO_OK, 'at the stretched depth: '//TRIM(em)) + CALL require(MAXVAL(ABS(u1 - u2)) <= tol .AND. MAXVAL(ABS(a1 - a2)) <= tol, & + 'stretched short-crested sum at the depth of the total elevation') + END DO + ! one direction everywhere + cb = COS(PAIRS(2, ip)*PI/180.0_wp) + sb = SIN(PAIRS(2, ip)*PI/180.0_wp) + CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, & + comp_cos=cb, comp_sin=sb) + CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, PAIRS(2, ip), .TRUE., om, k, amp, ph, e2, u2, a2, es, em) + CALL require(ABS(e1 - e2) <= tol .AND. MAXVAL(ABS(u1 - u2)) <= tol .AND. MAXVAL(ABS(a1 - a2)) <= tol, & + 'one direction per component equals the long-crested sum') + END DO + CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, & + comp_cos=cb) + CALL require(es /= CD_HYDRO_OK, 'comp_cos without comp_sin fails closed') + CALL CD_Component_Wave_Kinematics(X, Y, -4.0_wp, T, D, 0.0_wp, .TRUE., om, k, amp, ph, e1, u1, a1, es, em, & + comp_cos=cb(1:NH), comp_sin=sb(1:NH)) + CALL require(es /= CD_HYDRO_OK, 'direction tables of the wrong length fail closed') + END SUBROUTINE check_short_crested + + SUBROUTINE long_crested(x, y, z, t, depth, direction, table, eta, u, a, p, es, em) + !! The unstretched sum of the table rows omega, k, amplitude and phase along one direction. + REAL(wp), INTENT(IN) :: x, y, z, t, depth, direction, table(:, :) + REAL(wp), INTENT(OUT) :: eta, u(3), a(3), p + INTEGER, INTENT(OUT) :: es + CHARACTER(*), INTENT(OUT) :: em + REAL(wp) :: om(SIZE(table, 2)), k(SIZE(table, 2)), amp(SIZE(table, 2)), ph(SIZE(table, 2)) + om = table(1, :) + k = table(2, :) + amp = table(3, :) + ph = table(4, :) + CALL CD_Component_Wave_Kinematics(x, y, z, t, depth, direction, .FALSE., om, k, amp, ph, eta, u, a, es, em, pdyn=p) + END SUBROUTINE long_crested + SUBROUTINE sample(p, t, om, k, amp, ph, e2kd, tail, out) REAL(wp), INTENT(IN) :: p(3), t, om(:), k(:), amp(:), ph(:), e2kd(:), tail(:, :) REAL(wp), INTENT(OUT) :: out(7) diff --git a/tests/test_wavekin_modes.f90 b/tests/test_wavekin_modes.f90 index 896e9b30..2cfcf46b 100644 --- a/tests/test_wavekin_modes.f90 +++ b/tests/test_wavekin_modes.f90 @@ -17,7 +17,13 @@ !> - WvHiCOff: 1 rad/s keeps only the record's 8 s harmonic, exactly; none keeps both, as the !> 0.5 Hz default does (the record has nothing above it); !> - Currents 1: identical to 1e-12 m/s; -!> - the MoorDyn-C modes CableDyn does not support fail closed by name. +!> - WaveKin 7 rows of different directions, across the quadrants: the elevation is the sum of +!> the two single-direction seas and the velocity, acceleration likewise at the depth stretched +!> by the total elevation, to round-off; the beta of the omega = 0 row has no effect; +!> - WaveRefX, WaveRefY: an elevation record or a component file referred to a point equals the +!> same sea with the deck translated by it, for WaveKin 3 and 7 and WaterKin WaveKinMod 1; +!> - the MoorDyn-C modes CableDyn does not support, and a reference point with no record to +!> place, fail closed by name. PROGRAM test_wavekin_modes USE CableDyn_Precision, ONLY: wp USE CableDyn_DeckDriver, ONLY: CD_Deck_Fluid_Probe, CD_DECKDRV_OK @@ -32,6 +38,8 @@ PROGRAM test_wavekin_modes CALL check_cutoff() CALL check_mode('wavekin_mdc/currents1.dat', 'wavekin_mdc/mdc_currents1.csv', 'Currents 1', 1.0e-12_wp, & 1.0e-12_wp, 1.0e-12_wp, absolute=.TRUE.) + CALL check_short_crested() + CALL check_reference_point() CALL check_rejections() CALL check_still_water_ramp() IF (n_fail > 0) THEN @@ -234,6 +242,205 @@ SUBROUTINE expect_failure(wavekin, currents, needle) ' fails closed naming "'//needle//'": '//TRIM(em)) END SUBROUTINE expect_failure + SUBROUTINE expect_row_failure(wavekin, row, needle) + CHARACTER(*), INTENT(IN) :: wavekin, row, needle + REAL(wp) :: pos(3, 1), t(1), eta(1), vel(3, 1), acc(3, 1) + INTEGER :: es + CHARACTER(512) :: em + pos = 0.0_wp + pos(3, 1) = -5.0_wp + t = 0.0_wp + CALL write_fail_deck('wavekin_fail.dat', wavekin, '0', row) + CALL CD_Deck_Fluid_Probe('wavekin_fail.dat', pos, t, eta, vel, acc, es, em) + CALL require(es /= CD_DECKDRV_OK .AND. INDEX(em, needle) > 0, 'WaveKin '//wavekin//' with '//row// & + ' fails closed naming "'//needle//'": '//TRIM(em)) + END SUBROUTINE expect_row_failure + + SUBROUTINE write_components(beta_a, beta_b, which, beta_mean) + !! wave_frequencies.txt of two component sets, set a at direction beta_a and set b at beta_b + !! (radians), in file order a then b, so that the frequencies do not ascend; which: 0 both, + !! 1 a alone, 2 b alone. beta_mean: the direction of the omega = 0 row (default that of the + !! first component). + REAL(wp), INTENT(IN) :: beta_a, beta_b + INTEGER, INTENT(IN) :: which + REAL(wp), INTENT(IN), OPTIONAL :: beta_mean + INTEGER :: u + REAL(wp) :: b0 + b0 = MERGE(beta_b, beta_a, which == 2) + IF (PRESENT(beta_mean)) b0 = beta_mean + OPEN (NEWUNIT=u, FILE='wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE') + WRITE (u, '(A,ES24.16)') '0.0 0.0 0.0 ', b0 + IF (which /= 2) THEN + WRITE (u, '(A,ES24.16)') '0.6 0.7 0.2 ', beta_a + WRITE (u, '(A,ES24.16)') '1.1 0.25 -0.3 ', beta_a + END IF + IF (which /= 1) THEN + WRITE (u, '(A,ES24.16)') '0.6 -0.3 0.5 ', beta_b + WRITE (u, '(A,ES24.16)') '0.9 0.4 0.1 ', beta_b + END IF + CLOSE (u) + END SUBROUTINE write_components + + SUBROUTINE probe_components(beta_a, beta_b, which, pos, t, eta, vel, acc) + REAL(wp), INTENT(IN) :: beta_a, beta_b, pos(:, :), t(:) + INTEGER, INTENT(IN) :: which + REAL(wp), INTENT(OUT) :: eta(:), vel(:, :), acc(:, :) + INTEGER :: es + CHARACTER(512) :: em + CALL write_components(beta_a, beta_b, which) + CALL write_fail_deck('wavekin_short.dat', '7', '0') + CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, eta, vel, acc, es, em) + CALL require(es == CD_DECKDRV_OK, 'WaveKin 7 component deck evaluated: '//TRIM(em)) + END SUBROUTINE probe_components + + SUBROUTINE check_short_crested() + !! Rows of different directions against the two single-direction seas they are made of. The + !! deck stretches each sea by its own elevation, so a single sea is probed at the depth that + !! maps to the depth the total elevation stretches the point to. + INTEGER, PARAMETER :: NP = 6, NPAIR = 6 + REAL(wp), PARAMETER :: D = 50.0_wp + ! radians, within +-2 pi: the quadrants, the axes, and directions that are no multiple of 90 deg + REAL(wp), PARAMETER :: PAIRS(2, NPAIR) = RESHAPE([0.0_wp, 1.5707963267948966_wp, 0.3_wp, 2.4_wp, -1.0_wp, 3.7_wp, & + 5.1_wp, -2.2_wp, 3.141592653589793_wp, 4.71238898038469_wp, & + 0.65_wp, -6.0_wp], [2, NPAIR]) + REAL(wp) :: pos(3, NP), t(NP), eta(NP), vel(3, NP), acc(3, NP), ea(NP), va(3, NP), aa(3, NP) + REAL(wp) :: eb(NP), vb(3, NP), ab(3, NP), ze(NP), pa(3, NP), pb(3, NP), err_e, err_v, err_a, tmp(NP) + REAL(wp) :: tmpv(3, NP), tmpa(3, NP), e0(NP), v0(3, NP), a0(3, NP) + INTEGER :: ip, es + CHARACTER(512) :: em + pos(1, :) = [0.0_wp, 13.0_wp, -41.5_wp, 77.0_wp, -8.25_wp, 150.0_wp] + pos(2, :) = [0.0_wp, -7.0_wp, 22.0_wp, 90.5_wp, -33.0_wp, -120.0_wp] + pos(3, :) = [-2.0_wp, -5.0_wp, -11.0_wp, -20.0_wp, -33.0_wp, -46.0_wp] + t = [0.0_wp, 1.7_wp, 3.3_wp, 5.9_wp, 7.1_wp, 9.8_wp] + err_e = 0.0_wp + err_v = 0.0_wp + err_a = 0.0_wp + DO ip = 1, NPAIR + CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 0, pos, t, eta, vel, acc) + CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 1, pos, t, ea, tmpv, tmpa) + CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 2, pos, t, eb, tmpv, tmpa) + err_e = MAX(err_e, MAXVAL(ABS(eta - ea - eb))) + ze = (pos(3, :) - eta)*D/(D + eta) + pa = pos + pb = pos + pa(3, :) = ze*(D + ea)/D + ea + pb(3, :) = ze*(D + eb)/D + eb + CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 1, pa, t, tmp, va, aa) + CALL probe_components(PAIRS(1, ip), PAIRS(2, ip), 2, pb, t, tmp, vb, ab) + err_v = MAX(err_v, MAXVAL(ABS(vel - va - vb))) + err_a = MAX(err_a, MAXVAL(ABS(acc - aa - ab))) + END DO + WRITE (*, '(A,3ES10.3)') 'WaveKin 7 two directions vs two seas: elevation, velocity, acceleration ', err_e, & + err_v, err_a + CALL require(err_e <= 1.0e-13_wp, 'two directions: the elevation is the sum of the two seas') + CALL require(err_v <= 1.0e-13_wp, 'two directions: the velocity is the sum of the two seas') + CALL require(err_a <= 1.0e-13_wp, 'two directions: the acceleration is the sum of the two seas') + ! the mean row carries no direction: a different beta there leaves a one-direction table + CALL write_components(0.4_wp, 0.4_wp, 0) + CALL write_fail_deck('wavekin_short.dat', '7', '0') + CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, e0, v0, a0, es, em) + CALL require(es == CD_DECKDRV_OK, 'one direction evaluated: '//TRIM(em)) + CALL write_components(0.4_wp, 0.4_wp, 0, beta_mean=1.3_wp) + CALL CD_Deck_Fluid_Probe('wavekin_short.dat', pos, t, eta, vel, acc, es, em) + CALL require(es == CD_DECKDRV_OK, 'mean-row direction evaluated: '//TRIM(em)) + CALL require(ALL(eta == e0) .AND. ALL(vel == v0) .AND. ALL(acc == a0), & + 'the direction of the omega = 0 row has no effect') + END SUBROUTINE check_short_crested + + SUBROUTINE copy_deck_with(source, target, match, row, replace) + !! source with row after the line that holds match, or in its place (replace). + CHARACTER(*), INTENT(IN) :: source, target, match, row + LOGICAL, INTENT(IN) :: replace + CHARACTER(1024) :: line + INTEGER :: u_in, u_out, ios + OPEN (NEWUNIT=u_in, FILE=source, STATUS='OLD', ACTION='READ') + OPEN (NEWUNIT=u_out, FILE=target, STATUS='REPLACE', ACTION='WRITE') + DO + READ (u_in, '(A)', IOSTAT=ios) line + IF (ios /= 0) EXIT + IF (INDEX(line, match) > 0) THEN + IF (.NOT. replace) WRITE (u_out, '(A)') TRIM(line) + WRITE (u_out, '(A)') row + ELSE + WRITE (u_out, '(A)') TRIM(line) + END IF + END DO + CLOSE (u_in) + CLOSE (u_out) + END SUBROUTINE copy_deck_with + + SUBROUTINE compare_translated(plain, shifted, shift, label) + !! The deck shifted, referred to the point shift, at positions p + shift, against the deck + !! plain at p. + CHARACTER(*), INTENT(IN) :: plain, shifted, label + REAL(wp), INTENT(IN) :: shift(2) + INTEGER, PARAMETER :: NP = 5 + REAL(wp) :: pos(3, NP), moved(3, NP), t(NP), e0(NP), v0(3, NP), a0(3, NP), e1(NP), v1(3, NP), a1(3, NP), err + INTEGER :: es + CHARACTER(512) :: em + pos(1, :) = [0.0_wp, 11.0_wp, -37.0_wp, 60.5_wp, -5.0_wp] + pos(2, :) = [0.0_wp, -3.5_wp, 19.0_wp, 44.0_wp, -71.0_wp] + pos(3, :) = [-1.5_wp, -6.0_wp, -14.0_wp, -27.0_wp, -41.0_wp] + t = [0.3_wp, 1.9_wp, 4.4_wp, 6.2_wp, 9.1_wp] + moved = pos + moved(1, :) = pos(1, :) + shift(1) + moved(2, :) = pos(2, :) + shift(2) + CALL CD_Deck_Fluid_Probe(plain, pos, t, e0, v0, a0, es, em) + CALL require(es == CD_DECKDRV_OK, label//' plain deck evaluated: '//TRIM(em)) + IF (es /= CD_DECKDRV_OK) RETURN + CALL CD_Deck_Fluid_Probe(shifted, moved, t, e1, v1, a1, es, em) + CALL require(es == CD_DECKDRV_OK, label//' referred deck evaluated: '//TRIM(em)) + IF (es /= CD_DECKDRV_OK) RETURN + err = MAX(MAXVAL(ABS(e1 - e0)), MAXVAL(ABS(v1 - v0)), MAXVAL(ABS(a1 - a0))) + WRITE (*, '(A,A,ES10.3)') label, ' referred to a point vs the deck translated: ', err + CALL require(err <= 1.0e-12_wp, label//': a shifted reference point equals the deck translated') + ! and the point matters: the same positions without the shift are another sea + CALL CD_Deck_Fluid_Probe(shifted, pos, t, e1, v1, a1, es, em) + CALL require(MAXVAL(ABS(e1 - e0)) > 1.0e-3_wp, label//': the reference point moves the sea') + END SUBROUTINE compare_translated + + SUBROUTINE check_reference_point() + REAL(wp), PARAMETER :: SHIFT(2) = [250.37_wp, -70.19_wp] + CHARACTER(*), PARAMETER :: ROWS = '250.37 WaveRefX'//NEW_LINE('a')//'-70.19 WaveRefY' + INTEGER :: u + ! WaveKin 7, one direction and two + CALL write_components(0.9_wp, 0.9_wp, 0) + CALL write_fail_deck('wavekin_ref0.dat', '7', '0') + CALL write_fail_deck('wavekin_ref1.dat', '7', '0', ROWS) + CALL compare_translated('wavekin_ref0.dat', 'wavekin_ref1.dat', SHIFT, 'WaveKin 7') + CALL write_components(0.65_wp, 2.9_wp, 0) + CALL compare_translated('wavekin_ref0.dat', 'wavekin_ref1.dat', SHIFT, 'WaveKin 7, two directions') + ! WaveKin 3 + CALL copy_deck_with('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/wavekin3_ref.dat', 'dtWave', ROWS, .FALSE.) + CALL compare_translated('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/wavekin3_ref.dat', SHIFT, 'WaveKin 3') + ! WaterKin WaveKinMod 1, along a direction off the axes + OPEN (NEWUNIT=u, FILE='wavekin_mdc/waterkin_ref.dat', STATUS='REPLACE', ACTION='WRITE') + WRITE (u, '(A)') 'MoorDyn v2 water kinematics file' + WRITE (u, '(A)') '(header line 2)' + WRITE (u, '(A)') '--------------------------- WAVES -------------------------------------' + WRITE (u, '(A)') '1 WaveKinMod' + WRITE (u, '(A)') '"wave_elevation.txt" WaveKinFile' + WRITE (u, '(A)') '0.1 dtWave' + WRITE (u, '(A)') '37.0 WaveDir' + WRITE (u, '(A)') '0 X wave input' + WRITE (u, '(A)') '0.0' + WRITE (u, '(A)') '0 Y wave input' + WRITE (u, '(A)') '0.0' + WRITE (u, '(A)') '0 Z wave input' + WRITE (u, '(A)') '0.0' + WRITE (u, '(A)') '--------------------------- CURRENT -----------------------------------' + WRITE (u, '(A)') '0 CurrentMod' + WRITE (u, '(A)') '0 Z current' + WRITE (u, '(A)') '0.0' + WRITE (u, '(A)') '--------------------- need this line ----------------------------------' + CLOSE (u) + CALL copy_deck_with('wavekin_mdc/wavekin3.dat', 'wavekin_mdc/waterkin0.dat', 'WaveKin ', & + 'waterkin_ref.dat WaterKin', .TRUE.) + CALL copy_deck_with('wavekin_mdc/waterkin0.dat', 'wavekin_mdc/waterkin1.dat', 'WaterKin', & + 'waterkin_ref.dat WaterKin'//NEW_LINE('a')//ROWS, .TRUE.) + CALL compare_translated('wavekin_mdc/waterkin0.dat', 'wavekin_mdc/waterkin1.dat', SHIFT, 'WaterKin WaveKinMod 1') + END SUBROUTINE check_reference_point + SUBROUTINE check_still_water_ramp() !! rampTime is valid without waves (it also ramps prescribed motion). The still-water !! kinematics must then be exactly zero, including at t = 0 where the ramp is closed @@ -252,19 +459,17 @@ SUBROUTINE check_still_water_ramp() END SUBROUTINE check_still_water_ramp SUBROUTINE check_rejections() - INTEGER :: u CALL expect_failure('1', '0', 'MoorDyn-C API') CALL expect_failure('2', '0', 'FFT_GRID') CALL expect_failure('4', '0', 'no kinematics source in MoorDyn-C') CALL expect_failure('6', '0', 'no kinematics source in MoorDyn-C') CALL expect_failure('0', '3', 'no current source in MoorDyn-C') CALL expect_failure('0', '4', 'no current source in MoorDyn-C') - ! short-crested components (rows with different directions) - OPEN (NEWUNIT=u, FILE='wave_frequencies.txt', STATUS='REPLACE', ACTION='WRITE') - WRITE (u, '(A)') '0.0 0.0 0.0 0.0' - WRITE (u, '(A)') '0.6 0.7 0.2 0.0' - WRITE (u, '(A)') '0.9 -0.3 0.5 1.57' - CLOSE (u) - CALL expect_failure('7', '0', 'different directions') + ! a reference point needs a record or component file to place, and a bounded position + CALL write_components(0.4_wp, 0.4_wp, 0) + CALL expect_row_failure('0', '1.0 WaveRefX', 'reference point') + CALL expect_row_failure('0', '1.0 WaveRefY'//NEW_LINE('a')//'airy 2.0 8.0 0.0 waves', 'reference point') + CALL expect_row_failure('7', '2.0e6 WaveRefX', 'at most 1e6 m') + CALL expect_row_failure('7', 'nan WaveRefY', 'WaveRefY') END SUBROUTINE check_rejections END PROGRAM test_wavekin_modes From 4e91861f9ccf7f500e0a7c08cd697d93d836ada3 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sat, 7 Nov 2026 22:07:29 +0900 Subject: [PATCH 2/5] fix(python): judge wave records as the native reader does DeckFile accepted WaveRefX and WaveRefY with a WaterKin file of WaveKinMod 0, which the solver refuses, and did not read the rows of wave_frequencies.txt. It now reads the WaveKinMod of a WaterKin file and the rows of wave_frequencies.txt from the deck's folder, with the rules of the native reader: directions within 2 pi, the omega = 0 first row, positive frequencies, component bounds and the effective depth. A vesselRAO deck is judged by the same WaveKinMod. A file that is not there passes this check, and the solver, which must open it, refuses it by name while it is missing. The bound on WaveRefX and WaveRefY applies to the value in force, as in the native reader: a later row replaces an earlier one. A new parity test runs 41 decks through both DeckFile and the driver and requires the same verdict. The deck format notes record how a non-zero omega = 0 row differs from MoorDyn-C. Signed-off-by: Jae Hoon Seo --- doc/driver_format.md | 5 + python/cabledyn/deck_file.py | 136 +++++++++++++++++++-- python/tests/test_deck_native_parity.py | 156 ++++++++++++++++++++++++ 3 files changed, 286 insertions(+), 11 deletions(-) diff --git a/doc/driver_format.md b/doc/driver_format.md index 5dcc1a2e..b49f153b 100644 --- a/doc/driver_format.md +++ b/doc/driver_format.md @@ -1382,6 +1382,11 @@ time-varying and of a sheared current to MoorDyn-C's within 0.9 % on a 70 m rope ([options](options.rst)) are an extension of CableDyn and no part of MoorDyn-C or MoorDyn-F: a deck with them is not a plain MoorDyn deck, and a MoorDyn solver would not place the phases there. +A non-zero `omega = 0` row raises the still-water level. CableDyn then takes every wave number at +the depth `WtrDpth` plus that level, as it does for the mean of an elevation record, where +MoorDyn-C keeps `WtrDpth` and adds the level to the elevation alone. The two differ by the order +of the level over the depth in the long components. + (waterkin-file-modes)= ### MoorDyn-F `WaterKin` file modes diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index 107a242e..0ba6f8c5 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -663,6 +663,16 @@ def _is_subnormal(value: float) -> bool: return math.isfinite(value) and value != 0.0 and abs(value) < sys.float_info.min +def _plain_real(token: str) -> float | None: + """Mirror the native ``parse_real_token``: the finite, normal real a token reads as.""" + if any(character in " \"'/,;*\t" for character in token): + return None + value = _fortran_real_value(token) + if value is None or not math.isfinite(value) or _is_subnormal(value): + return None + return value + + def _real_token_defect(token: str) -> str | None: """Mirror the native ``real_token_defect``: why a numeric token is inadmissible. @@ -2428,13 +2438,17 @@ def validate(self) -> None: between rows. A deck that passes can still be rejected by the solver on physical grounds. - The ancillary files a deck names (motion, bathymetry, WaterKin, Syrope, - and the MoorDyn-C kinematics files) are not read here, so the checks - that depend on their contents stay with the solver: for example the - anchor height on a bathymetry surface, whether a WaterKin file carries - waves (``WaveKinMod``) for ``vesselRAO``, ``WaveRefX`` or ``WaveRefY``, - and whether its current doubles a ``current`` row or ``Currents 1``, - ``2`` or ``5``. + A WaterKin file's ``WaveKinMod`` (which decides whether ``vesselRAO``, + ``WaveRefX`` and ``WaveRefY`` have waves to act on) and the rows of + ``wave_frequencies.txt`` are read here, from the deck's folder. A file + that is not there, as for a deck held in memory under its nominal path, + is not an error at this stage: the deck passes these checks and the + solver, which must open the file, refuses it by name if it is still + missing. The other ancillary files (motion, bathymetry, Syrope and the + other MoorDyn-C kinematics files) are not read here, so the checks that + depend on their contents stay with the solver: for example the anchor + height on a bathymetry surface, and whether a WaterKin current doubles + a ``current`` row or ``Currents 1``, ``2`` or ``5``. Raises ------ @@ -2588,6 +2602,8 @@ def effective(key: str, default: float) -> float: has_waterkin_file = waterkin_value.lower() not in {"", "none", "seastate"} and any( character in _WATERKIN_FILENAME_LETTERS for character in waterkin_value ) + if has_waterkin_file and motion_source == "vesselrao": + has_waterkin_file = self._waterkin_wave_mode(waterkin_value) in {None, 1} has_trains = any(o.keyword.lower() in _WAVETRAIN_KEYWORDS for o in options) # WaveKin 3/7 installs linear component waves (native has_wave) if motion_source == "vesselrao" and not ( @@ -2680,11 +2696,19 @@ def effective(key: str, default: float) -> float: raise DeckFormatError( f"{self._label}: WaveKin {mdc_wave_kin} waves require dtM, TMax and WtrDpth" ) - # native validation: WaveRefX and WaveRefY place the phases of component waves alone. A - # WaterKin file may carry them (its WaveKinMod is read at run time). - wave_ref_rows = [o for o in options if o.keyword.lower() in {"waverefx", "waverefy"}] + if mdc_wave_kin == 7: + self._check_wave_frequencies(_native_float(option_by_key["water_depth"].values[0])) + # native validation: WaveRefX and WaveRefY place the phases of component waves alone, + # which a WaterKin file carries under WaveKinMod 1. + # The bound applies to the value in force: a later row of an option replaces an earlier one. + wave_ref_rows = { + o.keyword.lower(): o for o in options if o.keyword.lower() in {"waverefx", "waverefy"} + } if wave_ref_rows: - if any(abs(_native_float(o.values[0])) > _NATIVE_MAX_COORDINATE for o in wave_ref_rows): + if any( + abs(_native_float(o.values[0])) > _NATIVE_MAX_COORDINATE + for o in wave_ref_rows.values() + ): raise DeckFormatError( f"{self._label}: WaveRefX and WaveRefY must be at most 1e6 m in magnitude" ) @@ -2693,6 +2717,7 @@ def effective(key: str, default: float) -> float: has_record = bool(mdc_wave_kin) or ( waterkin_text.lower() not in {"", "none", "seastate"} and any(character in _WATERKIN_FILENAME_LETTERS for character in waterkin_text) + and self._waterkin_wave_mode(waterkin_text) in {None, 1} ) if self.caller_driven or has_waves or trains or not has_record: raise DeckFormatError( @@ -4337,6 +4362,95 @@ def is_body_kind(kind: str, pinned: bool) -> bool: continue raise DeckFormatError(f"{self._label}: unsupported OUTPUT channel {channel!r}") + def _ancillary_records(self, value: str) -> list[str] | None: + """Return the records of the ancillary file ``value`` names, ``None`` if it is not there. + + As the native ``path_join``: a relative name is taken from the deck's folder. + """ + name = _strip_quotes(value) + absolute = name.startswith(("/", "\\")) or name[1:2] == ":" + source = Path(name) if absolute else self.path.parent / name + try: + payload = source.read_bytes() + except OSError: + return None + records, _ = _split_records(payload.decode("utf-8", "surrogateescape")) + if records and records[0].startswith(_UTF8_BOM): + records[0] = records[0][len(_UTF8_BOM) :] + return records + + def _waterkin_wave_mode(self, value: str) -> int | None: + """Return the ``WaveKinMod`` of WaterKin file ``value``, ``None`` if it is not there. + + Mirrors the head of the native ``read_waterkin_records``: two header lines, a waves + header, then the WaveKinMod line. + """ + records = self._ancillary_records(value) + if records is None: + return None + if len(records) < 3: + raise DeckFormatError(f"{self._label}: WaterKin file ended inside its header") + if len(records) < 4: + raise DeckFormatError(f"{self._label}: WaterKin file has no WaveKinMod line") + tokens = records[3].split() + if not tokens: + raise DeckFormatError(f"{self._label}: WaterKin WaveKinMod line has no value") + if _strip_quotes(tokens[0]).lower() == "seastate": + return 2 + if re.fullmatch(r"[+-]?[0-9]{1,18}", tokens[0]) is None: + raise DeckFormatError( + f"{self._label}: WaterKin WaveKinMod must be 0, 1, 2, or SEASTATE" + ) + mode = int(tokens[0]) + if not 0 <= mode <= 2: + raise DeckFormatError( + f"{self._label}: WaterKin WaveKinMod supports 0 (none), 1 (elevation history), " + "or 2 (coupled SeaState hybrid)" + ) + return mode + + def _check_wave_frequencies(self, depth: float) -> None: + """Check ``wave_frequencies.txt`` in the deck's folder, if it is there. + + Mirrors the ``WaveKin 7`` branch of the native ``read_moordyn_c_kinematics``, the + component bounds of its validation and the effective depth of the wave precompute. + A row takes four numbers, or three and a direction of zero. + """ + records = self._ancillary_records(_MOORDYN_C_WAVE_FILES[7]) + if records is None: + return + label = f"{self._label}: WaveKin 7: wave_frequencies.txt" + rows: list[tuple[float, float, float, float]] = [] + for number, record in enumerate(records, start=1): + tokens = record.split() + if not tokens: + continue + values = [_plain_real(token) for token in tokens[:4]] + if len(values) == 4 and values[3] is None: + values = values[:3] + if len(values) < 3 or any(value is None for value in values): + raise DeckFormatError(f"{label} line {number}: row is not 3 to 4 plain numbers") + omega, real_part, imaginary, *beta = values + rows.append((omega, real_part, imaginary, beta[0] if beta else 0.0)) + if len(rows) < 2: + raise DeckFormatError(f"{label} needs at least two `omega Re Im` rows") + if rows[0][0] != 0.0: + raise DeckFormatError(f"{label}: the first row must be at omega = 0") + for index, (omega, real_part, imaginary, beta) in enumerate(rows): + if abs(beta) > 2.0 * math.pi: + raise DeckFormatError(f"{label} directions are radians within +-2 pi") + if index == 0: + continue + if omega <= 0.0: + raise DeckFormatError(f"{label} frequencies after the first row must be positive") + if math.hypot(real_part, imaginary) > 1.0e3 or omega > 1.0e4: + raise DeckFormatError( + f"{label} components need amplitudes of at most 1e3 m and frequencies of " + "at most 1e4 rad/s" + ) + if depth + rows[0][2] <= 0.0: + raise DeckFormatError(f"{label}: the mean level leaves a non-positive effective depth") + def _check_torsion_normal( self, where: str, diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index 66c529ff..9059e3da 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -86,6 +86,162 @@ def _rejects(text: str, message: str, *, caller_driven: bool = False) -> None: # --------------------------------------------------------------------------- options +_WATERKIN = ( + "MoorDyn v2 water kinematics file\n(header line 2)\n--- WAVES ---\n" + '{mode} WaveKinMod\n"eta.txt" WaveKinFile\n0.25 dtWave\n37.0 WaveDir\n' + "0 X\n0.0\n0 Y\n0.0\n0 Z\n0.0\n--- CURRENT ---\n0 CurrentMod\n0 Z\n0.0\n" + "--- need this line ---\n" +) +_ETA = "time elevation\n" + "".join( + f"{0.5 * i:g} {0.4 * (-1.0) ** (i // 4):g}\n" for i in range(41) +) +_RAO = (EXAMPLES / "data" / "vessel" / "sample_rao.txt").read_text(encoding="utf-8") +_ROWS = "0 0 0 0\n0.5 0.5 0 0.3\n0.6 0.3 0.1 2.4\n0.6 0.2 -0.1 -1.0\n" + + +def _waterkin(mode: str) -> dict[str, str]: + return {"wk.dat": _WATERKIN.format(mode=mode), "eta.txt": _ETA} + + +# (rows of the deck, files next to it): every case is judged by DeckFile and by the driver +_WAVE_RECORD_CASES = [ + # a reference point needs component waves + (("12.5 WaveRefX",), {}), + (("0 WaveRefY",), {}), + (("airy 2 8 0 waves", "1 WaveRefX"), {}), + (("airy 2 8 0 wavetrain", "1 WaveRefY"), {}), + (("7 WaterKin", "12.5 WaveRefX", "-3 WaveRefY"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "1.0e6 WaveRefX"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "1.0000001e6 WaveRefX"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "-2e6 WaveRefY"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "nan WaveRefY"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "north WaveRefY"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "1e-320 WaveRefX"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "2e6 WaveRefX", "1 WaveRefX"), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin", "1 WAVEREFX"), {"wave_frequencies.txt": _ROWS}), + # the WaveKinMod of a WaterKin file decides + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("1")), + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("0")), + (("wk.dat WaterKin", "0 WaveRefY"), _waterkin("0")), + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("2")), + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("SEASTATE")), + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("3")), + (("wk.dat WaterKin", "12.5 WaveRefX"), _waterkin("one")), + (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\n"}), + (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\nwaves\n"}), + (("wk.dat WaterKin", "12.5 WaveRefX"), {"wk.dat": "header\nheader\nwaves\n\n"}), + (("wk.dat WaterKin", "rao.txt vesselRAO"), {**_waterkin("0"), "rao.txt": _RAO}), + (("wk.dat WaterKin", "rao.txt vesselRAO"), {**_waterkin("1"), "rao.txt": _RAO}), + # the rows of wave_frequencies.txt + (("7 WaterKin",), {"wave_frequencies.txt": _ROWS}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 1.0\n0.5 0.5 0 0.3\n\n0.6 0.3 0.1 0.3\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0\n0.5 0.5 0\n0.6 0.3 0.1 -6.2\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 6.3\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 -6.3\n0.5 0.5 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 east\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 nan\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 nan 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "omega Re Im\n0 0 0 0\n0.5 0.5 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0.1 0 0 0\n0.5 0.5 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 0\n0 0.5 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 2000 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n2e4 0.1 0 0\n"}), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 -60 0\n0.5 0.5 0 0\n"}), +] + + +def _stage(folder: Path, rows: tuple[str, ...], files: dict[str, str]) -> Path: + folder.mkdir(parents=True, exist_ok=True) + for name, text in files.items(): + (folder / name).write_text(text, encoding="utf-8") + deck = folder / "deck.dat" + deck.write_text(_options(_HELD, *rows), encoding="utf-8") + return deck + + +def _accepted(deck: Path) -> bool: + try: + DeckFile.read(deck) + except DeckFormatError: + return False + return True + + +@pytest.mark.skipif( + not os.environ.get("CABLEDYN_DRIVER", "").strip(), + reason="needs the native driver (CABLEDYN_DRIVER)", +) +@pytest.mark.parametrize(("rows", "files"), _WAVE_RECORD_CASES) +def test_wave_record_rules_match_the_driver(tmp_path, rows, files): + deck = _stage(tmp_path, rows, files) + run = subprocess.run( + [os.environ["CABLEDYN_DRIVER"], str(deck), str(tmp_path / "out")], + capture_output=True, + text=True, + errors="replace", + env={**os.environ, "OMP_NUM_THREADS": "1", "OPENBLAS_NUM_THREADS": "1"}, + check=False, + ) + assert _accepted(deck) == (run.returncode == 0), run.stdout[-400:] + run.stderr[-400:] + + +@pytest.mark.parametrize( + ("rows", "files", "message"), + [ + ( + ("wk.dat WaterKin", "1 WaveRefX"), + _waterkin("0"), + "reference point of an elevation record", + ), + ( + ("wk.dat WaterKin", "1 WaveRefX"), + _waterkin("SEASTATE"), + "reference point of an elevation", + ), + (("wk.dat WaterKin", "1 WaveRefX"), _waterkin("3"), "WaveKinMod supports 0"), + (("wk.dat WaterKin", "1 WaveRefX"), _waterkin("one"), "must be 0, 1, 2, or SEASTATE"), + (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\n"}, "ended inside its header"), + (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\nc\n"}, "no WaveKinMod line"), + (("wk.dat WaterKin", "1 WaveRefX"), {"wk.dat": "a\nb\nc\n \n"}, "line has no value"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 0.5 0 6.3\n"}, "within \\+-2 pi"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n"}, "at least two"), + (("7 WaterKin",), {"wave_frequencies.txt": "1 0 0 0\n2 1 0 0\n"}, "at omega = 0"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0 1 0 0\n"}, "must be positive"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 x 0 0\n"}, "plain numbers"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 0 0\n0.5 2000 0 0\n"}, "at most 1e3 m"), + (("7 WaterKin",), {"wave_frequencies.txt": "0 0 -60 0\n0.5 1 0 0\n"}, "effective depth"), + ], +) +def test_wave_files_next_to_the_deck_are_read(tmp_path, rows, files, message): + with pytest.raises(DeckFormatError, match=message): + DeckFile.read(_stage(tmp_path, rows, files)) + + +def test_wave_files_that_are_not_there_are_left_to_the_solver(tmp_path): + # In memory under the nominal path, or on disk ahead of its files: the deck passes here, + # and the driver, which must open the file, refuses it by name while it is missing. + for rows in (("wk.dat WaterKin", "1 WaveRefX"), ("7 WaterKin", "1 WaveRefX")): + DeckFile.from_text(_options(_HELD, *rows)) + DeckFile.from_text(_options(_HELD, *rows), path=tmp_path / "deck.dat") + # once the file is there the same text is judged by it + (tmp_path / "wk.dat").write_text(_WATERKIN.format(mode="0"), encoding="utf-8") + with pytest.raises(DeckFormatError, match="reference point of an elevation record"): + DeckFile.from_text( + _options(_HELD, "wk.dat WaterKin", "1 WaveRefX"), path=tmp_path / "deck.dat" + ) + + +def test_a_waterkin_file_is_found_as_the_native_reader_finds_it(tmp_path): + # quoted, in a subfolder, and by an absolute path + (tmp_path / "sub").mkdir() + (tmp_path / "sub" / "wk.dat").write_text(_WATERKIN.format(mode="0"), encoding="utf-8") + for value in ('"sub/wk.dat"', "sub/wk.dat", (tmp_path / "sub" / "wk.dat").as_posix()): + with pytest.raises(DeckFormatError, match="reference point of an elevation record"): + DeckFile.read(_stage(tmp_path, (f"{value} WaterKin", "1 WaveRefX"), {})) + + @pytest.mark.parametrize("value", ["0", "-0.1"]) def test_dtwave_must_be_positive(value): _rejects(_options(_HELD, f"{value} dtWave"), "OPTION dtWave must be positive") From 7674f46e8e16b9076951cf7626ae4893fa02aed7 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sat, 7 Nov 2026 22:33:58 +0900 Subject: [PATCH 3/5] fix(python): read the ancillary files a deck names DeckFile left the contents of the files a deck names to the solver, so it certified decks the solver refuses. It now reads them from the deck's folder with the rules of their native readers: the WaterKin file and its wave record and current table, including a current that doubles a current row or Currents 1, 2 or 5; the Currents 1, 2 and 5 tables; motionFile rows on the route the deck takes, with TURBINES records and the roll column; vesselMotion rows; the vesselRAO table and the headings of the deck's sea; the bathymetry grid and the Fixed points on it; and the Syrope settings, OWC table and working-curve admissibility. A file that is not there passes this check, and the solver, which must open it, refuses it by name while it is missing. The parity test now runs 160 decks through DeckFile and the driver and requires the same verdict. Against the earlier validator, 28 of 1943 decks from the examples and the test suite are newly refused, and the driver refuses all 28. Signed-off-by: Jae Hoon Seo --- python/cabledyn/deck_file.py | 1257 ++++++++++++++++++++++- python/tests/test_deck_file.py | 7 +- python/tests/test_deck_native_parity.py | 323 +++++- 3 files changed, 1539 insertions(+), 48 deletions(-) diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index 0ba6f8c5..06cd88fd 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -12,6 +12,7 @@ from __future__ import annotations +import bisect import contextlib import copy import hashlib @@ -22,10 +23,11 @@ import re import sys import tempfile -from collections.abc import Callable, Mapping +from collections.abc import Callable, Iterable, Iterator, Mapping, Sequence from dataclasses import dataclass, field from pathlib import Path, PureWindowsPath from types import MappingProxyType +from typing import Any from cabledyn._paths import windows_device_component from cabledyn.errors import DeckFormatError @@ -663,16 +665,6 @@ def _is_subnormal(value: float) -> bool: return math.isfinite(value) and value != 0.0 and abs(value) < sys.float_info.min -def _plain_real(token: str) -> float | None: - """Mirror the native ``parse_real_token``: the finite, normal real a token reads as.""" - if any(character in " \"'/,;*\t" for character in token): - return None - value = _fortran_real_value(token) - if value is None or not math.isfinite(value) or _is_subnormal(value): - return None - return value - - def _real_token_defect(token: str) -> str | None: """Mirror the native ``real_token_defect``: why a numeric token is inadmissible. @@ -1455,6 +1447,872 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None: return None +# --------------------------------------------------------------------------- ancillary files +# The native readers of the files a deck names, as far as they refuse a deck. Each function +# follows the subroutine it names in src/CableDyn_DeckDriver.f90 (or the module given) and raises +# _AncillaryError with the native message; DeckFile adds the deck's label. + +_NATIVE_RECORD_LENGTH = 512 # DECK_RECLEN: the longest non-comment record of next_record +_NATIVE_TABLE_ROW_LENGTH = 65536 # LONG_ROW of read_moordyn_c_kinematics +_NATIVE_EPSILON = sys.float_info.epsilon +_NUMERIC_SHAPE = re.compile( + r"[+-]?(?:(?i:inf|infinity|nan)|(?:[0-9]+\.?[0-9]*|\.[0-9]+)(?:[eEdD][+-]?[0-9]+)?)" +) +_MAX_GRID_STEPS = float(2**31 - 2) +_MAX_MOTION_SPEED = 1.0e4 +_MAX_MOTION_ACCEL = 1.0e6 +_MAX_FLUID_SPEED = 1.0e3 +_MAX_WAVE_SAMPLES = 1.0e6 +_MAX_CURRENT_SAMPLES = 2000000 +_SYROPE_CURVE_SAMPLES = 30 # CD_SYROPE_NWC +_SYROPE_SCAN_SAMPLES = 200 # CD_SYROPE_NSCAN +_BATHY_CREASE_FRACTION = 0.05 +_BATHY_CREASE_MIN_HEIGHT = 1.0e-6 +_SEA_DEFAULT_DIRECTIONS = 9 + + +class _AncillaryError(ValueError): + """A refusal of a native ancillary-file reader.""" + + +def _token_real(token: str) -> float | None: + """Mirror the native ``parse_real_token``: the finite, normal real of a whole token.""" + if _NUMERIC_SHAPE.fullmatch(token) is None: + return None + try: + value = float(re.sub("[dD]", "e", token)) + except ValueError: + return None + if not math.isfinite(value) or _is_subnormal(value): + return None + return value + + +def _real_row(tokens: list[str], count: int) -> tuple[int, list[float]]: + """Mirror the native ``parse_real_row``: status 0, 1 (no number) or 2 (malformed), values.""" + if not tokens or _NUMERIC_SHAPE.fullmatch(tokens[0]) is None: + return 1, [] + values: list[float] = [] + for token in tokens[:count]: + value = _token_real(token) + if value is None: + return 2, [] + values.append(value) + return (0, values) if len(values) == count else (2, []) + + +def _data_rows(records: list[str], what: str) -> Iterator[tuple[int, list[str]]]: + """Mirror the native ``next_record``: non-blank rows without comments, by file line.""" + for number, raw in enumerate(records, start=1): + visible = _strip_inline_comment(raw) + if not visible.strip(_WHITESPACE): + continue + defect = _character_error(visible) + if defect is not None: + raise _AncillaryError(f"{what} line {number} {defect}") + if len(visible) > _NATIVE_RECORD_LENGTH: + raise _AncillaryError(f"{what} line {number} is longer than the record limit") + yield number, visible.split() + + +def _increasing(values: Sequence[float]) -> bool: + return all(later > earlier for earlier, later in itertools.pairwise(values)) + + +def _padded_fft_size(wanted: int) -> int: + """Mirror ``padded_fft_size``: twice the next 7-smooth half length.""" + half = max((wanted + 1) // 2, 1) + while True: + rest = half + for prime in (2, 3, 5, 7): + while rest % prime == 0: + rest //= prime + if rest == 1: + return 2 * half + half += 1 + + +def _check_elevation_record( + records: list[str], dt_wave: float, window: float, cutoff: float, depth: float | None +) -> None: + """Mirror ``read_wave_elevation_components`` and the component bounds of validation.""" + import numpy as np + + times: list[float] = [] + heights: list[float] = [] + for raw in records: + tokens = raw.split() + if not tokens: + continue + status, values = _real_row(tokens, 2) + if status == 2: + raise _AncillaryError(f'WaterKin WaveKinFile row "{raw}" is not "time elevation"') + if status != 0: + if times: + raise _AncillaryError("WaterKin WaveKinFile contains a non-data row in its history") + continue + times.append(values[0]) + heights.append(values[1]) + if len(times) < 4: + raise _AncillaryError("WaterKin WaveKinFile needs at least four `time elevation` rows") + tolerance = 1000.0 * _NATIVE_EPSILON * max(1.0, abs(times[-1])) + if abs(times[0]) > tolerance: + raise _AncillaryError("WaterKin WaveKinFile time history must start at t = 0") + if not _increasing(times): + raise _AncillaryError("WaterKin WaveKinFile times must be strictly increasing") + span = window if window > 0.0 else times[-1] + if not span / dt_wave <= _MAX_WAVE_SAMPLES: + raise _AncillaryError("WaterKin TMax/dtWave needs more than 1000000 wave samples") + wanted = math.ceil(span / dt_wave) + samples = _padded_fft_size(wanted) + series = np.zeros(samples) + series[:wanted] = np.interp(np.arange(wanted) * dt_wave, times, heights, right=0.0) + kept = samples // 2 + if cutoff >= 0.0: + step = 2.0 * math.pi / (samples * dt_wave) + kept = min(kept, math.floor(cutoff / step * (1.0 + 1.0e-15))) + if kept < 1: + raise _AncillaryError("WvHiCOff keeps no component of the WaterKin WaveKinFile record") + mean = float(series.sum() / samples) + amplitudes = np.abs(np.fft.rfft(series - mean))[1 : kept + 1] * (2.0 / samples) + if samples % 2 == 0 and kept == samples // 2: + amplitudes[-1] *= 0.5 + top = 2.0 * math.pi * kept / (samples * dt_wave) + if float(amplitudes.max()) > 1.0e3 or top > 1.0e4: + raise _AncillaryError( + "the wave components (WaterKin WaveKinMod 1 or WaveKin 3/7) need amplitudes of at " + "most 1e3 m and frequencies of at most 1e4 rad/s" + ) + if depth is not None and depth + mean <= 0.0: + raise _AncillaryError("WaterKin mean elevation leaves a non-positive effective depth") + + +def _check_waterkin_current(records: list[str]) -> None: + """Mirror the CurrentMod 1 table of ``read_waterkin_records`` (records after CurrentMod).""" + if len(records) < 2: + raise _AncillaryError("WaterKin file ended inside the current-table header") + malformed = "malformed WaterKin current row (expected z ux uy as plain numbers)" + rows: list[list[float]] = [] + start = 2 + for offset in range(4): + if start + offset >= len(records): + break + status, values = _real_row(records[start + offset].split(), 3) + if status == 2: + raise _AncillaryError(malformed) + if status == 0: + rows.append(values) + start += offset + 1 + break + if not rows: + raise _AncillaryError("could not find the WaterKin current profile table") + for record in records[start:]: + status, values = _real_row(record.split(), 3) + if status == 2: + raise _AncillaryError(malformed) + if status != 0: + break + rows.append(values) + if len(rows) < 2: + raise _AncillaryError("WaterKin CurrentMod 1 needs at least two `z ux uy` depth rows") + depths = [row[0] for row in rows] + if any(value > 0.0 for value in depths): + if not all(value >= 0.0 for value in depths): + raise _AncillaryError("WaterKin current depth column mixes signs") + depths = [-value for value in depths] + if depths[0] > depths[-1]: + depths.reverse() + if any(b - a <= 100.0 * _NATIVE_EPSILON for a, b in itertools.pairwise(depths)): + raise _AncillaryError("WaterKin current depths must be strictly monotonic") + _check_current_bounds(depths, [value for row in rows for value in row[1:]]) + + +def _check_current_bounds(depths: Sequence[float], velocities: Iterable[float]) -> None: + if any(abs(value) > _NATIVE_MAX_COORDINATE for value in depths) or any( + abs(value) > _MAX_FLUID_SPEED for value in velocities + ): + raise _AncillaryError( + "current profile depths must be at most 1e6 m and velocities at most 1e3 m/s in " + "magnitude" + ) + + +def _wave_frequency_rows(records: list[str], label: str) -> list[list[float]]: + """Mirror ``read_table``: rows `omega Re Im beta`, the direction zero when left out.""" + rows: list[list[float]] = [] + for number, record in enumerate(records, start=1): + tokens = record.split() + if not tokens: + continue + status, values = _real_row(tokens, 4) + if status != 0: + status, values = _real_row(tokens, 3) + values.append(0.0) + if status != 0: + raise DeckFormatError(f"{label} line {number}: row is not 3 to 4 plain numbers") + rows.append(values) + return rows + + +def _table_tokens(record: str) -> list[str]: + """Mirror ``row_tokens``: the tokens up to the first one holding a ``#``.""" + tokens: list[str] = [] + for token in record.split(): + if "#" in token: + break + tokens.append(token) + return tokens + + +def _table_number(token: str, name: str, number: int) -> float: + value = _token_real(token) + if value is None: + raise _AncillaryError(f'{name} line {number}: value "{token}" is not a plain finite number') + return value + + +def _check_current_profile(records: list[str], name: str) -> None: + """Mirror ``Currents 1``: `z ux uy uz` rows after their header rows.""" + rows: list[list[float]] = [] + for number, raw in enumerate(records, start=1): + tokens = raw.split() + if not tokens: + continue + status, values = _real_row(tokens, 4) + if status != 0: + status, values = _real_row(tokens, 3) + if status == 1 and not rows: + continue + if status != 0: + raise _AncillaryError(f"{name} line {number}: row is not 3 to 4 plain numbers") + rows.append(values) + if not rows: + raise _AncillaryError(f"Currents 1: {name} has no `z ux uy uz` rows") + depths = [row[0] for row in rows] + if not _increasing(depths): + raise _AncillaryError(f"Currents 1: {name} depths must be strictly increasing") + if any(b - a <= 100.0 * _NATIVE_EPSILON for a, b in itertools.pairwise(depths)): + raise _AncillaryError("current profile depths must be strictly increasing") + _check_current_bounds(depths, [value for row in rows for value in row[1:]]) + + +def _check_current_resample( + times: list[float], spatial: int, velocities: Any, coordinates: Iterable[float] +) -> None: + """Mirror ``CD_Current_Grid_Resample`` and the grid bounds; velocities is (nt, points).""" + import numpy as np + + step = min(b - a for a, b in itertools.pairwise(times)) + if times[-1] < 0.0: + raise _AncillaryError("current grid: the last time must not be negative") + if not times[-1] / _MAX_CURRENT_SAMPLES < step: + raise _AncillaryError( + "current grid: the last time over the smallest time interval exceeds the sample limit" + ) + samples = math.floor(times[-1] / step) + 1 + if samples * spatial > _MAX_CURRENT_SAMPLES: + raise _AncillaryError("current grid: more than 2e6 resampled space-time samples") + at = np.arange(samples) * step + upper = np.clip(np.searchsorted(times, at, side="right"), 1, len(times) - 1) + lower = upper - 1 + knots = np.asarray(times) + weight = np.clip((at - knots[lower]) / (knots[upper] - knots[lower]), 0.0, 1.0)[:, None] + resampled = (1.0 - weight) * velocities[lower] + weight * velocities[upper] + if any(abs(value) > _NATIVE_MAX_COORDINATE for value in coordinates) or ( + float(np.abs(resampled).max()) > _MAX_FLUID_SPEED + ): + raise _AncillaryError( + "current grid coordinates must be at most 1e6 m and velocities at most 1e3 m/s in " + "magnitude" + ) + + +def _check_dynamic_current(records: list[str], name: str) -> None: + """Mirror ``read_dynamic_profile`` (Currents 2).""" + import numpy as np + + depths: list[float] = [] + times: list[float] = [] + rows: list[list[float]] = [] + for number, raw in enumerate(records, start=1): + if len(raw) > _NATIVE_TABLE_ROW_LENGTH: + raise _AncillaryError(f"{name} line {number}: row is too long") + if number <= 4 or number == 6: + continue + tokens = _table_tokens(raw) + if number == 5: + if not tokens: + raise _AncillaryError(f"Currents 2: line 5 of {name} must list the depths") + for token in tokens: + depth = float(np.float32(_table_number(token, name, number))) + if not math.isfinite(depth): + raise _AncillaryError(f"{name} depth is not a single-precision number") + depths.append(depth) + continue + if not tokens: + continue + count = len(depths) + if len(tokens) not in (1 + count, 1 + 2 * count, 1 + 3 * count): + raise _AncillaryError(f"Currents 2: {name} line {number}: wrong number of values") + values = [_table_number(token, name, number) for token in tokens] + times.append(values[0]) + rows.append(values[1:] + [0.0] * (3 * count - len(values) + 1)) + if not depths or len(records) < 7: + raise _AncillaryError(f"Currents 2: {name} needs at least 7 lines") + if len(times) < 2: + raise _AncillaryError(f"Currents 2: {name} needs at least two time rows") + if not _increasing(depths): + raise _AncillaryError(f"Currents 2: {name} depths must be strictly increasing") + if not _increasing(times): + raise _AncillaryError(f"Currents 2: {name} times must be strictly increasing") + _check_current_resample(times, len(depths), np.asarray(rows), [*depths, 0.0]) + + +def _check_grid_current(records: list[str], name: str) -> None: + """Mirror ``read_4d_grid`` (Currents 5).""" + import numpy as np + + sizes: list[int] = [] + axes: list[list[float]] = [] + points: dict[tuple[int, int, int, int], list[float]] = {} + for number, raw in enumerate(records, start=1): + if len(raw) > _NATIVE_TABLE_ROW_LENGTH: + raise _AncillaryError(f"{name} line {number}: row is too long") + tokens = _table_tokens(raw) + if number == 1: + good = len(tokens) == 4 and all(re.fullmatch(r"[+-]?[0-9]+", t) for t in tokens) + sizes = [int(token) for token in tokens] if good else [] + if not good or any(size < 1 for size in sizes): + raise _AncillaryError(f"Currents 5: line 1 of {name} must be `nx ny nz nt`") + if math.prod(sizes) > _MAX_CURRENT_SAMPLES: + raise _AncillaryError(f"Currents 5: {name} declares more than 2e6 grid points") + elif number <= 5: + if len(tokens) != sizes[number - 2]: + raise _AncillaryError(f"Currents 5: line {number} of {name} has the wrong length") + axes.append([_table_number(token, name, number) for token in tokens]) + else: + if not tokens: + continue + if number == 6 and _NUMERIC_SHAPE.fullmatch(tokens[0]) is None: + continue + if len(tokens) != 7: + raise _AncillaryError(f"Currents 5: {name} line {number}: a row is 7 numbers") + row = [_table_number(token, name, number) for token in tokens] + where: list[int] = [] + for axis, value in zip(axes, row, strict=False): + if value not in axis: + raise _AncillaryError(f"Currents 5: {name} line {number}: off the grid") + where.append(axis.index(value)) + key = (where[0], where[1], where[2], where[3]) + if key in points: + raise _AncillaryError(f"Currents 5: {name} line {number}: repeats a grid point") + points[key] = row[4:] + if len(axes) < 4: + raise _AncillaryError(f"Currents 5: {name} needs the grid sizes and coordinates") + if not all(_increasing(axis) for axis in axes): + raise _AncillaryError(f"Currents 5: {name} coordinates must be strictly increasing") + if len(axes[3]) < 2: + raise _AncillaryError(f"Currents 5: {name} needs at least two times") + if len(points) < math.prod(sizes): + raise _AncillaryError(f"Currents 5: {name} has no row for a grid point") + spatial = math.prod(sizes[:3]) + field = np.zeros((sizes[3], spatial * 3)) + for (ix, iy, iz, it), velocity in points.items(): + column = 3 * ((ix * sizes[1] + iy) * sizes[2] + iz) + field[it, column : column + 3] = velocity + _check_current_resample(axes[3], spatial, field, [*axes[0], *axes[1], *axes[2]]) + + +def _motion_time_index(time: float, dt: float, n_time: int, tolerance: float) -> int | None: + """Mirror ``motion_grid_index``: the 0-based step, -1 past TMax, ``None`` off the grid.""" + if not abs(time / dt) < _MAX_GRID_STEPS: + return None + step = round(time / dt) + if step < 0: + return None + if step < n_time: + return step if abs(step * dt - time) <= tolerance else None + if abs(step * dt - time) <= 100.0 * _NATIVE_EPSILON * max(1.0, abs(time)): + return -1 + return None + + +def _check_motion_rows( + records: list[str], + dt: float, + tmax: float, + eligible: set[int] | None, + roll_ends: set[int] | None, + turbines: dict[int, list[int]] | None, +) -> None: + """Mirror ``read_motion_file``. + + eligible: the ids the file must move at every time (``None``: on a route that takes + other points as well, their rows are only read). roll_ends: with a line solved in torsion, + the points whose moving end is restrained in torsion. turbines: the point ids of each + TURBINES row, on a deck that has them. + """ + n_time = round(tmax / dt) + 1 + tolerance = 100.0 * _NATIVE_EPSILON * max(1.0, tmax) + what = "motionFile" + seen: set[tuple[int, int]] = set() + roll_rows: dict[int, dict[int, float]] = {} + for number, tokens in _data_rows(records, what): + where = f"{what} line {number}" + if turbines is not None: + status, row = _real_row(tokens, 21) + if status != 0: + raise _AncillaryError(f"{where}: TURBINES motion rows must be 21 numbers") + if ( + any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[2:9]) + or any(abs(value) > _MAX_MOTION_SPEED for value in row[9:15]) + or any(abs(value) > _MAX_MOTION_ACCEL for value in row[15:21]) + ): + raise _AncillaryError(f"{where}: value outside the admissible range") + turbine = int(row[1]) if row[1] == int(row[1]) else 0 + if turbine not in turbines: + raise _AncillaryError(f"{where}: J must be the integer id of a TURBINES row") + if abs(math.sqrt(sum(value * value for value in row[5:9])) - 1.0) > 1.0e-6: + raise _AncillaryError(f"{where}: the quaternion must have unit norm (within 1e-6)") + step = _motion_time_index(row[0], dt, n_time, tolerance) + if step is None: + raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") + if step == -1: + continue + if (-turbine, step) in seen: + raise _AncillaryError(f"{where}: duplicate row for turbine {turbine}") + seen.add((-turbine, step)) + seen.update((point, step) for point in turbines[turbine]) + continue + status, row = _real_row(tokens, 11) + if status != 0: + raise _AncillaryError(f"{where}: rows must be: time point_id x y z vx vy vz ax ay az") + if ( + any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[2:5]) + or any(abs(value) > _MAX_MOTION_SPEED for value in row[5:8]) + or any(abs(value) > _MAX_MOTION_ACCEL for value in row[8:11]) + ): + raise _AncillaryError(f"{where}: value outside the admissible range") + point = int(row[1]) if abs(row[1]) <= 2**31 - 1 and row[1] == int(row[1]) else 0 + if point < 1: + raise _AncillaryError(f"{where}: point_id must be a positive integer") + step = _motion_time_index(row[0], dt, n_time, tolerance) + if step is None: + raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") + if eligible is not None and point not in eligible: + raise _AncillaryError(f"{where}: references point id {point}, which is not eligible") + if step == -1: + continue + if (point, step) in seen: + raise _AncillaryError(f"{where}: duplicate row for point {point} at this time") + seen.add((point, step)) + if roll_ends is not None and len(tokens) >= 12 and _NUMERIC_SHAPE.fullmatch(tokens[11]): + roll = _token_real(tokens[11]) + if roll is None: + raise _AncillaryError(f"{where}: the roll column must be a finite number") + roll_rows.setdefault(point, {})[step] = roll + if eligible is not None and any( + (point, step) not in seen for point in eligible for step in range(n_time) + ): + raise _AncillaryError( + "motionFile must provide every Coupled/Vessel point at every dtM time" + ) + if roll_ends is None: + return + for point, rolls in roll_rows.items(): + if len(rolls) != n_time: + raise _AncillaryError( + f"motionFile: point {point} gives the roll column on some rows only" + ) + if any(rolls.values()): + if point not in roll_ends: + raise _AncillaryError( + f"motionFile: point {point} has a non-zero roll column, but no line restrained " + "in torsion has its moving end there" + ) + if rolls[0] != 0.0: + raise _AncillaryError( + f"motionFile: the roll column of point {point} must be 0 at t = 0" + ) + + +def _check_vessel_rows(records: list[str], dt: float, tmax: float) -> None: + """Mirror ``read_vessel_record``.""" + n_time = round(tmax / dt) + 1 + tolerance = 100.0 * _NATIVE_EPSILON * max(1.0, tmax) + what = "vesselMotion file" + form = 0 + seen: set[int] = set() + for number, tokens in _data_rows(records, what): + where = f"{what} line {number}" + if form == 0 and len(tokens) in (19, 20): + form = len(tokens) + status, row = (1, []) if len(tokens) != form else _real_row(tokens, form) + if status != 0: + raise _AncillaryError(f"{where}: rows must hold 19 or 20 values, one form for the file") + if any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[1:]): + raise _AncillaryError(f"{where}: value outside the admissible range") + step = _motion_time_index(row[0], dt, n_time, tolerance) + if step is None: + raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") + if step == -1: + continue + if step in seen: + raise _AncillaryError(f"{where}: duplicate row at this time") + if form == 20 and abs(math.sqrt(sum(v * v for v in row[4:8])) - 1.0) > 1.0e-6: + raise _AncillaryError(f"{where}: the quaternion must have unit norm (within 1e-6)") + seen.add(step) + if len(seen) < n_time: + raise _AncillaryError("vesselMotion file must provide one row at every dtM time") + + +def _check_rao_table(records: list[str], headings_in_use: Iterable[float]) -> None: + """Mirror ``read_rao_table``, ``CD_Vessel_RAO_Set`` and the heading test of ``_RAO_Eval``.""" + what = "vesselRAO file" + heads: list[float] = [] + blocks: list[list[float]] = [] + amplitudes: list[float] = [] + for number, tokens in _data_rows(records, what): + where = f"{what} line {number}" + if tokens[0].lower() == "heading": + value = _token_real(tokens[1]) if len(tokens) == 2 else None + if value is None or len(heads) >= 361: + raise _AncillaryError(f'{where}: a heading row must be "HEADING "') + heads.append(value) + blocks.append([]) + continue + status, row = (1, []) if len(tokens) != 13 else _real_row(tokens, 13) + if status != 0 or not heads or len(blocks[-1]) >= 2000: + raise _AncillaryError(f'{where}: data rows follow a "HEADING " row with 13 values') + blocks[-1].append(row[0]) + amplitudes.extend(row[1::2]) + if not heads: + raise _AncillaryError('vesselRAO file holds no "HEADING " block') + if any(not block for block in blocks): + raise _AncillaryError("a vesselRAO heading block holds no period rows") + if any(block != blocks[0] for block in blocks): + raise _AncillaryError("every vesselRAO heading block must list the same periods") + periods = sorted(blocks[0]) + heads = sorted(heads) + if periods[0] <= 0.0: + raise _AncillaryError("vesselRAO file: RAO periods must be positive") + if any(value < 0.0 for value in amplitudes): + raise _AncillaryError("vesselRAO file: RAO amplitudes must be non-negative") + if not _increasing(periods): + raise _AncillaryError("vesselRAO file: RAO periods must be strictly ascending") + if not _increasing(heads): + raise _AncillaryError("vesselRAO file: RAO headings must be strictly ascending") + if heads[-1] - heads[0] >= 360.0: + raise _AncillaryError("vesselRAO file: RAO headings must lie within one 360 deg turn") + for heading in headings_in_use: + turned = heads[0] + (heading - heads[0]) % 360.0 + if turned > heads[-1] + 1.0e-9 and abs(turned - 360.0 - heads[0]) <= 1.0e-9: + turned = heads[0] + if len(heads) == 1: + if abs(turned - heads[0]) > 1.0e-9: + raise _AncillaryError("relative wave heading is not in the one-heading RAO table") + elif turned > heads[-1] + 1.0e-9: + raise _AncillaryError("relative wave heading lies outside the RAO table heading range") + + +def _floor_axis(axis: list[float], q: float) -> tuple[int, float, bool]: + """Mirror ``locate_axis`` (0-based): the cell, the fraction and whether q is off the grid.""" + if q <= axis[0]: + return 0, 0.0, q < axis[0] + if q >= axis[-1]: + return len(axis) - 2, 1.0, q > axis[-1] + cell = bisect.bisect_left(axis, q) - 1 + return cell, (q - axis[cell]) / (axis[cell + 1] - axis[cell]), False + + +def _crease_band(axis: list[float], cell: int, q: float) -> tuple[int, float, float]: + """Mirror ``crease_band`` (0-based; -1: none).""" + if cell > 0: + band = _BATHY_CREASE_FRACTION * min( + axis[cell] - axis[cell - 1], axis[cell + 1] - axis[cell] + ) + if q - axis[cell] < band: + return cell, q - axis[cell], band + if cell + 2 < len(axis): + band = _BATHY_CREASE_FRACTION * min( + axis[cell + 1] - axis[cell], axis[cell + 2] - axis[cell + 1] + ) + if axis[cell + 1] - q < band: + return cell + 1, q - axis[cell + 1], band + return -1, 0.0, 0.0 + + +def _bathymetry_floor( + xs: list[float], ys: list[float], depth: dict[tuple[int, int], float], x: float, y: float +) -> tuple[float, float]: + """Mirror ``floor_eval``: the rounded and the bilinear seabed elevation at (x, y).""" + ix, tx, clamp_x = _floor_axis(xs, x) + iy, ty, clamp_y = _floor_axis(ys, y) + grid = ( + -(1.0 - tx) * (1.0 - ty) * depth[ix, iy] + - tx * (1.0 - ty) * depth[ix + 1, iy] + - (1.0 - tx) * ty * depth[ix, iy + 1] + - tx * ty * depth[ix + 1, iy + 1] + ) + if len(xs) <= 2 and len(ys) <= 2: + return grid, grid + + def x_jump(k: int, j: int, band: float) -> tuple[float, bool]: + jump = (depth[k, j] - depth[k + 1, j]) / (xs[k + 1] - xs[k]) - ( + depth[k - 1, j] - depth[k, j] + ) / (xs[k] - xs[k - 1]) + active = abs(jump) * band / 6.0 > _BATHY_CREASE_MIN_HEIGHT + return (jump if active else 0.0), active + + def y_jump(i: int, k: int, band: float) -> tuple[float, bool]: + jump = (depth[i, k] - depth[i, k + 1]) / (ys[k + 1] - ys[k]) - ( + depth[i, k - 1] - depth[i, k] + ) / (ys[k] - ys[k - 1]) + active = abs(jump) * band / 6.0 > _BATHY_CREASE_MIN_HEIGHT + return (jump if active else 0.0), active + + kx, u, band_x = (-1, 0.0, 0.0) if clamp_x else _crease_band(xs, ix, x) + ky, v, band_y = (-1, 0.0, 0.0) if clamp_y else _crease_band(ys, iy, y) + if kx < 0 and ky < 0: + return grid, grid + floor = grid + hx = 0.0 + if kx >= 0: + hx = (band_x - abs(u)) ** 3 / (6.0 * band_x * band_x) + jx0, on0 = x_jump(kx, iy, band_x) + jx1, on1 = x_jump(kx, iy + 1, band_x) + if on0 or on1: + floor += hx * ((1.0 - ty) * jx0 + ty * jx1) + if ky >= 0: + hy = (band_y - abs(v)) ** 3 / (6.0 * band_y * band_y) + jy0, on0 = y_jump(ix, ky, band_y) + jy1, on1 = y_jump(ix + 1, ky, band_y) + mixed = 0.0 + on_mixed = False + if kx >= 0: + below, _ = x_jump(kx, ky - 1, band_x) + centre, _ = x_jump(kx, ky, band_x) + above, _ = x_jump(kx, ky + 1, band_x) + mixed = (above - centre) / (ys[ky + 1] - ys[ky]) - (centre - below) / ( + ys[ky] - ys[ky - 1] + ) + on_mixed = abs(mixed) * band_x * band_y / 36.0 > _BATHY_CREASE_MIN_HEIGHT + if not on_mixed: + mixed = 0.0 + if on0 or on1 or on_mixed: + floor += hy * ((1.0 - tx) * jy0 + tx * jy1 + hx * mixed) + return floor, grid + + +def _check_bathymetry( + records: list[str], fixed_points: list[tuple[int, float, float, float]] +) -> None: + """Mirror ``load_bathymetry_file`` and ``check_fixed_points_on_bathymetry``.""" + what = "bathymetry file" + rows: list[list[float]] = [] + for number, tokens in _data_rows(records, what): + if len(tokens) != 3: + raise _AncillaryError(f'{what} line {number}: rows must be "x y depth"') + status, row = _real_row(tokens, 3) + if status != 0: + raise _AncillaryError(f"{what} line {number}: malformed row") + rows.append(row) + if len(rows) < 4: + raise _AncillaryError("bathymetry file needs at least a 2x2 grid") + if any(row[2] <= 0.0 for row in rows): + raise _AncillaryError("bathymetry file depths must be positive") + + def unique(values: list[float]) -> list[float]: + kept: list[float] = [] + for value in sorted(values): + if not kept or abs(value - kept[-1]) > 16.0 * _NATIVE_EPSILON * max( + 1.0, abs(value), abs(kept[-1]) + ): + kept.append(value) + return kept + + def nearest(axis: list[float], value: float) -> int: + at = bisect.bisect_left(axis, value) + return min( + (i for i in (at - 1, at) if 0 <= i < len(axis)), key=lambda i: abs(axis[i] - value) + ) + + xs = unique([row[0] for row in rows]) + ys = unique([row[1] for row in rows]) + if len(xs) < 2 or len(ys) < 2 or len(rows) != len(xs) * len(ys): + raise _AncillaryError("bathymetry file must define one complete rectangular grid") + depth: dict[tuple[int, int], float] = {} + for x, y, value in rows: + key = (nearest(xs, x), nearest(ys, y)) + if key in depth: + raise _AncillaryError("bathymetry file has duplicate x/y entries") + depth[key] = value + for point, x, y, z in fixed_points: + floor, grid = _bathymetry_floor(xs, ys, depth, x, y) + tolerance = max(1.0e-3, 1.0e-5 * abs(floor)) * (1.0 + 1.0e-9) + if abs(z - floor) <= tolerance or abs(z - grid) <= tolerance or z > floor: + continue + raise _AncillaryError(f"POINT {point} (Fixed) lies below the seabed: z = {z:g}") + + +def _read_syrope_settings(records: list[str]) -> tuple[str, str, float, float]: + """Mirror the settings rows of ``read_syrope_settings``: OWC path, WCType, k1, k2.""" + what = "Syrope settings file" + found: dict[str, str] = {} + for number, tokens in _data_rows(records, what): + if len(tokens) < 2: + continue + name = tokens[1].lower() + value = _strip_quotes(tokens[0]) + if name not in {"owc", "wctype", "k1", "k2"}: + continue + if name in found: + raise _AncillaryError(f"{what} line {number}: duplicate Syrope {tokens[1]} row") + if name == "wctype" and value.lower() not in {"linear", "quadratic", "exp"}: + raise _AncillaryError( + f"{what} line {number}: Syrope WCType must be LINEAR, QUADRATIC, or EXP" + ) + if name in {"k1", "k2"} and _token_real(value) is None: + raise _AncillaryError(f"{what} line {number}: Syrope {tokens[1]} is not a plain number") + found[name] = value + if len(found) < 4: + raise _AncillaryError("the Syrope settings file needs OWC, WCType, k1, and k2 rows") + k1 = _token_real(found["k1"]) + k2 = _token_real(found["k2"]) + assert k1 is not None and k2 is not None + return found["owc"], found["wctype"].lower(), k1, k2 + + +def _read_syrope_table(records: list[str]) -> list[list[float]]: + """Mirror the OWC table rows of ``read_syrope_settings``.""" + what = "Syrope OWC table" + rows: list[list[float]] = [] + for number, tokens in _data_rows(records, what): + if not any(_NUMERIC_SHAPE.fullmatch(token) for token in tokens): + continue + status, row = _real_row(tokens, 2) + if status != 0: + raise _AncillaryError(f"{what} line {number}: malformed row") + rows.append(row) + if len(rows) < 2: + raise _AncillaryError("the Syrope OWC table needs at least two rows") + return rows + + +def _syrope_curve_defect( + strain: list[float], + tension: list[float], + shape: str, + k1: float, + k2: float, + alpha: float, + beta: float, + t_max: float, +) -> bool: + """True when ``CD_Syrope_Working_Curve`` fails at the running maximum ``t_max``.""" + + def fast(value: float) -> float: + return math.log(1.0 + beta / alpha * value) / beta + + if not t_max > 0.0 or t_max < tension[0] or t_max > tension[-1]: + return True + at = bisect.bisect_left(tension, t_max) + if at == 0: + eps_max = strain[0] + elif at == len(tension): + eps_max = strain[-1] + else: + eps_max = strain[at - 1] + (strain[at] - strain[at - 1]) * (t_max - tension[at - 1]) / ( + tension[at] - tension[at - 1] + ) + eps_min = strain[0] + k1 * (eps_max - strain[0]) + if shape == "linear" and k1 >= 1.0: + eps_min = eps_max - t_max / k1 + if not (eps_min >= 0.0 and eps_max - eps_min > 0.0): + return True + previous: tuple[float, float] | None = None + for index in range(_SYROPE_CURVE_SAMPLES): + at_strain = eps_min + (eps_max - eps_min) * index / (_SYROPE_CURVE_SAMPLES - 1) + xi = (at_strain - eps_min) / (eps_max - eps_min) + if shape == "linear": + value = t_max * xi + elif shape == "quadratic": + value = t_max * xi * (k2 * xi + (1.0 - k2)) + else: + value = t_max * (1.0 - math.exp(k2 * xi)) / (1.0 - math.exp(k2)) + argument = 1.0 + beta / alpha * value + if argument <= 0.0: + return True + slow = at_strain - fast(value) + if previous is not None and (value <= previous[0] or slow <= previous[1]): + return True + previous = (value, slow) + return False + + +def _check_syrope_curve( + rows: list[list[float]], + shape: str, + k1: float, + k2: float, + alpha: float, + beta: float, +) -> None: + """Mirror ``CD_Syrope_Init`` (the constants BA_s and BA_d are checked with the row).""" + strain = [row[0] for row in rows] + tension = [row[1] for row in rows] + if not (_increasing(strain) and _increasing(tension)): + raise _AncillaryError("the OWC strain and tension columns must be strictly increasing") + if shape == "exp" and not k2 > 0.0: + raise _AncillaryError("the EXP working curve needs a positive shape parameter p2 (k2 > 0)") + if shape == "linear" and k1 < 0.0: + raise _AncillaryError("the LINEAR working curve requires k1 >= 0") + if shape == "quadratic" and (k1 < 0.0 or k1 >= 1.0 or k2 <= 0.0 or k2 > 1.0): + raise _AncillaryError("the QUADRATIC working curve requires 0 <= k1 < 1 and 0 < k2 <= 1") + if shape == "exp" and (k1 < 0.0 or k1 >= 1.0 or k2 < 1.0e-6 or k2 > 700.0): + raise _AncillaryError("the EXP working curve requires 0 <= k1 < 1 and 1e-6 <= k2 <= 700") + if any(1.0 + beta / alpha * value <= 0.0 for value in tension): + raise _AncillaryError("the OWC slow-spring static strain is non-finite") + slow = [ + e - math.log(1.0 + beta / alpha * t) / beta for e, t in zip(strain, tension, strict=True) + ] + if not _increasing(slow): + raise _AncillaryError("the OWC slow-spring static strain must be strictly increasing") + low = max(tension[0], 0.0) + high = tension[-1] + if all( + _syrope_curve_defect( + strain, + tension, + shape, + k1, + k2, + alpha, + beta, + low + (high - low) * k / _SYROPE_SCAN_SAMPLES, + ) + for k in range(1, _SYROPE_SCAN_SAMPLES + 1) + ): + raise _AncillaryError("no running maximum in the OWC table range admits a working curve") + + +def _train_direction(values: tuple[str, ...]) -> tuple[float, float]: + """Mirror ``CD_Sea_Parse_Train``: the direction and spreading of a wave row's values.""" + kind = values[0].lower() + if kind == "airy": + return _native_float(values[3]), 0.0 + index = {"jonswap": 4}.get(kind, 7 if kind.replace("-", "") in {"ochihubble", "ochi"} else 3) + spreading = _native_float(values[index + 1]) if len(values) > index + 1 else 0.0 + return _native_float(values[index]), spreading + + +def _spread_headings(direction: float, spreading: float, bins: int) -> list[float]: + """The headings of a train (``CD_Sea_Synthesise``): the bin centres of a spread one.""" + if spreading <= 0.0: + return [direction] + return [direction - 90.0 + 180.0 * (j + 0.5) / bins for j in range(bins)] + + @dataclass(frozen=True) class DeckRecord: """One parsed table row and its zero-based source-line location. @@ -2438,17 +3296,13 @@ def validate(self) -> None: between rows. A deck that passes can still be rejected by the solver on physical grounds. - A WaterKin file's ``WaveKinMod`` (which decides whether ``vesselRAO``, - ``WaveRefX`` and ``WaveRefY`` have waves to act on) and the rows of - ``wave_frequencies.txt`` are read here, from the deck's folder. A file - that is not there, as for a deck held in memory under its nominal path, - is not an error at this stage: the deck passes these checks and the - solver, which must open the file, refuses it by name if it is still - missing. The other ancillary files (motion, bathymetry, Syrope and the - other MoorDyn-C kinematics files) are not read here, so the checks that - depend on their contents stay with the solver: for example the anchor - height on a bathymetry surface, and whether a WaterKin current doubles - a ``current`` row or ``Currents 1``, ``2`` or ``5``. + The ancillary files a deck names (WaterKin and its wave record, the + MoorDyn-C kinematics files, motion and vessel records, the RAO table, + bathymetry, and the Syrope settings and OWC table) are read from the + deck's folder, with the rules of their native readers. A file that is + not there, as for a deck held in memory under its nominal path, is not + an error at this stage: the deck passes these checks and the solver, + which must open the file, refuses it by name if it is still missing. Raises ------ @@ -4361,15 +5215,353 @@ def is_body_kind(kind: str, pinned: bool) -> bool: ) continue raise DeckFormatError(f"{self._label}: unsupported OUTPUT channel {channel!r}") + try: + self._check_ancillary_files( + option_by_key, + options, + motion_source, + mdc_wave_kin, + mdc_currents, + has_waves or bool(trains), + current is not None and current.values[0].lower() != "none", + ( + { + node + for line_id in torsion_lines + for node, end in zip( + normalized_line_endpoints[line_id], ("a", "b"), strict=True + ) + if torsion_ends.get(line_id, {}).get(end, False) + and point_types[node] != "fixed" + } + if torsion_lines + else None + ), + bool(finite_length_rods), + has_rigid6, + ) + except _AncillaryError as exc: + raise DeckFormatError(f"{self._label}: {exc}") from None - def _ancillary_records(self, value: str) -> list[str] | None: - """Return the records of the ancillary file ``value`` names, ``None`` if it is not there. + def _check_ancillary_files( + self, + option_by_key: dict[str, OptionRecord], + options: tuple[OptionRecord, ...], + motion_source: str | None, + mdc_wave_kin: int, + mdc_currents: int, + has_sea: bool, + has_current: bool, + roll_ends: set[str] | None, + has_rods: bool, + has_rigid6: bool, + ) -> None: + """Run the native rules of the ancillary files the deck names, those that are there.""" + + def number(key: str) -> float | None: + row = option_by_key.get(key) + return None if row is None else _native_float(row.values[0]) + + dt = number("dtm") + tmax = number("tmax") or 0.0 + depth = number("water_depth") + record_dir: float | None = None + component_headings: list[float] | None = None + waterkin_row = option_by_key.get("waterkin") + if waterkin_row is not None and not mdc_wave_kin: + value = waterkin_row.values[0] + if value.lower() not in {"", "none", "seastate"} and any( + character in _WATERKIN_FILENAME_LETTERS for character in value + ): + record_dir = self._check_waterkin_file( + value, tmax, dt, depth, has_sea, has_current, mdc_currents, option_by_key + ) + if mdc_wave_kin == 3: + component_headings = [0.0] + elif mdc_wave_kin == 7: + component_headings = self._wave_frequency_headings() + if mdc_currents: + name = _MOORDYN_C_CURRENT_FILES[mdc_currents] + records = self._ancillary_records(name) + if records is not None: + check = { + 1: _check_current_profile, + 2: _check_dynamic_current, + 5: _check_grid_current, + }[mdc_currents] + check(records, name) + rows = self._rows("POINTS") + for key in ("bathymetry",): + row = option_by_key.get(key) + if row is None: + continue + records = self._ancillary_records(row.values[0]) + if records is not None: + _check_bathymetry( + records, + [ + ( + _native_int(point.tokens[0]), + _native_float(point.tokens[2]), + _native_float(point.tokens[3]), + _native_float(point.tokens[4]), + ) + for point in rows + if point.tokens[1].lower() == "fixed" + ], + ) + self._check_syrope_files() + if motion_source is None or dt is None or self.caller_driven: + return + source = next(o for o in reversed(options) if o.keyword.lower() == motion_source) + records = self._ancillary_records(source.values[0]) + if records is None: + return + if motion_source == "vesselmotion": + if has_rods or has_rigid6: + raise _AncillaryError( + "vesselMotion/vesselRAO move Coupled/Vessel line points; decks with prescribed " + "rods or Rigid6 bodies take motionFile rows" + ) + _check_vessel_rows(records, dt, tmax) + elif motion_source == "vesselrao": + if has_rods or has_rigid6: + raise _AncillaryError( + "vesselMotion/vesselRAO move Coupled/Vessel line points; decks with prescribed " + "rods or Rigid6 bodies take motionFile rows" + ) + headings = component_headings + if record_dir is not None: + headings = [record_dir] + if headings is None: + headings = self._sea_headings(option_by_key, options) + _check_rao_table(records, headings) + else: + self._check_motion_file(records, dt, tmax, roll_ends, has_rods, has_rigid6) + + def _check_waterkin_file( + self, + value: str, + tmax: float, + dt: float | None, + depth: float | None, + has_sea: bool, + has_current: bool, + mdc_currents: int, + option_by_key: dict[str, OptionRecord], + ) -> float | None: + """Mirror ``read_waterkin_records``; returns the WaveDir of a WaveKinMod 1 record.""" + records = self._ancillary_records(value) + mode = self._waterkin_wave_mode(value) + if records is None or mode is None: + return None - As the native ``path_join``: a relative name is taken from the deck's folder. + def line(index: int, name: str) -> str: + if index >= len(records): + raise _AncillaryError(f"WaterKin file has no {name} line") + return records[index] + + tokens = line(4, "WaveKinFile").split() + wave_file = _strip_quotes(tokens[0]) if tokens else "" + tokens = line(5, "dtWave").split() + dt_wave = _token_real(tokens[0]) if tokens else None + if dt_wave is None or dt_wave < 0.0: + raise _AncillaryError("WaterKin dtWave must be finite and non-negative") + tokens = line(6, "WaveDir").split() + wave_dir = _token_real(tokens[0]) if tokens else None + if wave_dir is None: + raise _AncillaryError("WaterKin WaveDir must be finite") + if len(records) < 13: + raise _AncillaryError("WaterKin file ended inside the waves block") + if mode == 1: + if self.caller_driven: + raise _AncillaryError( + "WaterKin WaveKinMod 1 is supported only by the standalone driver; coupled " + "aggregate runs must use WaveKinMod 2/SEASTATE through external_fluid" + ) + if not wave_file: + raise _AncillaryError("WaterKin WaveKinMod 1 requires WaveKinFile") + if dt_wave <= 0.0: + raise _AncillaryError("WaterKin WaveKinMod 1 requires positive dtWave") + if has_sea: + raise _AncillaryError( + "the deck declares BOTH a waves OPTION and WaterKin WaveKinMod 1 " + "(double-counting); keep one" + ) + if dt is None or depth is None: + raise _AncillaryError("deck waves require dtM, TMax and WtrDpth") + history = self._ancillary_records(wave_file) + if history is None: + history = self._ancillary_records( + str(self._ancillary_path(value).parent / wave_file) + ) + if history is not None: + cutoff_row = option_by_key.get("wvhicoff") + cutoff = math.pi + if cutoff_row is not None: + cutoff = ( + -1.0 + if cutoff_row.values[0].lower() == "none" + else _native_float(cutoff_row.values[0]) + ) + _check_elevation_record(history, dt_wave, tmax, cutoff, depth) + elif mode == 2 and not self.caller_driven: + raise _AncillaryError( + "WaterKin WaveKinMod 2/SEASTATE is coupled-only; the standalone driver has no host " + "SeaState field" + ) + line(13, "current header") + tokens = line(14, "CurrentMod").split() + if not tokens or re.fullmatch(r"[+-]?[0-9]{1,18}", tokens[0]) is None: + raise _AncillaryError("WaterKin CurrentMod must be an integer (0, 1, or 2)") + current_mode = int(tokens[0]) + if not 0 <= current_mode <= 2: + raise _AncillaryError( + "WaterKin CurrentMod supports 0 (none), 1 (depth profile), or 2 (SeaState hybrid)" + ) + if current_mode == 2: + if mode == 1: + raise _AncillaryError( + "WaterKin WaveKinMod 1 plus CurrentMod 2 mixes self-driven waves with a host " + "SeaState current, which the nodewise host field cannot separate" + ) + if not self.caller_driven: + raise _AncillaryError( + "WaterKin WaveKinMod 2/SEASTATE is coupled-only; the standalone driver has no " + "host SeaState field" + ) + if current_mode == 1: + _check_waterkin_current(records[15:]) + if has_current: + raise _AncillaryError( + "the deck declares BOTH a current OPTION and a WaterKin current " + "(double-counting); keep one" + ) + if mdc_currents: + raise _AncillaryError( + "the deck declares BOTH a current OPTION (or a WaterKin current) and " + f"Currents {mdc_currents} (double-counting); keep one" + ) + return wave_dir if mode == 1 else None + + def _wave_frequency_headings(self) -> list[float] | None: + """The directions [deg] of the components of ``wave_frequencies.txt``, if it is there.""" + records = self._ancillary_records(_MOORDYN_C_WAVE_FILES[7]) + if records is None: + return None + headings = [ + math.degrees(math.atan2(math.sin(beta), math.cos(beta))) + for omega, real_part, imaginary, beta in _wave_frequency_rows(records, "") + if omega > 0.0 and math.hypot(real_part, imaginary) > 0.0 + ] + return headings + + def _sea_headings( + self, option_by_key: dict[str, OptionRecord], options: tuple[OptionRecord, ...] + ) -> list[float]: + """The headings of a deck's own sea, as ``CD_Sea_Synthesise`` lays them out.""" + bins_row = option_by_key.get("wave_directions") + bins = _SEA_DEFAULT_DIRECTIONS + if bins_row is not None and round(_native_float(bins_row.values[0])) > 0: + bins = round(_native_float(bins_row.values[0])) + spreading_row = option_by_key.get("wave_spreading") + spreading = 0.0 if spreading_row is None else _native_float(spreading_row.values[0]) + headings: list[float] = [] + waves = option_by_key.get("waves") + if waves is not None and waves.values[0].lower() != "none": + direction, _ = _train_direction(waves.values) + headings.extend(_spread_headings(direction, spreading, bins)) + for option in options: + if option.keyword.lower() in _WAVETRAIN_KEYWORDS: + direction, own = _train_direction(option.values) + headings.extend(_spread_headings(direction, own, bins)) + return headings + + def _check_syrope_files(self) -> None: + """Mirror ``read_syrope_settings``, and ``CD_Syrope_Init`` for the types lines use.""" + used = {row.tokens[1].lower() for row in self.sections} + for row in self.line_types: + parts = row.tokens[3].split("|") + if not parts[0].lower().startswith("syrope:"): + continue + records = self._ancillary_records(parts[0][len("syrope:") :]) + if records is None: + continue + owc, shape, k1, k2 = _read_syrope_settings(records) + settings = self._ancillary_path(parts[0][len("syrope:") :]) + name = _strip_quotes(owc) + absolute = name.startswith(("/", "\\")) or name[1:2] == ":" + table = self._ancillary_records(name if absolute else str(settings.parent / name)) + if table is None: + continue + rows = _read_syrope_table(table) + if row.tokens[0].lower() in used: + _check_syrope_curve( + rows, shape, k1, k2, _native_float(parts[1]), _native_float(parts[2]) + ) + + def _check_motion_file( + self, + records: list[str], + dt: float, + tmax: float, + roll_ends: set[str] | None, + has_rods: bool, + has_rigid6: bool, + ) -> None: + """Mirror ``read_motion_file`` on the route the native dispatch takes. + + Only the routes with lines read the file, and only those of lines alone its roll column. """ + if not self.lines and (has_rods or has_rigid6): + return + moved = {"coupled", "vessel"} + ids: dict[int, str] = { + _native_int(row.tokens[0]): row.tokens[1].lower() for row in self.points + } + turbine_rows = self._rows("TURBINES") + eligible: set[int] | None = { + point for point, kind in ids.items() if kind in moved or _is_turbine_point_type(kind) + } + if has_rods and has_rigid6: + eligible = None + elif has_rods: + assert eligible is not None + eligible |= {point for point, kind in ids.items() if _rod_point(kind) is not None} + elif has_rigid6: + assert eligible is not None + eligible |= {point for point, kind in ids.items() if kind.startswith("body")} + if eligible is not None and not eligible: + raise _AncillaryError("motionFile requires at least one Coupled/Vessel point") + turbines: dict[int, list[int]] | None = None + if turbine_rows: + turbines = {_native_int(row.tokens[0]): [] for row in turbine_rows} + for point, kind in ids.items(): + if _is_turbine_point_type(kind): + turbines.setdefault(_turbine_id(kind), []).append(point) + _check_motion_rows( + records, + dt, + tmax, + eligible, + None + if roll_ends is None or has_rods or has_rigid6 + else {int(point) for point in roll_ends}, + turbines, + ) + + def _ancillary_path(self, value: str) -> Path: + """Return the file an ancillary value names, as the native ``path_join``.""" name = _strip_quotes(value) absolute = name.startswith(("/", "\\")) or name[1:2] == ":" - source = Path(name) if absolute else self.path.parent / name + return Path(name) if absolute else self.path.parent / name + + def _ancillary_records(self, value: str) -> list[str] | None: + """Return the records of the ancillary file ``value`` names, ``None`` if it is not there. + + A relative name is taken from the deck's folder. + """ + source = self._ancillary_path(value) try: payload = source.read_bytes() except OSError: @@ -4420,18 +5612,7 @@ def _check_wave_frequencies(self, depth: float) -> None: if records is None: return label = f"{self._label}: WaveKin 7: wave_frequencies.txt" - rows: list[tuple[float, float, float, float]] = [] - for number, record in enumerate(records, start=1): - tokens = record.split() - if not tokens: - continue - values = [_plain_real(token) for token in tokens[:4]] - if len(values) == 4 and values[3] is None: - values = values[:3] - if len(values) < 3 or any(value is None for value in values): - raise DeckFormatError(f"{label} line {number}: row is not 3 to 4 plain numbers") - omega, real_part, imaginary, *beta = values - rows.append((omega, real_part, imaginary, beta[0] if beta else 0.0)) + rows = _wave_frequency_rows(records, label) if len(rows) < 2: raise DeckFormatError(f"{label} needs at least two `omega Re Im` rows") if rows[0][0] != 0.0: diff --git a/python/tests/test_deck_file.py b/python/tests/test_deck_file.py index 543d24c1..25fea080 100644 --- a/python/tests/test_deck_file.py +++ b/python/tests/test_deck_file.py @@ -834,7 +834,7 @@ def _vessel_source(rows: str) -> str: def test_vessel_motion_options_parse_and_rebase(tmp_path): source = _vessel_source("vessel.dat vesselMotion\n80.0|0.0|-10.0 vesselRef") - path = tmp_path / "vessel.dat" + path = tmp_path / "deck.dat" path.write_text(source, encoding="utf-8") deck = DeckFile.read(path) @@ -2210,7 +2210,10 @@ def test_case_generation_rebases_relative_ancillary_paths(tmp_path): bathymetry = source_directory / "data" / "bathy_profile.dat" bathymetry.parent.mkdir() - bathymetry.write_text("placeholder\n", encoding="utf-8") + bathymetry.write_text( + "".join(f"{x} {y} 100.0\n" for x in (-2000.0, 2000.0) for y in (-2000.0, 2000.0)), + encoding="utf-8", + ) deck_text = (ROOT / "examples/chain_catenary_r3_100m.dat").read_text(encoding="utf-8") deck_text = deck_text.replace("100.0 WtrDpth", "data/bathy_profile.dat bathymetryFile") deck_path = source_directory / "base.dat" diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index 9059e3da..7f48943f 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -152,13 +152,320 @@ def _waterkin(mode: str) -> dict[str, str]: ] -def _stage(folder: Path, rows: tuple[str, ...], files: dict[str, str]) -> Path: +_CURRENT_TABLES = { + "2": "a\nb\nc\nd\n-40 -10\ne\n0 0.1 0.2\n5 0.2 0.3\n10 0.2 0.3\n", + "5": "1 1 2 2\n0\n0\n-40 -10\n0 10\n" + + "".join(f"0 0 {z} {t} 0.1 0 0\n" for z in (-40, -10) for t in (0, 10)), +} +_PROFILE = "z ux uy uz\n-40 0.1 0 0\n-10 0.2 0.1 0\n" +_TABLE = "-40 0.1 0.0\n-10 0.2 0.1\n" + + +def _kin(wave="0", current="0", table="", dt="0.25", direction="37.0", record='"eta.txt"') -> str: + return ( + "MoorDyn v2 water kinematics file\n(header line 2)\n--- WAVES ---\n" + f"{wave} WaveKinMod\n{record} WaveKinFile\n{dt} dtWave\n{direction} WaveDir\n" + "0 X\n0.0\n0 Y\n0.0\n0 Z\n0.0\n--- CURRENT ---\n" + f"{current} CurrentMod\nz ux uy\n(m) (m/s) (m/s)\n{table}--- need this line ---\n" + ) + + +def _hold(point=2, *, tail=(), steps=101, last=None) -> str: + rows = [f"{0.1 * i:.1f} {point} 0 0 0 0 0 0 0 0 0\n" for i in range(steps)] + return "".join(rows[:last]) + "".join(tail) + + +def _vessel(*, form=19, steps=101, tail=()) -> str: + pose = "0 0 0 0 0 0" if form == 19 else "0 0 0 1 0 0 0" + rows = [f"{0.1 * i:.1f} {pose} {'0 ' * 12}\n" for i in range(steps)] + return "".join(rows) + "".join(tail) + + +def _rao( + heads=(-180.0, 179.0), period_row="0.5 0 0 0 0 1 0 0 0 0 0 0 0\n100 0 0 0 0 1 0 0 0 0 0 0 0\n" +): + return "".join(f"HEADING {head}\n{period_row}" for head in heads) + + +def _sea_bed(depth=50.0, **changes) -> str: + rows = {(x, y): depth for x in (-1000.0, 1000.0) for y in (-1000.0, 1000.0)} + rows.update(changes.get("rows", {})) + return "".join(f"{x} {y} {value}\n" for (x, y), value in rows.items()) + changes.get("tail", "") + + +_SETTINGS = (EXAMPLES / "data" / "syrope" / "syrope_settings.dat").read_text(encoding="utf-8") +_OWC = (EXAMPLES / "data" / "syrope" / "syrope_owc.dat").read_text(encoding="utf-8") + + +def _syrope(settings=_SETTINGS, owc=_OWC) -> dict[str, str]: + return {"data/syrope/syrope_settings.dat": settings, "data/syrope/syrope_owc.dat": owc} + + +# a deck of objects and no lines: the only route that does not read a motion file +_OBJECTS = """Dry rod pinned to a heavy body +--- ROD TYPES --- +Name Diam Mass Cd Ca CdEnd CaEnd +(-) (m) (kg/m) (-) (-) (-) (-) +pend 0.2 50 0 0 0 0 +--- BODIES --- +ID Attachment X0 Y0 Z0 r0 p0 y0 Mass CG* I* Volume CdA* Ca* +(#) (-) (m) (m) (m) (deg) (deg) (deg) (kg) (m) (kg-m^2) (m^3) (m^2) (-) +1 Free 5 0 20 0 0 0 5.0e10 0 1.0e13 4.87804882926829E+07 0 0 +--- RODS --- +ID RodType Attachment XA YA ZA XB YB ZB NumSegs Outputs +(-) (-) (-) (m) (m) (m) (m) (m) (m) (-) (-) +1 pend Body1Pinned 0 0 0 5.0 0 -8.66025403784439 4 - +--- OPTIONS --- +100 WtrDpth +deck bodyIC +moordyn bodyWetting +0.05 dtM +1.0 TMax +--- OUTPUTS --- +Rod1Px Body1Pz +--- end --- +""" +_FREE = "wave_frequencies.txt" +_ONE_BETA = "0 0 0 0.5\n0.5 0.5 0 0.5\n" +_WK = ("wk.dat WaterKin",) +_C1 = "current_profile.txt" +_C2 = "current_profile_dynamic.txt" +_C5 = "current_profile_4d.txt" +_RAO_ROWS = ("rao.txt vesselRAO", "0|0|0 vesselRef") +# (deck, rows, files): every case is judged by DeckFile and by the driver +_ANCILLARY_CASES = [ + # WaterKin: the rest of the file + ("held", _WK, {"wk.dat": _kin()}), + ("held", _WK, {"wk.dat": _kin(dt="-1")}), + ("held", _WK, {"wk.dat": _kin(dt="x")}), + ("held", _WK, {"wk.dat": _kin(direction="north")}), + ("held", _WK, {"wk.dat": 'a\nb\nc\n0\n""\n0\n0\n0\n'}), + ("held", _WK, {"wk.dat": _kin(current="x")}), + ("held", _WK, {"wk.dat": _kin(current="3")}), + ("held", _WK, {"wk.dat": _kin(current="2")}), + ("held", _WK, {"wk.dat": _kin(wave="1", current="2"), "eta.txt": _ETA}), + ("held", _WK, {"wk.dat": _kin(wave="1", record='""')}), + ("held", _WK, {"wk.dat": _kin(wave="1", dt="0"), "eta.txt": _ETA}), + ("held", ("wk.dat WaterKin", "airy 2 8 0 waves"), {"wk.dat": _kin("1"), "eta.txt": _ETA}), + # the WaterKin current table, and a current doubled + ("held", _WK, {"wk.dat": _kin(current="1", table=_TABLE)}), + ("held", _WK, {"wk.dat": _kin(current="1", table="10 0.2 0.1\n40 0.1 0.0\n")}), + ( + "held", + ("wk.dat WaterKin", "uniform 0.1 0 0 current"), + {"wk.dat": _kin(current="1", table=_TABLE)}, + ), + ("held", ("wk.dat WaterKin", "uniform 0.1 0 0 current"), {"wk.dat": _kin()}), + ( + "held", + ("wk.dat WaterKin", "1 Currents"), + {"wk.dat": _kin(current="1", table=_TABLE), _C1: _PROFILE}, + ), + ("held", ("wk.dat WaterKin", "1 Currents"), {"wk.dat": _kin(), _C1: _PROFILE}), + ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n")}), + ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n10 0.2 0\n")}), + ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1 0\n-40 0.2 0\n")}), + ("held", _WK, {"wk.dat": _kin(current="1", table="-40 0.1\n-10 0.2 0\n")}), + ("held", _WK, {"wk.dat": _kin(current="1", table="-40 2000 0\n-10 0.2 0\n")}), + # the record of WaveKinMod 1 + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\n2 0\n"}), + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "1 0\n2 0\n3 0\n4 0\n"}), + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n2 0\n1 0\n3 0\n"}), + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\nend\n2 0\n3 0\n"}), + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 x\n2 0\n3 0\n"}), + ("held", _WK, {"wk.dat": _kin("1"), "eta.txt": _ETA.replace(" 0.4\n", " 4000\n")}), + ( + "held", + _WK, + { + "wk.dat": _kin("1"), + "eta.txt": _ETA.replace(" 0.4\n", " -70\n").replace(" -0.4\n", " -70\n"), + }, + ), + ("held", ("wk.dat WaterKin", "0.1 WvHiCOff"), {"wk.dat": _kin("1"), "eta.txt": _ETA}), + ("held", ("wk.dat WaterKin", "none WvHiCOff"), {"wk.dat": _kin("1"), "eta.txt": _ETA}), + # Currents 1, 2 and 5 + ("held", ("1 Currents",), {_C1: _PROFILE}), + ("held", ("1 Currents",), {_C1: "z ux uy uz\n"}), + ("held", ("1 Currents",), {_C1: "-10 0 0 0\n-40 0 0 0\n"}), + ("held", ("1 Currents",), {_C1: "-40 2000 0 0\n-10 0 0 0\n"}), + ("held", ("1 Currents",), {_C1: "-40 0 0\n-10 0 0 x\n"}), + ("held", ("1 Currents",), {_C1: "-40 0\n"}), + ("held", ("2 Currents",), {_C2: _CURRENT_TABLES["2"]}), + ("held", ("2 Currents",), {_C2: "a\nb\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 0.1 0.2\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-10 -40\ne\n0 1 1\n5 1 1\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 1 1\n0 1 1\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 1\n5 1 1\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 x 1\n5 1 1\n"}), + ("held", ("2 Currents",), {_C2: "a\nb\nc\nd\n-40 -10\ne\n0 2000 1\n5 1 1\n"}), + ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"]}), + ("held", ("5 Currents",), {_C5: "1 1 x 2\n"}), + ("held", ("5 Currents",), {_C5: "1 1 2 2\n0\n0\n-40\n0 10\n"}), + ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"] + "0 0 -40 0 0.1 0 0\n"}), + ( + "held", + ("5 Currents",), + {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1", "0 0 -41 0 0.1")}, + ), + ("held", ("5 Currents",), {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1 0 0\n", "")}), + ("held", ("5 Currents",), {_C5: "1 1 1 1\n0\n0\n-40\n0\n0 0 -40 0 0.1 0 0\n"}), + ( + "held", + ("5 Currents",), + {_C5: _CURRENT_TABLES["5"].replace("0 0 -40 0 0.1", "0 0 -40 0 2000")}, + ), + # motionFile + ("held", ("m.txt motionFile",), {"m.txt": _hold()}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("10.1 2 0 0 0 0 0 0 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(last=100)}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2 0 0 0 0 0 0 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.05 2 0 0 0 0 0 0 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 1 0 0 0 0 0 0 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2.5 0 0 0 0 0 0 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(tail=("0.0 2 0 0 0\n",))}), + ("held", ("m.txt motionFile",), {"m.txt": _hold(point=1)}), + ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 2e6 0 0 0 0 0 0 0 0\n"}), + ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 2e4 0 0 0 0 0\n"}), + ("held", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 0 0 0 2e6 0 0\n"}), + # vesselMotion + ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel()}), + ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(form=20)}), + ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(steps=100)}), + ( + "held", + ("v.txt vesselMotion", "0|0|0 vesselRef"), + {"v.txt": _vessel(form=19, tail=("0.0 " + "0 " * 19 + "\n",))}, + ), + ( + "held", + ("v.txt vesselMotion", "0|0|0 vesselRef"), + {"v.txt": _vessel(form=20).replace("0 0 0 1 0 0 0", "0 0 0 2 0 0 0", 1)}, + ), + ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": "0.05 " + "0 " * 18 + "\n"}), + ("held", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": "0.0 2e6 " + "0 " * 17 + "\n"}), + # vesselRAO, its rows and the headings of the sea + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao()}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(30.0,))}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0,))}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(40.0, 90.0))}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 90.0))}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 360.0))}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "1 0 0 0 0 1 0 0 0 0 0 0 0\n"}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n"}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30 40\n"}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n1 0 0\n"}), + ("held", ("airy 2 8 30 waves", *_RAO_ROWS), {"rao.txt": "HEADING 30\n-1 " + "0 " * 12 + "\n"}), + ( + "held", + ("airy 2 8 30 waves", *_RAO_ROWS), + {"rao.txt": "HEADING 30\n1 -1 " + "0 " * 11 + "\n"}, + ), + ( + "held", + ("airy 2 8 30 waves", *_RAO_ROWS), + {"rao.txt": "HEADING 30\n1 " + "0 " * 12 + "\nHEADING 40\n2 " + "0 " * 12 + "\n"}, + ), + ("held", ("jonswap 2 9 3.3 30 waves", *_RAO_ROWS), {"rao.txt": _rao(heads=(30.0,))}), + ( + "held", + ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(30.0,))}, + ), + ( + "held", + ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-60.0, 120.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", "5 WaveDirections", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-42.0, 102.0))}, + ), + ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(), _FREE: _ROWS}), + ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 100.0)), _FREE: _ROWS}), + ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(28.0, 29.0)), _FREE: _ONE_BETA}), + ( + "held", + ("wk.dat WaterKin", *_RAO_ROWS), + {"rao.txt": _rao(heads=(37.0,)), "wk.dat": _kin("1"), "eta.txt": _ETA}, + ), + ( + "held", + ("wk.dat WaterKin", *_RAO_ROWS), + {"rao.txt": _rao(heads=(38.0,)), "wk.dat": _kin("1"), "eta.txt": _ETA}, + ), + # bathymetry + ("bed", (), {"b.txt": _sea_bed()}), + ("bed", (), {"b.txt": _sea_bed(depth=49.0)}), + ("bed", (), {"b.txt": _sea_bed(depth=50.0004)}), + ("bed", (), {"b.txt": _sea_bed(depth=0.0)}), + ("bed", (), {"b.txt": "0 0 50\n1 0 50\n"}), + ("bed", (), {"b.txt": _sea_bed(tail="0 0 50\n")}), + ("bed", (), {"b.txt": _sea_bed(tail="-1000.0 -1000.0 50\n")}), + ("bed", (), {"b.txt": _sea_bed(tail="1 2\n")}), + ("bed", (), {"b.txt": _sea_bed(tail="1 x 3\n")}), + ("bed", (), {"b.txt": _sea_bed(tail="-1000.0 0.0 50\n1000.0 0.0 50\n")}), + # a crease at the anchor: the rounded and the tabulated surfaces differ there + ( + "bed", + (), + { + "b.txt": "".join( + f"{x} {y} {d}\n" + for x, d in ((-1000.0, 40.0), (400.0, 50.0), (1000.0, 40.0)) + for y in (-1000.0, 1000.0) + ) + }, + ), + # Syrope settings and OWC table + ("syrope", (), _syrope()), + ("syrope", (), _syrope(settings=_SETTINGS + "1.0 k1\n")), + ("syrope", (), _syrope(settings="data/x.dat OWC\n")), + ( + "syrope", + (), + _syrope( + settings=_SETTINGS.replace("LINEAR", "CUBIC") + .replace("QUADRATIC", "CUBIC") + .replace("EXP", "CUBIC") + ), + ), + ("syrope", (), _syrope(owc="strain tension\n0.0 0.0\n")), + ("syrope", (), _syrope(owc="0.0 0.0\n0.01 x\n")), + ("syrope", (), _syrope(owc="0.02 1e6\n0.01 2e6\n0.03 3e6\n")), + ("syrope", (), _syrope(owc="0.0 2e6\n0.01 1e6\n0.03 3e6\n")), + # the routes of the objects: which of them read a motion file + ("rigid6_buoy.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), + ("rod_moored_spar.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), + ("spar_pinned_rods.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), + ("mixed_body_rods_points.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), + ("objects", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), +] +_ALL_FILE_CASES = [("held", rows, files) for rows, files in _WAVE_RECORD_CASES] + _ANCILLARY_CASES + + +def _deck_text(deck: str, rows: tuple[str, ...]) -> str: + if deck == "bed": + return _options(_HELD, "b.txt bathymetryFile", *rows).replace("50.0 WtrDpth\n", "") + if deck == "syrope": + return (EXAMPLES / "syrope_polyester_mooring.dat").read_text(encoding="utf-8") + if deck == "objects": + return _OBJECTS.replace("1.0 TMax\n", "1.0 TMax\n" + "".join(f"{r}\n" for r in rows)) + if deck.endswith(".dat"): + text = (EXAMPLES / deck).read_text(encoding="utf-8") + return re.sub(r"(?m)^(\S+\s+TMax\b.*)$", "\\1\n" + "\n".join(rows), text, count=1) + return _options(_HELD, *rows) + + +def _stage(folder: Path, rows: tuple[str, ...], files: dict[str, str], deck: str = "held") -> Path: folder.mkdir(parents=True, exist_ok=True) for name, text in files.items(): + (folder / name).parent.mkdir(parents=True, exist_ok=True) (folder / name).write_text(text, encoding="utf-8") - deck = folder / "deck.dat" - deck.write_text(_options(_HELD, *rows), encoding="utf-8") - return deck + path = folder / "deck.dat" + path.write_text(_deck_text(deck, rows), encoding="utf-8") + return path def _accepted(deck: Path) -> bool: @@ -173,9 +480,9 @@ def _accepted(deck: Path) -> bool: not os.environ.get("CABLEDYN_DRIVER", "").strip(), reason="needs the native driver (CABLEDYN_DRIVER)", ) -@pytest.mark.parametrize(("rows", "files"), _WAVE_RECORD_CASES) -def test_wave_record_rules_match_the_driver(tmp_path, rows, files): - deck = _stage(tmp_path, rows, files) +@pytest.mark.parametrize(("base", "rows", "files"), _ALL_FILE_CASES) +def test_wave_record_rules_match_the_driver(tmp_path, base, rows, files): + deck = _stage(tmp_path, rows, files, base) run = subprocess.run( [os.environ["CABLEDYN_DRIVER"], str(deck), str(tmp_path / "out")], capture_output=True, @@ -799,7 +1106,7 @@ def test_generated_cases_carry_the_fixed_name_kinematics_files(tmp_path): ("mode", "name"), [("2", "current_profile_dynamic.txt"), ("5", "current_profile_4d.txt")] ) def test_generated_cases_carry_the_moordyn_c_current_tables(tmp_path, mode, name): - table = "a MoorDyn-C current table\n" + table = _CURRENT_TABLES[mode] base = _kinematic_base(tmp_path / "src", f"{mode} Currents", files={name: table}) out = tmp_path / "cases" generate_deck_cases(base, out, {"a": {"option.dtM": 0.05}}) From 708caab46a13fea81927c89ad29c655a76ec7f14 Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 8 Nov 2026 01:01:02 +0900 Subject: [PATCH 4/5] fix(python): honour wave kind spellings and the static check's clock The heading of a wave row missed the ochi_hubble spelling, so it read lambda1 as the direction and a one-heading RAO table at the sea's heading was refused. All the spellings the native reader takes now select the same column. A standalone static check, TMax = 0 without dtM, was refused with a WaterKin WaveKinMod 1 file, because the file rules took dtM from the explicit options. They now run on the nominal clock the native parser installs, and motion and vessel records follow its start-row rule: the grid is the t = 0 row, and a later row is not used, whatever its time step. The parity test now also takes every wave kind and keyword spelling on both clocks against a one-heading RAO table, every shipped deck as a static check without dtM, and every shipped wave deck under each spelling of the options these rules read: 507 decks through DeckFile and the driver, with the same verdict on each. Signed-off-by: Jae Hoon Seo --- python/cabledyn/deck_file.py | 41 ++++-- python/tests/test_deck_native_parity.py | 174 +++++++++++++++++++++++- 2 files changed, 203 insertions(+), 12 deletions(-) diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index 06cd88fd..8f2be9ee 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -1469,6 +1469,7 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None: _BATHY_CREASE_FRACTION = 0.05 _BATHY_CREASE_MIN_HEIGHT = 1.0e-6 _SEA_DEFAULT_DIRECTIONS = 9 +_STATIC_CHECK_DT = 1.0 # the nominal dtM of a standalone TMax = 0 without one (parse_deck) class _AncillaryError(ValueError): @@ -1823,8 +1824,16 @@ def _check_grid_current(records: list[str], name: str) -> None: _check_current_resample(axes[3], spatial, field, [*axes[0], *axes[1], *axes[2]]) -def _motion_time_index(time: float, dt: float, n_time: int, tolerance: float) -> int | None: - """Mirror ``motion_grid_index``: the 0-based step, -1 past TMax, ``None`` off the grid.""" +def _motion_time_index( + time: float, dt: float, n_time: int, tolerance: float, start_row_only: bool +) -> int | None: + """Mirror ``motion_grid_index``: the 0-based step, -1 past TMax, ``None`` off the grid. + + start_row_only: a static check without dtM, whose grid is the t = 0 row: a later row is not + used, whatever its time step. + """ + if start_row_only and time > tolerance: + return -1 if not abs(time / dt) < _MAX_GRID_STEPS: return None step = round(time / dt) @@ -1844,6 +1853,7 @@ def _check_motion_rows( eligible: set[int] | None, roll_ends: set[int] | None, turbines: dict[int, list[int]] | None, + start_row_only: bool, ) -> None: """Mirror ``read_motion_file``. @@ -1874,7 +1884,7 @@ def _check_motion_rows( raise _AncillaryError(f"{where}: J must be the integer id of a TURBINES row") if abs(math.sqrt(sum(value * value for value in row[5:9])) - 1.0) > 1.0e-6: raise _AncillaryError(f"{where}: the quaternion must have unit norm (within 1e-6)") - step = _motion_time_index(row[0], dt, n_time, tolerance) + step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only) if step is None: raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") if step == -1: @@ -1896,7 +1906,7 @@ def _check_motion_rows( point = int(row[1]) if abs(row[1]) <= 2**31 - 1 and row[1] == int(row[1]) else 0 if point < 1: raise _AncillaryError(f"{where}: point_id must be a positive integer") - step = _motion_time_index(row[0], dt, n_time, tolerance) + step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only) if step is None: raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") if eligible is not None and point not in eligible: @@ -1936,7 +1946,7 @@ def _check_motion_rows( ) -def _check_vessel_rows(records: list[str], dt: float, tmax: float) -> None: +def _check_vessel_rows(records: list[str], dt: float, tmax: float, start_row_only: bool) -> None: """Mirror ``read_vessel_record``.""" n_time = round(tmax / dt) + 1 tolerance = 100.0 * _NATIVE_EPSILON * max(1.0, tmax) @@ -1952,7 +1962,7 @@ def _check_vessel_rows(records: list[str], dt: float, tmax: float) -> None: raise _AncillaryError(f"{where}: rows must hold 19 or 20 values, one form for the file") if any(abs(value) > _NATIVE_MAX_COORDINATE for value in row[1:]): raise _AncillaryError(f"{where}: value outside the admissible range") - step = _motion_time_index(row[0], dt, n_time, tolerance) + step = _motion_time_index(row[0], dt, n_time, tolerance, start_row_only) if step is None: raise _AncillaryError(f"{where}: time is outside the dtM/TMax grid") if step == -1: @@ -2301,7 +2311,7 @@ def _train_direction(values: tuple[str, ...]) -> tuple[float, float]: kind = values[0].lower() if kind == "airy": return _native_float(values[3]), 0.0 - index = {"jonswap": 4}.get(kind, 7 if kind.replace("-", "") in {"ochihubble", "ochi"} else 3) + index = 4 if kind == "jonswap" else 7 if kind.startswith("ochi") else 3 spreading = _native_float(values[index + 1]) if len(values) > index + 1 else 0.0 return _native_float(values[index]), spreading @@ -5239,6 +5249,7 @@ def is_body_kind(kind: str, pinned: bool) -> bool: ), bool(finite_length_rods), has_rigid6, + static_check, ) except _AncillaryError as exc: raise DeckFormatError(f"{self._label}: {exc}") from None @@ -5255,14 +5266,18 @@ def _check_ancillary_files( roll_ends: set[str] | None, has_rods: bool, has_rigid6: bool, + static_check: bool, ) -> None: - """Run the native rules of the ancillary files the deck names, those that are there.""" + """Run the native rules of the ancillary files the deck names, those that are there. + + A static check runs on the nominal clock the native parser installs before it validates. + """ def number(key: str) -> float | None: row = option_by_key.get(key) return None if row is None else _native_float(row.values[0]) - dt = number("dtm") + dt = _STATIC_CHECK_DT if static_check else number("dtm") tmax = number("tmax") or 0.0 depth = number("water_depth") record_dir: float | None = None @@ -5323,7 +5338,7 @@ def number(key: str) -> float | None: "vesselMotion/vesselRAO move Coupled/Vessel line points; decks with prescribed " "rods or Rigid6 bodies take motionFile rows" ) - _check_vessel_rows(records, dt, tmax) + _check_vessel_rows(records, dt, tmax, static_check) elif motion_source == "vesselrao": if has_rods or has_rigid6: raise _AncillaryError( @@ -5337,7 +5352,9 @@ def number(key: str) -> float | None: headings = self._sea_headings(option_by_key, options) _check_rao_table(records, headings) else: - self._check_motion_file(records, dt, tmax, roll_ends, has_rods, has_rigid6) + self._check_motion_file( + records, dt, tmax, roll_ends, has_rods, has_rigid6, static_check + ) def _check_waterkin_file( self, @@ -5508,6 +5525,7 @@ def _check_motion_file( roll_ends: set[str] | None, has_rods: bool, has_rigid6: bool, + start_row_only: bool, ) -> None: """Mirror ``read_motion_file`` on the route the native dispatch takes. @@ -5548,6 +5566,7 @@ def _check_motion_file( if roll_ends is None or has_rods or has_rigid6 else {int(point) for point in roll_ends}, turbines, + start_row_only, ) def _ancillary_path(self, value: str) -> Path: diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index 7f48943f..d6cfd777 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -12,6 +12,7 @@ import json import os import re +import shutil import subprocess from pathlib import Path @@ -442,10 +443,158 @@ def _syrope(settings=_SETTINGS, owc=_OWC) -> dict[str, str]: ("mixed_body_rods_points.dat", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), ("objects", ("m.txt motionFile",), {"m.txt": "not a motion row\n"}), ] -_ALL_FILE_CASES = [("held", rows, files) for rows, files in _WAVE_RECORD_CASES] + _ANCILLARY_CASES +# every spelling the reader takes of the options and kinds these rules read +_KEYWORD_SPELLINGS = ( + ("waves", "wave"), + ("wavetrain", "wave_train"), + ("WaveSpreading", "wave_spreading"), + ("WaveDirections", "wave_directions"), + ("WaveComponents", "wave_components"), + ("WtrDpth", "water_depth"), + ("dtM", "dt"), + ("WaterKin", "WaveKin"), + ("bathymetryFile", "seafloor_file"), + ("rampTime", "tramp"), +) +_KIND_VALUES = { + "airy": ("2 8 30", "2 8 30"), + "Airy": ("2 8 30", "2 8 30"), + "jonswap": ("2 9 3.3 30", "2 9 3.3 30 0"), + "JONSWAP": ("2 9 3.3 30", "2 9 3.3 30 0"), + "pm": ("2 9 30", "2 9 30 0"), + "issc": ("2 9 30", "2 9 30 0"), + "bretschneider": ("2 9 30", "2 9 30 0"), + "torsethaugen": ("6 12 30", "6 12 30 0"), + "ochihubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"), + "ochi-hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"), + "ochi_hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"), + "Ochi_Hubble": ("2 8 1 1 12 1 30", "2 8 1 1 12 1 30 0"), +} +# every kind and keyword against a one-heading RAO table: at the sea's heading and off it, on +# the dynamic clock and on the static check's +_SPELLED_CASES = [ + (base, (f"{kind} {values[train]} {keyword}", *_RAO_ROWS), {"rao.txt": _rao(heads=(heading,))}) + for kind, values in _KIND_VALUES.items() + for train, keyword in ((0, "waves"), (0, "wave"), (1, "wavetrain"), (1, "wave_train")) + for heading in (30.0, 31.0) + for base in ("held", "held_static") +] +_SPELLED_CASES += [ + # the kinds whose rows spread, and the spellings of the options that spread them + (base, (f"pm 2 9 30 {keyword}", f"4 {spread}", *_RAO_ROWS), {"rao.txt": _rao(heads=heads)}) + for base in ("held", "held_static") + for keyword in ("waves", "wave") + for spread in ("WaveSpreading", "wave_spreading", "WAVESPREADING") + for heads in ((30.0,), (-180.0, 179.0)) +] + [ + (base, (f"ochi_hubble 2 8 1 1 12 1 30 4 {keyword}", *_RAO_ROWS), {"rao.txt": _rao(heads=heads)}) + for base in ("held", "held_static") + for keyword in ("wavetrain", "wave_train") + for heads in ((30.0,), (-180.0, 179.0)) +] +# a record or a kinematics file on the static check's clock, under each spelling of its selector +_STATIC_CASES = [ + ("held_static", (f"{value} {keyword}", *rows), files) + for keyword in ("WaterKin", "WaveKin", "WATERKIN") + for value, rows, files in ( + ("wk.dat", (), {"wk.dat": _kin("1"), "eta.txt": _ETA}), + ("wk.dat", ("12.5 WaveRefX",), {"wk.dat": _kin("1"), "eta.txt": _ETA}), + ("wk.dat", ("12.5 WaveRefX",), {"wk.dat": _kin("0"), "eta.txt": _ETA}), + ("wk.dat", (), {"wk.dat": _kin("1"), "eta.txt": "0 0\n1 0\n"}), + ("wk.dat", (), {"wk.dat": _kin("1", dt="100"), "eta.txt": _ETA}), + ( + "wk.dat", + _RAO_ROWS, + {"wk.dat": _kin("1"), "eta.txt": _ETA, "rao.txt": _rao(heads=(37.0,))}, + ), + ( + "wk.dat", + _RAO_ROWS, + {"wk.dat": _kin("1"), "eta.txt": _ETA, "rao.txt": _rao(heads=(38.0,))}, + ), + ("7", (), {_FREE: _ROWS}), + ("7", ("12.5 WaveRefY",), {_FREE: _ROWS}), + ("7", (), {_FREE: "0 0 0 0\n0.5 0.5 0 6.3\n"}), + ("3", ("0.1 dtWave",), {"wave_elevation.txt": _ETA.replace("time elevation\n", "")}), + ) +] + [ + ("held_static", ("m.txt motionFile",), {"m.txt": "0.0 2 0 0 0 0 0 0 0 0 0\n"}), + ( + "held_static", + ("m.txt motionFile",), + {"m.txt": "0.0 2 0 0 0 0 0 0 0 0 0\n0.05 2 0 0 0 0 0 0 0 0 0\n"}, + ), + ("held_static", ("m.txt motionFile",), {"m.txt": _hold()}), + ("held_static", ("m.txt motionFile",), {"m.txt": "1.0 2 0 0 0 0 0 0 0 0 0\n"}), + ("held_static", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel(steps=1)}), + ("held_static", ("v.txt vesselMotion", "0|0|0 vesselRef"), {"v.txt": _vessel()}), + ("bed_static", (), {"b.txt": _sea_bed()}), + ("held", ("airy 2 8 0 waves", "4 WaveSpreading"), {}), +] + + +def _shipped_decks() -> list[str]: + decks = [path.name for path in sorted(EXAMPLES.glob("*.dat"))] + decks.remove("cabledyn_options_reference.dat") + return [*decks, "moordynC_wavekin/chain_wavekin7_currents1.dat"] + + +def _carries_waves(name: str) -> bool: + text = (EXAMPLES / name).read_text(encoding="utf-8") + return re.search(r"(?mi)^(?!none\b|--)\S.*\s(?:waves?|wave_?train)\b", text) is not None or ( + re.search(r"(?mi)^(?!0\s|none\s)\S+\s+(?:WaterKin|WaveKin)\b", text) is not None + ) + + +# every shipped deck as a static check without dtM; every shipped wave deck under each spelling +_SHIPPED_CASES = [(f"shipped:static:{name}", (), {}) for name in _shipped_decks()] + [ + (f"shipped:{variant}:{name}", (), {}) + for name in _shipped_decks() + if _carries_waves(name) + for variant in ("spelled", "upper") +] +_ALL_FILE_CASES = [ + *[("held", rows, files) for rows, files in _WAVE_RECORD_CASES], + *_ANCILLARY_CASES, + *_SPELLED_CASES, + *_STATIC_CASES, + *_SHIPPED_CASES, +] + + +def _respell(text: str, old: str, new: str) -> str: + """``text`` with the first token ``old`` of each row (any case) spelled ``new``.""" + token = re.compile(rf"(?i)(? str: + """``text`` as a static check: TMax 0, and without dtM unless it is kept.""" + text = re.sub(r"(?mi)^\s*\S+(\s+TMax\b)", r"0.0\1", text, count=1) + return text if keep_dtm else re.sub(r"(?mi)^\s*\S+\s+dtM\b.*\n", "", text, count=1) + + +def _shipped_text(variant: str, name: str) -> str: + text = (EXAMPLES / name).read_text(encoding="utf-8") + if variant == "static": + return _static_check(text) + text = _static_check(text, keep_dtm=True) + for old, new in _KEYWORD_SPELLINGS: + text = _respell(text, old, new if variant == "spelled" else old.upper()) + if variant == "upper": + for kind in ("airy", "jonswap", "pm", "issc", "torsethaugen", "ochihubble", "stream"): + text = _respell(text, kind, kind.upper()) + return text def _deck_text(deck: str, rows: tuple[str, ...]) -> str: + if deck.startswith("shipped:"): + return _shipped_text(*deck.split(":")[1:]) + if deck.endswith("_static"): + return _static_check(_deck_text(deck[: -len("_static")], rows)) if deck == "bed": return _options(_HELD, "b.txt bathymetryFile", *rows).replace("50.0 WtrDpth\n", "") if deck == "syrope": @@ -463,6 +612,13 @@ def _stage(folder: Path, rows: tuple[str, ...], files: dict[str, str], deck: str for name, text in files.items(): (folder / name).parent.mkdir(parents=True, exist_ok=True) (folder / name).write_text(text, encoding="utf-8") + if deck.startswith("shipped:"): + source = (EXAMPLES / deck.split(":")[2]).parent + for item in source.iterdir(): + if item.is_dir(): + shutil.copytree(item, folder / item.name, dirs_exist_ok=True) + elif item.suffix != ".dat": + shutil.copy2(item, folder / item.name) path = folder / "deck.dat" path.write_text(_deck_text(deck, rows), encoding="utf-8") return path @@ -526,6 +682,22 @@ def test_wave_files_next_to_the_deck_are_read(tmp_path, rows, files, message): DeckFile.read(_stage(tmp_path, rows, files)) +@pytest.mark.parametrize("kind", ["ochihubble", "ochi-hubble", "ochi_hubble", "OCHI_HUBBLE"]) +def test_every_spelling_of_a_kind_takes_its_heading_from_its_own_column(tmp_path, kind): + rows = (f"{kind} 2 8 1 1 12 1 30 waves", *_RAO_ROWS) + DeckFile.read(_stage(tmp_path, rows, {"rao.txt": _rao(heads=(30.0,))})) + with pytest.raises(DeckFormatError, match="one-heading RAO table"): + DeckFile.read(_stage(tmp_path, rows, {"rao.txt": _rao(heads=(31.0,))})) + + +def test_a_static_check_reads_a_record_on_its_own_clock(tmp_path): + # TMax = 0 without dtM: the record resamples on its dtWave, whatever the nominal step + files = {"wk.dat": _kin("1"), "eta.txt": _ETA} + DeckFile.read(_stage(tmp_path, _WK, files, "held_static")) + with pytest.raises(DeckFormatError, match="four `time elevation` rows"): + DeckFile.read(_stage(tmp_path, _WK, {"eta.txt": "0 0\n1 0\n"}, "held_static")) + + def test_wave_files_that_are_not_there_are_left_to_the_solver(tmp_path): # In memory under the nominal path, or on disk ahead of its files: the deck passes here, # and the driver, which must open the file, refuses it by name while it is missing. From e913cf68db7c61e5c9ff93d21c9c93be790cb31f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Sun, 8 Nov 2026 02:34:18 +0900 Subject: [PATCH 5/5] fix(waves): judge RAO headings by the bins that carry energy A narrow spreading leaves the outer direction bins without energy, and the RAO synthesis skips their components, so DeckFile refused tables the solver accepts. It now mirrors the bin weights. A weight below the smallest normal number is exactly zero natively too, so which bins carry the sea depends on the spreading and the bin count alone, not on where the product with the spectrum underflows. Signed-off-by: Jae Hoon Seo --- python/cabledyn/deck_file.py | 33 +++++++++++++++++++++++-- python/tests/test_deck_native_parity.py | 31 +++++++++++++++++++++++ src/CableDyn_WaveSpectra.f90 | 5 +++- 3 files changed, 66 insertions(+), 3 deletions(-) diff --git a/python/cabledyn/deck_file.py b/python/cabledyn/deck_file.py index 8f2be9ee..5ba451b2 100644 --- a/python/cabledyn/deck_file.py +++ b/python/cabledyn/deck_file.py @@ -1469,6 +1469,7 @@ def _scalar_option_value_error(keyword: str, value: str) -> str | None: _BATHY_CREASE_FRACTION = 0.05 _BATHY_CREASE_MIN_HEIGHT = 1.0e-6 _SEA_DEFAULT_DIRECTIONS = 9 +_SPREADING_QUADRATURE = 128 _STATIC_CHECK_DT = 1.0 # the nominal dtM of a standalone TMax = 0 without one (parse_deck) @@ -2316,11 +2317,39 @@ def _train_direction(values: tuple[str, ...]) -> tuple[float, float]: return _native_float(values[index]), spreading +def _spreading_weights(spreading: float, bins: int) -> list[float]: + """Mirror ``CD_Spreading_Weights``: the cos-2s weight of each equal-angle bin.""" + constant = math.exp(math.lgamma(spreading + 1.0) - math.lgamma(spreading + 0.5)) + constant /= math.sqrt(math.pi) + weights = [] + for j in range(bins): + start = -0.5 * math.pi + math.pi * j / bins + step = (-0.5 * math.pi + math.pi * (j + 1) / bins - start) / _SPREADING_QUADRATURE + weight = 0.0 + for q in range(_SPREADING_QUADRATURE + 1): + density = constant * max(0.0, math.cos(start + step * q)) ** (2.0 * spreading) + end = q in (0, _SPREADING_QUADRATURE) + weight += density * (1.0 if end else 4.0 if q % 2 == 1 else 2.0) + weights.append(weight * step / 3.0) + total = sum(weights) + if total > 0.0: + weights = [weight / total for weight in weights] + return [weight if weight >= sys.float_info.min else 0.0 for weight in weights] + + def _spread_headings(direction: float, spreading: float, bins: int) -> list[float]: - """The headings of a train (``CD_Sea_Synthesise``): the bin centres of a spread one.""" + """The headings of a train (``CD_Sea_Synthesise``) that carry energy. + + A spread train has one heading per bin centre. A narrow spreading leaves the outer bins with + no weight, and ``synthesize_vessel_rao`` skips their zero-amplitude components. + """ if spreading <= 0.0: return [direction] - return [direction - 90.0 + 180.0 * (j + 0.5) / bins for j in range(bins)] + return [ + direction - 90.0 + 180.0 * (j + 0.5) / bins + for j, weight in enumerate(_spreading_weights(spreading, bins)) + if weight > 0.0 + ] @dataclass(frozen=True) diff --git a/python/tests/test_deck_native_parity.py b/python/tests/test_deck_native_parity.py index d6cfd777..f968114a 100644 --- a/python/tests/test_deck_native_parity.py +++ b/python/tests/test_deck_native_parity.py @@ -383,6 +383,37 @@ def _syrope(settings=_SETTINGS, owc=_OWC) -> dict[str, str]: ("jonswap 2 9 3.3 30 waves", "4 WaveSpreading", "5 WaveDirections", *_RAO_ROWS), {"rao.txt": _rao(heads=(-42.0, 102.0))}, ), + # a narrow spreading underflows the outer bins: only the bins with energy need a heading + ( + "held", + ("jonswap 2 9 3.3 0 waves", "1000 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-40.0, 40.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 0 waves", "1000 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-20.0, 20.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 0 waves", "838 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-40.0, 40.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 0 waves", "700 WaveSpreading", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-40.0, 40.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 0 1000 wavetrain", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-40.0, 40.0))}, + ), + ( + "held", + ("jonswap 2 9 3.3 0 1000 wavetrain", "5 WaveDirections", *_RAO_ROWS), + {"rao.txt": _rao(heads=(-36.0, 36.0))}, + ), ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(), _FREE: _ROWS}), ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(0.0, 100.0)), _FREE: _ROWS}), ("held", ("7 WaterKin", *_RAO_ROWS), {"rao.txt": _rao(heads=(28.0, 29.0)), _FREE: _ONE_BETA}), diff --git a/src/CableDyn_WaveSpectra.f90 b/src/CableDyn_WaveSpectra.f90 index f2b45c08..5ff84613 100644 --- a/src/CableDyn_WaveSpectra.f90 +++ b/src/CableDyn_WaveSpectra.f90 @@ -233,7 +233,9 @@ END FUNCTION CD_Spreading_Constant PURE SUBROUTINE CD_Spreading_Weights(s, n, offset, weight) !! Discretise D(theta) = K(s) cos^2s(theta) on [-pi/2, pi/2] into n equal-angle bins: !! offset(j) [rad] is the bin centre and weight(j) the bin integral of D (composite - !! Simpson quadrature, 128 intervals per bin), normalised to sum 1. + !! Simpson quadrature, 128 intervals per bin), normalised to sum 1. A bin below the smallest + !! normal number holds no energy: its weight is exactly zero, so which bins carry the sea + !! depends on s and n alone, not on where the spectrum's product with it underflows. REAL(wp), INTENT(IN) :: s INTEGER, INTENT(IN) :: n REAL(wp), INTENT(OUT) :: offset(n), weight(n) @@ -263,6 +265,7 @@ PURE SUBROUTINE CD_Spreading_Weights(s, n, offset, weight) END DO total = SUM(weight) IF (total > CD_ZERO) weight = weight/total + WHERE (weight < TINY(CD_ONE)) weight = CD_ZERO END SUBROUTINE CD_Spreading_Weights SUBROUTINE CD_Sea_Reset(sea)