Repository: FloatingArrayDesign/MoorDyn (language: python)

Issue:
Static line tension disagreement
Static lines initialize with different tensions between MDC_v2 and MDF_v2/MDC_v1. Differences also emerge with dynamic relaxation. Using [MoorPy](https://github.com/NREL/MoorPy), NREL's quasi-static model, the equilibrium fairlead tension is 4808 kN. This indicates that MDC_v1 and MDF_v2 are getting a solution closer to the correct answer. See the figures below for a single line resting partially on the seabed:

 Line profile:
![case_1_shape](https://github.com/FloatingArrayDesign/MoorDyn/assets/101124339/5c5bfe21-b3e8-467f-9181-c540eba7704f)

Fairlead tension with TmaxIC = 0:
![case1_comparison](https://github.com/FloatingArrayDesign/MoorDyn/assets/101124339/c2d75686-be37-47f8-8cd2-127303193d2c)

Fairlead tension with TmaxIC = 200:
![200s_dyrelax](https://github.com/FloatingArrayDesign/MoorDyn/assets/101124339/2785b163-220b-47e9-833a-422c9e00b936)




Discussion / hints:
If you have a MDC_v2 model file, I would definitely be interested to take a look. It looks like the disagreement is percentage wise pretty small, but it would definitely be good to track it down.
@AlexWKinley Sure thing, here it is: 
[case1.txt](https://github.com/FloatingArrayDesign/MoorDyn/files/12886501/case1.txt)

I've noticed dynamic relaxation differences popping up elsewhere too. I am currently working on tracking them down to put together a fix. @AlexWKinley if you see something that pops out to you that would be great to know. 
4808 x 1.025 = 4928

Might it be just a different default water density?

On Thu, 12 Oct 2023, 21:39 RyanDavies19, ***@***.***> wrote:

> I've noticed dynamic relaxation differences popping up elsewhere too. I am
> currently working on tracking them down to put together a fix.
> @AlexWKinley <https://github.com/AlexWKinley> if you see something that
> pops out to you that would be great to know.
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/FloatingArrayDesign/MoorDyn/issues/160#issuecomment-1760259401>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AAMXKKD3YIT2ONYDIH5PKPTX7BBOLAVCNFSM6AAAAAA55447RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONRQGI2TSNBQGE>
> .
> You are receiving this because you are subscribed to this thread.Message
> ID: ***@***.***>
>

@sanguinariojoe I don't believe so, at least not from the default value point. For MDC v2: https://github.com/FloatingArrayDesign/MoorDyn/blob/99cd2b6e6afa7e5d58e2ac7ac3a4b05e127f1f72/source/MoorDyn2.cpp#L127

For MDC v1: https://github.com/FloatingArrayDesign/MoorDyn/blob/ea5ccf8c050f359cdd13a0a8b6f5cc5f3158405c/MoorDyn.cpp#L307

And MDF requires a manual input of rho via a driver input file, it is set to 1025 there. 
I do believe though it is a default value or scale issue of some kind because of the constant difference 
It was just an idea. I let you guys deal with this. I have some other
projects to attend

On Thu, 12 Oct 2023, 22:47 RyanDavies19, ***@***.***> wrote:

> @sanguinariojoe <https://github.com/sanguinariojoe> I don't believe so,
> at least not from the default value point. For MDC v2:
> https://github.com/FloatingArrayDesign/MoorDyn/blob/99cd2b6e6afa7e5d58e2ac7ac3a4b05e127f1f72/source/MoorDyn2.cpp#L127
>
> For MDC v1:
> https://github.com/FloatingArrayDesign/MoorDyn/blob/ea5ccf8c050f359cdd13a0a8b6f5cc5f3158405c/MoorDyn.cpp#L307
>
> And MDF requires a manual input of rho via a driver input file, it is set
> to 1025 there.
>
> —
> Reply to this email directly, view it on GitHub
> <https://github.com/FloatingArrayDesign/MoorDyn/issues/160#issuecomment-1760343731>,
> or unsubscribe
> <https://github.com/notifications/unsubscribe-auth/AAMXKKBQ7OGDO7IRMZCSJYTX7BJOBAVCNFSM6AAAAAA55447RCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMYTONRQGM2DGNZTGE>
> .
> You are receiving this because you were mentioned.Message ID:
> ***@***.***>
>

@RyanDavies19 

Edit 1:
Actually I lied, I had made a change to the model that was having a more significant impact (increasing the number of segments). Rerunning without the change shows no significant improvement. 

Edit 2:
I found something that might actually matter, although I'm still not sure sure how important it is for getting agreement with the other version.

In Line.hhp Line;:getNodeTen
```c++
	inline vec getNodeTen(unsigned int i) const
	{
		if (i > N) {
			LOGERR << "Asking node " << i << " of line " << number
			       << ", which only has " << N + 1 << " nodes" << std::endl;
			throw moordyn::invalid_value_error("Invalid node index");
		}
		if ((i == 0) || (i == N))
			return (
			    Fnet[i] +
			    vec(0.0, 0.0, M[i](0, 0) * (-env->g))); // <<< update to use W

		// take average of tension in adjacent segments
		return (0.5 * (T[i] + T[i - 1]));
	};
```
It looks like we are adding an additional weight term to the end force. As far as I know, Fnet[N] should already contain that weight term. For reference the Line::getEndStuff function just returns Fnet[N] for ENDPOINT_TOP. 
A quick look at the v1 branch suggests to me this is a change carried over from there, but perhaps it was more justified in that case?

I'm not entirely convinced by this, but it would pretty accurately match the idea of a constant offset, and if you change it to just return Fnet[i] for the ends, the tensions decrease to more closely match the others.
![image](https://github.com/FloatingArrayDesign/MoorDyn/assets/129765925/39c9ff22-8e32-4f84-8a97-101f3d2b58e0)

I'll leave it up to you to determine if that behavior is present or correct in the other version of MoorDyn.

## Left because it still is a difference but doesn't actually make a difference
One thing that seems to have an effect is the fact that MoorDyn v2 does submergence calculations that make a decent approximation for cylinders. 

This causes the top segment to have slightly less buoyancy because it's like top corner is sticking out of the water.
Crude illustration of what part of the cylinder is out of the water
![image](https://github.com/FloatingArrayDesign/MoorDyn/assets/129765925/a6ba49f3-2a44-4e8a-ba37-f3dcf4f72122)

 If I remove the partial submergence calculation in Line::getStateDeriv
```c++
// F[i] = calcSubSeg(i, i + 1, surface_height);
F[i] = 1.0;
```
the fairlead tension does decrease to be a little closer to the others.
![image](https://github.com/FloatingArrayDesign/MoorDyn/assets/129765925/1ba96a4f-e095-4399-a15d-b77b58266bf0)

But it doesn't look like that is the only factor. You might want to validate my findings with your setup, I'm using a slightly old and slightly modified version of MoorDyn v2 (although my numbers generally agree with yours before the submergence change).
@AlexWKinley I think it was the first one you noted, with regards to adding an extra weight. V1 did have the added weight found here: https://github.com/FloatingArrayDesign/MoorDyn/blob/ea5ccf8c050f359cdd13a0a8b6f5cc5f3158405c/Line.cpp#L324
Removing the added weight in v2 causes all the versions to line up (under the TmaxIC = 0 case).  Here is the comparison with that weight removed:
![Mod_get_ten](https://github.com/FloatingArrayDesign/MoorDyn/assets/101124339/848f29fd-2e59-441c-82c9-4a12ea5d3701)

Another note though is that it seems like the dynamic relaxation methods have some differences. I will leave this issue open until we can resolve those as well. With the removal of that added weight, here are the fairtens with TmaxIC = 200:
![get_ten_mod_200s](https://github.com/FloatingArrayDesign/MoorDyn/assets/101124339/7c993c21-e112-4b32-83da-eb9d01da74b6)

@mattEhall do you have any inputs on either of these?



Hi all,
Yeah, I agree this line https://github.com/FloatingArrayDesign/MoorDyn/blob/d1010520a4285f650c696a0ba0545c3ce58fe756/source/Line.hpp#L434 is causing a double counting of the weight contribution to the tension of line end nodes. If we remove it the code will be consistent with v2 Fortran.

Regarding differences in the initial conditions between C and F, the first thing I'd wonder is if the codes hit the convergence tolerance at different times (assuming they don't go through all of TMaxIC=200 seconds). I don't recall if the convergence criterion is the same in both codes. 

@mattEhall That's likely it, I believe the Fortran code reaches it first. I am working on that now and will have a PR up soon with both the weight double counting fix and that relaxation criteria standardized. 

Fix patch:
diff --git a/source/Line.hpp b/source/Line.hpp
index 7df36f9c..e278ebb0 100644
--- a/source/Line.hpp
+++ b/source/Line.hpp
@@ -429,9 +429,11 @@ class Line final : public io::IO
 			throw moordyn::invalid_value_error("Invalid node index");
 		}
 		if ((i == 0) || (i == N))
+			// return (
+			//     Fnet[i] +
+			//     vec(0.0, 0.0, M[i](0, 0) * (-env->g))); // <<< update to use W
 			return (
-			    Fnet[i] +
-			    vec(0.0, 0.0, M[i](0, 0) * (-env->g))); // <<< update to use W
+			    Fnet[i]); // <<< update to use W
 
 		// take average of tension in adjacent segments
 		return (0.5 * (T[i] + T[i - 1]));
diff --git a/source/MoorDyn2.cpp b/source/MoorDyn2.cpp
index 682a4f14..9e801455 100644
--- a/source/MoorDyn2.cpp
+++ b/source/MoorDyn2.cpp
@@ -192,7 +192,7 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 
 	// Allocate past line fairlead tension array, which is used for convergence
 	// test during IC gen
-	const unsigned int convergence_iters = 10;
+	const unsigned int convergence_iters = 9; // 10 iterations, indexed 0-9
 	vector<real> FairTensLast_col(convergence_iters, 0.0);
 	for (unsigned int i = 0; i < convergence_iters; i++)
 		FairTensLast_col[i] = 1.0 * i;
@@ -296,7 +296,7 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 	// vector to store tensions for analyzing convergence
 	vector<real> FairTens(LineList.size(), 0.0);
 
-	unsigned int iic = 0;
+	unsigned int iic = 1; // To match MDF indexing
 	real t = 0;
 	bool converged = true;
 	real max_error = 0.0;
@@ -308,7 +308,10 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 	real best_score = (std::numeric_limits<real>::max)();
 	real best_score_t = 0.0;
 	unsigned int best_score_line = 0;
-	while ((t < ICTmax) && (!skip_ic)) {
+
+	// //dtIC set to fraction of input so convergence is over dtIC
+	ICdt = ICdt / (convergence_iters+1);
+	while (((ICTmax-t) > 0.00000001) && (!skip_ic)) { // tol of 0.00000001 should be smaller than anything anyone puts in as a ICdt
 		// Integrate one ICD timestep (ICdt)
 		real t_target = ICdt;
 		real dt;
@@ -382,7 +385,7 @@ moordyn::MoorDyn::Init(const double* x, const double* xd, bool skip_ic)
 		if (converged) {
 			LOGMSG << "Fairlead tensions converged" << endl;
 		} else {
-			LOGWRN << "Fairlead tensions did not converged" << endl;
+			LOGWRN << "Fairlead tensions did not converge" << endl;
 		}
 		LOGMSG << "Remaining error after " << t << " s = " << 100.0 * max_error
 		       << "% on line " << max_error_line << endl;
diff --git a/source/Rod.hpp b/source/Rod.hpp
index ef5ff3c2..8223cbc4 100644
--- a/source/Rod.hpp
+++ b/source/Rod.hpp
@@ -412,6 +412,7 @@ class Rod final : public io::IO
 	{
 		Cdn = Cdn * scaler;
 		Cdt = Cdt * scaler;
+		CdEnd = CdEnd * scaler;
 	}
 
 	/** @brief Set the line  simulation time


Test patch:
diff --git a/tests/pendulum.cpp b/tests/pendulum.cpp
index e0ae1b6f..8593ef3b 100644
--- a/tests/pendulum.cpp
+++ b/tests/pendulum.cpp
@@ -147,7 +147,7 @@ pendulum()
 	}
 
 	for (unsigned int i = 0; i < t_peaks.size(); i++) {
-		CHECK_VALUE("T", t_peaks[i], i * 0.5 * T, 0.02 * T, t_peaks[i]);
+		CHECK_VALUE("T", t_peaks[i], i * 0.5 * T, 0.025 * T, t_peaks[i]);
 		CHECK_VALUE("X", fabs(x_peaks[i]), x0, 0.01 * x0, t_peaks[i]);
 	}