BUG: WEC from impedance with more than one DOF (sandialabs#440)

---------

Co-authored-by: jtgrasb <[email protected]>

Author: cmichelenstrofer

tests/test_integration.py

@@ -24,14 +24,27 @@ def nfreq():
     """Number of frequencies in the array"""
     return 50
 
+
 @pytest.fixture(scope='module')
 def pto():
     """Basic PTO: unstructured, 1 DOF, mechanical power."""
     ndof = 1
     kinematics = np.eye(ndof)
     controller = wot.controllers.unstructured_controller()
-    pto = wot.pto.PTO(ndof=ndof, 
-                      kinematics=kinematics, 
+    pto = wot.pto.PTO(ndof=ndof,
+                      kinematics=kinematics,
+                      controller=controller)
+    return pto
+
+
+@pytest.fixture(scope='module')
+def pto_2d():
+    """Basic PTO: unstructured, 1 DOF, mechanical power."""
+    ndof = 1
+    kinematics = np.array([[1,0],])
+    controller = wot.controllers.unstructured_controller()
+    pto = wot.pto.PTO(ndof=ndof,
+                      kinematics=kinematics,
                       controller=controller)
     return pto
 
@@ -78,6 +91,28 @@ def fb():
     return fb
 
 
+@pytest.fixture(scope="module")
+def fb_2d():
+    """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.2
+    wb = geom.WaveBot()
+    mesh = wb.mesh(mesh_size_factor)
+    mesh_obj = load_from_meshio(mesh, 'WaveBot')
+    lid_mesh = mesh_obj.generate_lid(-2e-2)
+    fb = cpy.FloatingBody(mesh=mesh_obj, lid_mesh=lid_mesh, name="WaveBot")
+    fb.add_translation_dof(name="Heave")
+    fb.add_translation_dof(name="Surge")
+    return fb
+
+
 @pytest.fixture(scope="module")
 def hydro_data(f1, nfreq, fb):
     """Boundary element model (Capytaine) results (with friction added)"""
@@ -88,6 +123,16 @@ def hydro_data(f1, nfreq, fb):
     return hd
 
 
+@pytest.fixture(scope="module")
+def hydro_data_2d(f1, nfreq, fb_2d):
+    """Boundary element model (Capytaine) results (with friction added)"""
+    freq = wot.frequency(f1, nfreq, False)
+    hydro_data = wot.run_bem(fb_2d, freq)
+    hd = wot.add_linear_friction(hydro_data)
+    hd = wot.check_radiation_damping(hd, min_damping=min_damping)
+    return hd
+
+
 @pytest.fixture(scope='module')
 def regular_wave(f1, nfreq):
     """Single frequency wave"""
@@ -111,18 +156,35 @@ def long_crested_wave(f1, nfreq):
 
 
 @pytest.fixture(scope='module')
-def wec_from_bem(f1, nfreq, hydro_data, fb, pto):
+def wec_from_bem(hydro_data, pto):
     """Simple WEC: 1 DOF, no constraints."""
     f_add = {"PTO": pto.force_on_wec}
     wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
     return wec
 
 
+@pytest.fixture(scope='module')
+def wec_from_bem_2d(hydro_data_2d, pto_2d):
+    """Simple WEC: 2 DOF, no constraints."""
+    f_add = {"PTO": pto_2d.force_on_wec}
+    wec = wot.WEC.from_bem(hydro_data_2d, f_add=f_add)
+    return wec
+
+
 @pytest.fixture(scope='module')
 def wec_from_floatingbody(f1, nfreq, fb, pto):
     """Simple WEC: 1 DOF, no constraints."""
     f_add = {"PTO": pto.force_on_wec}
-    wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add, 
+    wec = wot.WEC.from_floating_body(fb, f1, nfreq, f_add=f_add,
+                                     min_damping=min_damping)
+    return wec
+
+
+@pytest.fixture(scope='module')
+def wec_from_floatingbody_2d(f1, nfreq, fb_2d, pto_2d):
+    """Simple WEC: 2 DOF, no constraints."""
+    f_add = {"PTO": pto_2d.force_on_wec}
+    wec = wot.WEC.from_floating_body(fb_2d, f1, nfreq, f_add=f_add,
                                      min_damping=min_damping)
     return wec
 
@@ -140,14 +202,39 @@ def wec_from_impedance(hydro_data, pto, fb):
     fb = fb.immersed_part()
     mass = bemc['inertia_matrix'].values
     hstiff = bemc['hydrostatic_stiffness'].values
-    K = np.expand_dims(hstiff, 2)
-    
+    K = np.expand_dims(hstiff, 0)
+
     freqs = omega / (2 * np.pi)
     impedance = (A + mass)*(1j*w) + B + K/(1j*w)
     exc_coeff = hydro_data['Froude_Krylov_force'] + hydro_data['diffraction_force']
     f_add = {"PTO": pto.force_on_wec}
 
-    wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add, 
+    wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add,
+                                 min_damping=min_damping)
+    return wec
+
+
+@pytest.fixture(scope='module')
+def wec_from_impedance_2d(hydro_data_2d, pto_2d, fb_2d):
+    """Simple WEC: 2 DOF, no constraints."""
+    bemc = hydro_data_2d.copy()
+    omega = bemc['omega'].values
+    w = np.expand_dims(omega, [1,2])
+    A = bemc['added_mass'].values
+    B = bemc['radiation_damping'].values
+    fb_2d.center_of_mass = [0, 0, 0]
+    fb_2d.rotation_center = fb_2d.center_of_mass
+    fb_2d = fb_2d.immersed_part()
+    mass = bemc['inertia_matrix'].values
+    hstiff = bemc['hydrostatic_stiffness'].values
+    K = np.expand_dims(hstiff, 0)
+
+    freqs = omega / (2 * np.pi)
+    impedance = (A + mass)*(1j*w) + B + K/(1j*w)
+    exc_coeff = hydro_data_2d['Froude_Krylov_force'] + hydro_data_2d['diffraction_force']
+    f_add = {"PTO": pto_2d.force_on_wec}
+
+    wec = wot.WEC.from_impedance(freqs, impedance, exc_coeff, hstiff, f_add,
                                  min_damping=min_damping)
     return wec
 
@@ -222,7 +309,28 @@ def test_same_wec_init(wec_from_bem,
     bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
     fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
     imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
-    
+
+    assert fb_res == approx(bem_res, rel=0.01)
+    assert imp_res == approx(bem_res, rel=0.01)
+
+
+def test_same_wec_init_2d(wec_from_bem_2d,
+                       wec_from_floatingbody_2d,
+                       wec_from_impedance_2d,
+                       f1,
+                       nfreq):
+    """Test that different init methods for WEC class produce the same object
+    """
+
+    waves = wot.waves.regular_wave(f1, nfreq, 0.3, 0.0625)
+    np.random.seed(0)
+    x_wec_0 = np.random.randn(wec_from_bem_2d.nstate_wec)
+    np.random.seed(1)
+    x_opt_0 = np.random.randn(wec_from_bem_2d.nstate_wec)
+    bem_res = wec_from_bem_2d.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
+    fb_res = wec_from_floatingbody_2d.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
+    imp_res = wec_from_impedance_2d.residual(x_wec_0, x_opt_0, waves.sel(realization=0))
+
     assert fb_res == approx(bem_res, rel=0.01)
     assert imp_res == approx(bem_res, rel=0.01)
 
@@ -246,7 +354,7 @@ def hydro_impedance(self, hydro_data):
         """Intrinsic hydrodynamic impedance"""
         Zi = wot.hydrodynamic_impedance(hydro_data)
         return Zi
-    
+
     @pytest.fixture(scope='class')
     def unstruct_wec(self,
                      hydro_data,
@@ -265,7 +373,7 @@ def long_crested_wave_unstruct_res(self,
                                        pto,
                                        hydro_data,
                                        nfreq):
-        """Solution for an unstructured controller with multiple long crested 
+        """Solution for an unstructured controller with multiple long crested
         waves"""
 
         f_add = {"PTO": pto.force_on_wec}
@@ -289,7 +397,7 @@ def test_p_controller_resonant_wave(self,
                                         hydro_impedance):
         """Proportional controller should match optimum for natural resonant
         wave"""
-        
+
         f_add = {"PTO": p_controller_pto.force_on_wec}
         wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
 
@@ -346,7 +454,7 @@ def test_pi_controller_regular_wave(self,
                          ).squeeze().sum('omega').item()
 
         assert power_sol == approx(power_optimal, rel=1e-4)
-        
+
     def test_unstructured_controller_long_crested_wave(self,
                                                        unstruct_wec,
                                                        long_crested_wave,
@@ -358,8 +466,8 @@ def test_unstructured_controller_long_crested_wave(self,
 
         power_sol = -1*long_crested_wave_unstruct_res[0]['fun']
 
-        res_fd, _ = unstruct_wec.post_process(unstruct_wec, long_crested_wave_unstruct_res, 
-                                        long_crested_wave, 
+        res_fd, _ = unstruct_wec.post_process(unstruct_wec, long_crested_wave_unstruct_res,
+                                        long_crested_wave,
                                         nsubsteps=1)
         Fex = res_fd.force.sel(realization=0,
                                type=['Froude_Krylov', 'diffraction']).sum('type')
@@ -370,7 +478,7 @@ def test_unstructured_controller_long_crested_wave(self,
 
     def test_unconstrained_solutions_multiple_phase_realizations(self,
                                                                  long_crested_wave_unstruct_res):
-        """Solutions for average power with an unstructured controller 
+        """Solutions for average power with an unstructured controller
         (no constraints) match for different phase realizations"""
 
         pow = [res['fun'] for res in long_crested_wave_unstruct_res]
@@ -390,27 +498,27 @@ def test_saturated_pi_controller(self,
         pto = {}
         wec = {}
         nstate_opt = {}
-        
+
         # Constraint
         f_max = 2000.0
 
         nstate_opt['us'] = 2*nfreq
         pto['us'] = pto_tmp
         def const_f_pto(wec, x_wec, x_opt, wave):
-            f = pto['us'].force_on_wec(wec, x_wec, x_opt, wave, 
+            f = pto['us'].force_on_wec(wec, x_wec, x_opt, wave,
                                        nsubsteps=4)
             return f_max - np.abs(f.flatten())
         wec['us'] = wot.WEC.from_bem(hydro_data,
                                      f_add={"PTO": pto['us'].force_on_wec},
                                      constraints=[{'type': 'ineq',
                                                    'fun': const_f_pto, }])
-        
-        
+
+
         ndof = 1
         nstate_opt['pi'] = 2
         # def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1):
-        #     return wot.pto.controller_pi(pto, wec, x_wec, x_opt, waves, 
-        #                                  nsubsteps, 
+        #     return wot.pto.controller_pi(pto, wec, x_wec, x_opt, waves,
+        #                                  nsubsteps,
         #                                  saturation=[-f_max, f_max])
         saturated_pi = wot.controllers.pid_controller(1,True,True,False,
                                                       saturation=[-f_max, f_max])
@@ -420,7 +528,7 @@ def const_f_pto(wec, x_wec, x_opt, wave):
         wec['pi'] = wot.WEC.from_bem(hydro_data,
                                      f_add={"PTO": pto['pi'].force_on_wec},
                                      constraints=[])
-        
+
         x_opt_0 = {'us': np.ones(nstate_opt['us'])*0.1,
                    'pi': [-1e3, 1e4]}
         scale_x_wec = {'us': 1e1,
@@ -431,7 +539,7 @@ def const_f_pto(wec, x_wec, x_opt, wave):
                      'pi': 1e-2}
         bounds_opt = {'us': None,
                       'pi': ((-1e4, 0), (0, 2e4),)}
-        
+
         res = {}
         pto_fdom = {}
         pto_tdom = {}
@@ -447,13 +555,13 @@ def const_f_pto(wec, x_wec, x_opt, wave):
                             optim_options={'maxiter': 200},
                             bounds_opt=bounds_opt[key]
                             )
-            
+
             nsubstep_postprocess = 4
-            pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key], 
-                                                                 res[key], 
-                                                                 regular_wave, 
+            pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key],
+                                                                 res[key],
+                                                                 regular_wave,
                                                                  nsubstep_postprocess)
-        
+
         xr.testing.assert_allclose(
             pto_tdom['pi'].sel(realization=0).power.squeeze().mean('time'),
             pto_tdom['us'].sel(realization=0).power.squeeze().mean('time'),

wecopttool/core.py

@@ -507,7 +507,7 @@ def from_impedance(
             Complex impedance of size :python:`(nfreq, ndof, ndof)`.
         exc_coeff
             Complex excitation transfer function of size
-            :python:`(ndof, nfreq)`.
+            :python:`(nfreq, nwavedir, ndof)`.
         hydrostatic_stiffness
             Linear hydrostatic restoring coefficient of size
             :python:`(ndof, ndof)`.
@@ -553,7 +553,7 @@ def from_impedance(
         omega = freqs * 2*np.pi
         omega0 = np.expand_dims(omega, [1,2])
         transfer_func = impedance * (1j*omega0)
-        transfer_func0 = np.expand_dims(hydrostatic_stiffness, 2)
+        transfer_func0 = np.expand_dims(hydrostatic_stiffness, 0)
         transfer_func = np.concatenate([transfer_func0, transfer_func], axis=0)
         transfer_func = -1 * transfer_func  # RHS of equation: ma = Σf
         force_impedance = force_from_rao_transfer_function(transfer_func)
@@ -586,9 +586,9 @@ def residual(self, x_wec: ndarray, x_opt: ndarray, wave: DataArray,
         x_opt
             Optimization (control) state.
         wave
-            2D :py:class:`xarray.DataArray` containing the wave's complex 
-            amplitude for a single realization as a function of wave 
-            angular frequency :python:`omega` (rad/s) and direction 
+            2D :py:class:`xarray.DataArray` containing the wave's complex
+            amplitude for a single realization as a function of wave
+            angular frequency :python:`omega` (rad/s) and direction
             :python:`wave_direction` (rad).
         """
         if not self.inertia_in_forces:
@@ -2032,8 +2032,8 @@ def force_from_wave(force_coeff: DataArray,
     Parameters
     ----------
     force_coeff
-        Complex excitation coefficients indexed by frequency and
-        direction angle.
+        Complex excitation coefficients indexed by frequency,
+        direction angle, and degree of freedom.
     """
     def force(wec, x_wec, x_opt, wave):
         force_fd = complex_to_real(wave_excitation(force_coeff, wave), False)
@@ -2197,7 +2197,7 @@ def run_bem(
                      'will be calculated assuming a solid and constant ' +
                      'density body.')
         else:
-            _log.warning('FloatingBody has no inertia_matrix field. ' + 
+            _log.warning('FloatingBody has no inertia_matrix field. ' +
                      'The FloatingBody mass is defined and will be ' +
                      'used for calculating the inertia matrix.')
         wec_im.inertia_matrix = wec_im.compute_rigid_body_inertia(rho=rho)
@@ -2289,12 +2289,12 @@ def wave_excitation(exc_coeff: DataArray, wave: DataArray) -> ndarray:
     Parameters
     ----------
     exc_coeff
-        Complex excitation coefficients indexed by frequency and
-        direction angle.
+        Complex excitation coefficients indexed by frequency,
+        direction angle, and degree of freedom.
     wave
-        2D :py:class:`xarray.DataArray` containing the wave's complex 
-        amplitude for a single realization as a function of wave 
-        angular frequency :python:`omega` (rad/s) and direction 
+        2D :py:class:`xarray.DataArray` containing the wave's complex
+        amplitude for a single realization as a function of wave
+        angular frequency :python:`omega` (rad/s) and direction
         :python:`wave_direction` (rad).
 
     Raises
@@ -2607,8 +2607,8 @@ def set_fb_centers(
                 log_str = (
                     "Using the center of gravity (COG) as the rotation center " +
                     "for hydrostatics. Note that the hydrostatics do not use the " +
-                    "axes defined by the FloatingBody degrees of freedom, and the " + 
-                    "rotation center should be set manually when using Capytaine to " + 
+                    "axes defined by the FloatingBody degrees of freedom, and the " +
+                    "rotation center should be set manually when using Capytaine to " +
                     "calculate hydrostatics about an axis other than the COG.")
             setattr(fb, property, def_val)
             _log.warning(log_str)