diff --git a/CHANGELOG.md b/CHANGELOG.md index 3d2fb7e..8ea3bb5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -68,6 +68,13 @@ All notable changes to CableDyn are recorded in this file. The format follows ### Fixed +- Finite-EI cables with seabed friction: a restored state (the C API `CableDyn_RestoreState`, + OpenFAST correction iterations) could continue differently in the last bits from the step taken + before the restore, breaking the documented bit-identity. A sliding node's moved friction + anchor reproduces the capped friction force only to rounding, so the committed force cached + from the step and a fresh evaluation differ slightly; the restore discarded the cache and took + the fresh path. The snapshot now keeps the cached force with the state, and the restore brings + it back unless a load setter invalidated it since. Runs without restores do not change. - OpenFAST with `dtM` larger than `DT`: an `openfast -restart` re-sampled the SeaState fields of the committed step at its end positions instead of restoring the fields the step held (sampled at its start positions), so the loads served until the next solve, and the run after them, diff --git a/src/CableDyn_HermiteCableDynamic.f90 b/src/CableDyn_HermiteCableDynamic.f90 index d239da1..3e5b24b 100644 --- a/src/CableDyn_HermiteCableDynamic.f90 +++ b/src/CableDyn_HermiteCableDynamic.f90 @@ -281,12 +281,22 @@ MODULE CableDyn_HermiteCableDynamic ! the current iterate (ws_Rk) and of a line-search trial (ws_Rt), and the force at the ! committed state F_n = R(q_n, v_n, t_n). F_n is cached from the previous step's converged ! evaluation and keyed on the exact state it was evaluated at (fc_q, fc_v, fc_t, fc_d0) - ! plus fc_valid, which every load-configuration setter and Restore clears. A key mismatch - ! (host-updated boundary, restored or externally written state) recomputes it, so the - ! cache only ever replaces a bitwise-identical evaluation. + ! plus fc_valid, which every load-configuration setter clears (counting the clears in + ! fc_clears). A key mismatch (host-updated boundary, externally written state) recomputes + ! it. The cache equals a fresh evaluation in value but not always in its last bits: the + ! friction return mapping of the commit moves a sliding node's anchor, and the spring force + ! of the moved anchor reproduces the capped force only to rounding. A run is therefore + ! bit-reproducible along one path: Snapshot keeps the cache with the state, and Restore + ! brings it back unless the cache was invalidated since (fc_clears), so a restored model + ! takes the same path as one that was never rolled back. REAL(wp), ALLOCATABLE :: ws_Rk(:), ws_Rt(:), fc_R(:), fc_q(:), fc_v(:) REAL(wp) :: fc_t = 0.0_wp, fc_d0(3, 2) = 0.0_wp LOGICAL :: fc_valid = .FALSE. + INTEGER :: fc_clears = 0 + REAL(wp), ALLOCATABLE :: snap_fc_R(:), snap_fc_q(:), snap_fc_v(:) + REAL(wp) :: snap_fc_t = 0.0_wp, snap_fc_d0(3, 2) = 0.0_wp + LOGICAL :: snap_fc_valid = .FALSE. + INTEGER :: snap_fc_clears = 0 ! Generalised-alpha blend: .TRUE. blends the forces (the default), .FALSE. evaluates the ! forces at the blended configuration q_alpha = (1-af) q_{n+1} + af q_n. LOGICAL :: force_blend = CD_HC_FORCE_BLEND_DEFAULT @@ -713,6 +723,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Init(model, l0, EA, EI, rho_a, w, seed, fixed_dof model%ws_recovery_q0(ndof), model%ws_recovery_v0(ndof), model%ws_recovery_a0(ndof)) ALLOCATE (model%ws_Rk(ndof), model%ws_Rt(ndof), model%fc_R(ndof), model%fc_q(ndof), model%fc_v(ndof)) model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. model%force_blend = CD_HC_FORCE_BLEND_DEFAULT model%max_step_rotation = CD_ZERO @@ -778,6 +789,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Contact(model, kn_node, cn_node, mu, ErrStat, model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN CALL fail_contact('model not initialized'); RETURN @@ -926,6 +938,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Friction_Axial(model, mu_axial, ErrStat, ErrM model%contact_fr_aniso = ABS(mu_axial - model%contact_mu) > CD_ZERO ! A friction-law change invalidates the cached committed force F_n and the carried factor. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. END SUBROUTINE CD_HermiteCable_Dyn_Set_Friction_Axial @@ -971,6 +984,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Friction_Anchors(model, anchors, ErrStat, Err model%fr_anchor = anchors model%fr_active = .TRUE. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. CALL solve_consistent_acceleration(model, es, em) IF (es /= CD_HCDYN_OK) THEN @@ -1049,6 +1063,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Attachments(model, node, mass, volume, cda, c RETURN END IF model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. a_old = model%a mgb_old = model%Mgb @@ -1701,6 +1716,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_ForceBlend(model, enabled, ErrStat, ErrMsg) END IF model%force_blend = enabled model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. END SUBROUTINE CD_HermiteCable_Dyn_Set_ForceBlend @@ -1924,6 +1940,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_EndConnection(model, k_rot, d0, ErrStat, ErrM ErrStat = CD_HCDYN_OK; ErrMsg = '' ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT @@ -2294,6 +2311,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Drag(model, rho_w, diam, cdn, cdt, waterline_ model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT; ErrMsg = 'CD_HermiteCable_Dyn_Set_Drag: model not initialised'; RETURN @@ -2393,6 +2411,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Axial_Damping(model, BA, ErrStat, ErrMsg) model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. ErrMsg = '' IF (.NOT. model%initialized) THEN @@ -2464,6 +2483,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_AddedMass(model, rho_w, diam, can, cat, water model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT; ErrMsg = 'CD_HermiteCable_Dyn_Set_AddedMass: model not initialised'; RETURN @@ -2560,6 +2580,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Waves(model, height, period, direction_deg, d model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT; ErrMsg = 'CD_HermiteCable_Dyn_Set_Waves: model not initialised'; RETURN @@ -2680,6 +2701,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Irregular_Waves(model, amplitude, period, pha model%rc_epoch = model%rc_epoch + 1 ! A load-configuration change invalidates the cached committed force F_n. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT @@ -2819,6 +2841,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Held_Fluid(model, fluid_velocity, fluid_accel ! factor stays usable (a per-step field update only perturbs the drag and inertia tangent; ! the first field's consistent-acceleration solve drops it). model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 IF (.NOT. model%initialized) THEN ErrStat = CD_HCDYN_BADINPUT ErrMsg = 'CD_HermiteCable_Dyn_Set_Held_Fluid: model not initialised'; RETURN @@ -3086,6 +3109,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Torsion_Frames(model, ends, ErrStat, ErrMsg, model%torsion%theta = th model%torsion%torque = (model%torsion%phi - th)/CD_HermiteTorsion_Compliance(model%torsion, model%l0) model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. END SUBROUTINE CD_HermiteCable_Dyn_Set_Torsion_Frames @@ -6141,6 +6165,7 @@ SUBROUTINE CD_HermiteCable_Dyn_Set_Torsion(model, torsion, ErrStat, ErrMsg) ! the snapshot predates the new load set (as for CD_HermiteCable_Dyn_Set_EndConnection) model%snap_valid = .FALSE. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. CALL CD_HermiteCable_Dyn_Recompute_Acceleration(model, es, em) IF (es /= CD_HCDYN_OK) THEN @@ -6332,6 +6357,16 @@ SUBROUTINE CD_HermiteCable_Dyn_Snapshot(model, ErrStat, ErrMsg) model%snap_torsion_phi = model%torsion%phi model%snap_torsion_ends = model%torsion%ends IF (model%fr_active) model%snap_fr_anchor = model%fr_anchor + ! the cached committed force goes with the state it describes (see fc_clears) + model%snap_fc_valid = model%fc_valid .AND. ALLOCATED(model%fc_R) + model%snap_fc_clears = model%fc_clears + IF (model%snap_fc_valid) THEN + model%snap_fc_R = model%fc_R + model%snap_fc_q = model%fc_q + model%snap_fc_v = model%fc_v + model%snap_fc_t = model%fc_t + model%snap_fc_d0 = model%fc_d0 + END IF model%snap_na_steps_done = model%na_steps_done model%snap_na_iter_sum = model%na_iter_sum model%snap_na_mn_active = model%na_mn_active @@ -6411,7 +6446,21 @@ SUBROUTINE CD_HermiteCable_Dyn_Restore(model, ErrStat, ErrMsg) model%max_step_rotation = model%snap_max_step_rotation model%recovery_event_count = model%snap_recovery_event_count model%recovery_max_substeps_used = model%snap_recovery_max_substeps_used - model%fc_valid = .FALSE. + ! The snapshot's cached committed force comes back with its state unless a setter has + ! invalidated the cache since the snapshot, so the next step takes the path it took after + ! the snapshot (the cache differs from a fresh evaluation in the last bits where a + ! friction anchor slid). Otherwise the next step evaluates it afresh. + IF (model%snap_fc_valid .AND. model%snap_fc_clears == model%fc_clears) THEN + model%fc_R = model%snap_fc_R + model%fc_q = model%snap_fc_q + model%fc_v = model%snap_fc_v + model%fc_t = model%snap_fc_t + model%fc_d0 = model%snap_fc_d0 + model%fc_valid = .TRUE. + ELSE + model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 + END IF model%tr_valid = .FALSE. END SUBROUTINE CD_HermiteCable_Dyn_Restore @@ -6493,7 +6542,10 @@ SUBROUTINE CD_HermiteCable_Dyn_End(model) model%endconn_k = CD_ZERO model%endconn_d0 = CD_ZERO IF (ALLOCATED(model%snap_q)) DEALLOCATE (model%snap_q, model%snap_v, model%snap_a) + IF (ALLOCATED(model%snap_fc_R)) DEALLOCATE (model%snap_fc_R, model%snap_fc_q, model%snap_fc_v) + model%snap_fc_valid = .FALSE. model%fc_valid = .FALSE. + model%fc_clears = model%fc_clears + 1 model%tr_valid = .FALSE. model%snap_t = CD_ZERO model%snap_endconn_d0 = CD_ZERO diff --git a/src/CableDyn_OpenFAST_HermiteFMF.f90 b/src/CableDyn_OpenFAST_HermiteFMF.f90 index 3fde501..ee0a94f 100644 --- a/src/CableDyn_OpenFAST_HermiteFMF.f90 +++ b/src/CableDyn_OpenFAST_HermiteFMF.f90 @@ -1694,6 +1694,7 @@ SUBROUTINE CD_HFMF_UnpackMirror(self, buf, t, ErrStat, ErrMsg) END IF ! the force-blend committed force reloads with the state it was evaluated at self%line%fc_valid = .FALSE. + self%line%fc_clears = self%line%fc_clears + 1 IF (self%line%force_blend) THEN IF (buf(k + 1) > 0.5_wp) THEN self%line%fc_R = buf(k + 2:k + 1 + self%line%ndof) diff --git a/tests/test_hermite_cable_control.f90 b/tests/test_hermite_cable_control.f90 index ba87404..590a74f 100644 --- a/tests/test_hermite_cable_control.f90 +++ b/tests/test_hermite_cable_control.f90 @@ -244,14 +244,15 @@ SUBROUTINE check_model_contract() CALL require(unchanged(a, ref) .AND. same_bits(a%l0_dot, ref%l0_dot) .AND. same_bits(a%q, ref%q) .AND. & same_bits(a%v, ref%v) .AND. same_bits(a%a, ref%a) .AND. same_bits2(a%fr_anchor, ref%fr_anchor), & 'A: restore brings back the lengths, rates, caches and state bit for bit') - ! Restore clears the carried Newton factor and the cached committed force, which ref still - ! holds. Where a step stops depends on the factor's age, so a restart is exact only with - ! tangent reuse off; and ref's next step would take its committed force from the cache - ! while a evaluates it afresh. Both models march with reuse off from an uncached - ! committed force, so they take the same path. + ! Restore clears the carried Newton factor, which ref still holds; where a step stops + ! depends on the factor's age, so a restart is exact only with tangent reuse off. The + ! cached committed force comes back with the snapshot (the length command after it does + ! not invalidate it), so a takes the same cached path as ref: on a friction line the cache + ! differs from a fresh evaluation in its last bits. a%tangent_reuse = .FALSE. ref%tangent_reuse = .FALSE. - ref%fc_valid = .FALSE. + CALL require(a%fc_valid .AND. ref%fc_valid .AND. same_bits(a%fc_R, ref%fc_R), & + 'A: restore brings back the cached committed force of the snapshot') DO k = 1, 3 CALL CD_HermiteCable_Dyn_Step(a, DT, 50, 1.0e-10_wp, es, em) CALL CD_HermiteCable_Dyn_Step(ref, DT, 50, 1.0e-10_wp, es, em)