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..cbed0f98 100644 --- a/source/Line.cpp +++ b/source/Line.cpp @@ -340,6 +340,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 +404,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 +499,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++) @@ -1660,9 +1677,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 +1759,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..a499ef96 100644 --- a/source/Rod.cpp +++ b/source/Rod.cpp @@ -317,6 +317,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 +391,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 +1568,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/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/Mooring/viv/plot_viv_timeseries.py b/tests/Mooring/viv/plot_viv_timeseries.py new file mode 100644 index 00000000..ce00b5e4 --- /dev/null +++ b/tests/Mooring/viv/plot_viv_timeseries.py @@ -0,0 +1,340 @@ +# @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") + + 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() + window = np.hanning(centered_values.size) + windowed_values = centered_values * window + spectrum = np.fft.rfft(windowed_values) + frequencies = np.fft.rfftfreq(windowed_values.size, d=dt) + 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/viv.cpp b/tests/viv.cpp index a325f76b..c8c3ebb1 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 @@ -49,6 +49,11 @@ 1.0e-2 // absolute tolerance. In setting this up, max error was 0.008069, // 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);