diff --git a/.github/workflows/build-mingw.yml b/.github/workflows/build-mingw.yml index 07bb7a16..b5364bc8 100644 --- a/.github/workflows/build-mingw.yml +++ b/.github/workflows/build-mingw.yml @@ -21,7 +21,7 @@ jobs: os: [windows-latest] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Create folders run: | diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 51787002..44ddf61d 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -21,7 +21,7 @@ jobs: release_id: ${{ steps.create_release.outputs.id }} upload_url: ${{ steps.create_release.outputs.upload_url }} steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 if: github.event_name == 'push' - name: moordyn_tag_name @@ -71,7 +71,7 @@ jobs: os: [ubuntu-latest, windows-latest, macos-latest] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: moordyn_version id: moordyn_version diff --git a/.github/workflows/check-version.yml b/.github/workflows/check-version.yml index d8e4a110..be8a5050 100644 --- a/.github/workflows/check-version.yml +++ b/.github/workflows/check-version.yml @@ -9,7 +9,7 @@ jobs: name: Check that the version is valid runs-on: ubuntu-latest steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: moordyn_version id: moordyn_version diff --git a/.github/workflows/matlab.yml b/.github/workflows/matlab.yml index 4c9447a6..edd3990f 100644 --- a/.github/workflows/matlab.yml +++ b/.github/workflows/matlab.yml @@ -25,7 +25,7 @@ jobs: steps: - name: Check out repository - uses: actions/checkout@v6 + uses: actions/checkout@v7 - name: Create install folder run: | diff --git a/.github/workflows/mdf_verification.yml b/.github/workflows/mdf_verification.yml index eedb7ef4..fa954237 100644 --- a/.github/workflows/mdf_verification.yml +++ b/.github/workflows/mdf_verification.yml @@ -21,7 +21,7 @@ jobs: os: [ubuntu-latest] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Install deps run: | @@ -30,7 +30,7 @@ jobs: if: runner.os == 'Linux' - name: Setup Python - uses: actions/setup-python@v6 + uses: actions/setup-python@v7 id: setup-python - name: Install Python dependencies diff --git a/.github/workflows/memcheck.yml b/.github/workflows/memcheck.yml index 99999289..616ce98d 100644 --- a/.github/workflows/memcheck.yml +++ b/.github/workflows/memcheck.yml @@ -19,7 +19,7 @@ jobs: os: [ubuntu-latest] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Install valgrind run: | diff --git a/.github/workflows/python-wheels-emulated.yml b/.github/workflows/python-wheels-emulated.yml index 3b168859..c201833a 100644 --- a/.github/workflows/python-wheels-emulated.yml +++ b/.github/workflows/python-wheels-emulated.yml @@ -47,10 +47,10 @@ jobs: rm -rf docs extern source tests steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 # Used to host cibuildwheel - - uses: actions/setup-python@v6 + - uses: actions/setup-python@v7 with: python-version: '3.x' @@ -75,7 +75,7 @@ jobs: if: runner.os == 'Linux' - name: Build wheels - uses: pypa/cibuildwheel@v3.4.1 + uses: pypa/cibuildwheel@v4.2.1 with: output-dir: dist diff --git a/.github/workflows/python-wheels-test.yml b/.github/workflows/python-wheels-test.yml index f693dc6d..8832dc02 100644 --- a/.github/workflows/python-wheels-test.yml +++ b/.github/workflows/python-wheels-test.yml @@ -24,9 +24,9 @@ jobs: os: ${{ fromJSON(inputs.os) }} steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - - uses: actions/setup-python@v6 + - uses: actions/setup-python@v7 with: python-version: '3.12' diff --git a/.github/workflows/python-wheels.yml b/.github/workflows/python-wheels.yml index 908f6200..5b9f98b2 100644 --- a/.github/workflows/python-wheels.yml +++ b/.github/workflows/python-wheels.yml @@ -15,10 +15,10 @@ jobs: os: [ubuntu-22.04] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 # Used to host cibuildwheel - - uses: actions/setup-python@v6 + - uses: actions/setup-python@v7 with: python-version: '3.x' @@ -133,7 +133,7 @@ jobs: os: [ubuntu-22.04] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Download artifacts uses: actions/download-artifact@v8 diff --git a/.github/workflows/python-wrapper.yml b/.github/workflows/python-wrapper.yml index e8af0d1b..3ed4e003 100644 --- a/.github/workflows/python-wrapper.yml +++ b/.github/workflows/python-wrapper.yml @@ -21,10 +21,10 @@ jobs: os: [ubuntu-latest] steps: - - uses: actions/checkout@v6 + - uses: actions/checkout@v7 - name: Setup Python - uses: actions/setup-python@v6 + uses: actions/setup-python@v7 id: setup-python - name: Install Python dependencies diff --git a/CHANGELOG.md b/CHANGELOG.md index 64583714..136c907f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,40 @@ +## [2.7.2] - 2026-10-06 + +### πŸš€ Features + +- *(lines,waves)* Implement improvements from OpenFAST +- *(io)* Added `!` and `%` as comment characters in input files (to match OpenFAST) +- *(core)* Add water velocity outputs to rods, disable currents if water-kinematic ramping enabled + +### πŸ› Bug Fixes + +- *(tests)* Improve VIV test peak finding, add VIV test plotting script +- *(tests)* Adjust VIV test sensitivity for windows +- *(core)* Do not raise floating point exceptions with CFL sentinel values +- *(core)* Do not divide by zero when setting up zero-length rods +- *(core)* Discard floating point exceptions raised by the catenary IC solver + +### πŸ’Ό Other + +- *(deps)* Bump actions/checkout from 6 to 7 +- *(deps)* Bump pypa/cibuildwheel from 3.4.1 to 4.1.0 +- *(deps)* Bump actions/setup-python from 6 to 7 +- *(deps)* Bump pypa/cibuildwheel from 4.1.0 to 4.1.1 +- *(deps)* Bump pypa/cibuildwheel from 4.1.1 to 4.2.1 + +### πŸ“š Documentation + +- Clarify difference between per-node tension outputs and FairTen outputs +- Add VIV paper to references, clarify line acc output and WEC-Sim errors +- Update core object docstrings to reflect input files + +### πŸ§ͺ Testing + +- Check that no floating point exceptions are raised + +### βš™οΈ Miscellaneous Tasks + +- Increased the version to 2.7.2 ## [2.7.1] - 2026-06-10 ### πŸ› Bug Fixes @@ -80,6 +117,9 @@ ### πŸš€ Features +- *(body)* Add centripetal forces for rotating bodies +- Cleaned up the hack-ish state to generalize it for visco stuff. Made # a comment character in input files. Readded diable VIV during ICgen +- Viscoelastic model with constant and load dependent dynamic stiffness is now live! - *(rust)* First test using the Rust wrapper - *(core)* New flexible states - *(core)* Create a base class for the instances which provides a unique identifier @@ -93,6 +133,27 @@ ### πŸ› Bug Fixes +- Read first the writelog option, and then anything else +- The quaternions shall be renormalized to get the rotation matrix +- Freeze when writeLog is not the first option +- Centripetal force for parallel axes shall be null +- Rebranding to include centripetal forces on getNetForceAndMass +- Add a centripetal force to bodies with a excentric COG +- EulerZYX -> EulerXYZ on moordyn::Euler2Quat() +- Accelerations of Coupled/fixed bodies/rods when there are several isntances of them +- Odd treatment was meant for indexes from 1 to 3, not 0 to 2, and the matrix indexes were transposed +- EulerXYZ intrinsic angles instead of extrinsic +- Drop the patch to move from extrinsic to intrinsic Euler angles +- Make rod submergence calcs match what is in MDF (verified code) +- Make rod submergence calcs match what is in MDF (verified code) +- Some small fixes for viscoelastic stuff +- Removing files +- Removing one more file +- Cleaned up some time scheme stuff (added in more spots where Misc states should be accounted for) +- Some more stiffness fixes. Noteably before when reading a non-linear look up table it assumed stress strain. This is not correct and doesnt match theory paper or docs. CHanges to stress tension. +- Dynamic current inflile reading fix +- Some more notes explaining the synchronization model +- Major clean up to the synchronization model - Added dummy AllOutput to old API - Old AllOutput back to voids, cleans it up a bit - Rename AllOutput -> WriteOutputs (To match MD-F) @@ -143,12 +204,25 @@ ### πŸ’Ό Other +- MinGW needs the DECLDIR on Body::setState - *(deps)* Bump pypa/cibuildwheel from 2.22.0 to 2.23.0 - *(deps)* Bump pypa/cibuildwheel from 2.23.0 to 2.23.1 ### πŸ“š Documentation +- Clarify that intrinsic angles are considered, and link to the external resources +- Added # comment character instructions to docs +- Minor change to update viscoelastic docs +- Clarifying dtout - Updated docs for OpenFAST/openfast#2597 + +### πŸ§ͺ Testing + +- Excentrical body test +- Test the centripetal force on a simple case +- Strip the VTK from the test and attach two points to the body to get a more stable orbit +- Test the rotations +- Typo on the config file description ## [2.3.8] - 2024-12-12 ### πŸš€ Features @@ -178,14 +252,6 @@ - *(python)* Autonomous MoorPy IC exporter - *(docs and tests)* Initialization notes with a practical application - *(python)* Absolute paths are not accepted any longer -- Some small fixes for viscoelastic stuff -- Removing files -- Removing one more file -- Cleaned up some time scheme stuff (added in more spots where Misc states should be accounted for) -- Some more stiffness fixes. Noteably before when reading a non-linear look up table it assumed stress strain. This is not correct and doesnt match theory paper or docs. CHanges to stress tension. -- Dynamic current inflile reading fix -- Some more notes explaining the synchronization model -- Major clean up to the synchronization model - *(python)* Added the missing time scheme manipulation API entries - *(python)* Missing GetXXXNumberNodes missing entries - *(fortran)* Added the time scheme manipulation entries @@ -213,8 +279,6 @@ - Updates and typo fixes - References file rename - Fixing broken links -- Minor change to update viscoelastic docs -- Clarifying dtout - Add disable output to the docs and options list ### πŸ§ͺ Testing @@ -284,46 +348,6 @@ ### βš™οΈ Miscellaneous Tasks - Disable memcheck on lowe_and_langley, which is way too slow -## [2.3.3] - 2024-06-24 - -### πŸš€ Features - -- *(body)* Add centripetal forces for rotating bodies -- Cleaned up the hack-ish state to generalize it for visco stuff. Made # a comment character in input files. Readded diable VIV during ICgen -- Viscoelastic model with constant and load dependent dynamic stiffness is now live! - -### πŸ› Bug Fixes - -- Read first the writelog option, and then anything else -- The quaternions shall be renormalized to get the rotation matrix -- Freeze when writeLog is not the first option -- Centripetal force for parallel axes shall be null -- Rebranding to include centripetal forces on getNetForceAndMass -- Add a centripetal force to bodies with a excentric COG -- EulerZYX -> EulerXYZ on moordyn::Euler2Quat() -- Accelerations of Coupled/fixed bodies/rods when there are several isntances of them -- Odd treatment was meant for indexes from 1 to 3, not 0 to 2, and the matrix indexes were transposed -- EulerXYZ intrinsic angles instead of extrinsic -- Drop the patch to move from extrinsic to intrinsic Euler angles -- Make rod submergence calcs match what is in MDF (verified code) -- Make rod submergence calcs match what is in MDF (verified code) - -### πŸ’Ό Other - -- MinGW needs the DECLDIR on Body::setState - -### πŸ“š Documentation - -- Clarify that intrinsic angles are considered, and link to the external resources -- Added # comment character instructions to docs - -### πŸ§ͺ Testing - -- Excentrical body test -- Test the centripetal force on a simple case -- Strip the VTK from the test and attach two points to the body to get a more stable orbit -- Test the rotations -- Typo on the config file description ## [2.3.2] - 2024-06-24 ### πŸ’Ό Other diff --git a/CMakeLists.txt b/CMakeLists.txt index 7a8b8b10..447b7df9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,7 +1,7 @@ cmake_minimum_required(VERSION 3.10) set(MOORDYN_MAJOR_VERSION 2) set(MOORDYN_MINOR_VERSION 7) -set(MOORDYN_PATCH_VERSION 1) +set(MOORDYN_PATCH_VERSION 2) set(MOORDYN_VERSION ${MOORDYN_MAJOR_VERSION}.${MOORDYN_MINOR_VERSION}) project(Moordyn VERSION ${MOORDYN_VERSION}) diff --git a/docs/drivers.rst b/docs/drivers.rst index c5115593..58f4625c 100644 --- a/docs/drivers.rst +++ b/docs/drivers.rst @@ -806,7 +806,10 @@ WEC-Sim WEC-Sim is coupled with MoorDyn-C v2 as of Spring 2024 with `WEC-Sim PR #1212 `_. The original coupling was with MoorDyn v1 and the the publication can be found in the -:ref:`theory section `. The WEC-Sim source code can be found `here `_. +:ref:`theory section `. The WEC-Sim source code can be found `here `_. Note that the current +library-based coupling does not handle MoorDyn errors, causing WEC-Sim/MATLAB to crash if MoorDyn fails. This is a known issue and will be +fixed in the future. For debugging your simulations, you can test your MoorDyn model standalone using the example driver on the +`examples page `_. DualSPHysics ^^^^^^^^^^^^ diff --git a/docs/inputs.rst b/docs/inputs.rst index 6d27bc8a..55e2dab0 100644 --- a/docs/inputs.rst +++ b/docs/inputs.rst @@ -247,8 +247,8 @@ two fixed points located far from where your system is located. Most of the sections are set up to contain a table of input information. These tables begin with two preset lines that contain the column names and the corresponding units. These lines are followed by any number of lines containing -the entries in that section's table of inputs. # is the general comment character. If you are adding notes -to self after any of the lines, # will prevent MoorDyn from reading them. +the entries in that section's table of inputs. !, #, and % are comment characters. If you are adding notes +to self after any of the lines, any of these characters will prevent MoorDyn from reading the rest of that line. Examples of input files for MoorDyn-C can be found in the `test directory `_ (note that these do not include outputs because they are for tests). @@ -491,13 +491,14 @@ The columns are as follows: This `LineOutputs` entry expects a string of one or more characters without spaces, each character activating a given output property. A placeholder character such as β€œ-” should be used if no -outputs are wanted. Ten output properties are currently possible: +outputs are wanted. Eleven output properties are currently possible: - p – node positions - v – node velocities + - a – node accelerations (0 at end nodes, acceleration is never computed for end nodes) - U – wave/current velocities at each node - D – hydrodynamic drag force at each node - - t – tension force at each segment + - t – internal tension force at each segment - c – internal damping force at each segment - V - the cross-flow VIV lift force at each node - K - the curvature at each node @@ -661,6 +662,8 @@ The list of possible options is: grid, 3 = kinematics in a regular grid, 7 = Wave Component Summing. Details on these flags can be found :ref:`here `. - dtWave (0.25): The time step to evaluate the waves, only for wave grid (WaveKin = 3) (s) + - waveKin_rampT (0.0): Ramp time for water kinematics (s). The wave kinematics are linearly ramped + from zero at t=0 to their full values at t=waveKin_rampT. A value of 0 disables the ramp. - Currents (0): The currents model to use. 0 = none, 1 = steady in a regular grid, 2 = dynamic in a regular grid, 3 = WIP, 4 = WIP, 5 = 4D Current Grid. Details on these flags can be found :ref:`here `. @@ -801,8 +804,8 @@ Reference Points: Footnotes: -- The tension on the Line n fairlead can be output with the FAIRTEN[n] flag (see examples above) -- The tension on the Line n anchor can be output with the ANCHTEN[n] flag (see examples above) +- The total force on the Line n fairlead can be output with the FAIRTEN[n] flag (see examples above) +- The total force on the Line n anchor can be output with the ANCHTEN[n] flag (see examples above) - Object indicates output is for whole object, Node indicates output is for node of object - Coupled/fixed bodies and points will output acceleration 0 because no forces are calculated - When looking at the rotational outputs of coupled pinned rods that are hanging near vertical, @@ -810,8 +813,13 @@ Footnotes: rotations. Hanging pinned rods are inverted (end A over end B) in MoorDyn and the output range for roll/pitch of rods is +/- 180 degrees. - There are a couple additional outputs left over from OpenFAST conventions that don’t follow the - same format: FairTen and AnchTen. FairTen[n] is the same as Line[n]TenB. For example, the - fairlead tension of line 1 would be FAIRTEN1 or LINE1TENB. + same format: FairTen and AnchTen. FairTen[n] is the same as Line[n]TenB. For example, the total + fairlead force of line 1 would be FAIRTEN1 or LINE1TENB. +- Lines can also output the per-node tensions. These are the internal tensions in the line segment. This + differs from the global outputs from Line[n]TenB / FairTen[n], which is the net force on the end node of + the line. The net force includes the tension, AND the per-node weight, buoyancy, seabed-contact, and + hydrodynamic forces. In general, this difference is small as tension dominates the top node's net force + but slight differences may be observed when comparing the two. - The output list is not case sensitive, however all MoorDyn-F outputs will be printed to the output file in all caps. When searching OpenFAST output channels, users will need to search for MoorDyn channels in all caps. Example: the channel fairten1 would appear in the output file as FAIRTEN1. diff --git a/docs/references.rst b/docs/references.rst index c4f59200..4e2a2ca0 100644 --- a/docs/references.rst +++ b/docs/references.rst @@ -64,6 +64,7 @@ Coupling with WEC-Sim or any Simulink code for wave energy converter simulation: Version 2 ^^^^^^^^^ +.. _version2: Version 2 builds upon the capabilities of Version 1. The theory behind the new features is described in the following references. @@ -117,7 +118,7 @@ Reef3D-MoorDyn Coupling: `Soydan, Ahmet, Widar Weizhi Wang, and Hans Bihs. "An Improved Direct Forcing Immersed Boundary Method With Integrated Mooring Algorithm for Floating Offshore Wind Turbines." American Society of Mechanical Engineers Digital Collection, 2024. `_ -Modeling of Bi-stable Nonlinear Energy Sinks in MoorDyn (most recent description of MoorDyn theory): +Modeling of Bi-stable Nonlinear Energy Sinks in MoorDyn (recent MoorDyn theory included in section 2.1): `Anargyros Michaloliakos, Wei-Ying Wong, Ryan Davies, Malakonda Reddy Lekkala, Matthew Hall, Lei Zuo, Alexander F. Vakakis, "Stabilizing dynamic subsea power cables using Bi-stable nonlinear energy sinks", Ocean Engineering, vol. 334, August 2025. `_ @@ -126,6 +127,11 @@ Syrope model for polyester ropes: `Wei, Zhilong, Harry B. Bingham, and Yanlin Shao. 2026. β€œESOMOOR D5.1: Extended Moordyn Solver and Validation Report”. Technical University of Denmark. `_ +Cross-flow VIV theory, implementation, and validation: + + `Ryan Davies, Malakonda Reddy Lekkala, Anargyros Michaloliakos, Wei-Ying Wong, Matthew Hall, Alexander Vakakis, Lei Zuo, "Time-domain vortex induced vibration modeling of + reference dynamic power cable for the Gulf of Maine", Journal of Fluids and Structures, vol. 146, September 2026. `_ + The Fortran version of MoorDyn is available as a module inside of OpenFAST: https://openfast.readthedocs.io/en/main/ diff --git a/source/Body.hpp b/source/Body.hpp index 72e721a3..3426d0bd 100644 --- a/source/Body.hpp +++ b/source/Body.hpp @@ -63,8 +63,7 @@ class Rod; * which can be extended through the usage of coupled ones. * * In the configuration file the options are: - * - * Name/ID, X0, Y0, Z0, Xcg, Ycg, Zcg, M, V, IX, IY, IZ, CdA-x,y,z Ca-x,y,z + * `Name/ID, Attachment, X0, Y0, Z0, R0, P0, Y0, Mass, (Xcg|Ycg|Zcg), (IX|IY|IZ), Volume, CdA(x|y|z|r|p|yw), Ca(x|y|z)` */ class DECLDIR Body final : public Instance diff --git a/source/Line.cpp b/source/Line.cpp index bcbcc684..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 @@ -340,6 +341,9 @@ Line::setup(int number_in, outfile = outfile_pointer.get(); // make outfile point to the right place channels = channels_in; // copy string of output channels to object + // set flag to store node accelerations if needed for VIV or acceleration output + store_rdd = (Cl > 0) || (channels.find("a") != string::npos); + LOGDBG << " Set up Line " << number << ". " << endl; }; @@ -401,6 +405,14 @@ Line::initialize() << setw(WIDTH) << right << ("Node" + to_string((int)i) + "vz"); } } + // output node accelerations + if (channels.find("a") != string::npos) { + for (unsigned int i = 0; i <= N; i++) { + *outfile << setw(WIDTH) << right << ("Node" + to_string((int)i) + "ax") + << setw(WIDTH) << right << ("Node" + to_string((int)i) + "ay") + << setw(WIDTH) << right << ("Node" + to_string((int)i) + "az"); + } + } // output wave velocities if (channels.find("U") != string::npos) { for (unsigned int i = 0; i <= N; i++) { @@ -488,6 +500,12 @@ Line::initialize() *outfile << setw(WIDTH) << right << "(m/s)"; } + // output node accelerations? + if (channels.find("a") != string::npos) { + for (unsigned int i = 0; i <= 3 * N + 2; i++) + *outfile << setw(WIDTH) << right + << "(m/s^2)"; + } // output wave velocities? if (channels.find("U") != string::npos) { for (unsigned int i = 0; i <= 3 * N + 2; i++) @@ -645,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, @@ -661,6 +687,7 @@ Line::initialize() Xl, Zl, Te); + std::fesetenv(&fenv); if (success >= 0) { @@ -1660,9 +1687,9 @@ Line::getStateDeriv(InstanceStateVarView drdt) drdt.row(i - 1).head<3>() = rd[i]; drdt.row(i - 1).segment<3>(3) = M[i].inverse() * Fnet[i]; - if (Cl > 0) + if (store_rdd) rdd_old[i] = drdt.row(i - 1).segment<3>( - 3); // saving the acceleration for VIV RMS calculation. End + 3); // store node acceleration for VIV and/or output. End // nodes are left at zero, VIV disabled for end nodes } @@ -1742,6 +1769,10 @@ Line::Output(real time) if (channels.find("v") != string::npos) { write_vec_array(rd, N + 1); } + // output node accelerations? + if (channels.find("a") != string::npos) { + write_vec_array(rdd_old, N + 1); + } // output wave velocities? if (channels.find("U") != string::npos) { auto [_z, U, _ud, _pdyn] = waves->getWaveKinLine(lineId); diff --git a/source/Line.hpp b/source/Line.hpp index 4691374c..5b1938a3 100644 --- a/source/Line.hpp +++ b/source/Line.hpp @@ -64,6 +64,14 @@ typedef std::shared_ptr WavesRef; * * The integration time step (moordyn::MoorDyn.dtM0) should be smaller than * this natural period to avoid numerical instabilities + * + * In the configuration file the types options are: + * + * `TypeName, Diam, Mass/m, EA, BA/-zeta, EI, Cd, Ca, CdAx, CaAx, Cl, dF, cF` + * + * and the list options are: + * + * `ID, LineType, AttachA, AttachB, UnstrLen, NumSegs, LineOutputs` */ class DECLDIR Line final : public Instance @@ -366,6 +374,8 @@ class DECLDIR Line final std::vector ydd_rms_old; /// node old accelerations std::vector rdd_old; + /// flag to store node accelerations (for VIV or output) + bool store_rdd; // end conditions /** @brief Types of end points diff --git a/source/Misc.cpp b/source/Misc.cpp index 11a323d8..c6897697 100644 --- a/source/Misc.cpp +++ b/source/Misc.cpp @@ -78,8 +78,10 @@ split(const string& str, const char sep) string token; vector words; while (std::getline(spliter, token, sep)) { - if (token.find("#") == string::npos) { - if (token.size()) // # is a comment in input files + if (token.find("#") == string::npos && + token.find("!") == string::npos && + token.find("%") == string::npos) { + if (token.size()) // !, #, % are comment characters in input files words.push_back(token); } else break; diff --git a/source/Misc.hpp b/source/Misc.hpp index d93013fd..9c9a49b8 100644 --- a/source/Misc.hpp +++ b/source/Misc.hpp @@ -1130,6 +1130,9 @@ struct EnvCond /// a ratio of static to dynamic friction ( = mu_static/mu_dynamic) double StatDynFricScale; + /// ramp time for water kinematics (s). 0 = no ramp + double waveKin_rampT; + /// a global switch for whether to show the units line in the output files /// (1, default), or skip it (0) int WriteUnits; diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp index 5c39c5da..abf2aa78 100644 --- a/source/MoorDyn2.cpp +++ b/source/MoorDyn2.cpp @@ -97,6 +97,7 @@ moordyn::MoorDyn::MoorDyn(const char* infilename, int log_level) , ICTmax(120.0) , ICthresh(0.001) , WaveKinTemp(waves::WAVES_NONE) + , CurrentModeTemp(waves::CURRENTS_NONE) , dtM0((std::numeric_limits::max)()) , cfl(0.5) , dtOut(0.0) @@ -142,6 +143,7 @@ moordyn::MoorDyn::MoorDyn(const char* infilename, int log_level) env->kb = 3.0e6; env->cb = 3.0e5; env->waterKinOptions = waves::WaterKinOptions(); + env->waveKin_rampT = 0.0; // no ramp by default env->WriteUnits = 1; // by default, write units line env->writeLog = 0; // by default, don't write out a log file env->FrictionCoefficient = 0.0; @@ -567,6 +569,12 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) // ------------------ do IC gen -------------------- if (!skip_ic) { + // If wave ramping is active, currents are disabled during IC gen and + // applied only at runtime (alongside the wave grid). Otherwise, + // currents are active during IC gen. + if (env->waveKin_rampT == 0.0) + env->waterKinOptions.currentMode = CurrentModeTemp; + for (unsigned int l = 0; l < LineList.size(); l++) { LineList[l]->IC_gen = true; // turn on IC_gen flag } @@ -602,9 +610,10 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic) } _t_integrator->SetTime(0.0); - // store passed WaveKin value to enable waves in simulation if applicable - // (they're not enabled during IC gen) + // store passed WaveKin and Currents values to enable waves/currents in + // simulation if applicable (they're not enabled during IC gen) env->waterKinOptions.waveMode = WaveKinTemp; + env->waterKinOptions.currentMode = CurrentModeTemp; try { // TODO - figure out how i want to do this better // because this is horrible. the solution is probably to move EnvCond @@ -1743,8 +1752,6 @@ moordyn::MoorDyn::ReadInFile() << "..." << endl << "'" << in_txt[i] << "'" << endl << "invalid output specifier: " << let1 - << ". Type must be oneof L/Line, P/Point, R/Rod, " - "or B/Body" << endl; dummy.OType = -1; continue; @@ -2564,9 +2571,11 @@ moordyn::MoorDyn::readOptionsLine(vector& in_txt, int i) LOGWRN << "Unknown WaveKin option value " << WaveKinTemp << endl; } else if (name == "dtwave") env->waterKinOptions.dtWave = stof(value); + else if (name == "wavekin_rampt") + env->waveKin_rampT = atof(value.c_str()); else if (name == "currents") { auto current_mode = (waves::currents_settings)stoi(value); - env->waterKinOptions.currentMode = current_mode; + CurrentModeTemp = current_mode; if ((current_mode < waves::CURRENTS_NONE) || (current_mode > waves::CURRENTS_4D)) LOGWRN << "Unknown Currents option value " << current_mode << endl; diff --git a/source/MoorDyn2.hpp b/source/MoorDyn2.hpp index 325e6c9c..d9841d8e 100644 --- a/source/MoorDyn2.hpp +++ b/source/MoorDyn2.hpp @@ -565,6 +565,9 @@ class MoorDyn final : public io::IO // temporary wave kinematics flag used to store input value while keeping // env.WaveKin=0 for IC gen moordyn::waves::waves_settings WaveKinTemp; + // temporary currents flag used to store input value; currents may be + // disabled during IC gen when wave ramping is active + moordyn::waves::currents_settings CurrentModeTemp; /// (s) desired mooring line model maximum time step real dtM0; /// desired mooring line model maximum CFL factor diff --git a/source/Point.hpp b/source/Point.hpp index d984ecb6..6fcd9e0d 100644 --- a/source/Point.hpp +++ b/source/Point.hpp @@ -62,6 +62,10 @@ typedef std::shared_ptr WavesRef; * unconnected termination point of a Line, which could have a clump * weight or float via the point's mass and volume parameters * - Coupled: The point position and velocity is externally imposed + * + * In the configuration file the options are: + * + * `Name/ID, Attachment, X, Y, Z, Mass, Volume, CdA, Ca` */ class DECLDIR Point final : public Instance diff --git a/source/Rod.cpp b/source/Rod.cpp index d1bc2344..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(); } @@ -317,6 +318,13 @@ Rod::openoutput() << setw(WIDTH) << right << ("Node" + to_string((int)i) + "Pdz"); } + // output water velocities + if (channels.find("U") != string::npos) { + for (unsigned int i = 0; i <= N; i++) + *outfile << setw(WIDTH) << right << ("Node" + to_string((int)i) + "Ux") + << setw(WIDTH) << right << ("Node" + to_string((int)i) + "Uy") + << setw(WIDTH) << right << ("Node" + to_string((int)i) + "Uz"); + } // output bottom contact force if (channels.find("b") != string::npos) { for (unsigned int i = 0; i <= N; i++) @@ -384,6 +392,11 @@ Rod::openoutput() for (unsigned int i = 0; i <= 3 * N + 2; i++) *outfile << setw(WIDTH) << right << "(Pa)"; } + // output wave velocities + if (channels.find("U") != string::npos) { + for (unsigned int i = 0; i <= 3 * N + 2; i++) + *outfile << setw(WIDTH) << right << "(m/s)"; + } // output bottom contact force if (channels.find("b") != string::npos) { for (unsigned int i = 0; i <= 3 * N + 2; i++) @@ -1556,6 +1569,10 @@ auto write_val = [&](real val) { if (channels.find("P") != string::npos) { write_vec_array(Pd); // dynamic pressure } + if (channels.find("U") != string::npos) { + auto [_z, U, _ud, _pdyn] = waves->getWaveKinRod(rodId); + write_vec_array(U); // wave fluid velocities + } if (channels.find("b") != string::npos) { write_vec_array(B); // seabed contact (bottom contact forces) } diff --git a/source/Rod.hpp b/source/Rod.hpp index b660bcb3..9a62ef65 100644 --- a/source/Rod.hpp +++ b/source/Rod.hpp @@ -58,6 +58,14 @@ class Line; * * Each end point of the rod can be fixed or pinned to another object, let free * or control it externally + * + * In the configuration file the types options are: + * + * `TypeName, Diam, Mass/m, Cd, Ca, CdEnd, CaEnd` + * + * and the list options are: + * + * `ID, RodType, Attachment, Xa, Ya, Za, Xb, Yb, Zb, NumSegs, RodOutputs` */ class DECLDIR Rod final : public Instance 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 diff --git a/source/Waves.cpp b/source/Waves.cpp index 20a208e3..8d0c4219 100644 --- a/source/Waves.cpp +++ b/source/Waves.cpp @@ -822,10 +822,9 @@ Waves::updateWaves() bodies.U[id][0] = waveKin.bodies.U[id][0] + curr_U; bodies.Ud[id][0] = waveKin.bodies.Ud[id][0] + curr_Ud; } - return; } // if there are both waves and currents, then we calculate and sum - if (waveKinematics && currentKinematics) { + else if (waveKinematics && currentKinematics) { kinematicsForAllNodes( nodeKin, [&](vec pos, vec& U, vec& Ud, real& zeta, real& pdyn) { vec3 wave_U{}, wave_Ud{}; @@ -845,10 +844,9 @@ Waves::updateWaves() U = wave_U + curr_U; Ud = wave_Ud + curr_Ud; }); - return; } // if there are just waves then we just do wave calculations - if (waveKinematics) { + else if (waveKinematics) { kinematicsForAllNodes( nodeKin, [&](vec pos, vec& U, vec& Ud, real& zeta, real& pdyn) { waveKinematics->getWaveKin(pos, @@ -859,17 +857,32 @@ Waves::updateWaves() &Ud, &pdyn); }); - return; } // if there are just currents then we just do current calculations - if (currentKinematics) { + else if (currentKinematics) { kinematicsForAllNodes( nodeKin, [&](vec pos, vec& U, vec& Ud, real& zeta, real& pdyn) { currentKinematics->getCurrentKin( pos, _t_integrator->GetTime(), floorProvider, &U, &Ud); }); - return; } + + // Apply ramp time for water kinematics if enabled + if (env->waveKin_rampT > 0.0) { + const real t = _t_integrator->GetTime(); + if (t < env->waveKin_rampT) { + const real ramp = t / env->waveKin_rampT; + kinematicsForAllNodes( + nodeKin, + [&](vec /*pos*/, vec& U, vec& Ud, real& zeta, real& pdyn) { + U *= ramp; + Ud *= ramp; + zeta *= ramp; + pdyn *= ramp; + }); + } + } + return; } } // ::moordyn 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/Mooring/viv/plot_viv_timeseries.py b/tests/Mooring/viv/plot_viv_timeseries.py new file mode 100644 index 00000000..22b24627 --- /dev/null +++ b/tests/Mooring/viv/plot_viv_timeseries.py @@ -0,0 +1,351 @@ +# @author: ryandavies19, Copilot (GPT-5.4 mini) + +"""Plot the VIV test timeseries and the expected period check. + +This script reads the generated test output in build/tests/Mooring/viv/ +and reproduces the peak-detection logic used by tests/viv.cpp. +""" + +from __future__ import annotations + +import argparse +from pathlib import Path + +import matplotlib.pyplot as plt +import numpy as np + + +EXPECTED_FREQUENCY_HZ = 1.04 +EXPECTED_PERIOD_S = 1.0 / (2.0 * EXPECTED_FREQUENCY_HZ) +PERIOD_TOLERANCE_S = 1.0e-2 +CHECK_START_TIME_S = 12.0 +MIN_PEAK_AMP = 0.5 # keep in sync with MIN_PEAK_AMP in tests/viv.cpp + + +def load_timeseries(path: Path) -> tuple[np.ndarray, np.ndarray]: + """Load the first two numeric columns from a MoorDyn output file.""" + + rows: list[tuple[float, float]] = [] + with path.open("r", encoding="utf-8", errors="surrogateescape") as stream: + for line in stream: + stripped = line.strip() + if not stripped: + continue + parts = stripped.split() + if len(parts) < 2: + continue + try: + time = float(parts[0]) + value = float(parts[1]) + except ValueError: + continue + rows.append((time, value)) + + if not rows: + raise ValueError(f"No numeric time series data found in {path}. Check that FairTen1 is an output in viv.txt") + + data = np.asarray(rows, dtype=float) + return data[:, 0], data[:, 1] + + +def detect_test_peaks(time: np.ndarray, tension: np.ndarray) -> tuple[np.ndarray, np.ndarray]: + """Reproduce the peak tracker used in tests/viv.cpp. + + The test tracks changes in the fairlead tension after 10 seconds and + records a new peak when the absolute difference from the last peak stops + growing. + """ + + ten_peaks = [0.0] + peak_times: list[float] = [] + d_peak = 0.0 + + for t, ten in zip(time, tension): + if t > CHECK_START_TIME_S: + delta = abs(ten - ten_peaks[-1]) + if delta > d_peak: + d_peak = delta + else: + if d_peak >= MIN_PEAK_AMP: + ten_peaks.append(float(ten)) + peak_times.append(float(t)) + d_peak = 0.0 + else: + # seed the reference from the latest pre-window sample so the swing + # is measured against the signal, not against the initial 0.0 + ten_peaks[-1] = float(ten) + d_peak = 0.0 + + return np.asarray(peak_times), np.asarray(ten_peaks[1:]) + + +def compute_fft_spectrum( + time: np.ndarray, + values: np.ndarray, + start_time: float, + end_time: float, +) -> tuple[np.ndarray, np.ndarray]: + """Compute a one-sided amplitude spectrum for the portion within a time window.""" + + mask = (time >= start_time) & (time <= end_time) + segment_time = time[mask] + segment_values = values[mask] + if segment_time.size < 4: + raise ValueError("Not enough samples after the FFT start time") + + # Set up time grid and interpolate, then center. + dt = float(np.median(np.diff(segment_time))) + uniform_time = np.arange(segment_time[0], segment_time[-1] + 0.5 * dt, dt) + uniform_values = np.interp(uniform_time, segment_time, segment_values) + centered_values = uniform_values - uniform_values.mean() + + # Taper the tails and apply Hanning window + tail_length = centered_values.size // 10 # Or something else + window_mask = np.zeros(centered_values.size, dtype=bool) + window_mask[:tail_length] = True + window_mask[-tail_length:] = True + windowed_values = np.copy(centered_values) + windowed_values[window_mask] *= np.hanning(2 * tail_length) + + # FFT and normalize + spectrum = np.fft.rfft(windowed_values) + frequencies = np.fft.rfftfreq(windowed_values.size, d=dt) + window = np.ones(centered_values.size) + window[window_mask] = np.hanning(2 * tail_length) + normalization = window.sum() / 2.0 + amplitudes = np.abs(spectrum) / normalization + return frequencies, amplitudes + + +def plot_fft_comparison( + time: np.ndarray, + tension: np.ndarray, + overlay_time: np.ndarray | None, + overlay_tension: np.ndarray | None, + fft_end_time: float, + output_path: Path, +) -> plt.Figure: + """Plot the FFT comparison for the post-startup portion of the simulation.""" + + fig, ax = plt.subplots(figsize=(12, 6)) + + frequencies, amplitudes = compute_fft_spectrum(time, tension, CHECK_START_TIME_S, fft_end_time) + ax.plot(frequencies, amplitudes, color="#1f4e79", lw=1.4, label="Fairlead tension") + + if overlay_time is not None and overlay_tension is not None: + overlay_frequencies, overlay_amplitudes = compute_fft_spectrum( + overlay_time, + overlay_tension, + CHECK_START_TIME_S, + fft_end_time, + ) + ax.plot( + overlay_frequencies, + overlay_amplitudes, + color="#7a3ea1", + lw=1.2, + alpha=0.85, + label="OpenFAST r-test", + ) + + ax.set_xlim(left=0.0) + ax.set_xlabel("Frequency [Hz]") + ax.set_ylabel("Amplitude") + ax.set_title(f"FFT comparison from {CHECK_START_TIME_S:.0f} to {fft_end_time:.3f} seconds") + ax.grid(True, alpha=0.25) + ax.legend(loc="best") + fig.tight_layout() + fig.savefig(output_path, dpi=200, bbox_inches="tight") + return fig + + +def plot_viv_timeseries( + time: np.ndarray, + tension: np.ndarray, + peak_times: np.ndarray, + overlay_time: np.ndarray | None, + overlay_tension: np.ndarray | None, + output_path: Path, + show: bool, +) -> plt.Figure: + """Plot the time series and the period check against the expected value.""" + + fig, (ax_tension, ax_period) = plt.subplots( + 2, + 1, + figsize=(12, 8), + sharex=True, + gridspec_kw={"height_ratios": [2.2, 1.0]}, + ) + + ax_tension.plot(time, tension, color="#1f4e79", lw=1.25, label="Fairlead tension") + if overlay_time is not None and overlay_tension is not None: + ax_tension.plot( + overlay_time, + overlay_tension, + color="#7a3ea1", + lw=1.0, + alpha=0.8, + label="OpenFAST r-test", + ) + if peak_times.size: + peak_indices = np.searchsorted(time, peak_times) + peak_indices = np.clip(peak_indices, 0, len(time) - 1) + ax_tension.scatter( + time[peak_indices], + tension[peak_indices], + s=18, + color="#c23b22", + zorder=3, + label="Detected peaks", + ) + + ax_tension.axvline(CHECK_START_TIME_S, color="#666666", ls="--", lw=1.0, alpha=0.7) + ax_tension.set_ylabel("Fairlead tension [N]") + ax_tension.set_title("VIV test output and expected period") + ax_tension.grid(True, alpha=0.25) + ax_tension.legend(loc="best") + + if peak_times.size >= 3: + observed_periods = peak_times[2:] - peak_times[:-2] + period_times = peak_times[2:] + ax_period.plot(period_times, observed_periods, "o", ms=4.5, color="#1f4e79", label="Observed period") + out_of_tolerance = np.abs(observed_periods - EXPECTED_PERIOD_S) > PERIOD_TOLERANCE_S + for period_time, period_value in zip(period_times[out_of_tolerance], observed_periods[out_of_tolerance]): + ax_period.axvline( + period_time, + color="#c23b22", + ls=":", + lw=1.0, + alpha=0.85, + ) + ax_period.annotate( + f"t={period_time:.3f} s", + xy=(period_time, 0.97), + xycoords=(ax_period.transData, ax_period.transAxes), + xytext=(0, 0), + textcoords="offset points", + rotation=90, + fontsize=8, + color="#c23b22", + ha="center", + va="top", + ) + else: + observed_periods = np.asarray([]) + ax_period.text(0.5, 0.5, "Not enough peaks detected to compute a period", ha="center", va="center", transform=ax_period.transAxes) + + ax_period.axhline(EXPECTED_PERIOD_S, color="#c23b22", lw=1.5, label=f"Expected period = {EXPECTED_PERIOD_S:.6f} s") + ax_period.axhspan( + EXPECTED_PERIOD_S - PERIOD_TOLERANCE_S, + EXPECTED_PERIOD_S + PERIOD_TOLERANCE_S, + color="#c23b22", + alpha=0.12, + label=f"Test tolerance Β±{PERIOD_TOLERANCE_S:.2e} s", + ) + ax_period.set_xlabel("Time [s]") + ax_period.set_ylabel("Period [s]") + ax_period.grid(True, alpha=0.25) + ax_period.legend(loc="best") + + if observed_periods.size: + mean_period = observed_periods.mean() + ax_period.set_title(f"Observed mean period = {mean_period:.6f} s") + + fig.tight_layout() + fig.savefig(output_path, dpi=200, bbox_inches="tight") + return fig + + +def main() -> int: + parser = argparse.ArgumentParser( + description="Plot the VIV test fairlead tension and the expected period check." + ) + parser.add_argument( + "input", + nargs="?", + default="../../../build/tests/Mooring/viv/viv.out", + help="Path to the viv.out file produced by the test run.", + ) + parser.add_argument( + "--output", + default=None, + help="Optional output image path. Defaults to the input path with .png appended.", + ) + parser.add_argument( + "--show", + action="store_true", + help="Display the plot interactively after saving it.", + ) + parser.add_argument( + "--overlay", + default="", + help="Optional second output file to overlay on the main timeseries plot. Use an empty string to disable.", + ) + parser.add_argument( + "--fft-output", + default=None, + help="Optional output image path for the FFT comparison. Defaults to the main output path with _fft.png appended.", + ) + args = parser.parse_args() + + input_path = Path(args.input).expanduser().resolve() + if args.output is None: + output_path = input_path.with_suffix(input_path.suffix + ".png") + else: + output_path = Path(args.output).expanduser().resolve() + + if args.fft_output is None: + fft_output_path = output_path.with_name(f"{output_path.stem}_fft{output_path.suffix}") + else: + fft_output_path = Path(args.fft_output).expanduser().resolve() + + time, tension = load_timeseries(input_path) + peak_times, _ = detect_test_peaks(time, tension) + fft_end_time = float(time.max()) + + overlay_time = None + overlay_tension = None + if args.overlay: + overlay_path = Path(args.overlay).expanduser().resolve() + overlay_time, overlay_tension = load_timeseries(overlay_path) + + print(f"Loaded {len(time)} samples from {input_path}") + print(f"Expected period: {EXPECTED_PERIOD_S:.6f} s") + if overlay_time is not None: + print(f"Loaded {len(overlay_time)} samples from {overlay_path}") + if peak_times.size >= 3: + observed_periods = peak_times[2:] - peak_times[:-2] + print(f"Detected {len(peak_times)} peaks after {CHECK_START_TIME_S:.1f} s") + print(f"Mean observed period: {observed_periods.mean():.6f} s") + else: + print("Not enough peaks detected to compute an observed period.") + + viv_fig = plot_viv_timeseries( + time, + tension, + peak_times, + overlay_time, + overlay_tension, + output_path, + args.show, + ) + fft_fig = plot_fft_comparison( + time, + tension, + overlay_time, + overlay_tension, + fft_end_time, + fft_output_path, + ) + if args.show: + plt.show() + plt.close(viv_fig) + plt.close(fft_fig) + print(f"Saved plot to {output_path}") + print(f"Saved FFT plot to {fft_output_path}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) \ No newline at end of file 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); +} diff --git a/tests/viv.cpp b/tests/viv.cpp index a325f76b..246afe4d 100644 --- a/tests/viv.cpp +++ b/tests/viv.cpp @@ -28,8 +28,8 @@ * POSSIBILITY OF SUCH DAMAGE. */ -/** @file pendulum.cpp - * Simple case of a pendulum motion +/** @file viv.cpp + * Simple case of a VIV simulation */ // Visual studio still uses this @@ -46,9 +46,14 @@ // NOTE: this is largely built on the pendulum test framework #define TOL \ - 1.0e-2 // absolute tolerance. In setting this up, max error was 0.008069, + 1.2e-2 // absolute tolerance. In setting this up, max error was 0.011 (OS dependent), // but this can vary slightly with every simulation +#define MIN_PEAK_AMP \ + 0.5 // minimum delta T (N) from the last peak before it gets counted as + // a peak. Tension oscillation is ~3.6 N pk-pk while numerical noise is + // <0.1 N, so this rejects duplicate detections around peaks. + bool compare(double v1, double v2, double tol) { @@ -96,13 +101,13 @@ TEST_CASE("VIV frequency simulation and check") double d_peak = 0.0; // last peak tracker. For checking if there is a change in the peak - while (t < 15) { // run for 15 seconds + while (t < 18) { // run for 18 seconds REQUIRE(MoorDyn_Step(system, NULL, NULL, NULL, &t, &dt) == MOORDYN_SUCCESS); REQUIRE(MoorDyn_GetLineFairTen(line1, &ten) == MOORDYN_SUCCESS); - if (t > 10) { // only check period after 10 seconds. Before then, there + if (t > 12) { // only check period after 12 seconds. Before then, there // isn't complete lock in // check for peaks in tension. These are used to determine the // period @@ -110,11 +115,17 @@ TEST_CASE("VIV frequency simulation and check") if (d > d_peak) d_peak = d; else { - ten_peaks.push_back(ten); - peak_times.push_back(t); + if (d_peak >= MIN_PEAK_AMP) { + ten_peaks.push_back(ten); + peak_times.push_back(t); + } d_peak = 0.0; } } + else { + // Save the current tension for comparison on the next timesteps + ten_peaks.push_back(ten); + } } REQUIRE(peak_times.size() != 0);