correct negative damping in BEM results (#134)

* add function to check nemoh results

currently checks for negative vals in diagonal FRFs

* don't issue warnings

may need to further consider whether/how to pass a warning up to higher
level functions and/or set debug/verbosity level across WecOptTool

* legend for debugging plot

* Allow fixing of negative damping values to be switched off

This is done through the parametric geometry option, which takes
an option called fixNegative, that defaults to true.

* add some inline comments to explain workflow

* Test for convergence if negative damping is fixed.

Doesn't seem to work.

* Check the value of the exitflag

* indent all functions

* comment out testPSNoConverge

our test case parameters for this
do not currently create a non-
convergent solution, so we must
wait to find a set of parameters
that result in such an error

* Update parametric test to PS case that might fail and organize
DeviceModel tests into new subpackage.

* Move tests

* Repose ParametricPSTest as a bug test

* Remove negative values on diagonal of damping matrix from NEMOH

* Test alternative formulation of CC equations

* Get meshes, where generated, from the device design routines and
simplify damping control equation

* Add mesh plotting function

* All tests passing. Needs a write up then conversion to active PR.

* Add description and test to plotMesh.

Use plotmesh in RM3/optimization.m example.

Fix other descriptions of the mesh output throughout.

* Remove incorrect docstring in WaveBot/Performance class

* Fix reference to removed types namespace and function

* Add copyright boilerplate

* Bug fix

* Fix Hydrodynamics docstring and add to API docs

* Minor docstring edit

* Update example docs

* Remove TODOs

* Reintroduce getAmplitudeSpectrum method into SeaState, but leave
interpolation to the user.

* Add getAmplitudeSpectrum to class docstring and give return value

* some minor typos

* use complex array for excitation

* explicitly set variable with load statement

Co-authored-by: Mathew Topper <damm_horse@yahoo.co.uk>

Author: Ryan Coe

docs/_static/example_meshes.svg

@@ -10,6 +10,9 @@ WecOptTool
 .. mat:autoclass:: +WecOptTool.AutoFolder
     :members:
 
+.. mat:autoclass:: +WecOptTool.Hydrodynamics
+    :members:
+
 .. mat:autoclass:: +WecOptTool.SeaState
     :members:
 

docs/user/api.rst

@@ -155,9 +155,9 @@ complete code is as follows:
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 63-91
+    :lines: 65-93
     :linenos:
-    :lineno-start: 63
+    :lineno-start: 65
 
 The following subsections will describe each stage of setting up the objective
 function.
@@ -202,9 +202,9 @@ object array. These inputs are combined and then added as arguments to the
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 69-72
+    :lines: 71-74
     :linenos:
-    :lineno-start: 69
+    :lineno-start: 71
 
 Calculate controlled device performance
 ---------------------------------------
@@ -249,9 +249,9 @@ as follows:
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 74-78
+    :lines: 76-80
     :linenos:
-    :lineno-start: 74
+    :lineno-start: 76
 
 Define the objective function value
 -----------------------------------
@@ -262,17 +262,17 @@ all spectra in the given seastate:
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 88-91
+    :lines: 92-95
     :linenos:
-    :lineno-start: 88
+    :lineno-start: 92
 
 Then, to make a minimisation problem, the negation of this value returned: 
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 80
+    :lines: 82
     :linenos:
-    :lineno-start: 80
+    :lineno-start: 82
 
 .. note::
     **Objective function:** The chosen objective function in |optimization.m|_ 
@@ -292,9 +292,9 @@ and then stash it for recovery later:
 
  .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 82-84
+    :lines: 84-88
     :linenos:
-    :lineno-start: 82
+    :lineno-start: 84
 
 Set optimization solver
 =======================
@@ -371,7 +371,7 @@ optimisation, the :mat:meth:`~+WecOptTool.AutoFolder.recoverVar` method of the
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 50-51
+    :lines: 50-52
     :linenos:
     :lineno-start: 50
 
@@ -381,24 +381,29 @@ match is found, as follows:
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 53-59
+    :lines: 54-60
     :linenos:
-    :lineno-start: 53
+    :lineno-start: 54
 
-Finally, by using the :mat:func:`+WecOptTool.+plot.powerPerFreq` function, the 
-spectral distribution of energy absorbed by the resulting design for each of 
-the eight sea states can be shown. 
+Finally, by using the :mat:func:`+WecOptTool.+plot.powerPerFreq` and
+:mat:func:`+WecOptTool.+plot.plotMesh` functions, the spectral distribution of 
+energy absorbed by the resulting design for each of the eight sea states can be 
+shown, alongside the simulated meshes.
 
 .. literalinclude:: /../examples/RM3/optimization.m
     :language: matlab
-    :lines: 61
+    :lines: 62-63
     :linenos:
-    :lineno-start: 61
+    :lineno-start: 62
 
 .. image:: /_static/example_spectralPowerPlot.svg
    :width: 400pt
    :alt: Absorbed Spectral power distribution for each sea state
 
+.. image:: /_static/example_meshes.svg
+   :width: 400pt
+   :alt: Absorbed Spectral power distribution for each sea state
+
 .. [Falnes] Falnes, Johannes. Ocean waves and oscillating systems: linear 
          interactions including wave-energy extraction. Cambridge University 
          Press, 2002.

docs/user/optimization.rst

@@ -1,13 +1,15 @@
-function hydro = designDevice(type, varargin) 
+function [hydro, meshes] = designDevice(type, varargin) 
 
+    meshes = [];
+
     switch type
 
         case 'existing'
             hydro = WecOptTool.geometry.existingNEMOH(varargin{:});
         case 'scalar'
             hydro = getHydroScalar(varargin{:});
         case 'parametric'
-            hydro = getHydroParametric(varargin{:});
+            [hydro, meshes] = getHydroParametric(varargin{:});
 
     end
 
@@ -20,7 +22,7 @@
     [filepath, ~, ~] = fileparts(p);
     dataPath = fullfile(filepath, 'RM3_BEM.mat');
 
-    load(dataPath, 'hydro');
+    hydro = struct2array(load(dataPath, 'hydro'));
 
     % dimensionalize w/ WEC-Sim built-in function
     hydro.rho = 1025;
@@ -34,15 +36,11 @@
     hydro.B = hydro.B .* lambda^2.5;
     hydro.A = hydro.A .* lambda^3;
     hydro.ex = complex(hydro.ex_re,hydro.ex_im) .* lambda^2;
-    hydro.ex_ma = abs(hydro.ex);
-    hydro.ex_ph = angle(hydro.ex);
-    hydro.ex_re = real(hydro.ex);
-    hydro.ex_im = imag(hydro.ex);
-           
+    hydro = rmfield(hydro,{'ex_re','ex_im'});
 end
 
 
-function hydro = getHydroParametric(folder, r1, r2, d1, d2, w)
+function [hydro, meshes] = getHydroParametric(folder, r1, r2, d1, d2, w)
 
     % Float
 

examples/RM3/designDevice.m

@@ -48,6 +48,7 @@
 [x, fval] = fmincon(objFun, x0, A, B, Aeq, Beq, lb, ub, NONLCON, opts);
 
 %% Recover device object of best simulation and plot its power per freq
+%% and mesh
 performances = folder.recoverVar("performances");
 
 for i = 1:length(performances)
@@ -59,6 +60,7 @@
 end
 
 WecOptTool.plot.powerPerFreq(bestPerformances);
+WecOptTool.plot.plotMesh(bestPerformances(1).meshes);
 
 %% Define objective function
 % This can take any form that complies with the requirements of the MATLAB
@@ -69,7 +71,7 @@
     w = seastate.getRegularFrequencies(0.5);
     geomParams = [folder.path num2cell(x) w];
 
-    deviceHydro = designDevice('parametric', geomParams{:});
+    [deviceHydro, meshes] = designDevice('parametric', geomParams{:});
 
     for j = 1:length(seastate)
         performances(j) = simulateDevice(deviceHydro,   ...
@@ -81,6 +83,8 @@
 
     [performances(:).w] = seastate.w;
     performances(1).x = x;
+    performances(1).meshes = meshes;
+    
     folder.stashVar(performances);
 
 end

examples/RM3/optimization.m

@@ -1,4 +1,4 @@
-function performance = simulateDevice(hydro, seastate, controlType, varargin)
+function [performance, motion] = simulateDevice(hydro, seastate, controlType, varargin)
 
     static = getStaticModel(hydro);
     motion = getDynamicModel(static, hydro, seastate);
@@ -48,9 +48,7 @@
 
     function result = interp_ex(hydro, dof, w)
 
-        h = complex(squeeze(hydro.ex_re(dof, 1, :)),   ...
-                    squeeze(hydro.ex_im(dof, 1, :)));
-
+        h = squeeze(hydro.ex(dof, 1, :));
         result = interp1(hydro.w, h ,w, 'linear', 0);
 
     end
@@ -89,6 +87,13 @@
 
     B33 = interp_rad(hydro, 3, 3, w) * hydro.rho .* w;
     A33 = interp_mass(hydro, 3, 3, w) * hydro.rho;
+    
+    % Check radiation damping
+    if max(abs(B39)) > max(max(B33), max(B99))
+        warning('WecOptTool:RM3:BadDamping',                            ...
+                ['Coupled radiation damping magnitude exceeds maximum ' ...
+                 'of single body values. Results may be unphysical.']);
+    end
 
     % Excitation
     H3 = interp_ex(hydro, 3, w) * hydro.g * hydro.rho;
@@ -162,22 +167,26 @@
 end
 
 function out = dampingControl(motion)
-            
-    % Max Power for a given Damping Coeffcient [Falnes 2002 
-    % (p.51-52)]
-    P_max = @(b) -0.5 * b *     ...
-                    sum(abs(motion.F0 ./ (motion.Zi + b)) .^ 2);
-
-    % Optimize the linear damping coeffcient(B)
-    B_opt = fminsearch(P_max, max(real(motion.Zi)));
 
-    % Power per frequency at optimial damping?
-    out.powPerFreq = 0.5 * B_opt * ...
-                    (abs(motion.F0 ./ (motion.Zi + B_opt)) .^ 2);
+    % Power per frequency at optimial damping
+    out.powPerFreq = 0.25 * abs(motion.F0) .^ 2 ./     ...
+                            (real(motion.Zi) + abs(motion.Zi));
 
 end
 
-function out = pseudoSpectralControl(motion, delta_Zmax, delta_Fmax)
+function out = pseudoSpectralControl(motion,        ...
+                                     delta_Zmax,    ...
+                                     delta_Fmax,    ...
+                                     display,       ...
+                                     OptimalityTolerance)
+                                 
+    arguments
+        motion
+        delta_Zmax
+        delta_Fmax
+        display = "off"
+        OptimalityTolerance = 1e-5
+    end
 
     % PSEUDOSPECTRAL Pseudo spectral control
     %   Returns power per frequency and frequency bins
@@ -202,7 +211,10 @@
     for ind_ph = 1 : n_ph
 
         ph = ph_mat(:, ind_ph);
-        [~, phasePowPerFreq] = getPSPhasePower(motion, ph);
+        [~, phasePowPerFreq] = getPSPhasePower(motion,  ...
+                                               ph,      ...
+                                               display, ...
+                                               OptimalityTolerance);
 
         for ind_freq = 1 : n_freqs
             powPerFreqMat(ind_ph, ind_freq) = phasePowPerFreq(ind_freq);
@@ -328,7 +340,10 @@
 
 end
 
-function [pow, powPerFreq] = getPSPhasePower(motion, ph)
+function [pow, powPerFreq] = getPSPhasePower(motion,    ...
+                                             ph,        ...
+                                             display,   ...
+                                             OptimalityTolerance)
     %Calculates power using the pseudospectral method given a phase and
     % a descrption of the body movement. Returns total phase power and 
     % power per frequency 
@@ -356,11 +371,10 @@
     % constrained optimiztion
     qp_options = optimoptions('fmincon',                        ...
                               'Algorithm', 'sqp',               ...
-                              'Display', 'off',                 ...
+                              'Display', display,               ...
                               'MaxIterations', 1e3,             ...
                               'MaxFunctionEvaluations', 1e5,    ...
-                              'OptimalityTolerance', 1e-8,      ...
-                              'StepTolerance', 1e-6);
+                              'OptimalityTolerance', OptimalityTolerance);
 
     siz = size(motion.A_ineq);
     X0 = zeros(siz(2),1);