MoorDyn_Init raises FE_OVERFLOW in CFL::cfl2dt (stationary IC solver); fatal when the host traps floating-point exceptions

[hgopalan]

Summary

MoorDyn_Init() raises FE_OVERFLOW inside the stationary initial-condition solver. The overflow is an intermediate value that is discarded, so the result is correct, but a host program that runs with floating-point traps enabled (as several CFD codes do to catch NaNs early) is killed inside MoorDyn_Init (SIGILL on macOS arm64, SIGFPE on Linux with feenableexcept).

Version: MoorDyn-C v2.7.1 (tag v2.7.1, commit 319ff99), built with CMake, default options (bundled Eigen, PYTHON_WRAPPER=OFF), AppleClang on macOS 26 arm64. The same expression is in master.

Where

moordyn::CFL::cfl2dt(const real cfl, const real v) in source/Util/CFL.hpp returns cfl * length() / v. The CFL base class initialises _l to std::numeric_limits<real>::max() (CFL.hpp line 62) and only lines set it (Line.cpp line 251), so for points, rods and bodies the product cfl * max() overflows to +inf.

StationaryScheme::Step() in source/Time.cpp (around line 206) calls it for every object to limit the next step:

real v = 0.5 * dt * _error;
for (auto obj : lines)  new_dt = (std::min)(new_dt, obj->cfl2dt(cfl, v));
for (auto obj : points) new_dt = (std::min)(new_dt, obj->cfl2dt(cfl, v));
...

The +inf from a point is discarded by std::min, which is why the solver still works, but the overflow flag is raised on every stationary step. A backtrace with the trap enabled (RelWithDebInfo build):

stop reason = EXC_BAD_INSTRUCTION
frame #0: libmoordyn.2.dylib`moordyn::CFL::cfl2dt(this=..., cfl=0.0453, v=0.00408) const at CFL.hpp:90 [inlined]
frame #1: moordyn::time::StationaryScheme::Step(double&)   (Time.cpp)
frame #2: moordyn::MoorDyn::icStationary()
frame #3: moordyn::MoorDyn::Init(double const*, double const*, bool)
frame #4: MoorDyn_Init

Reproducer

A single line between two fixed points in a light fluid (an overhead conductor in air), external wave kinematics on. No other dependency than libmoordyn.

span.txt:

MoorDyn-C input for a single fixed-fixed conductor span
----------------------- LINE TYPES ------------------------------------------
TypeName   Diam     Mass/m     EA         BA/-zeta    EI         Cd     Ca     CdAx    CaAx
(name)     (m)      (kg/m)     (N)        (N-s/-)     (N-m^2)    (-)    (-)    (-)     (-)
drake      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     drake      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)
1             WaveKin       the fluid kinematics are provided through the API (-)
0             ICgenDynamic  stationary initial-condition solver
1             disableOutput
1             disableOutTime
------------------------- need this line --------------------------------------

trap_repro.c:

#include <moordyn/MoorDyn2.h>
#include <fenv.h>
#include <stdio.h>
#include <string.h>
#pragma STDC FENV_ACCESS ON

static void enable_overflow_trap (void)
{
#if defined(__aarch64__) && defined(__APPLE__)
    fenv_t env; fegetenv(&env); env.__fpcr |= __fpcr_trap_overflow; fesetenv(&env);
#elif defined(__GLIBC__)
    feenableexcept(FE_OVERFLOW);
#endif
}

int main (int argc, char** argv)
{
    if (argc > 2 && strcmp(argv[2], "trap") == 0) { enable_overflow_trap(); }
    feclearexcept(FE_ALL_EXCEPT);
    MoorDyn s = MoorDyn_Create(argv[1]);
    MoorDyn_SetVerbosity(s, MOORDYN_ERR_LEVEL);
    const int rc = MoorDyn_Init(s, NULL, NULL);
    printf("MoorDyn_Init rc = %d; FE_OVERFLOW raised: %s; FE_DIVBYZERO %s; FE_INVALID %s\n", rc,
           fetestexcept(FE_OVERFLOW) ? "yes" : "no", fetestexcept(FE_DIVBYZERO) ? "yes" : "no",
           fetestexcept(FE_INVALID) ? "yes" : "no");
    MoorDyn_Close(s);
    return 0;
}

$ cc -I$PREFIX/include -L$PREFIX/lib -lmoordyn -Wl,-rpath,$PREFIX/lib trap_repro.c -o trap_repro
$ ./trap_repro span.txt
MoorDyn_Init rc = 0; FE_OVERFLOW raised: yes; FE_DIVBYZERO no; FE_INVALID no
$ ./trap_repro span.txt trap
Illegal instruction: 4

Only the overflow flag is raised; the invalid and divide-by-zero traps pass.

Suggested fix

Keep the sentinel out of the arithmetic, for example in CFL.hpp:

virtual inline real cfl2dt(const real cfl, const real v) const
{
    if (_l == (std::numeric_limits<real>::max)()) return _l;
    return cfl * length() / v;
}

or evaluate length() / v * cfl, which does not overflow for the sentinel when v >= cfl. Either way StationaryScheme::Step gets the same new_dt as today without raising the flag.

Happy to open a pull request with the one-line guard if you prefer.

[sanguinariojoe]

Can you please confirm that after patching everything works just fine, I mean no other floating point exceptions are happening during the simulation?

If so I can patch it and launch a minor release next week

Thanks!

[hgopalan]

Thanks @sanguinariojoe. Short answer: yes for this case, but your question was a good one, because other configurations had more of these.

This case: with the CFL guard, the span above runs MoorDyn_Init plus 60 s of simulated time (6000 steps with a gusty external cross-flow) with the invalid, div-by-zero and overflow traps all enabled. It raises no flags and gives the same fairlead tension as before.

Other configurations: I also ran the whole test suite with the three traps enabled in every test process. Even with the CFL guard, 12 of the 31 tests were still killed. The causes were:

• NatFreqCFL::cfl2dt(cfl) multiplies the cfl = max() sentinel that MoorDyn::Init() sets when dtM is prescribed. That overflows for any line whose segment natural period is > 1 s.
• CFL::cfl2dt(cfl, v) divides by v == 0 in the stationary solver.
• Rod::setup() computes i / N with N = 0 for zero-length rods (0/0).
• The Catenary() IC solver's Newton-Raphson iterations step out of domain for weightless or buoyant lines. The solver already detects the failure and the line falls back to the linear profile, but the flags leak to the host.

#405 fixes all of these, without changing results: the 76 *.out files written by the suite are byte-identical to dev. With the PR, all 32 tests (31 existing plus a new portable fpe test that checks fetestexcept after create/init/step) pass with all three traps on.