added tests and changes during testing

tests/test_utilities.py

@@ -0,0 +1,209 @@
+""" Unit tests for functions in the :python:`utilities.py` module.
+"""
+
+import pytest
+import numpy as np
+import xarray as xr
+from matplotlib.pyplot import Figure, Axes
+import wecopttool as wot
+from pytest import approx
+import capytaine as cpy
+
+
+
+# test function in the utilities.py
+
+
+
+@pytest.fixture(scope="module")
+def power_flows():
+    """Dictionary of power flows."""
+    pflows = {'Optimal Excitation': -100,
+                'Radiated': -20,
+                'Actual Excitation': -70,
+                'Electrical (solver)': -40,
+                'Mechanical (solver)': -50,
+                'Absorbed': -50,
+                'Unused Potential': -30,
+                'PTO Loss': -10
+    }
+    return pflows
+
+@pytest.fixture(scope="module")
+def f1():
+    """Fundamental frequency [Hz]."""
+    return 0.1
+
+
+@pytest.fixture(scope="module")
+def nfreq():
+    """Number of frequencies in frequency vector."""
+    return 5
+
+@pytest.fixture(scope="module")
+def ndof():
+    """Number of degrees of freedom."""
+    return 2
+
+@pytest.fixture(scope="module")
+def ndir():
+    """Number of wave directions."""
+    return 3
+
+@pytest.fixture(scope='module')
+def bem_data(f1, nfreq, ndof, ndir):
+    """Synthetic BEM data."""
+    # TODO - start using single BEM solution across entire test suite
+    coords = {
+        'omega': [2*np.pi*(ifreq+1)*f1 for ifreq in range(nfreq)],
+        'influenced_dof': [f'DOF_{idof+1}' for idof in range(ndof)],
+        'radiating_dof': [f'DOF_{idof+1}' for idof in range(ndof)],
+        'wave_direction': [2*np.pi/ndir*idir for idir in range(ndir)],
+    }
+    radiation_dims = ['omega', 'radiating_dof', 'influenced_dof']
+    excitation_dims = ['omega', 'influenced_dof', 'wave_direction']
+    hydrostatics_dims = ['radiating_dof', 'influenced_dof']
+
+    added_mass = np.ones([nfreq, ndof, ndof])
+    radiation_damping = np.ones([nfreq, ndof, ndof])
+    diffraction_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j
+    Froude_Krylov_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j
+    inertia_matrix = np.ones([ndof, ndof])
+    hydrostatic_stiffness = np.ones([ndof, ndof])
+
+    data_vars = {
+        'added_mass': (radiation_dims, added_mass),
+        'radiation_damping': (radiation_dims, radiation_damping),
+        'diffraction_force': (excitation_dims, diffraction_force),
+        'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force),
+        'inertia_matrix': (hydrostatics_dims, inertia_matrix),
+        'hydrostatic_stiffness': (hydrostatics_dims, hydrostatic_stiffness)
+    }
+    return xr.Dataset(data_vars=data_vars, coords=coords)
+
+@pytest.fixture(scope='module')
+def intrinsic_impedance(bem_data):
+    bem_data = wot.add_linear_friction(bem_data)
+    intrinsic_impedance = wot.hydrodynamic_impedance(bem_data)
+    return intrinsic_impedance
+
+@pytest.fixture(scope='module')
+def pi_controller_pto():
+    """Basic PTO: proportional-integral (PI) controller, 1DOF, mechanical
+    power."""
+    ndof = 1
+    pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof),
+                      controller=wot.pto.controller_pi,
+                      names=["PI controller PTO"])
+    return pto
+
+@pytest.fixture(scope='module')
+def regular_wave(f1, nfreq):
+    """Single frequency wave"""
+    wfreq = 0.3
+    wamp = 0.0625
+    wphase = 0
+    wdir = 0
+    waves = wot.waves.regular_wave(f1, nfreq, wfreq, wamp, wphase, wdir)
+    return waves
+
+@pytest.fixture(scope="module")
+def fb():
+    """Capytaine FloatingBody object"""
+    try:
+        import wecopttool.geom as geom
+    except ImportError:
+        pytest.skip(
+            'Skipping integration tests due to missing optional geometry ' +
+            'dependencies. Run `pip install wecopttool[geometry]` to run ' +
+            'these tests.'
+            )
+    mesh_size_factor = 0.5
+    wb = geom.WaveBot()
+    mesh = wb.mesh(mesh_size_factor)
+    fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot")
+    fb.add_translation_dof(name="Heave")
+    return fb
+
+
+@pytest.fixture(scope="module")
+def wb_bem(f1, nfreq, fb):
+    """Boundary elemement model (Capytaine) results"""
+    freq = wot.frequency(f1, nfreq, False)
+    return wot.run_bem(fb, freq)
+
+@pytest.fixture(scope='class')
+def wb_hydro_impedance(wb_bem):
+    """Intrinsic hydrodynamic impedance"""
+    hd = wot.add_linear_friction(wb_bem)
+    hd = wot.check_linear_damping(hd)
+    Zi = wot.hydrodynamic_impedance(hd)
+    return Zi
+
+
+
+
+def test_plot_hydrodynamic_coefficients(bem_data,ndof):
+    bem_figure_list = wot.utilities.plot_hydrodynamic_coefficients(bem_data)
+    correct_len = ndof*(ndof+1)/2   #using only the subdiagonal elements
+    #added mass
+    fig_am = bem_figure_list[0][0]
+    assert correct_len == len(fig_am.axes)
+    assert isinstance(fig_am,Figure)
+    #radiation damping
+    fig_rd = bem_figure_list[1][0]
+    assert correct_len == len(fig_rd.axes)
+    assert isinstance(fig_rd,Figure)
+    #radiation damping
+    fig_ex = bem_figure_list[2][0]
+    assert ndof == len(fig_ex.axes)
+    assert isinstance(fig_ex,Figure)
+
+def test_plot_bode_impedance(intrinsic_impedance, ndof):
+    fig_Zi, axes_Zi = wot.utilities.plot_bode_impedance(intrinsic_impedance)    
+
+    assert 2*ndof*ndof == len(fig_Zi.axes)
+    assert isinstance(fig_Zi,Figure)
+    assert all([isinstance(ax, Axes) for ax in np.reshape(axes_Zi,-1)])
+
+
+def test_plot_power_flow(power_flows):
+    fig_sankey, ax_sankey = wot.utilities.plot_power_flow(power_flows)    
+
+    assert isinstance(fig_sankey, Figure)
+    assert isinstance(ax_sankey, Axes) 
+
+def test_calculate_power_flow(wb_bem,
+                              regular_wave,
+                              pi_controller_pto,
+                              wb_hydro_impedance):
+    """PI controller matches optimal for any regular wave, 
+        thus we check if the radiated power is equal the absorber power
+        and if the Optimal excitation is equal the actual excitation"""
+
+    f_add = {"PTO": pi_controller_pto.force_on_wec}
+    wec = wot.WEC.from_bem(wb_bem, f_add=f_add)
+
+    res = wec.solve(waves=regular_wave,
+                    obj_fun=pi_controller_pto.average_power,
+                    nstate_opt=2,
+                    x_wec_0=1e-1*np.ones(wec.nstate_wec),
+                    x_opt_0=[-1e3, 1e4],
+                    scale_x_wec=1e2,
+                    scale_x_opt=1e-3,
+                    scale_obj=1e-2,
+                    optim_options={'maxiter': 50},
+                    bounds_opt=((-1e4, 0), (0, 2e4),)
+                    )
+
+    pflows = wot.utilities.calculate_power_flows(wec, 
+                          pi_controller_pto, 
+                          res[0], 
+                          regular_wave.sel(realization = 0), 
+                          wb_hydro_impedance)
+
+    assert pflows['Absorbed'] == approx(pflows['Radiated'], rel=1e-4)
+    assert pflows['Optimal Excitation'] == approx(pflows['Actual Excitation'], rel=1e-4)
+
+
+

wecopttool/utilities.py

@@ -28,6 +28,8 @@
 from xarray import DataArray, concat
 import matplotlib.pyplot as plt
 from matplotlib.figure import Figure
+from matplotlib.axes import Axes
+
 from matplotlib.sankey import Sankey
 
 
@@ -39,7 +41,8 @@
 
 
 def plot_hydrodynamic_coefficients(bem_data,
-                                   wave_dir: Optional[float] = 0.0)-> None:
+                                   wave_dir: Optional[float] = 0.0
+                                   )-> list(tuple(Figure, Axes)):
     """Plots hydrodynamic coefficients (added mass, radiation damping,
        and wave excitation) based on BEM data.
 
@@ -131,12 +134,12 @@ def plot_hydrodynamic_coefficients(bem_data,
                 fig_am.delaxes(ax_am[i, j])
                 fig_rd.delaxes(ax_rd[i, j])
     fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False)
-
+    return [(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)]
 
 def plot_bode_impedance(impedance: DataArray, 
-                        title: Optional[str]= None,
+                        title: Optional[str]= '',
                         #plot_natural_freq: Optional[bool] = False,
-)-> None:
+)-> tuple(Figure, Axes):
     """Plot Bode graph from wecoptool impedance data array.
 
     Parameters
@@ -254,19 +257,19 @@ def calculate_power_flows(wec,
         P_e.append((1/4)*(Fe_FD_t@np.conj(U_FD) + U_FD_t@np.conj(Fe_FD)))
 
     power_flows = {
-        'Optimal Excitation' : 2* np.sum(np.real(P_max)),#eq 6.68 positive inflow
-        'Radiated': -1*np.sum(np.real(P_r)), #negative because "out"flow
-        'Actual Excitation': -1*np.sum(np.real(P_e)), #negative because "out"flow
-        'Electrical (solver)': P_elec,  #solver determins sign
-        'Mechanical (solver)': P_mech, #solver determins sign
+        'Optimal Excitation' : -2* np.sum(np.real(P_max)),#eq 6.68 
+        'Radiated': -1*np.sum(np.real(P_r)), 
+        'Actual Excitation': -1*np.sum(np.real(P_e)), 
+        'Electrical (solver)': P_elec, 
+        'Mechanical (solver)': P_mech, 
                   }
 
     power_flows['Absorbed'] =  (
         power_flows['Actual Excitation'] 
         - power_flows['Radiated']
             )
     power_flows['Unused Potential'] =  (
-        -1*power_flows['Optimal Excitation'] 
+        power_flows['Optimal Excitation'] 
         - power_flows['Actual Excitation']
             )
     power_flows['PTO Loss'] = (
@@ -276,7 +279,7 @@ def calculate_power_flows(wec,
     return power_flows
 
 
-def plot_power_flow(power_flows: dict[str, float])-> Figure:
+def plot_power_flow(power_flows: dict[str, float])-> tuple(Figure, Axes):
     """Plot power flow through a WEC as Sankey diagram.
 
     Parameters
@@ -290,19 +293,24 @@ def plot_power_flow(power_flows: dict[str, float])-> Figure:
 
 
     """
-    fig = plt.figure(figsize = [8,4])
-    ax = fig.add_subplot(1, 1, 1,)
-    plt.viridis()
+    # fig = plt.figure(figsize = [8,4])
+    # ax = fig.add_subplot(1, 1, 1,)
+    fig, ax = plt.subplots(1, 1,
+        tight_layout=True, 
+        figsize=(8, 4), 
+        )
+
+    # plt.viridis()
     sankey = Sankey(ax=ax, 
-                    scale= 1/power_flows['Optimal Excitation'],
+                    scale= -1/power_flows['Optimal Excitation'],
                     offset= 0,
                     format = '%.1f',
                     shoulder = 0.02,
-                    tolerance=1e-03*power_flows['Optimal Excitation'],
+                    tolerance=-1e-03*power_flows['Optimal Excitation'],
                     unit = 'W'
     )
 
-    sankey.add(flows=[power_flows['Optimal Excitation'],
+    sankey.add(flows=[-1*power_flows['Optimal Excitation'],
                     power_flows['Unused Potential'],
                     power_flows['Actual Excitation']], 
             labels = ['Optimal Excitation', 
@@ -358,9 +366,9 @@ def plot_power_flow(power_flows: dict[str, float])-> Figure:
             diagram.texts[2].set_text('')
 
     plt.axis("off") 
-    plt.show()
+    # plt.show()
 
-    return fig
+    return fig, ax
 
 
 # def add_zerofreq_to_xr(data):