From 1e05c9c045319dab2fe21be507e07fd882151114 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20A=2E=20Michel=C3=A9n=20Str=C3=B6fer?= Date: Tue, 30 Sep 2025 13:09:57 -0600 Subject: [PATCH 1/5] fix bug: wrong index --- wecopttool/core.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/wecopttool/core.py b/wecopttool/core.py index 43e7b15dc..f0a03dc9b 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -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: @@ -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) @@ -2292,9 +2292,9 @@ def wave_excitation(exc_coeff: DataArray, wave: DataArray) -> ndarray: Complex excitation coefficients indexed by frequency and direction angle. 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) From 1baf54c61698f39bf2a0821d48c787ea864d49d6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20A=2E=20Michel=C3=A9n=20Str=C3=B6fer?= Date: Tue, 30 Sep 2025 13:18:09 -0600 Subject: [PATCH 2/5] docstrings --- wecopttool/core.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/wecopttool/core.py b/wecopttool/core.py index f0a03dc9b..9672ac761 100644 --- a/wecopttool/core.py +++ b/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)`. @@ -2033,7 +2033,7 @@ def force_from_wave(force_coeff: DataArray, ---------- force_coeff Complex excitation coefficients indexed by frequency and - direction angle. + 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) @@ -2289,8 +2289,8 @@ 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 From 4d798f805197df6816d90a30ef12b94476d36d36 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20A=2E=20Michel=C3=A9n=20Str=C3=B6fer?= Date: Tue, 30 Sep 2025 13:36:49 -0600 Subject: [PATCH 3/5] tests --- tests/test_integration.py | 164 +++++++++++++++++++++++++++++++------- 1 file changed, 136 insertions(+), 28 deletions(-) diff --git a/tests/test_integration.py b/tests/test_integration.py index 7d33f410a..c9bb152ed 100644 --- a/tests/test_integration.py +++ b/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['Froude_Krylov_force'] + hydro_data['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'), From c39f3058ea7c6f8bcc0cc058ecd19d7192172906 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Carlos=20A=2E=20Michel=C3=A9n=20Str=C3=B6fer?= Date: Tue, 30 Sep 2025 14:28:51 -0600 Subject: [PATCH 4/5] tests --- tests/test_integration.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_integration.py b/tests/test_integration.py index c9bb152ed..19a74f12f 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -231,7 +231,7 @@ def wec_from_impedance_2d(hydro_data_2d, pto_2d, fb_2d): 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'] + 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, From 88f6bef55611ffe1998b0d526b2d1e1749512d1d Mon Sep 17 00:00:00 2001 From: jtgrasb <87095491+jtgrasb@users.noreply.github.com> Date: Thu, 9 Oct 2025 11:08:04 -0600 Subject: [PATCH 5/5] Fix kinematics matrix --- tests/test_integration.py | 2 +- wecopttool/core.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/test_integration.py b/tests/test_integration.py index 19a74f12f..aba4cf58c 100644 --- a/tests/test_integration.py +++ b/tests/test_integration.py @@ -41,7 +41,7 @@ def pto(): def pto_2d(): """Basic PTO: unstructured, 1 DOF, mechanical power.""" ndof = 1 - kinematics = np.array([1,0]) + kinematics = np.array([[1,0],]) controller = wot.controllers.unstructured_controller() pto = wot.pto.PTO(ndof=ndof, kinematics=kinematics, diff --git a/wecopttool/core.py b/wecopttool/core.py index 9672ac761..f78c2dab5 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -2032,8 +2032,8 @@ def force_from_wave(force_coeff: DataArray, Parameters ---------- force_coeff - Complex excitation coefficients indexed by frequency and - direction angle and degree of freedom. + 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)