From 1a060c62d701a95cb80e7e9b78ff382dd5f29cb8 Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Fri, 2 Oct 2026 07:13:23 -0600 Subject: [PATCH 1/4] fix(core): Do not raise floating point exceptions with CFL sentinel values Points, rods and bodies keep length = max() as a sentinel, and MoorDyn::Init() sets cfl = max() when dtM is prescribed. Multiplying those sentinels raised FE_OVERFLOW, and a null velocity in the stationary solver raised FE_DIVBYZERO. The resulting inf/nan values were discarded by std::min(), so the time step is unchanged, but the host program is killed if it traps those exceptions. Fixes #404 Co-Authored-By: Claude Opus 5.5 --- source/Util/CFL.hpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/source/Util/CFL.hpp b/source/Util/CFL.hpp index 0bba6a4a..4a38c41b 100644 --- a/source/Util/CFL.hpp +++ b/source/Util/CFL.hpp @@ -87,6 +87,11 @@ class DECLDIR CFL */ virtual inline real cfl2dt(const real cfl, const real v) const { + // Objects without a characteristic length keep the sentinel value, + // and a null velocity does not limit the timestep either. Return the + // sentinel instead of raising FE_OVERFLOW or FE_DIVBYZERO + if ((length() == (std::numeric_limits::max)()) || (v <= 0.0)) + return (std::numeric_limits::max)(); return cfl * length() / v; } @@ -134,7 +139,14 @@ class DECLDIR NatFreqCFL : public CFL * @param cfl CFL factor * @return The timestUtilep */ - inline real cfl2dt(const real cfl) const { return cfl * period(); } + inline real cfl2dt(const real cfl) const + { + // cfl = max() is used to not limit the timestep, see + // MoorDyn::Init(). Do not multiply it, which raises FE_OVERFLOW + if (cfl == (std::numeric_limits::max)()) + return cfl; + return cfl * period(); + } /** @brief Get the CFL factor from a timestep * @param dt Timestep From 8a44649e6ba7933d39de642a7be7ae17d3b2935f Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Fri, 2 Oct 2026 07:13:23 -0600 Subject: [PATCH 2/4] fix(core): Do not divide by zero when setting up zero-length rods The interpolation factor i / N was 0 / 0 for zero-length rods (N = 0), raising FE_INVALID and setting the node position to nan until it was overwritten later on. Co-Authored-By: Claude Opus 5.5 --- source/Rod.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/source/Rod.cpp b/source/Rod.cpp index a499ef96..07516fba 100644 --- a/source/Rod.cpp +++ b/source/Rod.cpp @@ -165,7 +165,8 @@ Rod::setup(int number_in, const vec org = endCoords(Eigen::seqN(0, 3)); const vec dst = endCoords(Eigen::seqN(3, 3)); for (unsigned int i = 0; i <= N; i++) { - const real f = i / (real)N; + // N = 0 for zero-length rods, where i / N would be 0 / 0 + const real f = N ? i / (real)N : 0.0; r[i] = org + f * (dst - org); rd[i] = vec::Zero(); } From eec80259b5d03ff502a7318ced45079c5aac96d0 Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Fri, 2 Oct 2026 07:13:23 -0600 Subject: [PATCH 3/4] fix(core): Discard floating point exceptions raised by the catenary IC solver The Newton-Raphson iterations of Catenary() may visit points out of the domain of the involved functions (e.g. weightless or buoyant lines), raising FE_INVALID, FE_DIVBYZERO or FE_OVERFLOW. The solver detects those failures by itself, and Line::initialize() falls back to the linear profile, so the raised exceptions are now held while the solver runs and discarded afterwards. Co-Authored-By: Claude Opus 5.5 --- source/Line.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/source/Line.cpp b/source/Line.cpp index cbed0f98..be47313e 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -34,6 +34,7 @@ #include "QSlines.hpp" #include "Util/Interp.hpp" #include +#include // #include #include @@ -662,6 +663,14 @@ Line::initialize() COSPhi = (r[N][0] - r[0][0]) / XF; SINPhi = (r[N][1] - r[0][1]) / XF; + // The Newton-Raphson iterations of the catenary solver may visit + // points out of the domain of the involved functions (e.g. weightless + // or buoyant lines), raising floating point exceptions. The solver + // detects those failures by itself, so the raised exceptions are + // discarded instead of reaching (and eventually trapping in) the + // calling program + std::fenv_t fenv; + std::feholdexcept(&fenv); int success = Catenary(XF, ZF, UnstrLen, @@ -678,6 +687,7 @@ Line::initialize() Xl, Zl, Te); + std::fesetenv(&fenv); if (success >= 0) { From b9f00b1e425bb9d57f542a7c7ea9b0a61b7d452b Mon Sep 17 00:00:00 2001 From: "Harish @ NLR" Date: Fri, 2 Oct 2026 07:13:23 -0600 Subject: [PATCH 4/4] test: Check that no floating point exceptions are raised Co-Authored-By: Claude Opus 5.5 --- tests/CMakeLists.txt | 1 + tests/Mooring/fpe/span.txt | 24 +++++++++ tests/fpe.cpp | 108 +++++++++++++++++++++++++++++++++++++ 3 files changed, 133 insertions(+) create mode 100644 tests/Mooring/fpe/span.txt create mode 100644 tests/fpe.cpp diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 46dcc407..19cee474 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -57,6 +57,7 @@ set(CATCH2_TESTS aca wilson line_break + fpe ) function(make_executable test_name, extension) diff --git a/tests/Mooring/fpe/span.txt b/tests/Mooring/fpe/span.txt new file mode 100644 index 00000000..5dbfe6e9 --- /dev/null +++ b/tests/Mooring/fpe/span.txt @@ -0,0 +1,24 @@ +MoorDyn input file for a single line between two fixed points in a light fluid +----------------------- LINE TYPES ------------------------------------------ +TypeName Diam Mass/m EA BA/-zeta EI Cd Ca CdAx CaAx +(name) (m) (kg/m) (N) (N-s/-) (N-m^2) (-) (-) (-) (-) +cable 0.0281 1.628 3e+07 -0.5 0 1 1.0 0.0 0.0 +---------------------- POINT PROPERTIES -------------------------------- +ID Type X Y Z Mass Volume CdA Ca +(#) (-) (m) (m) (m) (kg) (m^3) (m^2) (-) +1 Fixed 0.0 0.0 -100 0 0 0 0 +2 Fixed 300 0.0 -100 0 0 0 0 +---------------------- LINES ---------------------------------------- +ID LineType AttachA AttachB UnstrLen NumSegs LineOutputs +(#) (name) (#) (#) (m) (-) (-) +1 cable 1 2 301.5 20 - +---------------------- OPTIONS ----------------------------------------- +0 writeLog Write a log file +0.001 dtM time step to use in the line integration (s) +9.81 g gravity (m/s^2) +1.2 WtrDnsty fluid density (kg/m^3) +1000 WtrDpth depth of the flat bottom (m) +0 ICgenDynamic stationary initial-condition solver +1 disableOutput +1 disableOutTime +------------------------- need this line -------------------------------------- diff --git a/tests/fpe.cpp b/tests/fpe.cpp new file mode 100644 index 00000000..cd9221d8 --- /dev/null +++ b/tests/fpe.cpp @@ -0,0 +1,108 @@ +/* + * Copyright (c) 2026 Harish Gopalan + * + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * 3. Neither the name of the copyright holder nor the names of its + * contributors may be used to endorse or promote products derived from + * this software without specific prior written permission. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + */ + +/** @file fpe.cpp + * Check that no floating point exceptions are raised while creating, + * initializing and integrating some systems. + * + * Programs that trap the floating point exceptions (e.g. with feenableexcept()) + * are killed by FE_INVALID, FE_DIVBYZERO or FE_OVERFLOW, even if MoorDyn + * discards the resulting value afterwards + */ + +#include "MoorDyn2.h" +#include +#include +#include +#include + +#define FPE_FLAGS (FE_INVALID | FE_DIVBYZERO | FE_OVERFLOW) + +/** @brief Create, initialize and integrate a system for some time steps + * @param path The input file + * @param cpld_point The coupled point index, 0 if there are no coupled DOFs + * @return The floating point exception flags raised + */ +int +raised_fpe(const std::string& path, unsigned int cpld_point = 0) +{ + std::feclearexcept(FE_ALL_EXCEPT); + + MoorDyn system = MoorDyn_Create(path.c_str()); + REQUIRE(system); + REQUIRE(MoorDyn_SetVerbosity(system, MOORDYN_ERR_LEVEL) == MOORDYN_SUCCESS); + unsigned int n_dof; + REQUIRE(MoorDyn_NCoupledDOF(system, &n_dof) == MOORDYN_SUCCESS); + REQUIRE(n_dof == (cpld_point ? 3 : 0)); + std::vector x(n_dof, 0.0), xd(n_dof, 0.0), f(n_dof, 0.0); + if (cpld_point) { + auto point = MoorDyn_GetPoint(system, cpld_point); + REQUIRE(point); + REQUIRE(MoorDyn_GetPointPos(point, x.data()) == MOORDYN_SUCCESS); + } + REQUIRE(MoorDyn_Init(system, x.data(), xd.data()) == MOORDYN_SUCCESS); + + double t = 0.0, dt = 0.01; + for (unsigned int i = 0; i < 10; i++) { + REQUIRE(MoorDyn_Step(system, x.data(), xd.data(), f.data(), &t, &dt) == + MOORDYN_SUCCESS); + } + REQUIRE(MoorDyn_Close(system) == MOORDYN_SUCCESS); + + return std::fetestexcept(FPE_FLAGS); +} + +TEST_CASE("Line between fixed points, stationary IC") +{ + // Points, rods and bodies have no characteristic length for the CFL + REQUIRE(raised_fpe("Mooring/fpe/span.txt") == 0); +} + +TEST_CASE("Weightless line catenary IC") +{ + // The catenary solver cannot find a solution, so it is discarded + REQUIRE(raised_fpe("Mooring/pendulum.txt") == 0); +} + +TEST_CASE("Body with a soft line and a prescribed time step") +{ + // cfl = max() if dtM is prescribed, and the line natural period > 1 s + REQUIRE(raised_fpe("Mooring/body_tests/bodyDrag.txt") == 0); +} + +TEST_CASE("Zero-length rods") +{ + REQUIRE(raised_fpe("Mooring/local_euler/complex_system.txt", 4) == 0); +} + +TEST_CASE("Coupled fairlead, stationary IC") +{ + REQUIRE(raised_fpe("Mooring/polyester/simple.txt", 2) == 0); +}