Enable Variable Mass (WEC-Sim#1465)

* refactor radiation CIC calculation

* add radiation conv integral surface calc for variable hydro

* add CIC surface tests

* increase time step to speed up cicTest.m

* add cicTest to wecSimTest.m

* update ci surface docstring

* rename new function to camelCase

* create IRF surfaces in the body class

* minor clean up

* add excitation IRF time to body.variableHydro

* clarify where simu.cicDt is used

* excitation convolution integral surface calc

* add surface calcs and variant subsystems to library

* update library tag

* remove check preventing VH with userDefinedExcitation

* add goto tags for excitation and total forces

* fix goto tag name and turn off a warning

* comment to explain why a warning is turned off

* update excitation CI surface calc to use persistent amplitude

* add test for excitationConvolutionIntegralSurface

* fix issue template typo

* add empty bodies.hydroForceIndex to responseClass

* move the adjusted mass used in simulation to body.hydroForce

* fix hydroForce.mass definition for all cases

* clear persistent functions before simulation

* progress on library

* working but fundamentally incorrect excitation IRF surface

* remove excitation IRF surface function

* fix one more bug removing the excitation irf surface

* library to r2020b

* update post processing added mass force with VH

* fix added mass force postprocessing again

* docs - comment on added mass treatment with variable hydro

* docs - explain IRF surface for variable hydro

* update explanation of our handling of added mass

* fix reference typo

* add figures of an irf surface

* updates to bodyClass

* remove comments

* initialize inertialMass

* clarify gravity force in added mass dev docs

* resolve some issues

* update docs on use of mass vs adjustedMass

* updates

* Fix indexing

* add inertia products option and warning message

* amend warning

* revise warning

* remove extraneous moordyn readme

* clarifying comment in docs

* add docstring to body.variableHydro

* expand warning message

* add error to prevent variable hydro and elevationImport

* disable nonlinear hydro and variable hydro

* update docs and error message

* typo

* use hydroforce.mass for nonlinear hydro

* Revert "use hydroforce.mass for nonlinear hydro"

This reverts commit fa7ca47.

* storage.mass to mass

* fix typos

---------

Co-authored-by: akeeste <akeeste@sandia.gov>

Author: jtgrasb

docs/_include/added_mass.rst

@@ -15,9 +15,9 @@ shown in the manipulation of the governing equation below:
 
 .. math::
 
-    M\ddot{X_i} &= \Sigma F(t,\omega) - A\ddot{X_i} \\
-    (M+A)\ddot{X_i} &= \Sigma F(t,\omega) \\
-    M_{adjusted}\ddot{X_i} &= \Sigma F(t,\omega)
+    M\ddot{X_i} &= \Sigma F(t,\omega) + m\overrightarrow{g} - A\ddot{X_i} \\
+    (M+A)\ddot{X_i} &= \Sigma F(t,\omega) + m\overrightarrow{g} \\
+    M_{adjusted}\ddot{X_i} &= \Sigma F(t,\omega) + m\overrightarrow{g} 
 
 where :math:`F_{am} = -A\ddot{X_i}` is the added mass component of the radiation force that is pulled out of the force summation for manipulation,
 :math:`M` is the mass matrix, :math:`A` is the single frequency or infinite frequency added mass, and subscript :math:`i` represents the timestep being solved for. 
@@ -42,6 +42,10 @@ The core issue with this combined mass formulation is that Simscape Multibody, a
 For example, a rigid body can't have one mass for surge motion and another mass for heave motion. 
 Simscape rigid bodies only have one translational mass, a 1x3 moment of inertia matrix, and 1x3 product of inertia matrix. 
 
+Note that the body mass used to define the gravitational force is not adjusted as gravity does not act on any added mass contributions. 
+The gravity force will always depend on the unadjusted body mass. 
+The gravity force is lumped with the sum of remaining forces in the rest of this section for simplicity.
+
 WEC-Sim's Implementation
 ^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -54,7 +58,7 @@ There is a 1-1 mapping between the body's inertia tensor and rotational added ma
 These added mass coefficients are entirely lumped with the body's inertia.
 Additionally, the surge-surge (1,1), sway-sway (2,2), heave-heave (3,3) added mass coefficients correspond to the translational mass of the body, but must be treated identically.
 
-WEC-Sim implements this added mass treatment by adding a change in mass matrix ($$dM$$) to both sides of the equation, creating both a modified added mass matrix and a modified body mass matrix:
+WEC-Sim implements this added mass treatment by adding a change in mass matrix (:math:`dM`) to both sides of the equation, creating both a modified added mass matrix and a modified body mass matrix:
 
 .. math::
 
@@ -121,12 +125,7 @@ This will convert the algebraic loop equation of motion to a solvable one:
     M_{adjusted}\ddot{X_i} &= \Sigma F(t,\omega) - A_{adjusted}\ddot{X}_{i - (1 - 10^{-7}/dt)} \\
 
 Body-to-body Interactions
-"""""""""""""""""""""""""""
-F = A * acc
-first dimension/index = down, 2nd = across
-non b2b: A = [6x6], acc = [6x1], 
-b2b: A = [6x12], acc = [12x1] in order of body numbers regardless of the current body number
-
+""""""""""""""""""""""""""
 The above implementation extends readily to the case where there are body interactions to account for.
 In this example, we assume there are two bodies with body interaction.
 Then the right hand side added mass and acceleration matrices above are (without generalized modes)
@@ -154,17 +153,45 @@ Working with the Added Mass Implementation
 WEC-Sim's added mass implementation should not affect a user's modeling workflow.
 WEC-Sim handles the manipulation and restoration of the mass and forces in the bodyClass functions ``bodyClass.adjustMassMatrix()`` called by ``initializeWecSim`` and ``bodyClass.restoreMassMatrix()``, ``bodyClass.storeForceAddedMass()`` called by ``postProcessWecSim``.
 However viewing ``body.hydroForce.hf*.fAddedMass`` between calls to ``initializeWecSim`` and ``postProcessWecSim`` will not show the values from the BEM dataset.
-Users can get the adjusted mass, moments of inertia, products of inertia, added mass coefficients, added mass force, and total force from 
-``body.hydroForce.hf*.mass``, ``body.hydroForce.hf*.inertia``, ``body.hydroForce.hf*.inertiaProducts``, ``body.hydroForce.hf*.storage.hydroForce_fAddedMass``, 
-``body.hydroForce.hf*.storage.output_forceAddedMass``, and ``body.hydroForce.hf*.storage.output_forceTotal`` respectively after the simulation.
-However, in the case that a user wants to view the added mass force *during* a simulation (typically during debugging), 
+See the following table for the variables containing both unadjusted and adjusted parameters after the simulation:
+
+.. list-table:: Nominal and adjusted mass parameters
+   :widths: 25 37 37
+   :header-rows: 1
+
+   * - Parameter
+     - Nominal value
+     - Adjusted value
+   * - Mass
+     - ``body.mass``, ``body.hydroForce.hf*.mass``, ``body.variableHydro.mass``
+     - ``body.hydroForce.hf*.adjustedMass``
+   * - Inertia
+     - ``body.inertia``, ``body.variableHydro.inertia``
+     - ``body.hydroForce.hf*.adjustedInertia``
+   * - Inertia Products
+     - ``body.inertiaProducts``, ``body.variableHydro.inertiaProducts``
+     - ``body.hydroForce.hf*.adjustedInertiaProducts``
+   * - Added mass coefficients
+     - ``body.hydroForce.hf*.hydroForce_fAddedMass``
+     - ``body.hydroForce.hf*.storage.hydroForce_fAddedMass``
+   * - Added mass force timeseries
+     - ``output.bodies(*).forceAddedMass``
+     - ``body.hydroForce.hf*.storage.output_forceAddedMass``
+   * - Total force timeseries
+     - ``output.bodies(*).forceTotal``
+     - ``body.hydroForce.hf*.storage.output_forceTotal``
+
+The nominal (*unadjusted*) body mass is stored in ``body.hydroForce.hf*.mass``. 
+This nominal value is retained during the simulation because it is required for the calculation of the gravitational force.
+However, in the case that a user wants to view the added mass force *during* a simulation (typically during debugging in a failed simulation), 
 the change in mass, :math:`dM` above, must be taken into account. 
 Refer to how ``body.calculateForceAddedMass()`` calculates the entire added mass force in WEC-Sim post-processing.
 
-When using variable hydrodynamics, the added mass matrix can change with the varying state.
-However, the body mass matrix cannot vary. 
-So, :math:`dM` is taken to be constant based on the hydrodynamic dataset specified by `body.variableHydro.hydroForceIndexInitial`.
-All added mass datasets within `body.hydroForce` are changed using this consistent definition of :math:`dM`.
+When using variable hydrodynamics, the added mass matrix and mass matrix can change with the varying state.
+In the above derivation, all values of :math:`dM`, mass matrix, inertia matrix, added mass force, and the resultant
+adjusted mass and added mass matrices are calculated for each hydrodynamic dataset.
+Note that varying any part of the body mass matrix with variable hydro requires using a special simscape block.
+Refer to the variable hydrodynamics applications.
 
 .. Note::
 	Depending on the wave formulation used, :math:`A` can either be a function of wave frequency :math:`A(\omega)`, or equal to the added mass at infinite wave frequency :math:`A_{\infty}`

docs/_include/variable_hydro.rst

@@ -18,6 +18,11 @@ User defined logic and user selected signals determine
 when the state changes. The Variable Hydrodynamics feature does not determine
 a specific scenario, state, signal, or discretization required.
 
+For a case which varies the mass of the body, the body library link in the Simulink 
+model will need to be broken and the default File Solid : Body Properties block 
+replaced with the General Variable Mass block. See the variable hydro WEC-Sim 
+Application.
+
 Example
 """"""""
 
@@ -108,7 +113,7 @@ every body in a simulation. To implement variable hydrodynamics for a given body
     * FIR Filter radiation calculations
     * Generalized body modes
     * Non-hydrodynamic and drag bodies
-    * Conditions that require a variable mass, center of gravity, or center of buoyancy
+    * Nonlinear hydrodynamics is not compatible when using a case with variable mass.
 
 Impulse Response Function with Variable Hydrodynamics
 """""""""""""""""""""""""""""""""""""""""""""""""""""

docs/developer/advanced_features.rst

@@ -51,8 +51,8 @@ structure:
 
 * ``h5File``
    * Instead of initializing the body class with one string 
-     (one H5 file name), a cell array of strings can be passed, e.g.
-   * body(1) = bodyClass({'H5FILE_1.h5','H5FILE_2.h5','H5FILE_3.h5');
+     (one H5 file name), a cell array of strings can be passed, e.g
+     ``body(1) = bodyClass({'H5FILE_1.h5','H5FILE_2.h5','H5FILE_3.h5');``
 * ``hydroData``
    * Each H5 file is processed into a hydroData structure. All structs are 
      concatenated into ``hydroData`` which is now an array of structures instead

source/functions/initializeWecSim.m

@@ -372,6 +372,15 @@
     end; clear iW
 end
 
+% Check for elevationImport with variable hydro
+for ii = 1:simu.numHydroBodies
+    for iW = 1:length(waves)
+        if strcmp(waves(iW).type,'elevationImport') && body(ii).variableHydro.option == 1
+            error('Cannot run WEC-Sim with Variable Hydrodynamics (body(ii).variableHydro.option>0) and "elevationImport" wave type')
+        end
+    end; clear iW
+end
+
 % Check for morisonElement inputs for body(ii).morisonElement.option == 1 || body(ii).morisonElement.option == 2
 for ii = 1:length(body(1,:))
     if body(ii).morisonElement.option == 1
@@ -562,7 +571,7 @@
 fprintf('\nSimulating the WEC device defined in the SimMechanics model %s...   \n',simu.simMechanicsFile)
 % Modify some stuff for simulation
 for iBod = 1:simu.numHydroBodies
-    body(iBod).adjustMassMatrix(simu.b2b);
+    body(iBod).adjustMassMatrix(simu.b2b, simu.rho);
 end; clear iBod
 
 % Create the buses for hydroForce