Merge pull request FloatingArrayDesign#140 from RyanDavies19/ryan_dev

Updates from dev before new version release

Author: RyanDavies19

docs/compiling.rst

@@ -165,7 +165,65 @@ the Git window select "Clone Existing Repository".
 
    The Git GUI to clone repositories
 
-We are starting with Eigen3, so in the first box of the window that pops up set
+As a default, MoorDyn uses the Egien package that is internal in the source code.
+If you would like to use an external copy of Eigen, please follow the instructions
+in the :ref:`installing with external eigen <external_eigen>` note at this point. 
+
+We will install MoorDyn following a very similar process.
+Launch CMake again, and set "C:\MoorDyn\MoorDyn" in the source box and
+"C:\MoorDyn\MoorDyn.build" in the binaries box, clicking "Configure" afterwards.
+Select again the "MinGW Makefiles" for the generator.
+When the configuration options appear, set CMAKE_BUILD_TYPE as "Release", and
+enable FORTRAN_WRAPPER and PYTHON_WRAPPER:
+
+.. figure:: win_cmake_moordyn.png
+   :alt: Configuration options for MoorDyn
+
+   Configuration options for MoorDyn
+
+You can also enable MATLAB_WRAPPER if you have MATLAB installed in your system.
+We are ready, click "Configure" once more and then "Generate".
+
+Now go back to your Command Prompt from earlier (which has administrative rights), and
+type the following commands:
+
+.. code-block:: bash
+
+  cd C:\MoorDyn\MoorDyn.build
+  mingw32-make
+  mingw32-make install
+
+This will generate three libraries in the MoorDyn/build/source directory labeled 
+libmoordyn, libmoordyn.2, and libmoordyn.2.2. The first two are symbolic links to the 
+latter, setup that way for development purposes. In your project you should use 
+libmoordyn.
+
+
+NOTE: If you want to generate a Windows installer, disable the PYTHON_WRAPPER
+option and type
+
+.. code-block:: bash
+
+  cd C:\MoorDyn\MoorDyn.build
+  mingw32-make
+  cpack -C Release
+
+NOTE: If you are working on a proxy serveryou may need to add the .crt file for your proxy 
+configuration to the folder ``C:/msys64/etc/pki/ca-trust/source/anchors`` or equivalent for your 
+system.
+
+NOTE: You may need to upgrade or install the build tool using pip
+
+.. code-block:: bash
+  
+  \<path-to-python>/python<version>.exe -m pip install --upgrade build
+
+NOTE: Installing External Eigen
+
+.. _external_eigen:
+
+To use an external copy of Eigen, ensure that the DEXTERNAL_EIGEN flag is turned on.
+In the first box of the window that pops up set
 "https://gitlab.com/libeigen/eigen.git", and in the second "C:\MoorDyn\eigen":
 
 .. figure:: win_git_eigen.png
@@ -235,55 +293,6 @@ Now you just need to type the following commands:
 
 We will need to use cmd with administrative rights later on, so do not close it.
 
-Now we will install MoorDyn following a very similar process.
-Launch CMake again, and set "C:\MoorDyn\MoorDyn" in the source box and
-"C:\MoorDyn\MoorDyn.build" in the binaries box, clicking "Configure" afterwards.
-Select again the "MinGW Makefiles" for the generator.
-When the configuration options appear, set CMAKE_BUILD_TYPE as "Release", and
-enable FORTRAN_WRAPPER and PYTHON_WRAPPER:
-
-.. figure:: win_cmake_moordyn.png
-   :alt: Configuration options for MoorDyn
-
-   Configuration options for MoorDyn
-
-You can also enable MATLAB_WRAPPER if you have MATLAB installed in your system.
-We are ready, click "Configure" once more and then "Generate".
-
-Now go back to your Command Prompt from earlier (which has administrative rights), and
-type the following commands:
-
-.. code-block:: bash
-
-  cd C:\MoorDyn\MoorDyn.build
-  mingw32-make
-  mingw32-make install
-
-This will generate three libraries in the MoorDyn/build/source directory labeled 
-libmoordyn, libmoordyn.2, and libmoordyn.2.2. The first two are symbolic links to the 
-latter, setup that way for development purposes. In your project you should use 
-libmoordyn.
-
-
-NOTE: If you want to generate a Windows installer, disable the PYTHON_WRAPPER
-option and type
-
-.. code-block:: bash
-
-  cd C:\MoorDyn\MoorDyn.build
-  mingw32-make
-  cpack -C Release
-
-NOTE: If you are working on a proxy serveryou may need to add the .crt file for your proxy 
-configuration to the folder ``C:/msys64/etc/pki/ca-trust/source/anchors`` or equivalent for your 
-system.
-
-NOTE: You may need to upgrade or install the build tool using pip
-
-.. code-block:: bash
-  
-  \<path-to-python>/python<version>.exe -m pip install --upgrade build
-
 Linux and Mac
 ^^^^^^^^^^^^^
 

source/Line.cpp

@@ -485,34 +485,27 @@ Line::initialize()
 		                       &VF,
 		                       &HA,
 		                       &VA,
-		                       N,
+		                       N+1,
 		                       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;
-				for (unsigned int i = 1; i < N; i++) {
-					vec l(Xl[i] * COSPhi, Xl[i] * SINPhi, Zl[i]);
-					r[i] = r[0] + l;
-				}
+			
+			// the catenary solve is successful, update the node positions
+			LOGDBG << "Catenary initial profile available for Line "
+					<< number << endl;
+			for (unsigned int i = 1; i < N; i++) {
+				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;
+			       << ", initalizing as linear " << endl;
 		}
 	} else {
-		LOGDBG << "Vertical initial profile for Line " << number << endl;
+		LOGDBG << "Vertical linear initial profile for Line " << number << endl;
 	}
 
 	LOGMSG << "Initialized Line " << number << endl;

source/Line.hpp

@@ -137,7 +137,7 @@ class Line final : public io::IO
 	moordyn::real UnstrLend;
 	/// line diameter
 	moordyn::real d;
-	/// line linear density
+	/// line density (kg/m^3)
 	moordyn::real rho;
 	/// line elasticity modulus (Young's modulus) [Pa]
 	moordyn::real E;

source/MoorDyn2.cpp

@@ -108,22 +108,21 @@ moordyn::MoorDyn::MoorDyn(const char* infilename, int log_level)
 		_basepath = _filepath.substr(0, lastSlash + 1);
 	}
 
-	LOGMSG << "\n Running MoorDyn (v2.a10, 2022-01-01)" << endl
+	LOGMSG << "\n Running MoorDyn (v2.0.0, 2023-09-18)" << endl
 	       << "         MoorDyn v2 has significant ongoing input file changes "
 	          "from v1."
 	       << endl
-	       << "   Copyright: (C) 2021 National Renewable Energy Laboratory, "
+	       << "   Copyright: (C) 2023 National Renewable Energy Laboratory, "
 	          "(C) 2014-2019 Matt Hall"
 	       << endl
-	       << "   This program is released under the GNU General Public "
-	          "License v3."
+	       << "   This program is released under the  BSD 3-Clause license."
 	       << endl;
 
 	LOGMSG << "The filename is " << _filepath << endl;
 	LOGDBG << "The basename is " << _basename << endl;
 	LOGDBG << "The basepath is " << _basepath << endl;
 
-	env->g = 9.8;
+	env->g = 9.80665;
 	env->WtrDpth = 0.;
 	env->rho_w = 1025.;
 	env->kb = 3.0e6;
@@ -956,7 +955,7 @@ moordyn::MoorDyn::ReadInFile()
 				LOGERR << "Error in " << _filepath << ":" << i + 1 << "..."
 				       << endl
 				       << "'" << in_txt[i] << "'" << endl
-				       << "7 fields are required, but just " << entries.size()
+				       << "7 fields are required, but only " << entries.size()
 				       << " are provided" << endl;
 				return MOORDYN_INVALID_INPUT;
 			}
@@ -2091,8 +2090,8 @@ moordyn::MoorDyn::detachLines(FailProps* failure)
 moordyn::error_id
 moordyn::MoorDyn::AllOutput(double t, double dt)
 {
-	if (env->writeLog == 0)
-		return MOORDYN_SUCCESS;
+	// if (env->writeLog == 0)
+	// 	return MOORDYN_SUCCESS;
 
 	if (dtOut > 0)
 		if (t < (floor((t - dt) / dtOut) + 1.0) * dtOut)

source/QSlines.hpp

@@ -95,7 +95,7 @@ Catenary(T XF,
          T* HAout,
          T* VAout,
          unsigned int Nnodes,
-         const vector<T>& s,
+         vector<T>& s,
          vector<T>& X,
          vector<T>& Z,
          vector<T>& Te)
@@ -109,6 +109,16 @@ Catenary(T XF,
 	if (longwinded == 1)
 		cout << "hi" << endl;
 
+	// Checking inverted points at line ends: from catenary.py in MoorPy
+	bool reverseFlag;
+	if ( ZF <  0.0 ){ // True if the fairlead has passed below its anchor
+        ZF = -ZF;
+        reverseFlag = true;
+		if (longwinded == 1)
+     		cout << " Warning from catenary: "
+			 	 << "Anchor point is above the fairlead point" << endl;
+    } else reverseFlag = false;
+
 	// Maximum stretched length of the line with seabed interaction beyond which
 	// the line would have to double-back on itself; here the line forms an "L"
 	// between the anchor and fairlead (i.e. it is horizontal along the seabed
@@ -226,8 +236,8 @@ Catenary(T XF,
 		// The current mooring line must be slack and not vertical
 		Lamda0 = sqrt(3.0 * ((L * L - ZF2) / XF2 - 1.0));
 	}
-	// ! As above, set the lower limit of the guess value of HF to the tolerance
-	HF = max((T)abs(0.5 * W * XF / Lamda0), Tol);
+
+	HF = abs(0.5 * W * XF / Lamda0);
 	VF = 0.5 * W * (ZF / tanh(Lamda0) + L);
 
 	/*
@@ -570,5 +580,25 @@ Catenary(T XF,
 	*HAout = HA;
 	*VAout = VA;
 
-	return 1;
+	if (reverseFlag) { // return values to normal
+		reverse(s.begin(), s.end());
+		reverse(X.begin(), X.end());
+		reverse(Z.begin(), Z.end());
+		reverse(Te.begin(), Te.end());
+		for (unsigned int I = 0; I < Nnodes; I++){
+			s[I] = L - s[I];
+			X[I] = XF - X[I];
+			Z[I] = Z[I] - ZF;
+		}
+		ZF = -ZF;
+	}
+
+	// It might happen that the output solution does not respect the
+	// queried final point. See the pendulum example
+	if (abs(Z[Nnodes-1] - ZF) > Tol) {
+		if (longwinded == 1)
+			cout << "Fairlead and anchor vertical seperation has changed"
+				 << ", aborting catenary solver ..." <<  endl;
+		return -1;
+	} else return 1;
 }

source/Rod.cpp

@@ -1132,7 +1132,7 @@ Rod::doRHS()
 			Mext += vec(Mtemp * sinBeta, -Mtemp * cosBeta, 0.0);
 
 			// axial drag
-			Dq[i] += VOF[i] * Area * env->rho_w * CdEnd * vq_mag * vq;
+			Dq[i] += 0.5 * VOF[i] * Area * env->rho_w * CdEnd * vq_mag * vq;
 
 			// long-wave diffraction force
 			const real V_temp = 2.0 / 3.0 * pi * d * d * d / 8.0;
@@ -1159,7 +1159,7 @@ Rod::doRHS()
 			Mext += vec(Mtemp * sinBeta, -Mtemp * cosBeta, 0.0);
 
 			// axial drag
-			Dq[i] += VOF[i] * Area * env->rho_w * CdEnd * vq_mag * vq;
+			Dq[i] += 0.5 * VOF[i] * Area * env->rho_w * CdEnd * vq_mag * vq;
 
 			// long-wave diffraction force
 			const real V_temp = 2.0 / 3.0 * pi * d * d * d / 8.0;

tests/beam.cpp

@@ -47,7 +47,7 @@
 #define W 100.0
 #define EA 1.5e9
 #define EI 1.0e7
-#define G 9.81
+#define G 9.80665
 
 bool
 compare(double v1, double v2)

tests/pendulum.cpp

@@ -113,7 +113,7 @@ pendulum()
 	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 w = sqrt(9.80665 / l0);
 	const double T = 2.0 * M_PI / w;
 	double dt = T / 100.0;
 	double t = 0.0;

tests/rods.cpp

@@ -57,7 +57,7 @@ added_mass()
 	double r = 0.25;
 	double V = length * M_PI * r * r;
 	double rho_w = 1025;
-	double g = 9.8;
+	double g = 9.80665;
 	double buoyancy = rho_w * V * g;
 	double weight = mass * g;
 	double Fnet = buoyancy - weight;

tests/seafloor.cpp

@@ -108,7 +108,7 @@ pendulum()
 	double z_max = 5.0;
 	double z_step = 1.0;
 
-	const double w = sqrt(9.81 / l0);
+	const double w = sqrt(9.80665 / l0);
 	const double T = 2.0 * M_PI / w;
 	double dt = T / 20.0;
 	double t = 0.0;