Frequency Domain Linear Solver (#419)

For single DOF PTO

Author: cmichelenstrofer

tests/test_utilities.py

@@ -70,7 +70,7 @@ def bem_data(f1, nfreq, ndof, ndir):
     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),
@@ -160,24 +160,24 @@ def test_plot_hydrodynamic_coefficients(bem_data,ndof):
     assert isinstance(fig_ex,Figure)
 
 def test_plot_bode_impedance(intrinsic_impedance, ndof):
-    fig_Zi, axes_Zi = wot.utilities.plot_bode_impedance(intrinsic_impedance)    
-  
+    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)    
-  
+    fig_sankey, ax_sankey = wot.utilities.plot_power_flow(power_flows)
+
     assert isinstance(fig_sankey, Figure)
-    assert isinstance(ax_sankey, Axes) 
+    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, 
+    """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"""
 
@@ -196,14 +196,39 @@ def test_calculate_power_flow(wb_bem,
                     bounds_opt=((-1e4, 0), (0, 2e4),)
                     )
 
-    pflows = wot.utilities.calculate_power_flows(wec, 
-                          pi_controller_pto, 
-                          res, 
-                          regular_wave, 
+    pflows = wot.utilities.calculate_power_flows(wec,
+                          pi_controller_pto,
+                          res,
+                          regular_wave,
                           wb_hydro_impedance)
 
     assert pflows['Absorbed'] == approx(pflows['Radiated'], rel=1e-4)
     assert pflows['Optimal Excitation'] == approx(pflows['Actual Excitation'], rel=1e-4)
 
-
-
+def test_linear_solve(wb_bem, regular_wave):
+    omega = wb_bem.omega.values
+    gear_ratio = 12.0
+    torque_constant = 6.7
+    winding_resistance = 0.5
+    winding_inductance = 0.0
+    drivetrain_inertia = 2.0
+    drivetrain_friction = 1.0
+    drivetrain_stiffness = 0.0
+
+    drivetrain_impedance = (1j*omega*drivetrain_inertia +
+                            drivetrain_friction +
+                            1/(1j*omega)*drivetrain_stiffness)
+
+    winding_impedance = winding_resistance + 1j*omega*winding_inductance
+
+
+    pto_impedance_11 = -1* gear_ratio**2 * drivetrain_impedance
+    off_diag = np.sqrt(3.0/2.0) * torque_constant * gear_ratio
+    pto_impedance_12 = -1*(off_diag+0j) * np.ones(omega.shape)
+    pto_impedance_21 = -1*(off_diag+0j) * np.ones(omega.shape)
+    pto_impedance_22 = winding_impedance
+    pto_impedance = np.array([[pto_impedance_11, pto_impedance_12],
+                                [pto_impedance_21, pto_impedance_22]])
+
+    power, _, _, _ = wot.utilities.linear_solve(wb_bem, pto_impedance, regular_wave.isel(realization=0), np.eye(1))
+    assert power == approx(-29.2, abs=0.05)

tests/test_waves.py

@@ -250,7 +250,7 @@ def test_time_series(self, pm_spectrum, pm_f1, pm_nfreq):
         direction = 0.0
         nrealizations = 1
         wave = wot.waves.long_crested_wave(pm_spectrum, nrealizations, direction)
-        wave_ts = wot.fd_to_td(wave.sel(realization=0).values, pm_f1, pm_nfreq, False)
+        wave_ts = wot.fd_to_td(wave.sel(realization=0).values, pm_f1, pm_nfreq, 1, False)
         # calculate the spectrum from the time-series
         t = wot.time(pm_f1, pm_nfreq)
         fs = 1/t[1]
@@ -279,7 +279,7 @@ def ndbc_spectrum(self,):
         files = [f'41013{i}2020.txt' for i in markers]
         spec = ws.read_ndbc_ascii([os.path.join(dir, file) for file in files])
         return spec.sel(time=time).interp(freq=freq)
-    
+
     @pytest.fixture(scope="class")
     def nrealizations(self):
         """Number of wave realizations."""
@@ -308,7 +308,7 @@ def test_shape(self, ndbc_spectrum, elevation, nrealizations):
     def test_type(self, elevation):
         """Test that the elevation dataArray has the correct type."""
         assert np.iscomplexobj(elevation)
-    
+
     def test_realizations(self, elevation):
         """Test that the number of realizations is correct."""
         realization_out = elevation.realization.values
@@ -344,4 +344,4 @@ def test_float(self,):
         matrix.
         """
         phase = wot.waves.random_phase()
-        assert (phase < np.pi) and (phase >= -np.pi)
+        assert (phase < np.pi) and (phase >= -np.pi)

wecopttool/core.py

@@ -1278,7 +1278,7 @@ def fd_to_td(self, fd: ndarray) -> ndarray:
         --------
         fd_to_td, WEC.td_to_fd
         """
-        return fd_to_td(fd, self.f1, self.nfreq, True)
+        return fd_to_td(fd, self.f1, self.nfreq, 1, True)
 
     def td_to_fd(
         self,
@@ -1718,6 +1718,7 @@ def fd_to_td(
     fd: ArrayLike,
     f1: Optional[float] = None,
     nfreq: Optional[int] = None,
+    nsubsteps: int = 1,
     zero_freq: Optional[bool] = True,
 ) -> ndarray:
     """Convert a complex array of Fourier coefficients to a real array
@@ -1753,6 +1754,9 @@ def fd_to_td(
         Fundamental frequency :python:`f1` [:math:`Hz`].
     nfreq
         Number of frequencies.
+    nsubsteps
+        Number of steps between the default (implied) time steps.
+        A value of :python:`1` corresponds to the default step length.
     zero_freq
         Whether the mean (DC) component is included.
 
@@ -1773,7 +1777,7 @@ def fd_to_td(
         assert np.allclose(np.imag(fd[0, :]), 0), msg
 
     if (f1 is not None) and (nfreq is not None):
-        tmat = time_mat(f1, nfreq, zero_freq=zero_freq)
+        tmat = time_mat(f1, nfreq, nsubsteps, zero_freq=zero_freq)
         td = tmat @ complex_to_real(fd, zero_freq)
     elif (f1 is None) and (nfreq is None):
         n = 2*(fd.shape[0]-1)