Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Oct 9, 2025
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
1e05c9c
fix bug: wrong index
cmichelenstrofer Sep 30, 2025
1baf54c
docstrings
cmichelenstrofer Sep 30, 2025
4d798f8
tests
cmichelenstrofer Sep 30, 2025
c39f305
tests
cmichelenstrofer Sep 30, 2025
88f6bef
Fix kinematics matrix
jtgrasb Oct 9, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
164 changes: 136 additions & 28 deletions tests/test_integration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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)"""
Expand All @@ -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"""
Expand All @@ -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

Expand All @@ -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

Expand Down Expand Up @@ -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)

Expand All @@ -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,
Expand All @@ -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}
Expand All @@ -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)

Expand Down Expand Up @@ -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,
Expand All @@ -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')
Expand All @@ -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]
Expand All @@ -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])
Expand All @@ -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,
Expand All @@ -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 = {}
Expand All @@ -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'),
Expand Down
30 changes: 15 additions & 15 deletions wecopttool/core.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)`.
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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:
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down
Loading