explicit buoyancy/gravity (#53)

* add mesh and volume properties

* make f_add a dict

* use dict to pass f_add

* use dict for f_add

* test for excess buoyancy - not working yet

* working example for surface device

* summarize changes

* with previous

* fix spelling/typos

* use function scope (default) for tests

* spelling typos

* add explanation in docs

* use value iterator

* hyperlink to class API

Co-authored-by: Ryan Coe <rcoe@sandia.gov>

CHANGES.md

@@ -4,6 +4,9 @@
 ## Version 1.0.3
 
 * minor updates to README
+* `f_add` should be passed as a `dict`, e.g., `{'my_name': my_func}`
+  * optionally treat buoyancy/gravity explicitly via user-defined functions passed to `f_add`
+  * time and freq domain results calculated for entries of `f_add` after `solve` completes
 * logging of decision vector and objective function 
   * controlled entirely via logging package config
   * move to `info` logging level (`debug` gives too much from other packages, e.g., matplotlib)

docs/source/theory.rst

@@ -91,7 +91,7 @@ Consider, for example, a general case without a controller structure, in which :
 For a wave tank scale device, one might expect velocities of :math:`\mathcal{O}(1e-1)`, but the forces could be :math:`\mathcal{O}(1e3)`.
 For larger WECs, this discrepancy in the orders of magnitude may be even worse.
 Scaling mismatches in the decision variable :math:`x` and with the objective function :math:`J(x)` can lead to problems with convergence.
-To alleviate this issue, WecOptTool allows users to set scale factors for the components of :math:`x` as well as the objective function (see :meth:`core.WEC.solve`).
+To alleviate this issue, WecOptTool allows users to set scale factors for the components of :math:`x` as well as the objective function (see :meth:`wecopttool.WEC.solve`).
 Additionally, you may set :code:`import logging, logging.basicConfig(level=logging.INFO)` to output the mean values of `x` and the objective function during the solution process.
 
 Constraints
@@ -100,7 +100,14 @@ Constraints, such as maximum PTO force, maximum piston force, or maintaining ten
 This functionality is well-illustrated in :doc:`_examples/tutorial_1_wavebot`.
 An important practical factor when using this functionality is to make sure that the constraint is evaluated at a sufficient number of collocation points.
 It may be required to enforce constraints at more points than the dynamics (as defined by the frequency array).
-In WecOptTool's example PTO module, this is controlled by the :code:`nsubsteps` argument (see, e.g., :meth:`pto.PseudoSpectralPTO.force`).
+In WecOptTool's example PTO module, this is controlled by the :code:`nsubsteps` argument (see, e.g., :py:meth:`pto.PseudoSpectralPTO.force`).
+
+Buoyancy/gravity
+^^^^^^^^^^^^^^^^
+As WecOptTool is intended primarily to utilize linear potential flow hydrodynamics, a linear hydrostatic stiffness is used.
+The implicit assumption of this approach is that the body is neutrally buoyant (i.e., gravitational and buoyancy forces are in balance at the zero position).
+However, some WECs achieve equilibrium through a pretension applied via mooring and/or the PTO.
+In this case, the device can still be modeled with the linear hydrostatic stiffness, but if you wish to solve for the pretension force in your simulations, you may explicitly include the buoyancy, gravity, and pretension forces via the :code:`f_add` argument to :py:class:`wecopttool.core.WEC`.
 
 .. _WEC-Sim: https://wec-sim.github.io/WEC-Sim/master/index.html
 .. _Autograd: https://github.com/HIPS/autograd

tests/test_wecopttool.py

@@ -14,7 +14,7 @@
 from wecopttool.core import power_limit
 
 
-@pytest.fixture(scope="module")
+@pytest.fixture()
 def wec():
     # water properties
     rho = 1000.0
@@ -58,7 +58,7 @@ def const_f_pto(wec, x_wec, x_opt):
     constraints = [ineq_cons]
 
     # WEC
-    f_add = pto.force_on_wec
+    f_add = {'PTO': pto.force_on_wec}
 
     wec = wot.WEC(fb, mass, stiffness, f0, nfreq,  rho=rho,
                     f_add=f_add, constraints=constraints)
@@ -69,7 +69,7 @@ def const_f_pto(wec, x_wec, x_opt):
     return wec
 
 
-@pytest.fixture(scope="module")
+@pytest.fixture()
 def regular_wave(wec):
     freq = 0.5
     amplitude = 0.25
@@ -78,7 +78,7 @@ def regular_wave(wec):
     return wave
 
 
-@pytest.fixture(scope="module")
+@pytest.fixture()
 def resonant_wave(wec):
     freq = wec.natural_frequency()[0].squeeze().item()
     amplitude = 0.25
@@ -87,7 +87,7 @@ def resonant_wave(wec):
     return wave
 
 
-@pytest.fixture(scope="module")
+@pytest.fixture()
 def pto(wec):
     kinematics = np.eye(wec.ndof)
     pto = wot.pto.PseudoSpectralPTO(wec.nfreq, kinematics)
@@ -266,18 +266,18 @@ def test_wavebot_ps_theoretical_limit(wec,regular_wave,pto):
 
 def test_wavebot_p_cc(wec,resonant_wave):
     """Check that power from proportional damping controller can match
-    theorectical limit at the natural resonance.
+    theoretical limit at the natural resonance.
     """
 
-    # remove contraints
+    # remove constraints
     wec.constraints = []
 
     # update PTO
     kinematics = np.eye(wec.ndof)
     pto = wot.pto.ProportionalPTO(kinematics)
     obj_fun = pto.average_power
     nstate_opt = pto.nstate
-    wec.f_add = pto.force_on_wec
+    wec.f_add = {'PTO': pto.force_on_wec}
 
     _, fdom, _, xopt, average_power, _ = wec.solve(resonant_wave, obj_fun, nstate_opt,
         optim_options={'maxiter': 1000, 'ftol': 1e-8}, scale_x_opt=1e3)
@@ -289,18 +289,18 @@ def test_wavebot_p_cc(wec,resonant_wave):
 
 def test_wavebot_pi_cc(wec,regular_wave):
     """Check that power from proportional integral (PI) controller can match
-    theorectical limit at any single wave frequency (i.e., regular wave).
+    theoretical limit at any single wave frequency (i.e., regular wave).
     """
 
-    # remove contraints
+    # remove constraints
     wec.constraints = []
 
     # update PTO
     kinematics = np.eye(wec.ndof)
     pto = wot.pto.ProportionalIntegralPTO(kinematics)
     obj_fun = pto.average_power
     nstate_opt = pto.nstate
-    wec.f_add = pto.force_on_wec
+    wec.f_add = {'PTO': pto.force_on_wec}
 
     tdom, fdom, xwec, xopt, average_power, res = wec.solve(regular_wave,
         obj_fun, nstate_opt,
@@ -320,3 +320,43 @@ def test_examples_device_wavebot_mesh():
 def test_examples_device_wavebot_plot_cross_section():
     wb = WaveBot()
     wb.plot_cross_section()
+
+
+def test_buoyancy_excess(wec, pto, regular_wave):
+    """Give too much buoyancy and check that equilibrium point found matches
+    that given by the hydrostatic stiffness"""
+
+    delta = np.random.randn() # excess buoyancy factor
+
+    # remove constraints
+    wec.constraints = []
+
+    def f_b(wec, x_wec, x_opt):
+        V = wec.fb.keep_immersed_part(inplace=False).mesh.volume
+        rho = wec.rho
+        g = 9.81
+        return (1+delta) * rho * V * g * np.ones([wec.ncomponents, wec.ndof])
+
+    def f_g(wec, x_wec, x_opt):
+        g = 9.81
+        m = wec.mass.item()
+        return -1 * m * g * np.ones([wec.ncomponents, wec.ndof])
+
+    wec.f_add = {**wec.f_add,
+                 'Fb':f_b,
+                 'Fg':f_g,
+                 }
+
+    tdom, fdom, x_wec, x_opt, avg_pow, _ = wec.solve(regular_wave, 
+                                                obj_fun = pto.average_power,
+                                                nstate_opt = pto.nstate,
+                                                scale_x_wec = 1.0,
+                                                scale_x_opt = 0.01,
+                                                scale_obj = 1e-1,
+                                                optim_options={})
+
+    mean_pos = tdom.pos.squeeze().mean().item()
+    expected = (wec.rho * wec.fb.mesh.volume * wec.g * delta) \
+        / wec.hydrostatic_stiffness.item()
+    assert pytest.approx (expected, 1e-1) == mean_pos
+

wecopttool/core.py

@@ -13,7 +13,7 @@
 
 
 import logging
-from typing import Iterable, Callable, Any
+from typing import Iterable, Callable, Any, Optional, Mapping
 from pathlib import Path
 
 import numpy.typing as npt
@@ -55,8 +55,8 @@ def __init__(self, fb: cpy.FloatingBody, mass: np.ndarray,
                  hydrostatic_stiffness: np.ndarray, f0: float, nfreq: int,
                  dissipation: np.ndarray | None = None,
                  stiffness: np.ndarray | None = None,
-                 f_add: Callable[[WEC, np.ndarray, np.ndarray], np.ndarray] |
-                 None = None, constraints: list[dict] = [],
+                 f_add: Optional[Mapping[str, Callable[[WEC, np.ndarray, np.ndarray], np.ndarray]]] = None,
+                 constraints: list[dict] = [],
                  rho: float = _default_parameters['rho'],
                  depth: float = _default_parameters['depth'],
                  g: float = _default_parameters['g']) -> None:
@@ -83,9 +83,10 @@ def __init__(self, fb: cpy.FloatingBody, mass: np.ndarray,
             Additional stiffness for the impedance calculation in
             ``capytaine.post_pro.impedance``. Shape:
             (``ndof`` x ``ndof`` x ``1``) or (``ndof`` x ``ndof`` x ``nfreq``).
-        f_add: function
-            Additional forcing terms (e.g. PTO, mooring, etc.) for the
-            WEC dynamics in the time-domain. Takes three inputs:
+        f_add: dict[str, Callable]
+            Additional forcing terms (e.g. buoyancy, gravity, PTO, mooring, 
+            etc.) for the WEC dynamics in the time-domain. Dictionary entries 
+            should be ``entry = {'name': function_handle}``. Takes three inputs:
             (1) the WEC object,
             (2) the WEC dynamics state (1D np.ndarray), and
             (3) the optimization state (1D np.ndarray)
@@ -115,7 +116,11 @@ def __init__(self, fb: cpy.FloatingBody, mass: np.ndarray,
         super().__setattr__('freq', (f0, nfreq))
 
         # additional WEC dynamics forces
-        super().__setattr__('f_add', f_add)
+        if callable(f_add):
+            log.debug(f"Assigning dictionary entry 'f_add'" + 
+                      "for Callable argument {f_add}")
+            f_add = {'f_add': f_add}
+        self.f_add = f_add
         if stiffness is None:
             stiffness = 0.0
         super().__setattr__('stiffness', stiffness)
@@ -133,8 +138,7 @@ def __init__(self, fb: cpy.FloatingBody, mass: np.ndarray,
         super().__setattr__('hydro', None)
 
     def __setattr__(self, name, value):
-        """Delete dependent attributes  when user manually modifies
-        an attribute.
+        """Delete dependent attributes when user manually modifies an attribute.
         """
         _attrs_delete_mass = ['fb']
         _attrs_delete_stiffness = ['fb', 'rho', 'g']
@@ -186,6 +190,20 @@ def __repr__(self):
         return str_info
 
     # PROPERTIES
+    # properties: f_add
+    @property
+    def f_add(self):
+        """Additonal forces on the WEC (e.g., PTO, mooring, buoyancy, gravity)"""
+        return self._f_add
+
+    @f_add.setter
+    def f_add(self, f_add):
+        if callable(f_add):
+            log.debug(f"Assigning dictionary entry 'f_add'" + 
+                      "for Callable argument {f_add}")
+            f_add = {'f_add': f_add}
+        super().__setattr__('_f_add', f_add)
+
     # properties: frequency
     @property
     def freq(self):
@@ -273,6 +291,23 @@ def derivative_mat(self):
         """Derivative matrix for the state vector."""
         return self._derivative_mat
 
+    # properties mesh
+    @property
+    def mesh(self):
+        return self.fb.mesh
+
+    @property
+    def volume(self):
+        return self.mesh.volume
+
+    @property
+    def submerged_mesh(self):
+        return self.mesh.keep_immersed_part()
+
+    @property
+    def submerged_volume(self):
+        return self.submerged_mesh.volume
+
     ## METHODS
     # methods: class I/O
     def to_file(self, fpath: str | Path) -> None:
@@ -736,7 +771,7 @@ def callback(x):
 
         # post-process
         x_wec, x_opt = self.decompose_decision_var(res.x)
-        fd_x, td_x = self._post_process_wec_dynamics(x_wec)
+        fd_x, td_x = self._post_process_wec_dynamics(x_wec, x_opt)
         fd_we = fd_we.reset_coords(drop=True)
         time_dom = xr.merge([td_x, td_we])
         freq_dom = xr.merge([fd_x, fd_we])
@@ -772,13 +807,15 @@ def _dynamic_residual(self, x: np.ndarray, f_exc: np.ndarray
         f_i = np.dot(self.time_mat, f_i)
 
         # additional forces
-        if self.f_add is not None:
-            f_add = self.f_add(self, x_wec, x_opt)
-        else:
-            f_add = 0.0
+        f_add = 0.0
+        for f_add_fun in self.f_add.values():
+            f_add = f_add + f_add_fun(self, x_wec, x_opt)
+
         return f_i - f_exc - f_add
 
-    def _post_process_wec_dynamics(self, x_wec: np.ndarray
+    def _post_process_wec_dynamics(self, 
+                                   x_wec: np.ndarray,
+                                   x_opt: np.ndarray
                                    ) -> tuple[xr.DataArray, xr.DataArray]:
         """Transform the results from optimization solution to a form
         that the user can work with directly.
@@ -831,6 +868,13 @@ def _post_process_wec_dynamics(self, x_wec: np.ndarray
         acc_fd = xr.DataArray(
             acc_fd, dims=dims_fd, coords=coords_fd, attrs=attrs_acc)
         freq_dom = xr.Dataset({'pos': pos_fd, 'vel': vel_fd, 'acc': acc_fd},)
+
+        # user-defined additional forces (in WEC DoFs)
+        for f_add_key, f_add_fun in self.f_add.items():
+            time_dom[f_add_key] = (('time', 'influenced_dof'), 
+                                   f_add_fun(self, x_wec, x_opt))
+            freq_dom[f_add_key] = (('omega', 'influenced_dof'),
+                                   self.td_to_fd(time_dom[f_add_key]))
 
         return freq_dom, time_dom