Update MD-F comparison to point to OpenFAST/dev

.github/workflows/mdf_verification.yml

@@ -28,13 +28,14 @@ jobs:
       id: setup-python
 
     - name: Install Python dependencies
-      run: pip install --upgrade build pytest Bokeh numpy
+      run: 
+        pip install --upgrade build pytest Bokeh numpy matplotlib
 
     - name: Download OpenFAST
       shell: bash
       run: |
         cd ${{github.workspace}}/
-        git clone --recursive https://github.com/OpenFAST/openfast.git
+        git clone --recursive --single-branch -b dev https://github.com/OpenFAST/openfast.git
 
     - name: Compile MoorDynF
       shell: bash
@@ -73,3 +74,11 @@ jobs:
 
     - name: Run the tests
       run: python ${{github.workspace}}/tests/.mdf_verification/verify.py ${{github.workspace}}/
+
+    - name: Upload the artifacts
+      uses: actions/upload-artifact@v4
+      if: failure()
+      with:
+        name: MD-F comparison plots
+        path: |
+          ${{github.workspace}}/*.png

.gitignore

@@ -47,6 +47,13 @@ compile/SO/test*
 
 #program output files
 *.out
+*.png
+*.log
+*.dat
+tests/.mdf_verification/*.out
+tests/.mdf_verification/*.log
+tests/.mdf_verification/*.png
+tests/.mdf_verification/*.dat
 
 # Python packaging #
 ####################

docs/inputs.rst

@@ -595,6 +595,7 @@ The list of possible options is:
    If this is enabled initial conditions are calculated with scaled drag according to CdScaleIC. 
    The new stationary solver in MoorDyn-C is more stable and more precise than the dynamic solver, 
    but it can take longer to reach equilibrium.
+ - print_time (1): MoorDyn-C switch for printing the current timestep to the console. 
 
 A note about time steps in MoorDyn-C: The internal time step is first taken from the dtM option. If
 no CFL factor is provided, then the user provided time step is used to calculate CFL and MoorDyn-C 
@@ -629,6 +630,7 @@ The following MoorDyn-C options are not supported by MoorDyn-F:
  - FricDamp: Same as CV in MoorDyn-F.
  - StatDynFricScale: Same as MC in MoorDyn-F.
  - ICgenDynamic: MoorDyn-F does not have a stationary solver for initial conditions
+ - print_time: Console printing in MoorDyn-F is controled by OpenFAST
 
 The following options from MoorDyn-F are not supported by MoorDyn-C: 
 

source/MoorDyn2.cpp

@@ -91,6 +91,7 @@ moordyn::MoorDyn::MoorDyn(const char* infilename, int log_level)
   , dtOut(0.0)
   , _t_integrator(NULL)
   , ICgenDynamic(false)
+  , print_time(true)
   , env(std::make_shared<EnvCond>())
   , GroundBody(NULL)
   , waves(nullptr)
@@ -611,12 +612,15 @@ moordyn::MoorDyn::Step(const double* x,
                        double& dt)
 {
 	// should check if wave kinematics have been set up if expected!
-	const auto default_precision{std::cout.precision()};
-	std::cout << std::fixed << setprecision(1);
-	LOGDBG << "t = " << t << "s     \r";
-	std::cout << std::defaultfloat << setprecision(default_precision);
+
+	if (print_time) {
+		const auto default_precision{std::cout.precision()};
+		std::cout << std::fixed << setprecision(1);
+		LOGDBG << "t = " << t << "s     \r";
+		std::cout << std::defaultfloat << setprecision(default_precision);
 
-	cout << "\rt = " << t << " " << flush;
+		cout << "\rt = " << t << " " << flush;
+	}
 
 	if (dt <= 0) {
 		// Nothing to do, just recover the forces if there are coupled DOFs
@@ -2183,6 +2187,8 @@ moordyn::MoorDyn::readOptionsLine(vector<string>& in_txt, int i)
 		this->seafloor->setup(env, filepath);
 	} else if (name == "ICgenDynamic")
 		ICgenDynamic = bool(atof(entries[0].c_str()));
+	else if (name == "print_time")
+		print_time = bool(atof(entries[0].c_str()));
 	else
 		LOGWRN << "Warning: Unrecognized option '" << name << "'" << endl;
 }

source/MoorDyn2.hpp

@@ -550,6 +550,8 @@ class MoorDyn final : public io::IO
 	real ICthresh;
 	// use dynamic (true) or stationary (false) inital condition solver
 	bool ICgenDynamic;
+	/// print the timesteps. Disable for MD_F reg tests
+	bool print_time;
 	// temporary wave kinematics flag used to store input value while keeping
 	// env.WaveKin=0 for IC gen
 	moordyn::waves::waves_settings WaveKinTemp;

tests/.mdf_verification/verify.py

@@ -85,17 +85,19 @@ def to_num(s):
 
         if line.startswith('"'):
             # The value is a string
-            end = line[1:].find('"') + 1 
-            value = line[1:end]
+            end = line[1:].find('"') + 2 
+            value = line[1:end-1]
         else:
-            end = line.find(' ')
+            indicies = [line.find(' '), line.find('  '), line.find('\t')]
+            if -1 in indicies:
+                indicies.remove(-1)
+            end = min(indicies)
             value = to_num(line[:end])
-
-        line = line[end:]
+        line = line[end:].strip()
         while line.find("  ") != -1:
             line = line.replace("  ", " ")
-        key = line.split(" ")[1]
 
+        key = line.split(" ")[0]
         data[key] = value
 
     return data
@@ -106,8 +108,11 @@ def read_driver(test):
     test_root = os.path.join(args.root,
                              'openfast/reg_tests/r-test/modules/moordyn/',
                              test)
-    with open(os.path.join(test_root, "md_driver.inp"), "r") as f:
-        lines = f.readlines()
+    try: 
+        with open(os.path.join(test_root, "md_driver.inp"), "r") as f:
+            lines = f.readlines()
+    except FileNotFoundError: 
+        return -1,-1 
 
     def get_section(lines, name):
         start, end = None, None
@@ -153,7 +158,7 @@ def create_input_file(env, md):
                 end = i
                 break
         opts = mdopts2dict(opts[:end])
-        for optin, optout in (('Gravity', 'gravity'),
+        for optin, optout in (('Gravity', 'gravity'), # todo: make case insensitive
                               ('rhoW', 'rho'),
                               ('WtrDpth', 'WtrDpth')):
             if optin not in env.keys():
@@ -168,6 +173,8 @@ def create_input_file(env, md):
                         f"{optout}={opts[optout]} in MoorDyn config file")
             lines.insert(start + end,
                          f"{env[optin]}  {optout}  option set by the driver\n")
+        lines.insert(start+end, "0   print_time      - console printing of timesteps disabled for GH actions" )
+
         for line in lines:
             fout.write(line)
     return fname
@@ -254,34 +261,79 @@ def to_num(s):
         return float(s)
 
     data = []
+    heads = []
     with open(fpath, "r") as fin:
-        lines = fin.readlines()[skiplines:]
-        for line in lines:
+        lines = fin.readlines()[skiplines-2:]
+        line = lines[0]
+        line = line.strip().replace("\t", " ")
+        while line.find("  ") != -1:
+            line = line.replace("  ", " ")
+        heads = [field for field in line.split()]
+
+        for line in lines[2:]:
             line = line.strip().replace("\t", " ")
             while line.find("  ") != -1:
                 line = line.replace("  ", " ")
             data.append([to_num(field) for field in line.split()])
-    return np.transpose(data)
+    return np.transpose(data), heads
 
 
-def plot(ref, data, fpath):
+def plot(ref, ref_heads, data, data_heads, passed, rtol_mag, atol_mag, fpath):
     if plt is None:
         return
-    colors = list(mcolors.XKCD_COLORS.values())
+    fig,axes = plt.subplots(ref.shape[0]-1, 2, sharex = True, figsize=(12.8,4*(ref.shape[0]-1)))
+    if ref.shape[0] <= 2:
+        ax = [axes]
+    else:
+        ax = axes
     for i in range(1, ref.shape[0]):
-        plt.plot(ref[0, :], ref[i, :], linestyle='dashed',
-                 color=colors[i - 1])
-        plt.plot(data[0, :], data[i, :], linestyle='solid',
-                 color=colors[i - 1], label=f'channel {i}')
-    plt.legend(loc='best')
-    plt.savefig(fpath)
+        # Plot the channels
+        ax[i-1][0].plot(data[0, :], data[i, :], linestyle='solid',
+                color='b', label=f"MD-C: {(data_heads[i])}")
+        ax[i-1][0].plot(ref[0, :], ref[i, :], linestyle='dashed',
+                color='r', label=f"MD-F: {(ref_heads[i])}")
+        ax[i-1][0].legend(loc='best')
+
+        # Plot the difference. First calculate threshold (from OpenFAST testing erorPlotting.py _plotError)
+        NUMEPS = 1e-12
+        ATOL_MIN = 1e-6
+        baseline_offset = ref[i, :] - np.min(ref[i, :])
+        b_order_of_magnitude = np.floor( np.log10( baseline_offset + NUMEPS ) )
+        rtol = 10**(-1 * rtol_mag)
+        atol = 10**(max(b_order_of_magnitude) - atol_mag)
+        atol = max(atol, ATOL_MIN)
+        passfail_line = atol + rtol * abs(ref[i, :])
+        ax[i-1][1].plot(data[0, :], passfail_line, color='g', 
+                        label=f"Threshold {(data_heads[i])}")
+        ax[i-1][1].plot(data[0, :], abs(ref[i, :] - data[i, :]), color='gray', 
+                        label=f"Error {(data_heads[i])}")
+        ax[i-1][1].legend(loc='best')
+
+        # titles
+        if i == 1:
+            if passed:
+                ax[i-1][0].set_title("PASSED")
+            else:
+                ax[i-1][0].set_title("FAILED")
+
+            ax[i-1][1].set_title(f"abs diff btwn MD-C and MD-F")
+
+
+
+    fig.tight_layout()
+    fig.savefig(fpath, dpi=400)
+    plt.close()
 
 
 # Run the tests...
 summary = {}
 for test in tests:
+    print("\n ------------------")
     print(f"Test {test}...")
     env, md = read_driver(test)
+    if env == md == -1:
+        print(f"WARNING: Error in loading {test}. Does not match md_driver.inp convention")
+        continue
     fname = create_input_file(env, md)
     system = moordyn.Create(fname)
     # Get the NDoFs and check if the motions are right
@@ -305,17 +357,18 @@ def plot(ref, data, fpath):
         moordyn.Step(system, rorg, u, t, dt)
     moordyn.Close(system)
     # Read the ouputs and compare
-    ref = read_outs(md["OutRootName"], skiplines=8)
-    new = read_outs(os.path.splitext(fname)[0] + ".out", skiplines=2)
-    # Drop the eventual points at the tail that ight come from precision errors
+    ref, ref_heads = read_outs(md["OutRootName"], skiplines=8)
+    new, new_heads = read_outs(os.path.splitext(fname)[0] + ".out", skiplines=2)
+    # Drop the eventual points at the tail that might come from precision errors
     # on the time
     n_samples = min(ref.shape[1], new.shape[1])
     ref = ref[:, :n_samples]
     new = new[:, :n_samples]
-    plot(ref, new, test + ".png")
     passing = np.all(
         pass_fail.passing_channels(ref, new, args.rtol, args.atol))
 
+    plot(ref, ref_heads, new, new_heads, passing, args.rtol, args.atol, test + ".png")
+
     summary[test] = passing
 
 print("")