Repository: FloatingArrayDesign/MoorDyn (language: python)

Issue:
Small typo in code Line::setup
Hi, 
There is a small mistake in Line::setup. I think you forgot to multiply with one more d when calculating volume. I don't think it causes any problem, but still. :)
for (int i=0; i<N; i++)	
	{	l[i] = UnstrLen/double(N);	// distribute line length evenly over segments
		V[i] = l[i]*0.25*pi*d;
	}


Discussion / hints:
Thanks for noting this! I've corrected it in the latest commit, ea5ccf8c050f359cdd13a0a8b6f5cc5f3158405c.

Fix patch:
diff --git a/source/Connection.hpp b/source/Connection.hpp
index 4f0311bc..deaaee7d 100644
--- a/source/Connection.hpp
+++ b/source/Connection.hpp
@@ -138,6 +138,10 @@ class Connection final : public io::IO
 	/// total force on node
 	vec Fnet;
 
+	// time
+	/// simulation time
+	moordyn::real t;
+
 	/// node mass + added mass matrices
 	mat M;
 
@@ -184,6 +188,11 @@ class Connection final : public io::IO
 	/// Connection type
 	types type;
 
+	/** @brief Set the line  simulation time
+	 * @param time Simulation time
+	 */
+	inline void setTime(real time) { t = time; }
+
 	/** @brief Setup the connection
 	 *
 	 * Always call this function after the construtor
diff --git a/source/Line.cpp b/source/Line.cpp
index 976344bc..1360f3a0 100644
--- a/source/Line.cpp
+++ b/source/Line.cpp
@@ -455,51 +455,58 @@ Line::initialize()
 
 	// if conditions are ideal, try to calculate initial line profile using
 	// catenary routine (from FAST v.7)
-	real waterDepth = getWaterDepth(r[0][0], r[0][1]);
-	if (r[0][2] == waterDepth) {
-		real XF = dir(Eigen::seqN(0, 2)).norm(); // horizontal spread
+
+	real XF = dir(Eigen::seqN(0, 2)).norm(); // horizontal spread
+	if (XF > 0.0) {
+
 		real ZF = dir[2];
 		real LW = ((rho - env->rho_w) * A) * env->g;
 		real CB = 0.;
 		real Tol = 0.00001;
 
-		if ((XF > 0.0) && (ZF > 0.0)) {
-			// locations of line nodes along line length - evenly distributed
-			// here
-			vector<real> snodes(N + 1, 0.0);
-			for (unsigned int i = 1; i <= N; i++)
-				snodes[i] = snodes[i - 1] + l[i - 1];
-			// double check to ensure the last node does not surpass the line
-			// length
-			snodes[N] = UnstrLen;
-
-			// output variables
-			real HF, VF, HA, VA, COSPhi, SINPhi;
-			vector<real> Xl(N + 1, 0.0);
-			vector<real> Zl(N + 1, 0.0);
-			vector<real> Te(N + 1, 0.0);
-
-			COSPhi = (r[N][0] - r[0][0]) / XF;
-			SINPhi = (r[N][1] - r[0][1]) / XF;
-
-			int success = Catenary(XF,
-			                       ZF,
-			                       UnstrLen,
-			                       E * A,
-			                       LW,
-			                       CB,
-			                       Tol,
-			                       &HF,
-			                       &VF,
-			                       &HA,
-			                       &VA,
-			                       N,
-			                       snodes,
-			                       Xl,
-			                       Zl,
-			                       Te);
-
-			if (success >= 0) {
+		// locations of line nodes along line length - evenly distributed
+		// here
+		vector<real> snodes(N + 1, 0.0);
+		for (unsigned int i = 1; i <= N; i++)
+			snodes[i] = snodes[i - 1] + l[i - 1];
+		// double check to ensure the last node does not surpass the line
+		// length
+		snodes[N] = UnstrLen;
+
+		// output variables
+		real HF, VF, HA, VA, COSPhi, SINPhi;
+		vector<real> Xl(N + 1, 0.0);
+		vector<real> Zl(N + 1, 0.0);
+		vector<real> Te(N + 1, 0.0);
+
+		COSPhi = (r[N][0] - r[0][0]) / XF;
+		SINPhi = (r[N][1] - r[0][1]) / XF;
+
+		int success = Catenary(XF,
+		                       ZF,
+		                       UnstrLen,
+		                       E * A,
+		                       LW,
+		                       CB,
+		                       Tol,
+		                       &HF,
+		                       &VF,
+		                       &HA,
+		                       &VA,
+		                       N,
+		                       snodes,
+		                       Xl,
+		                       Zl,
+		                       Te);
+
+		if (success >= 0) {
+			// It might happens that the output solution does not respect the
+			// queried final point. See the pendulum example
+			if (abs(Zl[N] - ZF) > Tol) {
+				LOGWRN << "Wrong catenary initial profile for Line, intializing as linear " << number
+				       << endl;
+			}
+			else {
 				// the catenary solve is successful, update the node positions
 				LOGDBG << "Catenary initial profile available for Line "
 				       << number << endl;
@@ -507,12 +514,18 @@ Line::initialize()
 					vec l(Xl[i] * COSPhi, Xl[i] * SINPhi, Zl[i]);
 					r[i] = r[0] + l;
 				}
-			} else {
-				LOGWRN << "Catenary initial profile failed for Line " << number
-				       << endl;
 			}
 		}
+		else {
+			LOGWRN << "Catenary initial profile failed for Line " << number
+					<< endl;
+		}
+	}
+	else {
+		LOGDBG << "Vertical initial profile for Line " << number
+		       << endl;
 	}
+	
 
 	LOGMSG << "Initialized Line " << number << endl;
 
@@ -1115,9 +1128,9 @@ Line::getStateDeriv()
 			Ap[i] = env->rho_w * (1. + Can) * 0.5 * (V[i] + V[i - 1]) * ap;
 		// tangential Froude-Krylov force
 		if (i == 0)
-			Aq[i] = 0.5 * env->rho_w * (1. + Cat) * 0.5 * (V[i]) * aq;
+			Aq[i] = env->rho_w * (1. + Cat) * 0.5 * (V[i]) * aq;
 		else if (i == N)
-			Aq[i] = 0.5 * env->rho_w * (1. + Cat) * 0.5 * (V[i - 1]) * aq;
+			Aq[i] = env->rho_w * (1. + Cat) * 0.5 * (V[i - 1]) * aq;
 		else
 			Aq[i] = env->rho_w * (1. + Cat) * 0.5 * (V[i] + V[i - 1]) * aq;
 
diff --git a/source/Time.cpp b/source/Time.cpp
index b2a44f58..4bf11f0c 100644
--- a/source/Time.cpp
+++ b/source/Time.cpp
@@ -198,9 +198,8 @@ RK2Scheme::Step(real& dt)
 	r[1] = r[0] + rd[0] * (0.5 * dt);
 	Update(0.5 * dt, 1);
 	// And so we can compute the new derivative and apply it
-	CalcStateDeriv(0);
-	r[0] = r[0] + rd[0] * dt;
-
+	CalcStateDeriv(1);
+	r[0] = r[0] + rd[1] * dt;
 	t += 0.5 * dt;
 	Update(dt, 0);
 	TimeSchemeBase::Step(dt);
diff --git a/source/Time.hpp b/source/Time.hpp
index b68bf3b8..63e4eb27 100644
--- a/source/Time.hpp
+++ b/source/Time.hpp
@@ -751,7 +751,7 @@ class HeunScheme : public TimeSchemeBase<1, 2>
  *
  * This was the traditionally applied time scheme in MoorDyn v1
  */
-class RK2Scheme : public TimeSchemeBase<2, 1>
+class RK2Scheme : public TimeSchemeBase<2, 2>
 {
   public:
 	/** @brief Costructor


Test patch:
diff --git a/tests/pendulum.cpp b/tests/pendulum.cpp
index 2cc08fdd..496cce0a 100644
--- a/tests/pendulum.cpp
+++ b/tests/pendulum.cpp
@@ -103,6 +103,15 @@ pendulum()
 		cerr << "Failure during the mooring initialization: " << err << endl;
 		return false;
 	}
+	err = MoorDyn_GetConnectPos(conn, pos);
+	if (err != MOORDYN_SUCCESS) {
+		cerr << "Failure getting the initial mass position: " << err << endl;
+		return false;
+	}
+	const double x = pos[0];
+	const double z = pos[2];
+	const double l = sqrt(x * x + z * z);
+	CHECK_VALUE("L0", l0, l, 0.01 * l0, 0);
 
 	const double w = sqrt(9.81 / l0);
 	const double T = 2.0 * M_PI / w;