fix(waves): a MoorDyn-C mean level is added to the elevation alone

CHANGELOG.md

@@ -96,6 +96,13 @@ All notable changes to CableDyn are recorded in this file. The format follows
 - Every initialisation note of a coupled line reaches the console; long lists were cut or dropped.
 - The OpenFAST module driver runs a rotating motion in a single-precision build.
 - A range-graph file that cannot be written completely is removed instead of left behind.
+- A mean level in a MoorDyn-C kinematics file, a non-zero `omega = 0` row of `WaveKin 7` or the
+  mean of a `WaveKin 3` record, is added to the elevation alone, as MoorDyn-C does: the wave
+  numbers stay those of `WtrDpth` and the stretching maps the water under the total elevation.
+  CableDyn took the wave numbers at `WtrDpth` plus the level. With a 0.8 m level in 50 m of water
+  and six components in six directions, `WaveKin 7` now matches MoorDyn-C v2.7.2 to 1e-8 instead
+  of 7.7e-3 m in elevation and 1 % in velocity. The mean of a WaterKin `WaveKinMod 1` record still
+  raises the still-water level; MoorDyn-F drops it.
 - The C API messages of `CableDyn_UpdatePointFluidFields` and `CableDyn_CalcOutputDerivatives` on
   the coupled aggregate name its real scope: decks with finite-EI lines, bodies or rods.
 
@@ -121,6 +128,10 @@ All notable changes to CableDyn are recorded in this file. The format follows
   `EI × Curv` only to within the strain. Tension, positions and forces are unchanged, and so are
   `OutputRecovery element` values on elements without a stiffener sleeve.
 - Azimuth values of the `EI = 0` lines named above read 360 degrees lower on those routes.
+- On a `WaveKin 3` or `WaveKin 7` deck whose mean level is not zero, the elevation away from the
+  origin, the fluid velocity and acceleration, and so the loads and every output, change by the
+  order of the level over the depth: up to 1 % in velocity at 0.8 m in 50 m. A deck with no mean
+  level is unchanged bit for bit.
 - OpenFAST checkpoints change format: they carry the range graphs and the input of the last
   commit, so a checkpoint written by an earlier CableDyn cannot be restarted with this one.
 

CMakeLists.txt

@@ -2321,6 +2321,12 @@ file(GLOB CD_WAVEKIN_MDC_FILES RELATIVE ${CMAKE_SOURCE_DIR}/tests/data/wavekin_m
 foreach(f IN LISTS CD_WAVEKIN_MDC_FILES)
   configure_file(tests/data/wavekin_mdc/${f} ${CMAKE_BINARY_DIR}/wavekin_mdc/${f} COPYONLY)
 endforeach()
+# the same with a mean level (an omega = 0 row, a raised record), against MoorDyn-C v2.7.2
+file(GLOB CD_WAVEKIN_MDC_MEAN_FILES RELATIVE ${CMAKE_SOURCE_DIR}/tests/data/wavekin_mdc_mean
+     ${CMAKE_SOURCE_DIR}/tests/data/wavekin_mdc_mean/*)
+foreach(f IN LISTS CD_WAVEKIN_MDC_MEAN_FILES)
+  configure_file(tests/data/wavekin_mdc_mean/${f} ${CMAKE_BINARY_DIR}/wavekin_mdc_mean/${f} COPYONLY)
+endforeach()
 add_executable(test_wavekin_modes tests/test_wavekin_modes.f90)
 target_link_libraries(test_wavekin_modes PRIVATE cabledyn_core)
 add_test(NAME wavekin_modes COMMAND test_wavekin_modes WORKING_DIRECTORY ${CMAKE_BINARY_DIR})

doc/driver_format.md

@@ -1385,10 +1385,13 @@ 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.
+A mean level, a non-zero `omega = 0` row or the mean of `wave_elevation.txt`, is added to the
+elevation alone, as MoorDyn-C does: every wave number stays that of `WtrDpth`, and the Wheeler
+stretching maps the water column under the total elevation, the level included. The
+`wavekin_modes` test holds both files with a 0.8 m level, under six `WaveKin 7` components in
+six directions, to MoorDyn-C v2.7.2: `WaveKin 7` to `1e-8` of its scale and `WaveKin 3` to the
+tolerances above. The mean of a WaterKin `WaveKinMod 1` record raises the still-water level
+instead ([below](#moordyn-f-waterkin-file-modes)).
 
 (waterkin-file-modes)=
 ### MoorDyn-F `WaterKin` file modes
@@ -1406,6 +1409,11 @@ A MoorDyn-F deck may select a WaterKin file (grammar in [Auxiliary file formats]
 | `CurrentMod = 1` | N-level WaterKin depth table, interpolated at structural nodes |
 | `CurrentMod = 2` | current from the coupled host SeaState field |
 
+The mean of a `WaveKinMod 1` record raises the still-water level: its wave numbers are those of
+`WtrDpth` plus the mean, and the stretching maps the water column under the raised level.
+MoorDyn-F instead drops the record's mean, and its Nyquist bin, so a record whose mean is not zero
+puts the surface of the two codes apart by that mean.
+
 **`WaveKinMod 1`** is standalone-only. It needs a `WaveKinFile` of at least four finite
 `time elevation` rows starting at `t = 0` with strictly increasing times, a positive `dtWave`, a
 positive `TMax`, and the ordinary `waves` requirements (`dtM`/`TMax` and `WtrDpth`). It cannot be

doc/theory.rst

@@ -1185,6 +1185,13 @@ has one direction :math:`\beta_i = \beta`. The rows of a ``WaveKin 7`` table may
 their own: a short-crested sea, whose elevation and kinematics are the sums of those of its
 directions, stretched by the one total elevation.
 
+The mean :math:`\bar\eta` of an elevation history, or the ``omega = 0`` row of a ``WaveKin 7``
+table, is a steady level. A WaterKin history raises the still-water level by it: the wave numbers
+are those of the depth :math:`h + \bar\eta`, and the stretching maps the water column under the
+raised level. A MoorDyn-C file (``WaveKin 3`` or ``7``) adds it to the elevation alone, as
+MoorDyn-C does: the wave numbers stay those of :math:`h`, and the stretching maps
+:math:`[-h,\ \bar\eta + \eta]` onto :math:`[-h,\ 0]`.
+
 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

src/CableDyn_DeckDriver.f90

@@ -675,7 +675,8 @@ MODULE CableDyn_DeckDriver
     ! irregular seas; wave_eta_mean carries the zero-frequency waterline offset.
     REAL(wp), ALLOCATABLE :: wave_component_omega(:), wave_component_k(:)
     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)
+    ! exp(-2 k d) at the depth d of the wave numbers (wtrdpth, plus wave_eta_mean for a WaterKin
+    ! record: component_level_shift) 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
@@ -14305,6 +14306,8 @@ SUBROUTINE read_moordyn_c_kinematics(dir, opts, ErrStat, ErrMsg)
     !!    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.
+    !!  The mean level of either file is added to the elevation alone, which the stretching then
+    !!  maps, and the wave numbers stay those of WtrDpth, as in MoorDyn-C (component_level_shift).
     !!  - 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
@@ -32751,7 +32754,7 @@ SUBROUTINE deck_wave_kinematics(opts, x, y, z, time, eta, velocity, acceleration
     INTEGER, INTENT(OUT) :: ErrStat
     CHARACTER(*), INTENT(OUT) :: ErrMsg
     REAL(wp), INTENT(OUT), OPTIONAL :: pdyn
-    REAL(wp) :: ramp, ramp_rate
+    REAL(wp) :: ramp, ramp_rate, shift
 
     IF (opts%has_wave .AND. .NOT. opts%wave_cache_ready) THEN
       CALL fail(ErrStat, ErrMsg, 'deck wave kinematics called before wave precompute')
@@ -32787,26 +32790,29 @@ SUBROUTINE deck_wave_kinematics(opts, x, y, z, time, eta, velocity, acceleration
       CALL CD_Stream_Kinematics(opts%stream, x, y, z, time, ramp, eta, velocity, acceleration, ErrStat, ErrMsg, &
                                 pdyn)
     CASE (WAVE_COMPONENTS)
+      ! the mean of a WaterKin record raises the still-water level (shift); that of a MoorDyn-C
+      ! file is added to the surface alone (eta_offset), the wave numbers staying at WtrDpth
+      shift = component_level_shift(opts)
       IF (ramp < CD_ONE) THEN
-        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., &
+        CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - shift, time, &
+                                          opts%wtrdpth + shift, 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, comp_cos=opts%wave_component_cos, &
-                                          comp_sin=opts%wave_component_sin)
+                                          comp_sin=opts%wave_component_sin, eta_offset=opts%wave_eta_mean - shift)
       ELSE
-        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., &
+        CALL CD_Component_Wave_Kinematics(x - opts%wave_ref_x, y - opts%wave_ref_y, z - shift, time, &
+                                          opts%wtrdpth + shift, 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, comp_cos=opts%wave_component_cos, &
-                                          comp_sin=opts%wave_component_sin)
+                                          comp_sin=opts%wave_component_sin, eta_offset=opts%wave_eta_mean - shift)
       END IF
-      eta = eta + opts%wave_eta_mean
+      eta = eta + shift
     CASE DEFAULT
       eta = CD_ZERO
       velocity = CD_ZERO
@@ -32818,6 +32824,15 @@ SUBROUTINE deck_wave_kinematics(opts, x, y, z, time, eta, velocity, acceleration
       acceleration = acceleration + (ramp_rate/ramp)*velocity
   END SUBROUTINE deck_wave_kinematics
 
+  PURE REAL(wp) FUNCTION component_level_shift(opts) RESULT(shift)
+    !! The rise of the still-water level under a component table: the mean of a WaterKin
+    !! WaveKinMod 1 record, so its wave numbers are those of WtrDpth plus the mean; zero for a
+    !! MoorDyn-C kinematics file (WaveKin 3, 7), whose mean MoorDyn-C adds to the elevation alone.
+    TYPE(DeckOptions), INTENT(IN) :: opts
+    shift = opts%wave_eta_mean
+    IF (opts%mdc_wave_kin > 0) shift = CD_ZERO
+  END FUNCTION component_level_shift
+
   PURE SUBROUTINE deck_wave_ramp(ramp_time, time, ramp, ramp_rate)
     !! Half-cosine start-up ramp r(t) = (1 - cos(pi*t/T))/2 on [0, T], 1 afterwards, and
     !! its rate dr/dt; T <= 0 disables the ramp (r = 1). r and dr/dt are continuous, so
@@ -33193,12 +33208,12 @@ SUBROUTINE precompute_deck_wave(opts, ErrStat, ErrMsg)
       END IF
       DO k = 1, SIZE(opts%wave_component_omega)
         CALL CD_Solve_Dispersion_Wavenumber(opts%wave_component_omega(k), &
-                                            opts%wtrdpth + opts%wave_eta_mean, opts%g, &
+                                            opts%wtrdpth + component_level_shift(opts), opts%g, &
                                             opts%wave_component_k(k), ErrStat, ErrMsg)
         IF (ErrStat /= CD_DECKDRV_OK) EXIT
       END DO
       IF (ErrStat == CD_DECKDRV_OK) THEN
-        opts%wave_component_e2kd = EXP(-2.0_wp*(opts%wave_component_k*(opts%wtrdpth + opts%wave_eta_mean)))
+        opts%wave_component_e2kd = EXP(-2.0_wp*(opts%wave_component_k*(opts%wtrdpth + component_level_shift(opts))))
         CALL set_component_tail(opts%wave_component_omega, opts%wave_component_k, opts%wave_component_amplitude, &
                                 opts%wave_component_tail)
       END IF

src/CableDyn_Hydro.f90

@@ -2239,7 +2239,8 @@ 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, comp_cos, comp_sin)
+                                          amplitude_scale, pdyn, exp_m2kd, tail_sums, comp_cos, comp_sin, &
+                                          eta_offset)
     !! 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 --
@@ -2278,6 +2279,10 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     !! 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.
+    !! eta_offset is a steady level added to the surface elevation the stretching maps and eta
+    !! returns, the components' wave numbers staying those of depth: the zero-frequency row of a
+    !! MoorDyn-C kinematics file, which MoorDyn-C adds to the elevation alone. It is not scaled
+    !! by amplitude_scale.
     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(:)
@@ -2290,8 +2295,9 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
     ! (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(:, :), comp_cos(:), comp_sin(:)
+    REAL(wp), INTENT(IN), OPTIONAL :: eta_offset
     INTEGER :: n
-    REAL(wp) :: cb, sb, scl
+    REAL(wp) :: cb, sb, scl, offset
     LOGICAL :: trust_inputs
 
     eta = CD_ZERO
@@ -2356,6 +2362,12 @@ SUBROUTINE CD_Component_Wave_Kinematics(x, y, z, t, depth, direction_deg, stretc
       CALL fail(ErrStat, ErrMsg, 'point, time, and direction must be finite')
       RETURN
     END IF
+    offset = CD_ZERO
+    IF (PRESENT(eta_offset)) offset = eta_offset
+    IF (.NOT. CD_Is_Finite(offset)) THEN
+      CALL fail(ErrStat, ErrMsg, 'eta_offset must be finite')
+      RETURN
+    END IF
     ! Preserve the regular-wave reduction as an executable identity, independent of
     ! compiler reassociation choices in the general summation loop. Besides making the
     ! documented one-component contract portable, this avoids loop setup in a common
@@ -2364,27 +2376,29 @@ 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 .AND. .NOT. PRESENT(comp_cos)) THEN
+    IF (n == 1 .AND. ABS(phase(1)) <= CD_ZERO .AND. scl*amplitude(1) > CD_ZERO .AND. .NOT. PRESENT(comp_cos) &
+        .AND. ABS(offset) <= CD_ZERO) 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)
       RETURN
     END IF
     cb = COS(direction_deg*DEG2RAD)
     sb = SIN(direction_deg*DEG2RAD)
-    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)
+    CALL component_sums(n, x, y, z, t, depth, stretch, scl, offset, 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, 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)
+  SUBROUTINE component_sums(n, x, y, z, t, depth, stretch, scl, offset, 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. 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.
+    !! offset (eta_offset) is added to the summed elevation before the stretching.
     INTEGER, INTENT(IN) :: n
-    REAL(wp), INTENT(IN) :: x, y, z, t, depth, scl, cb, sb
+    REAL(wp), INTENT(IN) :: x, y, z, t, depth, scl, offset, 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)
@@ -2432,7 +2446,7 @@ SUBROUTINE component_sums(n, x, y, z, t, depth, stretch, scl, cb, sb, omega_comp
         CALL add_elevation_lanes(nk, amplitude(i0:i1), bc, lane_eta)
       END IF
     END DO
-    eta = scl*CD_Wave_Sum_Lanes(lane_eta, lanes)
+    eta = scl*CD_Wave_Sum_Lanes(lane_eta, lanes) + offset
 
     IF (stretch) THEN
       IF (z > eta) RETURN