diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb index 81548bd5..065ae977 100644 --- a/examples/tutorial_3_LUPA.ipynb +++ b/examples/tutorial_3_LUPA.ipynb @@ -1000,7 +1000,7 @@ " fp = 1/tp\n", " spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n", " efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n", - " return wot.waves.long_crested_wave(efth,nrealizations=5)\n", + " return wot.waves.long_crested_wave(efth,nrealizations)\n", "\n", "for case, data in wave_cases.items():\n", " waves[case] = irregular_wave(data['Hs'], data['Tp'])" diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb index 3f0e7a8d..a817b7a7 100644 --- a/examples/tutorial_4_Pioneer.ipynb +++ b/examples/tutorial_4_Pioneer.ipynb @@ -90,7 +90,7 @@ "fp = 1/Tp\n", "spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, Hs)\n", "efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n", - "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations = nrealizations)" + "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations)" ] }, { diff --git a/tests/test_waves.py b/tests/test_waves.py index 87403c74..42d17303 100644 --- a/tests/test_waves.py +++ b/tests/test_waves.py @@ -190,7 +190,7 @@ def elevation(self, ndbc_omnidirectional, direction, nrealizations): """Complex sea state elevation amplitude [m] indexed by frequency and direction.""" elev = wot.waves.long_crested_wave( - ndbc_omnidirectional.efth, direction, nrealizations) + ndbc_omnidirectional.efth, nrealizations, direction) return elev @pytest.fixture(scope="class") @@ -244,7 +244,7 @@ def test_direction(self, elevation, direction): dir_out = elevation.wave_direction.values.item() assert np.isclose(dir_out, wot.degrees_to_radians(direction)) - def test_realizations(self, elevation, direction): + def test_realizations(self, elevation): """Test that the number of realizations is correct.""" realization_out = elevation.realization.values assert (realization_out == [0,1]).all() @@ -259,7 +259,7 @@ def test_time_series(self, pm_spectrum, pm_f1, pm_nfreq): # create time-series direction = 0.0 nrealizations = 1 - wave = wot.waves.long_crested_wave(pm_spectrum, direction, nrealizations) + wave = wot.waves.long_crested_wave(pm_spectrum, nrealizations, direction) wave_ts = wot.fd_to_td(wave.sel(realization=0).values, pm_f1, pm_nfreq, False) # calculate the spectrum from the time-series t = wot.time(pm_f1, pm_nfreq) @@ -289,12 +289,17 @@ def ndbc_spectrum(self,): files = [f'41013{i}2020.txt' for i in markers] spec = ws.read_ndbc_ascii([os.path.join(dir, file) for file in files]) return spec.sel(time=time).interp(freq=freq) + + @pytest.fixture(scope="class") + def nrealizations(self): + """Number of wave realizations.""" + return 2 @pytest.fixture(scope="class") - def elevation(self, ndbc_spectrum): + def elevation(self, ndbc_spectrum, nrealizations): """Complex sea state elevation amplitude [m] indexed by frequency and direction.""" - return wot.waves.irregular_wave(ndbc_spectrum.efth) + return wot.waves.irregular_wave(ndbc_spectrum.efth, nrealizations) def test_coordinates(self, elevation): """Test that the elevation dataArray has the correct @@ -304,15 +309,20 @@ def test_coordinates(self, elevation): for icoord in coordinates: assert icoord in elevation.coords, f'missing coordinate: {icoord}' - def test_shape(self, ndbc_spectrum, elevation): + def test_shape(self, ndbc_spectrum, elevation, nrealizations): """Test that the elevation dataArray has the correct shape.""" nfreq = len(ndbc_spectrum.freq) ndir = len(ndbc_spectrum.dir) - assert np.squeeze(elevation.values).shape == (nfreq, ndir) + assert np.squeeze(elevation.values).shape == (nfreq, ndir, nrealizations) def test_type(self, elevation): """Test that the elevation dataArray has the correct type.""" assert np.iscomplexobj(elevation) + + def test_realizations(self, elevation): + """Test that the number of realizations is correct.""" + realization_out = elevation.realization.values + assert (realization_out == [0,1]).all() class TestRandomPhase: @@ -323,7 +333,8 @@ def shape(self,): """Shape of the phase matrix, randomized each time the test is run. """ - return (np.random.randint(10, 100), np.random.randint(10, 100)) + return (np.random.randint(10, 100), np.random.randint(10, 100), + np.random.randint(10, 100)) @pytest.fixture(scope="class") def phase_mat(self, shape): diff --git a/wecopttool/core.py b/wecopttool/core.py index 11514813..890bcb5f 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -683,11 +683,13 @@ def solve(self, obj_fun=pto.average_power, nstate_opt=2*nfreq+1) - To get the post-processed results for the - :py:class:`wecopttool.pto.PTO`, you may call + To get the post-processed results for the :py:class:`wecopttool.WEC` + and :py:class:`wecopttool.pto.PTO` for a single realization, you + may call - >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt[0]) - >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[0]) + >>> realization = 0 # realization index + >>> res_wec_fd, res_wec_td = wec.post_process(wec,res_opt[realization]) + >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt[realization]) See Also -------- @@ -712,41 +714,13 @@ def solve(self, # composite scaling vector scale = np.concatenate([scale_x_wec, scale_x_opt]) - - # objective function - sign = -1.0 if maximize else 1.0 - - def obj_fun_scaled(x): - x_wec, x_opt = self.decompose_state(x/scale) - return obj_fun(self, x_wec, x_opt, wave)*scale_obj*sign - - # constraints - constraints = self.constraints.copy() - - for i, icons in enumerate(self.constraints): - icons_new = {"type": icons["type"]} - - def make_new_fun(icons): - def new_fun(x): - x_wec, x_opt = self.decompose_state(x/scale) - return icons["fun"](self, x_wec, x_opt, wave) - return new_fun - - icons_new["fun"] = make_new_fun(icons) - if use_grad: - icons_new['jac'] = jacobian(icons_new['fun']) - constraints[i] = icons_new - - # system dynamics through equality constraint, ma - Σf = 0 - def scaled_resid_fun(x): - x_s = x/scale - x_wec, x_opt = self.decompose_state(x_s) - return self.residual(x_wec, x_opt, wave) - - eq_cons = {'type': 'eq', 'fun': scaled_resid_fun} - if use_grad: - eq_cons['jac'] = jacobian(scaled_resid_fun) - constraints.append(eq_cons) + + # decision variable initial guess + if x_wec_0 is None: + x_wec_0 = np.random.randn(self.nstate_wec) + if x_opt_0 is None: + x_opt_0 = np.random.randn(nstate_opt) + x0 = np.concatenate([x_wec_0, x_opt_0])*scale # bounds if (bounds_wec is None) and (bounds_opt is None): @@ -771,32 +745,60 @@ def scaled_resid_fun(x): bounds = Bounds(lb=np.hstack([le.lb for le in bounds_list])*scale, ub=np.hstack([le.ub for le in bounds_list])*scale) - # callback - if callback is None: - def callback_scipy(x): - x_wec, x_opt = self.decompose_state(x) - max_x_opt = np.nan if np.size(x_opt)==0 else np.max(np.abs(x_opt)) - _log.info("Scaled [max(x_wec), max(x_opt), obj_fun(x)]: " - + f"[{np.max(np.abs(x_wec)):.2e}, " - + f"{max_x_opt:.2e}, " - + f"{obj_fun_scaled(x):.2e}]") - else: - def callback_scipy(x): - x_s = x/scale - x_wec, x_opt = self.decompose_state(x_s) - return callback(self, x_wec, x_opt, wave) - for realization, wave in waves.groupby('realization'): _log.info("Solving pseudo-spectral control problem " + f"for realization number {realization}.") + + # objective function + sign = -1.0 if maximize else 1.0 + + def obj_fun_scaled(x): + x_wec, x_opt = self.decompose_state(x/scale) + return obj_fun(self, x_wec, x_opt, wave)*scale_obj*sign + + # constraints + constraints = self.constraints.copy() + + for i, icons in enumerate(self.constraints): + icons_new = {"type": icons["type"]} + + def make_new_fun(icons): + def new_fun(x): + x_wec, x_opt = self.decompose_state(x/scale) + return icons["fun"](self, x_wec, x_opt, wave) + return new_fun - # decision variable initial guess - if x_wec_0 is None: - x_wec_0 = np.random.randn(self.nstate_wec) - if x_opt_0 is None: - x_opt_0 = np.random.randn(nstate_opt) - x0 = np.concatenate([x_wec_0, x_opt_0])*scale + icons_new["fun"] = make_new_fun(icons) + if use_grad: + icons_new['jac'] = jacobian(icons_new['fun']) + constraints[i] = icons_new + + # system dynamics through equality constraint, ma - Σf = 0 + def scaled_resid_fun(x): + x_s = x/scale + x_wec, x_opt = self.decompose_state(x_s) + return self.residual(x_wec, x_opt, wave) + + eq_cons = {'type': 'eq', 'fun': scaled_resid_fun} + if use_grad: + eq_cons['jac'] = jacobian(scaled_resid_fun) + constraints.append(eq_cons) + + # callback + if callback is None: + def callback_scipy(x): + x_wec, x_opt = self.decompose_state(x) + max_x_opt = np.nan if np.size(x_opt)==0 else np.max(np.abs(x_opt)) + _log.info("Scaled [max(x_wec), max(x_opt), obj_fun(x)]: " + + f"[{np.max(np.abs(x_wec)):.2e}, " + + f"{max_x_opt:.2e}, " + + f"{obj_fun_scaled(x):.2e}]") + else: + def callback_scipy(x): + x_s = x/scale + x_wec, x_opt = self.decompose_state(x_s) + return callback(self, x_wec, x_opt, wave) # optimization problem optim_options['disp'] = optim_options.get('disp', True) diff --git a/wecopttool/waves.py b/wecopttool/waves.py index 9f45e9ac..97ddc1a2 100644 --- a/wecopttool/waves.py +++ b/wecopttool/waves.py @@ -95,7 +95,7 @@ def elevation_fd( freq = frequency(f1, nfreq, False) omega = freq*2*np.pi - dims = ('omega', 'wave_direction','realization') + dims = ('omega', 'wave_direction', 'realization') omega_attr = {'long_name': 'Radial frequency', 'units': 'rad/s'} freq_attr = {'long_name': 'Frequency', 'units': 'Hz'} dir_attr = {'long_name': 'Wave direction', 'units': 'rad'} @@ -107,11 +107,17 @@ def elevation_fd( if amplitudes is None: amplitudes = np.zeros([nfreq, ndirections, nrealizations]) + else: + if amplitudes.shape == (nfreq, ndirections): + amplitudes = np.expand_dims(amplitudes,axis=2) + assert amplitudes.shape == (nfreq, ndirections, 1) or \ + amplitudes.shape == (nfreq, ndirections, nrealizations) if phases is None: - phases = random_phase([nfreq, ndirections, nrealizations],seed) + phases = random_phase([nfreq, ndirections, nrealizations], seed) else: phases = degrees_to_radians(phases, False) + assert phases.shape == (nfreq, ndirections, nrealizations) camplitude = amplitudes * np.exp(1j*phases) @@ -190,8 +196,8 @@ def regular_wave( def long_crested_wave( efth: DataArray, + nrealizations: int, direction: Optional[float] = 0.0, - nrealizations: Optional[float] = 1, seed: Optional[float] = None, ) -> DataArray: """Create a complex frequency-domain wave elevation from an @@ -210,11 +216,11 @@ def long_crested_wave( efth Omnidirection wave spectrum in units of m^2/Hz, in the format used by :py:class:`wavespectra.SpecArray`. - direction - Direction (in degrees) of the long-crested wave. nrealizations Number of wave phase realizations to be created for the long-crested wave. + direction + Direction (in degrees) of the long-crested wave. seed Seed for random number generator. Used for reproducibility. Generally should not be used except for testing. @@ -225,7 +231,6 @@ def long_crested_wave( values = efth.values values[values<0] = np.nan amplitudes = np.sqrt(2 * values * df) - amplitudes = np.expand_dims(amplitudes,axis=2) attr = { 'Wave type': 'Long-crested irregular', @@ -236,7 +241,7 @@ def long_crested_wave( def irregular_wave(efth: DataArray, - nrealizations: Optional[float] = 1, + nrealizations: int, seed: Optional[float] = None,) -> DataArray: """Create a complex frequency-domain wave elevation from a spectrum. @@ -270,7 +275,6 @@ def irregular_wave(efth: DataArray, values = efth.values values[values<0] = np.nan amplitudes = np.sqrt(2 * values * df * dd) - amplitudes = np.expand_dims(amplitudes,axis=2) attr = {'Wave type': 'Irregular'}