From 07c0b9e5c959380dd7b4423db3930cc5c463fe4f Mon Sep 17 00:00:00 2001 From: Jae Hoon Seo Date: Wed, 14 Oct 2026 19:10:28 +0900 Subject: [PATCH] perf: evaluate each directional-sea phase and depth profile once Profiling the release build showed the directional-sea kinematics (CD_Sea_Kinematics, the deck's wavetrain rows and spectral seas) taking 80 % of a chain run in a two-train sea. For every node, step and component it computed the phase cosine twice (for the elevation and again for the kinematics) and two exponentials for the depth profile. The cosines of the elevation pass are now kept for the kinematics pass (up to 4096 components; any further ones are recomputed), and the depth profile uses EXP(-k (z + 2h)) = e2kd / EXP(k z) with the table's e2kd = EXP(-2 k h), as the JONSWAP component path does. A component whose EXP(k z) or e2kd underflows, or a call at another depth than the synthesis, takes the second exponential instead; the sea now records the depth of its synthesis for that test. With the release build the two-train chain example runs in 1.57 s instead of 1.89 s (-17 %). The formulas are unchanged. With gfortran the outputs of the two directional-sea examples are bit-identical; with IFX the reuse removes a fused sin/cos evaluation of the second pass, and one channel of chain_torsethaugen_spread.dat differs by one unit in its eighth printed digit at isolated rows. A new unit test checks the kinematics against a direct per-component sum, with all cosines kept, with cosines past the kept length, with underflowing exponentials, and at another depth. Signed-off-by: Jae Hoon Seo --- CHANGELOG.md | 10 +++++ src/CableDyn_WaveSpectra.f90 | 53 +++++++++++++++++----- tests/test_wave_spectra.f90 | 86 +++++++++++++++++++++++++++++++++++- 3 files changed, 137 insertions(+), 12 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index efd8f9d..b67c720 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -8,6 +8,16 @@ All notable changes to CableDyn are recorded in this file. The format follows ## [Unreleased] +### Changed + +- Faster wave kinematics of directional and multi-train seas (`wavetrain` rows, spectral + `waves` rows, `WaveSpreading`). Each component's phase cosine is now computed once for the + surface elevation and the kinematics, and its depth profile with one exponential instead of + two. `examples/chain_two_train_sea.dat` runs 17 % faster with the release build. The formulas + are unchanged. With gfortran the results are bit-identical. With the Intel compiler of the + release, which had computed the second cosine through a fused sine/cosine evaluation, they + can differ in the last printed digit. + ### Fixed - The development environment (`environment.yml`) excludes OpenBLAS 0.3.34, which caused the diff --git a/src/CableDyn_WaveSpectra.f90 b/src/CableDyn_WaveSpectra.f90 index befe8b0..2b617f1 100644 --- a/src/CableDyn_WaveSpectra.f90 +++ b/src/CableDyn_WaveSpectra.f90 @@ -64,6 +64,7 @@ MODULE CableDyn_WaveSpectra ! the synthesised component table (all trains) INTEGER :: n_comp = 0 REAL(wp), ALLOCATABLE :: omega(:), k(:), amplitude(:), phase(:), cosb(:), sinb(:), e2kd(:) + REAL(wp) :: depth = CD_ZERO ! the water depth of the synthesis: k and e2kd = EXP(-2 k depth) LOGICAL :: ready = .FALSE. END TYPE CD_SeaType @@ -569,6 +570,7 @@ SUBROUTINE CD_Sea_Synthesise(sea, depth, gravity, seed, ErrStat, ErrMsg) sea%e2kd(i) = EXP(-2.0_wp*sea%k(i)*depth) END DO sea%n_comp = n_total + sea%depth = depth sea%ready = .TRUE. END SUBROUTINE CD_Sea_Synthesise @@ -626,8 +628,12 @@ SUBROUTINE CD_Sea_Kinematics(sea, x, y, z, t, depth, scale, eta, velocity, accel CHARACTER(*), INTENT(OUT) :: ErrMsg REAL(wp), INTENT(OUT), OPTIONAL :: pdyn INTEGER, PARAMETER :: NB = 256 - REAL(wp) :: cth(NB), sth(NB), ch(NB), sh(NB) - INTEGER :: i0, i1, i, m + ! The phase cosines of the elevation pass, kept for the kinematics pass (whose chunk loop + ! would compute the same values again), for the first NKEEP components; the rest are + ! recomputed there. A multiple of NB, so each chunk is either kept whole or recomputed. + INTEGER, PARAMETER :: NKEEP = 16*NB + REAL(wp) :: cth(NB), sth(NB), ch(NB), sh(NB), ckeep(NKEEP) + INTEGER :: i0, i1, i, m, n_small REAL(wp) :: z_eval, a, e1, e2, den, uh eta = CD_ZERO @@ -648,6 +654,7 @@ SUBROUTINE CD_Sea_Kinematics(sea, x, y, z, t, depth, scale, eta, velocity, accel DO i = i0, i1 eta = eta + (scale*sea%amplitude(i))*cth(i - i0 + 1) END DO + IF (i1 <= NKEEP) ckeep(i0:i1) = cth(1:i1 - i0 + 1) END DO IF (z > eta) RETURN IF (.NOT. (depth + eta > EPSILON(depth)*depth)) THEN @@ -659,20 +666,44 @@ SUBROUTINE CD_Sea_Kinematics(sea, x, y, z, t, depth, scale, eta, velocity, accel DO i0 = 1, sea%n_comp, NB i1 = MIN(sea%n_comp, i0 + NB - 1) m = i1 - i0 + 1 - DO i = i0, i1 - cth(i - i0 + 1) = COS(sea%k(i)*(x*sea%cosb(i) + y*sea%sinb(i)) - sea%omega(i)*t + sea%phase(i)) - END DO + IF (i1 <= NKEEP) THEN + cth(1:m) = ckeep(i0:i1) + ELSE + DO i = i0, i1 + cth(i - i0 + 1) = COS(sea%k(i)*(x*sea%cosb(i) + y*sea%sinb(i)) - sea%omega(i)*t + sea%phase(i)) + END DO + END IF DO i = i0, i1 sth(i - i0 + 1) = SIN(sea%k(i)*(x*sea%cosb(i) + y*sea%sinb(i)) - sea%omega(i)*t + sea%phase(i)) END DO + ! cosh(k(z+h))/sinh(kh) and sinh(k(z+h))/sinh(kh) without overflow, from e1 = EXP(k z) + ! and e2 = EXP(-k (z + 2h)) = e2kd/e1 at the synthesis depth: one exponential per + ! component. A component whose e1 or e2kd underflows, or a call at another depth than the + ! synthesis, takes the second exponential instead. DO i = i0, i1 - ! cosh(k(z+h))/sinh(kh) and sinh(k(z+h))/sinh(kh) without overflow - e1 = EXP(sea%k(i)*z_eval) - e2 = EXP(-sea%k(i)*(z_eval + 2.0_wp*depth)) - den = CD_ONE - sea%e2kd(i) - ch(i - i0 + 1) = (e1 + e2)/den - sh(i - i0 + 1) = (e1 - e2)/den + ch(i - i0 + 1) = EXP(sea%k(i)*z_eval) END DO + n_small = 0 + IF (ABS(depth - sea%depth) <= CD_ZERO) THEN + n_small = COUNT(ch(1:m) <= TINY(CD_ONE) .OR. sea%e2kd(i0:i1) <= TINY(CD_ONE)) + DO i = i0, i1 + e1 = ch(i - i0 + 1) + e2 = sea%e2kd(i)/MAX(e1, TINY(CD_ONE)) + den = CD_ONE - sea%e2kd(i) + ch(i - i0 + 1) = (e1 + e2)/den + sh(i - i0 + 1) = (e1 - e2)/den + END DO + END IF + IF (n_small > 0 .OR. ABS(depth - sea%depth) > CD_ZERO) THEN + DO i = i0, i1 + e1 = EXP(sea%k(i)*z_eval) + IF (ABS(depth - sea%depth) <= CD_ZERO .AND. e1 > TINY(CD_ONE) .AND. sea%e2kd(i) > TINY(CD_ONE)) CYCLE + e2 = EXP(-sea%k(i)*(z_eval + 2.0_wp*depth)) + den = CD_ONE - sea%e2kd(i) + ch(i - i0 + 1) = (e1 + e2)/den + sh(i - i0 + 1) = (e1 - e2)/den + END DO + END IF DO i = 1, m a = scale*sea%amplitude(i0 + i - 1)*sea%omega(i0 + i - 1) uh = a*ch(i)*cth(i) diff --git a/tests/test_wave_spectra.f90 b/tests/test_wave_spectra.f90 index 5b4f1bc..c1a3373 100644 --- a/tests/test_wave_spectra.f90 +++ b/tests/test_wave_spectra.f90 @@ -12,7 +12,9 @@ PROGRAM test_wave_spectra !! band-averaged periodogram matches S(omega); two trains add in energy; a spread !! sea gives the lateral/longitudinal velocity-variance ratio 1/(2s + 1). !! 4. A long-crested JONSWAP train reproduces the original JONSWAP synthesis and its - !! kinematics (same seed, same components, same Wheeler convention). + !! kinematics (same seed, same components, same Wheeler convention); the kinematics + !! of spread seas below and above the kept-phase length of CD_Sea_Kinematics equal a + !! direct per-component sum. !! 5. Deck rows: spectral waves rows, WaveSpreading, wavetrain rows run end to end, and !! malformed or conflicting rows are rejected. USE CableDyn_Precision, ONLY: wp @@ -36,6 +38,7 @@ PROGRAM test_wave_spectra CALL case_two_trains() CALL case_spread_velocity_ratio() CALL case_legacy_jonswap_identity() + CALL case_kinematics_reference() CALL case_parse_rows() CALL case_deck_runs() IF (nfail > 0) THEN @@ -493,6 +496,87 @@ SUBROUTINE case_legacy_jonswap_identity() END DO END SUBROUTINE case_legacy_jonswap_identity + SUBROUTINE case_kinematics_reference() + !! CD_Sea_Kinematics against a direct per-component sum of the same linear kinematics with + !! Wheeler stretching, at points below the surface (one near the seabed) and one above it: + !! 1. 9 x 200 components: every phase cosine kept between the two passes; + !! 2. 31 x 200 components: the components past the kept length recomputed; + !! 3. a short sea in 1000 m: EXP(-2 k h), and near the seabed EXP(k z), underflow for the + !! short components, which take the two-exponential form; + !! 4. kinematics asked at another depth than the synthesis: the two-exponential form. + ! allocatable: a default-initialised local with allocatable components draws spurious + ! -Wuninitialized warnings from gfortran 16 (DEVELOPMENT.md) + TYPE(CD_SeaType), ALLOCATABLE :: sea + TYPE(CD_WaveTrainType) :: tr + REAL(wp), PARAMETER :: SCL = 0.85_wp + REAL(wp), PARAMETER :: PX(4) = [0.0_wp, 37.5_wp, -210.0_wp, 4.0_wp] + REAL(wp), PARAMETER :: PY(4) = [0.0_wp, -12.0_wp, 95.0_wp, 1.0_wp] + REAL(wp), PARAMETER :: PT(4) = [0.0_wp, 13.7_wp, 401.25_wp, 7.0_wp] + INTEGER, PARAMETER :: NDIR(4) = [9, 31, 9, 9] + REAL(wp), PARAMETER :: SYNTH_DEPTH(4) = [120.0_wp, 120.0_wp, 1000.0_wp, 120.0_wp] + REAL(wp), PARAMETER :: CALL_DEPTH(4) = [120.0_wp, 120.0_wp, 1000.0_wp, 150.0_wp] + REAL(wp), PARAMETER :: PERIOD(4) = [10.0_wp, 10.0_wp, 4.0_wp, 10.0_wp] + REAL(wp) :: eta, v(3), a(3), p, eta_r, v_r(3), a_r(3), p_r, th, z_e, e1, e2, ch, sh, amp, depth, pz(4) + INTEGER :: es, js, ip, i + CHARACTER(256) :: em + CHARACTER(48) :: tag + ALLOCATE (sea) + DO js = 1, SIZE(NDIR) + CALL CD_Sea_Reset(sea) + tr%kind = CD_TRAIN_JONSWAP + tr%height = MERGE(1.0_wp, 6.0_wp, PERIOD(js) < 5.0_wp) + tr%period = PERIOD(js) + tr%gamma = 3.3_wp + tr%direction = 20.0_wp + tr%spreading = 4.0_wp + sea%n_dir = NDIR(js) + CALL CD_Sea_Add_Train(sea, tr, es, em) + CALL CD_Sea_Synthesise(sea, SYNTH_DEPTH(js), G, 3, es, em) + WRITE (tag, '(A,I0,A,I0)') 'sea-reference:case ', js, ', n_comp=', sea%n_comp + CALL require(es == CD_SEA_OK .AND. sea%n_comp == 200*NDIR(js), TRIM(tag)//':synth') + IF (es /= CD_SEA_OK) CYCLE + depth = CALL_DEPTH(js) + pz = [-1.5_wp, -20.0_wp, -0.95_wp*depth, 30.0_wp] + DO ip = 1, SIZE(PX) + CALL CD_Sea_Kinematics(sea, PX(ip), PY(ip), pz(ip), PT(ip), depth, SCL, eta, v, a, es, em, pdyn=p) + CALL require(es == CD_SEA_OK, TRIM(tag)//':kinematics: '//TRIM(em)) + eta_r = 0.0_wp + DO i = 1, sea%n_comp + th = sea%k(i)*(PX(ip)*sea%cosb(i) + PY(ip)*sea%sinb(i)) - sea%omega(i)*PT(ip) + sea%phase(i) + eta_r = eta_r + SCL*sea%amplitude(i)*COS(th) + END DO + v_r = 0.0_wp + a_r = 0.0_wp + p_r = 0.0_wp + IF (pz(ip) <= eta_r) THEN + z_e = MIN(0.0_wp, (pz(ip) - eta_r)*depth/(depth + eta_r)) + DO i = 1, sea%n_comp + th = sea%k(i)*(PX(ip)*sea%cosb(i) + PY(ip)*sea%sinb(i)) - sea%omega(i)*PT(ip) + sea%phase(i) + e1 = EXP(sea%k(i)*z_e) + e2 = EXP(-sea%k(i)*(z_e + 2.0_wp*depth)) + ch = (e1 + e2)/(1.0_wp - sea%e2kd(i)) + sh = (e1 - e2)/(1.0_wp - sea%e2kd(i)) + amp = SCL*sea%amplitude(i)*sea%omega(i) + v_r = v_r + amp*[ch*COS(th)*sea%cosb(i), ch*COS(th)*sea%sinb(i), sh*SIN(th)] + amp = amp*sea%omega(i) + a_r = a_r + amp*[ch*SIN(th)*sea%cosb(i), ch*SIN(th)*sea%sinb(i), -sh*COS(th)] + p_r = p_r + amp/sea%k(i)*ch*COS(th) + END DO + END IF + CALL require(ABS(eta - eta_r) <= 1.0e-11_wp*MAX(1.0_wp, ABS(eta_r)), TRIM(tag)//':eta') + CALL require(nan_max_abs(v - v_r) <= 1.0e-11_wp*MAX(1.0_wp, MAXVAL(ABS(v_r))), TRIM(tag)//':velocity') + CALL require(nan_max_abs(a - a_r) <= 1.0e-11_wp*MAX(1.0_wp, MAXVAL(ABS(a_r))), TRIM(tag)//':acceleration') + CALL require(ABS(p - p_r) <= 1.0e-11_wp*MAX(1.0_wp, ABS(p_r)), TRIM(tag)//':pressure') + END DO + ! the last point lies above the surface: the elevation only, no kinematics + CALL require(pz(SIZE(pz)) > eta_r .AND. MAXVAL(ABS(v)) <= 0.0_wp .AND. MAXVAL(ABS(a)) <= 0.0_wp, & + TRIM(tag)//':above-surface') + ! case 3 reaches both underflows that the two-exponential form handles + IF (js == 3) CALL require(ANY(sea%e2kd <= TINY(1.0_wp)) .AND. & + ANY(EXP(sea%k*(-0.95_wp*depth)) <= TINY(1.0_wp)), TRIM(tag)//':underflow-reached') + END DO + END SUBROUTINE case_kinematics_reference + SUBROUTINE case_parse_rows() TYPE(CD_WaveTrainType) :: tr CHARACTER(16) :: kw