feat: OpenFAST synchronization, docs updates, and VIV test fix

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 <https://github.com/WEC-Sim/WEC-Sim/pull/1212>`_. 
 The original coupling was with MoorDyn v1 and the the publication can be found in the 
-:ref:`theory section <theory>`. The WEC-Sim source code can be found `here <https://github.com/WEC-Sim/WEC-Sim>`_.
+:ref:`theory section <theory>`. The WEC-Sim source code can be found `here <https://github.com/WEC-Sim/WEC-Sim>`_. 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 <https://github.com/FloatingArrayDesign/MoorDyn/tree/master/example>`_.
 
 DualSPHysics
 ^^^^^^^^^^^^

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 <https://github.com/FloatingArrayDesign/MoorDyn/tree/master/tests/Mooring>`_ (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 <waterkinematics>`.
  - 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 <waterkinematics>`.
@@ -801,17 +804,22 @@ 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, 
   the best approach is to attach a rod to a zero-mass, zero-volume pinned body and output the body 
   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.

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. <https://doi.org/10.1115/1.4067117>`_
 
-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. <https://doi.org/10.1016/j.oceaneng.2025.121613>`_
@@ -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. <https://doi.org/10.11583/DTU.31408806>`_
 
+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. <https://doi.org/10.1016/j.jfluidstructs.2026.104651>`_
+
 The Fortran version of MoorDyn is available as a module inside of OpenFAST:
 
   https://openfast.readthedocs.io/en/main/

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

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);

source/Line.hpp

@@ -64,6 +64,14 @@ typedef std::shared_ptr<Waves> 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<moordyn::real> ydd_rms_old;
 	/// node old accelerations
 	std::vector<vec> rdd_old;
+	/// flag to store node accelerations (for VIV or output)
+	bool store_rdd;
 
 	// end conditions
 	/** @brief Types of end points

source/Misc.cpp

@@ -78,8 +78,10 @@ split(const string& str, const char sep)
 	string token;
 	vector<string> 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;

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;

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<real>::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<string>& 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;

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

source/Point.hpp

@@ -62,6 +62,10 @@ typedef std::shared_ptr<Waves> 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

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)
 	}

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