actually test power solution for irregular wave

tests/test_integration.py

@@ -10,7 +10,7 @@
 
 
 kplim = -1e1
-
+min_damping = 45
 
 @pytest.fixture(scope="module")
 def f1():
@@ -73,10 +73,13 @@ def fb():
 
 
 @pytest.fixture(scope="module")
-def bem(f1, nfreq, fb):
-    """Boundary elemement model (Capytaine) results"""
+def hydro_data(f1, nfreq, fb):
+    """Boundary element model (Capytaine) results (with friction added)"""
     freq = wot.frequency(f1, nfreq, False)
-    return wot.run_bem(fb, freq)
+    hydro_data = wot.run_bem(fb, freq)
+    hd = wot.add_linear_friction(hydro_data)
+    hd = wot.check_radiation_damping(hd, min_damping=min_damping)
+    return hd
 
 
 @pytest.fixture(scope='module')
@@ -91,36 +94,42 @@ def regular_wave(f1, nfreq):
 
 
 @pytest.fixture(scope='module')
-def irregular_wave(f1, nfreq):
+def long_crested_wave(f1, nfreq):
     """Idealized (Pierson-Moskowitz) spectrum wave"""
     freq = wot.frequency(f1, nfreq, False)
     fp = 0.3
     hs = 0.0625*1.9
-    spec = wot.waves.pierson_moskowitz_spectrum(freq, fp, hs)
-    waves = wot.waves.long_crested_wave(spec)
+    spec_fun = lambda f: wot.waves.pierson_moskowitz_spectrum(freq=f, 
+                                                              fp=fp, 
+                                                              hs=hs)
+    efth = wot.waves.omnidirectional_spectrum(f1=f1, nfreq=nfreq, 
+                                              spectrum_func=spec_fun,
+                                              )
+    waves = wot.waves.long_crested_wave(efth, nrealizations=1)
     return waves
 
 
 @pytest.fixture(scope='module')
-def wec_from_bem(f1, nfreq, bem, fb, pto):
+def wec_from_bem(f1, nfreq, hydro_data, fb, pto):
     """Simple WEC: 1 DOF, no constraints."""
     f_add = {"PTO": pto.force_on_wec}
-    wec = wot.WEC.from_bem(bem, f_add=f_add)
+    wec = wot.WEC.from_bem(hydro_data, 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_impedance(bem, pto, fb):
+def wec_from_impedance(hydro_data, pto, fb):
     """Simple WEC: 1 DOF, no constraints."""
-    bemc = bem.copy()
+    bemc = hydro_data.copy()
     omega = bemc['omega'].values
     w = np.expand_dims(omega, [1,2])
     A = bemc['added_mass'].values
@@ -134,17 +143,18 @@ def wec_from_impedance(bem, pto, fb):
 
     freqs = omega / (2 * np.pi)
     impedance = (A + mass)*(1j*w) + B + K/(1j*w)
-    exc_coeff = bem['Froude_Krylov_force'] + bem['diffraction_force']
+    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 resonant_wave(f1, nfreq, fb, bem):
+def resonant_wave(f1, nfreq, fb, hydro_data):
     """Regular wave at natural frequency of the WEC"""
-    hd = wot.add_linear_friction(bem)
+    hd = wot.add_linear_friction(hydro_data)
     Zi = wot.hydrodynamic_impedance(hd)
     wn = Zi['omega'][np.abs(Zi).argmin(dim='omega')].item()
     waves = wot.waves.regular_wave(f1, nfreq, freq=wn/2/np.pi, amplitude=0.1)
@@ -231,23 +241,21 @@ def hstiff(self, fb):
         return fb.compute_hydrostatic_stiffness()
 
     @pytest.fixture(scope='class')
-    def hydro_impedance(self, bem):
+    def hydro_impedance(self, hydro_data):
         """Intrinsic hydrodynamic impedance"""
-        hd = wot.add_linear_friction(bem)
-        hd = wot.check_radiation_damping(hd)
-        Zi = wot.hydrodynamic_impedance(hd)
+        Zi = wot.hydrodynamic_impedance(hydro_data)
         return Zi
 
     def test_p_controller_resonant_wave(self,
-                                        bem,
+                                        hydro_data,
                                         resonant_wave,
                                         p_controller_pto,
                                         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(bem, f_add=f_add)
+        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
 
         res = wec.solve(waves=resonant_wave,
                         obj_fun=p_controller_pto.average_power,
@@ -272,14 +280,14 @@ def test_p_controller_resonant_wave(self,
         assert power_sol == approx(power_optimal, rel=0.03)
 
     def test_pi_controller_regular_wave(self,
-                                        bem,
+                                        hydro_data,
                                         regular_wave,
                                         pi_controller_pto,
                                         hydro_impedance):
         """PI controller matches optimal for any regular wave"""
 
         f_add = {"PTO": pi_controller_pto.force_on_wec}
-        wec = wot.WEC.from_bem(bem, f_add=f_add)
+        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
 
         res = wec.solve(waves=regular_wave,
                         obj_fun=pi_controller_pto.average_power,
@@ -302,20 +310,20 @@ def test_pi_controller_regular_wave(self,
                          ).squeeze().sum('omega').item()
 
         assert power_sol == approx(power_optimal, rel=1e-4)
-    def test_unstructured_controller_irregular_wave(self,
-                                                    fb,
-                                                    bem,
-                                                    regular_wave,
-                                                    pto,
-                                                    nfreq,
-                                                    hydro_impedance):
+        
+    def test_unstructured_controller_long_crested_wave(self,
+                                                       hydro_data,
+                                                       long_crested_wave,
+                                                       pto,
+                                                       nfreq,
+                                                       hydro_impedance):
         """Unstructured (numerical optimal) controller matches optimal for any
-        irregular wave when unconstrained"""
+        irregular (long crested) wave when unconstrained"""
 
         f_add = {"PTO": pto.force_on_wec}
-        wec = wot.WEC.from_bem(bem, f_add=f_add)
+        wec = wot.WEC.from_bem(hydro_data, f_add=f_add)
 
-        res = wec.solve(waves=regular_wave,
+        res = wec.solve(waves=long_crested_wave,
                         obj_fun=pto.average_power,
                         nstate_opt=2*nfreq,
                         x_wec_0=1e-1*np.ones(wec.nstate_wec),
@@ -326,16 +334,18 @@ def test_unstructured_controller_irregular_wave(self,
 
         power_sol = -1*res[0]['fun']
 
-        res_fd, _ = wec.post_process(res[0], regular_wave.sel(realization=0), nsubsteps=1)
+        res_fd, _ = wec.post_process(res[0], 
+                                        long_crested_wave.sel(realization=0), 
+                                        nsubsteps=1)
         Fex = res_fd.force.sel(
             type=['Froude_Krylov', 'diffraction']).sum('type')
         power_optimal = (np.abs(Fex)**2/8 / np.real(hydro_impedance.squeeze())
-                         ).squeeze().sum('omega').item()
+                            ).squeeze().sum('omega').item()
 
-        assert power_sol == approx(power_optimal, rel=1e-3)
+        assert power_sol == approx(power_optimal, rel=1e-2)
 
     def test_saturated_pi_controller(self,
-                                    bem,
+                                    hydro_data,
                                     regular_wave,
                                     pto,
                                     nfreq):
@@ -357,7 +367,7 @@ def const_f_pto(wec, x_wec, x_opt, waves):
             f = pto['us'].force_on_wec(wec, x_wec, x_opt, waves, 
                                        nsubsteps=4)
             return f_max - np.abs(f.flatten())
-        wec['us'] = wot.WEC.from_bem(bem,
+        wec['us'] = wot.WEC.from_bem(hydro_data,
                                      f_add={"PTO": pto['us'].force_on_wec},
                                      constraints=[{'type': 'ineq',
                                                    'fun': const_f_pto, }])
@@ -372,7 +382,7 @@ def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1):
         pto['pi'] = wot.pto.PTO(ndof=ndof,
                                 kinematics=np.eye(ndof),
                                 controller=saturated_pi,)
-        wec['pi'] = wot.WEC.from_bem(bem,
+        wec['pi'] = wot.WEC.from_bem(hydro_data,
                                      f_add={"PTO": pto['pi'].force_on_wec},
                                      constraints=[])