Skip to content

Commit

Permalink
actually test power solution for irregular wave (sandialabs#327)
Browse files Browse the repository at this point in the history
  • Loading branch information
@ryancoe
ryancoe authored Apr 3, 2024
1 parent 775d273 commit a7bef5f
Showing 1 changed file with 50 additions and 40 deletions.
90 changes: 50 additions & 40 deletions tests/test_integration.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@


kplim = -1e1

min_damping = 45

@pytest.fixture(scope="module")
def f1():
Expand Down Expand Up @@ -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')
Expand All @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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,
Expand All @@ -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,
Expand All @@ -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),
Expand All @@ -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):
Expand All @@ -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, }])
Expand All @@ -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=[])

Expand Down

0 comments on commit a7bef5f

Please sign in to comment.