Additional line performance improvements

Author: AlexWKinley

bench/MDBench.cpp

@@ -21,8 +21,10 @@ LineGetStateDeriv(benchmark::State& state, std::string input_file)
 
 	auto line = system.GetLines().front();
 
+	std::vector<moordyn::vec> vel{};
+	std::vector<moordyn::vec> acc{};
 	for (auto _ : state) {
-		line->getStateDeriv();
+		line->getStateDeriv(vel, acc);
 	}
 }
 BENCHMARK_CAPTURE(LineGetStateDeriv,

source/Line.cpp

@@ -666,18 +666,16 @@ Line::setPin(std::vector<real> p)
 }
 
 void
-Line::setState(std::vector<vec> pos, std::vector<vec> vel)
+Line::setState(const std::vector<vec>& pos, const std::vector<vec>& vel)
 {
 	if ((pos.size() != N - 1) || (vel.size() != N - 1)) {
 		LOGERR << "Invalid input size" << endl;
 		throw moordyn::invalid_value_error("Invalid input size");
 	}
 
 	// set interior node positions and velocities based on state vector
-	for (unsigned int i = 1; i < N; i++) {
-		r[i] = pos[i - 1];
-		rd[i] = vel[i - 1];
-	}
+	std::copy(pos.begin(), pos.end(), r.begin() + 1);
+	std::copy(vel.begin(), vel.end(), rd.begin() + 1);
 }
 
 void
@@ -769,8 +767,8 @@ Line::getEndSegmentMoment(EndPoints end_point, EndPoints rod_end_point) const
 	return qEnd * EIEnd / dlEnd;
 }
 
-std::pair<std::vector<vec>, std::vector<vec>>
-Line::getStateDeriv()
+void
+Line::getStateDeriv(std::vector<vec>& vel, std::vector<vec>& acc)
 {
 	// NOTE:
 	// Jose Luis Cercos-Pita: This is by far the most consuming function of the
@@ -802,6 +800,22 @@ Line::getStateDeriv()
 		ldstr[i] = qs[i].dot(rd[i + 1] - rd[i]); // strain rate of segment
 
 		// V[i] = A * l[i]; // volume attributed to segment
+
+		if (nEApoints)
+			E = getNonlinearE(lstr[i], l[i]);
+
+		if (lstr[i] > l[i]) {
+			T[i] = E * A * (lstr[i] - l[i]) / l[i] * qs[i];
+		} else {
+			// cable can't "push" ...
+			// or can it, if bending stiffness is nonzero? <<<<<<<<<
+			T[i] = vec::Zero();
+		}
+
+		// line internal damping force
+		if (nCpoints)
+			c = getNonlinearC(ldstr[i], l[i]);
+		Td[i] = c * A * ldstr[i] / l[i] * qs[i];
 	}
 
 	// calculate unit tangent vectors (q) for each internal node. note: I've
@@ -876,51 +890,9 @@ Line::getStateDeriv()
 	}
 	//============================================================================================
 
-	// calculate mass matrix
-	for (unsigned int i = 0; i <= N; i++) {
-		real m_i; // node mass
-		real v_i; // node submerged volume
-
-		if (i == 0) {
-			m_i = pi / 8. * d * d * l[0] * rho;
-			v_i = 0.5 * F[0] * V[0];
-		} else if (i == N) {
-			m_i = pi / 8. * d * d * l[N - 1] * rho;
-			v_i = 0.5 * F[i - 1] * V[i - 1];
-		} else {
-			m_i = pi / 8. * (d * d * rho * (l[i] + l[i - 1]));
-			v_i = 0.5 * (F[i - 1] * V[i - 1] + F[i] * V[i]);
-		}
-
-		// Make node mass matrix
-		const mat I = mat::Identity();
-		const mat Q = q[i] * q[i].transpose();
-		M[i] = m_i * I + env->rho_w * v_i * (Can * (I - Q) + Cat * Q);
-	}
-
 	// ============  CALCULATE FORCES ON EACH NODE
 	// ===============================
 
-	// loop through the segments
-	for (unsigned int i = 0; i < N; i++) {
-		// line tension
-		if (nEApoints > 0)
-			E = getNonlinearE(lstr[i], l[i]);
-
-		if (lstr[i] / l[i] > 1.0) {
-			T[i] = E * A * (lstr[i] - l[i]) / l[i] * qs[i];
-		} else {
-			// cable can't "push" ...
-			// or can it, if bending stiffness is nonzero? <<<<<<<<<
-			T[i] = vec::Zero();
-		}
-
-		// line internal damping force
-		if (nCpoints > 0)
-			c = getNonlinearC(ldstr[i], l[i]);
-		Td[i] = c * A * ldstr[i] / l[i] * qs[i];
-	}
-
 	// Bending loads
 	// first zero out the forces from last run
 	for (unsigned int i = 0; i <= N; i++) {
@@ -1092,95 +1064,76 @@ Line::getStateDeriv()
 		}
 	}
 
+	// This mostly just simplifies the code, but also sometimes helps the
+	// compiler avoid repeated pointer lookup.
+	const real rho_w = env->rho_w;
+	const real g = env->g;
+
 	// loop through the nodes
 	for (unsigned int i = 0; i <= N; i++) {
-		W[i][0] = W[i][1] = 0.0;
+		// Nodes are considered to have a line length of half the length of the
+		// line on either side.
+		const real length_left = i == 0 ? 0.0 : l[i - 1];
+		const real length_right = i == N ? 0.0 : l[i];
+		const real length = 0.5 * (length_left + length_right);
+		const real submergence_left = i == 0 ? 0.0 : F[i - 1];
+		const real submergence_right = i == N ? 0.0 : F[i];
+		// Submerged length is half the submerged length of the neighboring
+		// segments
+		const real submerged_length = 0.5 * (length_left * submergence_left +
+		                                     length_right * submergence_right);
+
+		const real submerged_volume = submerged_length * A;
+
+		const real m_i = A * length * rho; // node mass
+		const real v_i = submerged_volume; // node submerged volume
+
+		// Make node mass matrix
+		const mat I = mat::Identity();
+		const mat Q = q[i] * q[i].transpose();
+		// calculate mass matrix
+		M[i] = m_i * I + rho_w * v_i * (Can * (I - Q) + Cat * Q);
+
 		// submerged weight (including buoyancy)
-		if (i == 0)
-			W[i][2] = 0.5 * A * (l[i] * (rho - F[i] * env->rho_w)) * (-env->g);
-		else if (i == N)
-			W[i][2] = 0.5 * A * (l[i - 1] * (rho - F[i - 1] * env->rho_w)) *
-			          (-env->g); // missing the "W[i][2] =" previously!
-		else
-			W[i][2] = 0.5 * A *
-			          (l[i] * (rho - F[i] * env->rho_w) +
-			           l[i - 1] * (rho - F[i - 1] * env->rho_w)) *
-			          (-env->g);
+		W[i][0] = W[i][1] = 0.0;
+		W[i][2] = -g * (m_i - v_i * rho_w);
 
 		// relative flow velocity over node
 		const vec vi = U[i] - rd[i];
 		// tangential relative flow component
 		// <<<<<<< check sign since I've reversed q
-		const moordyn::real vql = vi.dot(q[i]);
+		const real vql = vi.dot(q[i]);
 		const vec vq = vql * q[i];
 		// transverse relative flow component
 		const vec vp = vi - vq;
 
-		const moordyn::real vq_mag = vq.norm();
-		const moordyn::real vp_mag = vp.norm();
+		const real vq_mag = vq.norm();
+		const real vp_mag = vp.norm();
 
+		const real shared_part = 0.5 * rho_w * d * submerged_length;
 		// transverse drag
-		if (i == 0)
-			Dp[i] = 0.25 * vp_mag * env->rho_w * Cdn * d * (F[i] * l[i]) * vp;
-		else if (i == N)
-			Dp[i] = 0.25 * vp_mag * env->rho_w * Cdn * d *
-			        (F[i - 1] * l[i - 1]) * vp;
-		else
-			Dp[i] = 0.25 * vp_mag * env->rho_w * Cdn * d *
-			        (F[i] * l[i] + F[i - 1] * l[i - 1]) * vp;
-
+		Dp[i] = vp_mag * Cdn * shared_part * vp;
 		// tangential drag
-		if (i == 0)
-			Dq[i] =
-			    0.25 * vq_mag * env->rho_w * Cdt * pi * d * (F[i] * l[i]) * vq;
-		else if (i == N)
-			Dq[i] = 0.25 * vq_mag * env->rho_w * Cdt * pi * d *
-			        (F[i - 1] * l[i - 1]) * vq;
-		else
-			Dq[i] = 0.25 * vq_mag * env->rho_w * Cdt * pi * d *
-			        (F[i] * l[i] + F[i - 1] * l[i - 1]) * vq;
+		Dq[i] = vq_mag * Cdt * pi * shared_part * vq;
 
-		// tangential component of fluid acceleration
 		// <<<<<<< check sign since I've reversed q
-		const moordyn::real aql = Ud[i].dot(q[i]);
+		const real aql = Ud[i].dot(q[i]);
 		const vec aq = aql * q[i];
 		// normal component of fluid acceleration
 		const vec ap = Ud[i] - aq;
 
 		// transverse Froude-Krylov force
-		if (i == 0)
-			Ap[i] = env->rho_w * (1. + Can) * 0.5 * (F[i] * V[i]) * ap;
-		else if (i == N)
-			Ap[i] = env->rho_w * (1. + Can) * 0.5 * (F[i - 1] * V[i - 1]) * ap;
-		else
-			Ap[i] = env->rho_w * (1. + Can) * 0.5 *
-			        (F[i] * V[i] + F[i - 1] * V[i - 1]) * ap;
+		Ap[i] = rho_w * (1. + Can) * submerged_volume * ap;
 		// tangential Froude-Krylov force
-		if (i == 0)
-			Aq[i] = env->rho_w * (1. + Cat) * 0.5 * (F[i] * V[i]) * aq;
-		else if (i == N)
-			Aq[i] = env->rho_w * (1. + Cat) * 0.5 * (F[i - 1] * V[i - 1]) * aq;
-		else
-			Aq[i] = env->rho_w * (1. + Cat) * 0.5 *
-			        (F[i] * V[i] + F[i - 1] * V[i - 1]) * aq;
+		Aq[i] = rho_w * (1. + Cat) * submerged_volume * aq;
 
 		// bottom contact (stiffness and damping, vertical-only for now) -
 		// updated for general case of potentially anchor or fairlead end in
 		// contact
 		const real waterDepth = getWaterDepth(r[i][0], r[i][1]);
 		if (r[i][2] < waterDepth) {
-			if (i == 0)
-				B[i][2] =
-				    ((waterDepth - r[i][2]) * env->kb - rd[i][2] * env->cb) *
-				    0.5 * d * (l[i]);
-			else if (i == N)
-				B[i][2] =
-				    ((waterDepth - r[i][2]) * env->kb - rd[i][2] * env->cb) *
-				    0.5 * d * (l[i - 1]);
-			else
-				B[i][2] =
-				    ((waterDepth - r[i][2]) * env->kb - rd[i][2] * env->cb) *
-				    0.5 * d * (l[i - 1] + l[i]);
+			B[i][2] = ((waterDepth - r[i][2]) * env->kb - rd[i][2] * env->cb) *
+			          d * (length);
 
 			// new rough-draft addition of seabed friction
 			real FrictionMax =
@@ -1223,44 +1176,16 @@ Line::getStateDeriv()
 			Bs[i] + Pb[i];
 	}
 
-	//	if (t > 5)
-	//	{
-	//		cout << " in getStateDeriv of line " << number << endl;
-	//
-	//		B[0][0] = 0.001; // meaningless
-	//	}
-
 	// loop through internal nodes and compute the accelerations
-	std::vector<vec> u, a;
-	u.reserve(N - 1);
-	a.reserve(N - 1);
+	vel.reserve(N - 1);
+	acc.reserve(N - 1);
 	for (unsigned int i = 1; i < N; i++) {
-		//	double M_out[9];
-		//	double F_out[3];
-		//	for (int I=0; I<3; I++)
-		//	{	F_out[I] = Fnet[i][I];
-		//		for (int J=0; J<3; J++) M_out[3*I + J] = M[i][I][J];
-		//	}
-
-		// solve for accelerations in [M]{a}={f} using LU decomposition
-		//	double LU[9];                        // serialized matrix that will
-		// hold LU matrices combined 	Crout(3, M_out, LU);                  //
-		// perform LU decomposition on mass matrix 	double acc[3]; //
-		// acceleration vector to solve for 	solveCrout(3, LU, F_out, acc);
-		// // solve for acceleration vector
-
-		//	LUsolve3(M[i], acc, Fnet[i]);
-
-		//	Solve3(M[i], acc, (const double*)Fnet[i]);
-
 		// For small systems it is usually faster to compute the inverse
 		// of the matrix. See
 		// https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html
-		a.push_back(M[i].inverse() * Fnet[i]);
-		u.push_back(rd[i]);
+		acc[i - 1] = M[i].inverse() * Fnet[i];
+		vel[i - 1] = rd[i];
 	}
-
-	return make_pair(u, a);
 };
 
 // write output file for line  (accepts time parameter since retained time value

source/Line.hpp

@@ -816,7 +816,7 @@ class Line final : public io::IO, public NatFreqCFL
 	 * @see moordyn::Line::setEndState
 	 * @throws invalid_value_error If either @p r or @p u have wrong sizes
 	 */
-	void setState(std::vector<vec> r, std::vector<vec> u);
+	void setState(const std::vector<vec>& r, const std::vector<vec>& u);
 
 	/** @brief Set the position and velocity of an end point
 	 * @param r Position
@@ -856,11 +856,11 @@ class Line final : public io::IO, public NatFreqCFL
 	vec getEndSegmentMoment(EndPoints end_point, EndPoints rod_end_point) const;
 
 	/** @brief Calculate forces and get the derivative of the line's states
-	 * @return The velocities of the internal nodes (first) and the accelerations
-	 * of the internal nodes (second)
+	 * @param vel Where to store the velocities of the internal nodes
+	 * @param acc Where to store the accelerations of the internal nodes
 	 * @throws nan_error If nan values are detected in any node position
 	 */
-	std::pair<std::vector<vec>, std::vector<vec>> getStateDeriv();
+	void getStateDeriv(std::vector<vec>& vel, std::vector<vec>& acc);
 
 	// void initiateStep(vector<double> &rFairIn, vector<double> &rdFairIn,
 	// double time);

source/State.hpp

@@ -50,12 +50,6 @@ template<typename T, typename V = T>
 class StateVar
 {
   public:
-	/// @brief Costructor
-	StateVar() {}
-
-	/// @brief Destructor
-	~StateVar() {}
-
 	/// The position
 	T pos;
 	/// The velocity
@@ -108,12 +102,6 @@ template<class T, class V = T>
 class StateVarDeriv
 {
   public:
-	/// @brief Constructor
-	StateVarDeriv() {}
-
-	/// @brief Destructor
-	~StateVarDeriv() {}
-
 	/// The velocity
 	T vel;
 	/// The acceleration
@@ -261,12 +249,6 @@ typedef StateVarDeriv<XYZQuat, vec6> DBodyStateDt;
 class MoorDynState
 {
   public:
-	/// @brief Constructor
-	MoorDynState() {}
-
-	/// @brief Destructor
-	~MoorDynState() {}
-
 	/// The states of the lines
 	std::vector<LineState> lines;
 
@@ -321,12 +303,6 @@ class MoorDynState
 class DMoorDynStateDt
 {
   public:
-	/// @brief Constructor
-	DMoorDynStateDt() {}
-
-	/// @brief Destructor
-	~DMoorDynStateDt() {}
-
 	/// The state derivatives of the lines
 	std::vector<DLineStateDt> lines;
 

source/Time.cpp

@@ -98,8 +98,8 @@ TimeSchemeBase<NSTATE, NDERIV>::CalcStateDeriv(unsigned int substep)
 	for (unsigned int i = 0; i < lines.size(); i++) {
 		if (!_calc_mask.lines[i])
 			continue;
-		std::tie(rd[substep].lines[i].vel, rd[substep].lines[i].acc) =
-		    lines[i]->getStateDeriv();
+		lines[i]->getStateDeriv(rd[substep].lines[i].vel,
+		                        rd[substep].lines[i].acc);
 	}
 
 	for (unsigned int i = 0; i < points.size(); i++) {