Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions source/Line.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -34,6 +34,7 @@
#include "QSlines.hpp"
#include "Util/Interp.hpp"
#include <tuple>
#include <cfenv>
// #include <random>
#include <iomanip>

Expand Down Expand Up @@ -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,
Expand All @@ -678,6 +687,7 @@ Line::initialize()
Xl,
Zl,
Te);
std::fesetenv(&fenv);

if (success >= 0) {

Expand Down
3 changes: 2 additions & 1 deletion source/Rod.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
}
Expand Down
14 changes: 13 additions & 1 deletion source/Util/CFL.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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<real>::max)()) || (v <= 0.0))
return (std::numeric_limits<real>::max)();
return cfl * length() / v;
}

Expand Down Expand Up @@ -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<real>::max)())
return cfl;
return cfl * period();
}

/** @brief Get the CFL factor from a timestep
* @param dt Timestep
Expand Down
1 change: 1 addition & 0 deletions tests/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ set(CATCH2_TESTS
aca
wilson
line_break
fpe
)

function(make_executable test_name, extension)
Expand Down
24 changes: 24 additions & 0 deletions tests/Mooring/fpe/span.txt
Original file line number Diff line number Diff line change
@@ -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 --------------------------------------
108 changes: 108 additions & 0 deletions tests/fpe.cpp
Original file line number Diff line number Diff line change
@@ -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 <cfenv>
#include <string>
#include <vector>
#include <catch2/catch_test_macros.hpp>

#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<double> 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);
}