From 598e875cd48751ddb652657391205f6e208714f0 Mon Sep 17 00:00:00 2001
From: jtgrasb <87095491+jtgrasb@users.noreply.github.com>
Date: Thu, 7 Dec 2023 18:34:28 -0500
Subject: [PATCH] Add multiple phase realizations to waves (#297)
* Add realization to waves variable and adjust tutorials to handle it
* Update tutorial_3_LUPA.ipynb
* Update wave and core functionality
* Update tests for multiple phases
* Remove demonstration from pioneer
* Update tests
* Test
* Fix tutorials
* Reformat core.py
* Update AH tutorial
* Update LUPA and Pioneer nrealizations and explanation
* Update based on review
---
examples/tutorial_1_WaveBot.ipynb | 20 +--
examples/tutorial_2_AquaHarmonics.ipynb | 16 +--
examples/tutorial_3_LUPA.ipynb | 50 ++++---
examples/tutorial_4_Pioneer.ipynb | 49 ++++---
tests/test_core.py | 27 ++--
tests/test_integration.py | 24 ++--
tests/test_waves.py | 54 +++++---
wecopttool/core.py | 177 +++++++++++++-----------
wecopttool/waves.py | 40 ++++--
9 files changed, 265 insertions(+), 192 deletions(-)
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
index 2f82207e..030077ba 100644
--- a/examples/tutorial_1_WaveBot.ipynb
+++ b/examples/tutorial_1_WaveBot.ipynb
@@ -334,7 +334,7 @@
" scale_obj=scale_obj,\n",
" )\n",
"\n",
- "opt_mechanical_average_power = results.fun\n",
+ "opt_mechanical_average_power = results[0].fun\n",
"print(f'Optimal average mechanical power: {opt_mechanical_average_power} W')"
]
},
@@ -353,8 +353,8 @@
"outputs": [],
"source": [
"nsubsteps = 5\n",
- "pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps)\n",
- "wec_fdom, wec_tdom = wec.post_process(results, waves, nsubsteps=nsubsteps)"
+ "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves.sel(realization=0), nsubsteps=nsubsteps)\n",
+ "wec_fdom, wec_tdom = wec.post_process(results[0], waves.sel(realization=0), nsubsteps=nsubsteps)"
]
},
{
@@ -534,12 +534,12 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj,\n",
")\n",
- "opt_average_power = results_2.fun\n",
+ "opt_average_power = results_2[0].fun\n",
"print(f'Optimal average electrical power: {opt_average_power} W')\n",
"\n",
"# Post-process\n",
- "wec_fdom_2, wec_tdom_2 = wec_2.post_process(results_2, waves, nsubsteps)\n",
- "pto_fdom_2, pto_tdom_2 = pto_2.post_process(wec_2, results_2, waves, nsubsteps)"
+ "wec_fdom_2, wec_tdom_2 = wec_2.post_process(results_2[0], waves.sel(realization=0), nsubsteps)\n",
+ "pto_fdom_2, pto_tdom_2 = pto_2.post_process(wec_2, results_2[0], waves.sel(realization=0), nsubsteps)"
]
},
{
@@ -569,12 +569,12 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj,\n",
")\n",
- "opt_average_power = results_2_nocon.fun\n",
+ "opt_average_power = results_2_nocon[0].fun\n",
"print(f'Optimal average electrical power: {opt_average_power} W')\n",
"wec_fdom_2_nocon, wec_tdom_2_nocon = wec_2_nocon.post_process(\n",
- " results_2_nocon, waves, nsubsteps)\n",
+ " results_2_nocon[0], waves.sel(realization=0), nsubsteps)\n",
"pto_fdom_2_nocon, pto_tdom_2_nocon = pto_2.post_process(\n",
- " wec_2_nocon, results_2_nocon, waves, nsubsteps)"
+ " wec_2_nocon, results_2_nocon[0], waves.sel(realization=0), nsubsteps)"
]
},
{
@@ -784,7 +784,7 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj)\n",
"\n",
- " return res.fun"
+ " return res[0].fun"
]
},
{
diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
index 36a95c44..83df7e2d 100644
--- a/examples/tutorial_2_AquaHarmonics.ipynb
+++ b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -565,7 +565,7 @@
" scale_obj=scale_obj,\n",
")\n",
"\n",
- "print(f'Optimal average power: {results.fun/1e3:.2f} kW')"
+ "print(f'Optimal average power: {results[0].fun/1e3:.2f} kW')"
]
},
{
@@ -583,8 +583,8 @@
"metadata": {},
"outputs": [],
"source": [
- "wec_fdom, wec_tdom = wec.post_process(results, waves, nsubsteps=nsubsteps)\n",
- "pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps)"
+ "wec_fdom, wec_tdom = wec.post_process(results[0], waves.sel(realization=0), nsubsteps=nsubsteps)\n",
+ "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves.sel(realization=0), nsubsteps=nsubsteps)"
]
},
{
@@ -642,7 +642,7 @@
" ax=ax[1], \n",
" label='PTO force in WEC frame',\n",
")\n",
- "x_wec, x_opt = wot.decompose_state(results.x, ndof=ndof, nfreq=nfreq)\n",
+ "x_wec, x_opt = wot.decompose_state(results[0].x, ndof=ndof, nfreq=nfreq)\n",
"ax[1].plot(\n",
" wec.time_nsubsteps(nsubsteps),\n",
" f_pto_line(wec, x_wec, x_opt, waves, nsubsteps)/1e3,\n",
@@ -784,7 +784,7 @@
" scale_x_wec=scale_x_wec,\n",
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj)\n",
- " opt_average_power = res.fun\n",
+ " opt_average_power = res[0].fun\n",
" print(\n",
" f'* Mass: {mass_var:.2f} kg\\n' +\n",
" f'* Optimal average electrical power: {opt_average_power/1e3:.2f} kW' \n",
@@ -792,10 +792,10 @@
"\n",
" # wec_fdom, wec_tdom = wec_mass.post_process(res, waves, nsubsteps=nsubsteps)\n",
" global max_torque \n",
- " x_wec, x_opt = wot.decompose_state(res.x, ndof=ndof, nfreq=nfreq)\n",
- " max_torque.append(np.max(f_pto_line(wec_mass, x_wec, x_opt, waves, nsubsteps)))\n",
+ " x_wec, x_opt = wot.decompose_state(res[0].x, ndof=ndof, nfreq=nfreq)\n",
+ " max_torque.append(np.max(f_pto_line(wec_mass, x_wec, x_opt, waves.sel(realization=0), nsubsteps)))\n",
" \n",
- " return res.fun"
+ " return res[0].fun"
]
},
{
diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
index 10c9dc55..065ae977 100644
--- a/examples/tutorial_3_LUPA.ipynb
+++ b/examples/tutorial_3_LUPA.ipynb
@@ -956,7 +956,11 @@
"For each site/scale they provide four wave conditions to test at the Oregon State Large Wave Flume (LWF): the maximum 90th percentile, maximum percent annual energy, maximum occurrence, and minimum 10th percentile. \n",
"\n",
"In this tutorial we will use the PacWave South conditions scaled to the LWF and will design for the maximum occurrence wave. \n",
- "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum."
+ "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum. \n",
+ "\n",
+ "The irregular wave is created with multiple phase realizations. For this tutorial, the number of realizations has been overwritten to 2 to reduce the total run-time. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes.\n",
+ "\n",
+ "$t_{total} = \\frac{nrealizations}{f1} $"
]
},
{
@@ -974,6 +978,12 @@
"wavedir = 0\n",
"waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n",
"\n",
+ "# number of realizations to reach 20 minutes of total simulation time\n",
+ "minutes_needed = 20\n",
+ "nrealizations = minutes_needed*60*f1\n",
+ "print(f'Number of realizations for a 20 minute total simulation time: {nrealizations}')\n",
+ "nrealizations = 2 # overwrite nrealizations to reduce run-time\n",
+ "\n",
"# irregular wave cases from OSU\n",
"wave_cases = {\n",
" 'south_max_90': {'Hs': 0.21, 'Tp': 3.09}, \n",
@@ -990,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)\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'])"
@@ -1010,21 +1020,21 @@
"outputs": [],
"source": [
"fig, ax = plt.subplots()\n",
- "plt1 = np.abs(waves['south_max_90']).plot(\n",
+ "plt1 = np.abs(waves['south_max_90'].sel(realization=0)).plot(\n",
" ax=ax, color='C0', linestyle='solid', label='PW South, 90th percentile')\n",
- "plt2 = np.abs(waves['south_max_annual']).plot(\n",
+ "plt2 = np.abs(waves['south_max_annual'].sel(realization=0)).plot(\n",
" ax=ax, color='C0', linestyle='dotted', label='PW South, Max Annual')\n",
- "plt3 = np.abs(waves['south_max_occurrence']).plot(\n",
+ "plt3 = np.abs(waves['south_max_occurrence'].sel(realization=0)).plot(\n",
" ax=ax, color='C0', linestyle='dashed', label='PW South, Max Occurrence')\n",
- "plt4 = np.abs(waves['south_min_10']).plot(\n",
+ "plt4 = np.abs(waves['south_min_10'].sel(realization=0)).plot(\n",
" ax=ax, color='C0', linestyle='dashdot', label='PW South, 10th percentile')\n",
- "plt5 = np.abs(waves['north_max_90']).plot(\n",
+ "plt5 = np.abs(waves['north_max_90'].sel(realization=0)).plot(\n",
" ax=ax, color='C1', linestyle='solid', label='PW North, 90th percentile')\n",
- "plt6 = np.abs(waves['north_max_annual']).plot(\n",
+ "plt6 = np.abs(waves['north_max_annual'].sel(realization=0)).plot(\n",
" ax=ax, color='C1', linestyle='dotted', label='PW North, Max Annual')\n",
- "plt7 = np.abs(waves['north_max_occurrence']).plot(\n",
+ "plt7 = np.abs(waves['north_max_occurrence'].sel(realization=0)).plot(\n",
" ax=ax, color='C1', linestyle='dashed', label='PW North, Max Occurrence')\n",
- "plt8 = np.abs(waves['north_min_10']).plot(\n",
+ "plt8 = np.abs(waves['north_min_10'].sel(realization=0)).plot(\n",
" ax=ax, color='C1', linestyle='dashdot', label='PW North, 10th percentile')\n",
"\n",
"ax.set_title('PacWave wave spectra, LWF scale', fontweight='bold')\n",
@@ -1066,8 +1076,8 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj,\n",
")\n",
- "\n",
- "print(f'Optimal average electrical power: {results.fun} W')"
+ "power_results = [result.fun for result in results]\n",
+ "print(f'Optimal average electrical power: {np.mean(power_results)} W')"
]
},
{
@@ -1076,7 +1086,7 @@
"metadata": {},
"outputs": [],
"source": [
- "x_wec, x_opt = wec.decompose_state(results.x)\n",
+ "x_wec, x_opt = wec.decompose_state(results[0].x)\n",
"print(len(x_opt))"
]
},
@@ -1088,8 +1098,8 @@
"source": [
"# Post-process\n",
"nsubsteps = 5\n",
- "wec_fdom, wec_tdom = wec.post_process(results, waves['south_max_occurrence'], nsubsteps)\n",
- "pto_fdom, pto_tdom = pto.post_process(wec, results, waves['south_max_occurrence'], nsubsteps)"
+ "wec_fdom, wec_tdom = wec.post_process(results[0], waves['south_max_occurrence'].sel(realization=0), nsubsteps)\n",
+ "pto_fdom, pto_tdom = pto.post_process(wec, results[0], waves['south_max_occurrence'].sel(realization=0), nsubsteps)"
]
},
{
@@ -1348,21 +1358,21 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj,\n",
" )\n",
- " print(f'Optimal average electrical power: {results.fun:.2f} W')\n",
- " x_wec, x_opt = wec.decompose_state(results.x)\n",
- " x_wec_jac , x_opt_jac = wec.decompose_state(results.jac)\n",
+ " print(f'Optimal average electrical power: {results[0].fun:.2f} W')\n",
+ " x_wec, x_opt = wec.decompose_state(results[0].x)\n",
+ " x_wec_jac , x_opt_jac = wec.decompose_state(results[0].jac)\n",
" ds = xr.Dataset(data_vars=dict(\n",
" x_wec=('wec_state', x_wec),\n",
" x_opt=('opt_state', x_opt),\n",
" x_wec_jac=('wec_state', x_wec_jac),\n",
" x_opt_jac=('opt_state', x_opt_jac),\n",
- " fval=results.fun),\n",
+ " fval=results[0].fun),\n",
" coords=dict(\n",
" wec_state=range(wec.nstate_wec),\n",
" opt_state=range(nstate_opt))\n",
" )\n",
" wot.write_netcdf(os.path.join(dir, 'data', f'tutorial_3_results_{x}.nc'), ds)\n",
- " return results.fun"
+ " return results[0].fun"
]
},
{
diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
index 0fd03425..a817b7a7 100644
--- a/examples/tutorial_4_Pioneer.ipynb
+++ b/examples/tutorial_4_Pioneer.ipynb
@@ -49,7 +49,11 @@
"\n",
"### 1.1 Waves\n",
"We start with the setting up the different waves we want to model, as this will inform what to select for our frequency range, which we need throughout the rest of the problem setup. \n",
- "We will consider two waves: a regular wave and an irregular wave, both with typical characteristics of the deployment site. The regular wave is roughly at 0.35 Hz, the known pitch resonance frequency of the buoy. The irregular wave has a peak period of 5 seconds, matching that of the deployment site."
+ "We will consider two waves: a regular wave and an irregular wave, both with typical characteristics of the deployment site. The regular wave is roughly at 0.35 Hz, the known pitch resonance frequency of the buoy. The irregular wave has a peak period of 5 seconds, matching that of the deployment site. \n",
+ "\n",
+ "The irregular wave is created with multiple phase realizations. For this tutorial, the number of realizations has been overwritten to 3 to reduce the total run-time. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations and average the resultant optimal power. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes.\n",
+ "\n",
+ "$t_{total} = \\frac{nrealizations}{f1} $"
]
},
{
@@ -77,10 +81,16 @@
"Hs = 1.5\n",
"Tp = 5 \n",
"\n",
+ "# number of realizations to reach 20 minutes of total simulation time\n",
+ "minutes_needed = 20\n",
+ "nrealizations = minutes_needed*60*f1\n",
+ "print(f'Number of realizations for a 20 minute total simulation time: {nrealizations}')\n",
+ "nrealizations = 3 # overwrite nrealizations to reduce run-time\n",
+ "\n",
"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)"
+ "waves_irregular = wot.waves.long_crested_wave(efth, nrealizations)"
]
},
{
@@ -100,7 +110,7 @@
"#TODO: highlight the harmonics if wave freq and Tp with other markers+colors\n",
"fig, ax = plt.subplots()\n",
"np.abs(waves_regular).plot(marker = 'x', label=\"regular\")\n",
- "np.abs(waves_irregular).plot(marker = '*', label=\"irregular\")\n",
+ "np.abs(waves_irregular.sel(realization=0)).plot(marker = '*', label=\"irregular\")\n",
"ax.set_title('Wave elevation spectrum', fontweight='bold')\n",
"plt.legend()"
]
@@ -629,7 +639,7 @@
" scale_x_opt=1e-2,\n",
" scale_obj=1e-2,\n",
")\n",
- "print(f'Optimal average power: {results.fun:.2f} W')"
+ "print(f'Optimal average power: {results[0].fun:.2f} W')"
]
},
{
@@ -647,17 +657,17 @@
"outputs": [],
"source": [
"nsubsteps = 5\n",
- "wec_fdom, wec_tdom = wec.post_process(results, waves_regular, nsubsteps=nsubsteps)\n",
+ "wec_fdom, wec_tdom = wec.post_process(results[0], waves_regular.sel(realization=0), nsubsteps=nsubsteps)\n",
"\n",
"# Manually post-process PTO and flywheel outputs\n",
- "x_wec, x_opt = wot.decompose_state(results.x, 1, nfreq)\n",
+ "x_wec, x_opt = wot.decompose_state(results[0].x, 1, nfreq)\n",
"fw_pos = np.dot(wec.time_mat_nsubsteps(nsubsteps), x_opt[nstate_pto:])\n",
- "pto_pos = rel_position(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
- "pto_vel = rel_velocity(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
- "pto_force = f_motor(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
+ "pto_pos = rel_position(wec, x_wec, x_opt, waves_regular.sel(realization=0), nsubsteps)\n",
+ "pto_vel = rel_velocity(wec, x_wec, x_opt, waves_regular.sel(realization=0), nsubsteps)\n",
+ "pto_force = f_motor(wec, x_wec, x_opt, waves_regular.sel(realization=0), nsubsteps)\n",
"pto_force_fd = wec.td_to_fd(pto_force[::nsubsteps])\n",
- "pto_mech_power = mechanical_power(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
- "pto_elec_power = electrical_power(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
+ "pto_mech_power = mechanical_power(wec, x_wec, x_opt, waves_regular.sel(realization=0), nsubsteps)\n",
+ "pto_elec_power = electrical_power(wec, x_wec, x_opt, waves_regular.sel(realization=0), nsubsteps)\n",
"avg_mech_power = np.mean(pto_mech_power)\n",
"avg_elec_power = np.mean(pto_elec_power)\n"
]
@@ -862,7 +872,8 @@
" scale_x_opt=1e-2,\n",
" scale_obj=1e-2,\n",
")\n",
- "print(f'Optimal average power: {results.fun:.2f} W')"
+ "power_results = [result.fun for result in results]\n",
+ "print(f'Optimal average power: {np.mean(power_results):.2f} W')"
]
},
{
@@ -879,17 +890,17 @@
"outputs": [],
"source": [
"nsubsteps = 5\n",
- "wec_fdom, wec_tdom = wec.post_process(results, waves_irregular, nsubsteps=nsubsteps)\n",
+ "wec_fdom, wec_tdom = wec.post_process(results[0], waves_irregular.sel(realization=0), nsubsteps=nsubsteps)\n",
"\n",
"# Manually post-process PTO and flywheel outputs\n",
- "x_wec, x_opt = wot.decompose_state(results.x, 1, nfreq)\n",
+ "x_wec, x_opt = wot.decompose_state(results[0].x, 1, nfreq)\n",
"fw_pos = np.dot(wec.time_mat_nsubsteps(nsubsteps), x_opt[nstate_pto:])\n",
- "pto_pos = rel_position(wec, x_wec, x_opt, waves_irregular, nsubsteps)\n",
- "pto_vel = rel_velocity(wec, x_wec, x_opt, waves_irregular, nsubsteps)\n",
- "pto_force = f_motor(wec, x_wec, x_opt, waves_irregular, nsubsteps)\n",
+ "pto_pos = rel_position(wec, x_wec, x_opt, waves_irregular.sel(realization=0), nsubsteps)\n",
+ "pto_vel = rel_velocity(wec, x_wec, x_opt, waves_irregular.sel(realization=0), nsubsteps)\n",
+ "pto_force = f_motor(wec, x_wec, x_opt, waves_irregular.sel(realization=0), nsubsteps)\n",
"pto_force_fd = wec.td_to_fd(pto_force[::nsubsteps])\n",
- "pto_mech_power = mechanical_power(wec, x_wec, x_opt, waves_irregular, nsubsteps)\n",
- "pto_elec_power = electrical_power(wec, x_wec, x_opt, waves_irregular, nsubsteps)\n",
+ "pto_mech_power = mechanical_power(wec, x_wec, x_opt, waves_irregular.sel(realization=0), nsubsteps)\n",
+ "pto_elec_power = electrical_power(wec, x_wec, x_opt, waves_irregular.sel(realization=0), nsubsteps)\n",
"avg_mech_power = np.mean(pto_mech_power)\n",
"avg_elec_power = np.mean(pto_elec_power)"
]
diff --git a/tests/test_core.py b/tests/test_core.py
index 7773525a..b8beb4fe 100644
--- a/tests/test_core.py
+++ b/tests/test_core.py
@@ -97,7 +97,7 @@ def wave_regular(f1, nfreq):
phase = 24.5 # degrees
wave = wot.waves.regular_wave(f1, nfreq, freq, amp, phase)
params = {'n': n, 'amp': amp, 'phase': phase}
- return wave, params
+ return wave.sel(realization=0), params
@pytest.fixture(scope="module")
@@ -160,21 +160,22 @@ def waves_multi(f1, nfreq):
"""
n = np.random.randint(1, nfreq)
directions = [0.0, 30.0]
+ nrealizations = 1
ndir = len(directions)
- amplitudes = np.zeros([nfreq, ndir])
- phases = np.zeros([nfreq, ndir])
+ amplitudes = np.zeros([nfreq, ndir, nrealizations])
+ phases = np.zeros([nfreq, ndir, nrealizations])
amp0, amp1 = 1.2, 2.1
phase0, phase1 = 26, -13
- amplitudes[n-1, :] = [amp0, amp1]
- phases[n-1, :] = [phase0, phase1]
+ amplitudes[n-1, :, :] = [[amp0], [amp1]]
+ phases[n-1, :, :] = [[phase0], [phase1]]
- waves = wot.waves.elevation_fd(f1, nfreq, directions, amplitudes, phases)
+ waves = wot.waves.elevation_fd(f1, nfreq, directions, nrealizations, amplitudes, phases)
- params = {'n': n, 'directions': directions, 'amp0': amp0, 'amp1': amp1,
- 'phase0': phase0, 'phase1': phase1}
+ params = {'n': n, 'directions': directions, 'nrealizations': nrealizations, 'amp0': amp0,
+ 'amp1': amp1, 'phase0': phase0, 'phase1': phase1}
- return waves, params
+ return waves.sel(realization=0), params
@pytest.fixture(scope="module")
@@ -1218,8 +1219,8 @@ def forces(self, hydro_data):
radiation = forces['radiation'](wec, x_wec, None, None)
hydrostatics = forces['hydrostatics'](wec, x_wec, None, None)
friction = forces['friction'](wec, x_wec, None, None)
- fk = forces['Froude_Krylov'](wec, None, None, waves)
- diffraction = forces['diffraction'](wec, None, None, waves)
+ fk = forces['Froude_Krylov'](wec, None, None, waves.sel(realization=0))
+ diffraction = forces['diffraction'](wec, None, None, waves.sel(realization=0))
# true/expected inertia and forces
inertia_truth = mass * acc
radiation_truth = -(rad*vel + addmass*acc)
@@ -1342,7 +1343,7 @@ def test_error_directions_not_subset(self, f1, nfreq, exc_coeff):
a subset of the hydrodata directions.
"""
with pytest.raises(ValueError):
- waves = wot.waves.elevation_fd(f1, nfreq, [0.0, 25])
+ waves = wot.waves.elevation_fd(f1, nfreq, [0.0, 25], 1)
wot.wave_excitation(exc_coeff, waves)
def test_error_different_frequencies(
@@ -1354,7 +1355,7 @@ def test_error_different_frequencies(
_, params_multi = waves_multi
directions = np.deg2rad(params_multi['directions'])
with pytest.raises(ValueError):
- waves = wot.waves.elevation_fd(f1+0.01, nfreq, directions)
+ waves = wot.waves.elevation_fd(f1+0.01, nfreq, directions, 1)
wot.wave_excitation(exc_coeff, waves)
diff --git a/tests/test_integration.py b/tests/test_integration.py
index 4cef7702..ea116fa4 100644
--- a/tests/test_integration.py
+++ b/tests/test_integration.py
@@ -192,7 +192,7 @@ def test_solve_bounds(bounds_opt, wec_from_bem, regular_wave,
bounds_opt=bounds_opt,
)
- assert pytest.approx(kplim, 1e-5) == res['x'][-1]
+ assert pytest.approx(kplim, 1e-5) == res[0]['x'][-1]
def test_same_wec_init(wec_from_bem,
@@ -208,9 +208,9 @@ def test_same_wec_init(wec_from_bem,
x_wec_0 = np.random.randn(wec_from_bem.nstate_wec)
np.random.seed(1)
x_opt_0 = np.random.randn(wec_from_bem.nstate_wec)
- bem_res = wec_from_bem.residual(x_wec_0, x_opt_0, waves)
- fb_res = wec_from_floatingbody.residual(x_wec_0, x_opt_0, waves)
- imp_res = wec_from_impedance.residual(x_wec_0, x_opt_0, waves)
+ 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)
@@ -261,9 +261,9 @@ def test_p_controller_resonant_wave(self,
bounds_opt=((-1*np.infty, 0),),
)
- power_sol = -1*res['fun']
+ power_sol = -1*res[0]['fun']
- res_fd, _ = wec.post_process(res, resonant_wave,nsubsteps=1)
+ res_fd, _ = wec.post_process(res[0], resonant_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())
@@ -293,9 +293,9 @@ def test_pi_controller_regular_wave(self,
bounds_opt=((-1e4, 0), (0, 2e4),)
)
- power_sol = -1*res['fun']
+ power_sol = -1*res[0]['fun']
- res_fd, _ = wec.post_process(res, regular_wave, nsubsteps=1)
+ res_fd, _ = wec.post_process(res[0], regular_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())
@@ -324,9 +324,9 @@ def test_unstructured_controller_irregular_wave(self,
scale_obj=1e-2,
)
- power_sol = -1*res['fun']
+ power_sol = -1*res[0]['fun']
- res_fd, _ = wec.post_process(res, regular_wave, nsubsteps=1)
+ res_fd, _ = wec.post_process(res[0], regular_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())
@@ -405,8 +405,8 @@ def saturated_pi(pto, wec, x_wec, x_opt, waves=None, nsubsteps=1):
nsubstep_postprocess = 4
pto_fdom[key], pto_tdom[key] = pto[key].post_process(wec[key],
- res[key],
- regular_wave,
+ res[key][0],
+ regular_wave.sel(realization=0),
nsubstep_postprocess)
xr.testing.assert_allclose(pto_tdom['pi'].power.squeeze().mean('time'),
diff --git a/tests/test_waves.py b/tests/test_waves.py
index 67fd7ce7..42d17303 100644
--- a/tests/test_waves.py
+++ b/tests/test_waves.py
@@ -58,13 +58,13 @@ class TestElevationFD:
def elevation(self, f1, nfreq, directions):
"""Complex sea state elevation amplitude [m] indexed by
frequency and direction."""
- return wot.waves.elevation_fd(f1, nfreq, directions)
+ return wot.waves.elevation_fd(f1, nfreq, directions, 1)
def test_coordinates(self, elevation):
"""Test that the elevation dataArray has the correct
coordinates.
"""
- coordinates = ['wave_direction', 'omega', 'freq']
+ coordinates = ['wave_direction', 'omega', 'freq', 'realization']
for icoord in coordinates:
assert icoord in elevation.coords, f'missing coordinate: {icoord}'
@@ -115,7 +115,7 @@ def test_coordinates(self, elevation):
"""Test that the elevation dataArray has the correct
coordinates.
"""
- coordinates = ['wave_direction', 'omega', 'freq']
+ coordinates = ['wave_direction', 'omega', 'freq', 'realization']
for icoord in coordinates:
assert icoord in elevation.coords, f'missing coordinate: {icoord}'
@@ -181,11 +181,16 @@ def direction(self, ndbc_omnidirectional, ndir):
return ndbc_omnidirectional.dir.values[np.random.randint(0, ndir)]
@pytest.fixture(scope="class")
- def elevation(self, ndbc_omnidirectional, direction):
+ def nrealizations(self):
+ """Number of wave realizations."""
+ return 2
+
+ @pytest.fixture(scope="class")
+ 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)
+ ndbc_omnidirectional.efth, nrealizations, direction)
return elev
@pytest.fixture(scope="class")
@@ -222,13 +227,13 @@ def test_coordinates(self, elevation):
"""Test that the elevation dataArray has the correct
coordinates.
"""
- coordinates = ['wave_direction', 'omega', 'freq']
+ coordinates = ['wave_direction', 'omega', 'freq', 'realization']
for icoord in coordinates:
assert icoord in elevation.coords, f'missing coordinate: {icoord}'
- def test_shape(self, elevation, nfreq, ndir):
+ def test_shape(self, elevation, nfreq, ndir, nrealizations):
"""Test that the elevation dataArray has the correct shape."""
- assert np.squeeze(elevation.values).shape == (nfreq, )
+ assert np.squeeze(elevation.values).shape == (nfreq, nrealizations)
def test_type(self, elevation):
"""Test that the elevation dataArray has the correct type."""
@@ -239,6 +244,11 @@ 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):
+ """Test that the number of realizations is correct."""
+ realization_out = elevation.realization.values
+ assert (realization_out == [0,1]).all()
+
def test_spectrum(self, pm_spectrum, pm_hs):
"""Test that the constructed spectrum has the expected Hs."""
efth = ws.SpecArray(pm_spectrum)
@@ -248,8 +258,9 @@ def test_time_series(self, pm_spectrum, pm_f1, pm_nfreq):
"""Test that the created time series has the desired spectrum."""
# create time-series
direction = 0.0
- wave = wot.waves.long_crested_wave(pm_spectrum, direction)
- wave_ts = wot.fd_to_td(wave.values, pm_f1, pm_nfreq, False)
+ nrealizations = 1
+ 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)
fs = 1/t[1]
@@ -278,30 +289,40 @@ 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
coordinates.
"""
- coordinates = ['wave_direction', 'omega', 'freq']
+ coordinates = ['wave_direction', 'omega', 'freq', 'realization']
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:
@@ -312,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 a569c0da..890bcb5f 100644
--- a/wecopttool/core.py
+++ b/wecopttool/core.py
@@ -608,7 +608,7 @@ def solve(self,
bounds_wec: Optional[Bounds] = None,
bounds_opt: Optional[Bounds] = None,
callback: Optional[TStateFunction] = None,
- ) -> OptimizeResult:
+ ) -> list[OptimizeResult]:
"""Simulate WEC dynamics using a pseudo-spectral solution
method and returns the raw results dictionary produced by
:py:func:`scipy.optimize.minimize`.
@@ -683,18 +683,20 @@ 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)
- >>> res_pto_fd, res_pto_td = pto.post_process(wec,res_opt)
+ >>> 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
--------
wecopttool.waves,
"""
- _log.info("Solving pseudo-spectral control problem.")
+ results = []
# x_wec scaling vector
if scale_x_wec == None:
@@ -712,48 +714,13 @@ def solve(self,
# composite scaling vector
scale = np.concatenate([scale_x_wec, scale_x_opt])
-
+
# 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
-
- # 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, waves)*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, waves)
- 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, waves)
-
- eq_cons = {'type': 'eq', 'fun': scaled_resid_fun}
- if use_grad:
- eq_cons['jac'] = jacobian(scaled_resid_fun)
- constraints.append(eq_cons)
+ x0 = np.concatenate([x_wec_0, x_opt_0])*scale
# bounds
if (bounds_wec is None) and (bounds_opt is None):
@@ -778,51 +745,93 @@ 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):
+ 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
+
+ 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 callback(self, x_wec, x_opt, waves)
-
- # optimization problem
- optim_options['disp'] = optim_options.get('disp', True)
- problem = {'fun': obj_fun_scaled,
- 'x0': x0,
- 'method': 'SLSQP',
- 'constraints': constraints,
- 'options': optim_options,
- 'bounds': bounds,
- 'callback': callback_scipy,
- }
- if use_grad:
- problem['jac'] = grad(obj_fun_scaled)
-
- # minimize
- optim_res = minimize(**problem)
-
- msg = f'{optim_res.message} (Exit mode {optim_res.status})'
- if optim_res.status == 0:
- _log.info(msg)
- elif optim_res.status == 9:
- _log.warning(msg)
- else:
- raise Exception(msg)
+ return self.residual(x_wec, x_opt, wave)
- # unscale
- optim_res.x = optim_res.x / scale
- optim_res.fun = optim_res.fun / scale_obj
- optim_res.jac = optim_res.jac / scale_obj * scale
+ 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)
+ problem = {'fun': obj_fun_scaled,
+ 'x0': x0,
+ 'method': 'SLSQP',
+ 'constraints': constraints,
+ 'options': optim_options,
+ 'bounds': bounds,
+ 'callback': callback_scipy,
+ }
+ if use_grad:
+ problem['jac'] = grad(obj_fun_scaled)
- return optim_res
+ # minimize
+ optim_res = minimize(**problem)
+
+ msg = f'{optim_res.message} (Exit mode {optim_res.status})'
+ if optim_res.status == 0:
+ _log.info(msg)
+ elif optim_res.status == 9:
+ _log.warning(msg)
+ else:
+ raise Exception(msg)
+
+ # unscale
+ optim_res.x = optim_res.x / scale
+ optim_res.fun = optim_res.fun / scale_obj
+ optim_res.jac = optim_res.jac / scale_obj * scale
+
+ results.append(optim_res)
+
+ return results
def post_process(self,
res: OptimizeResult,
diff --git a/wecopttool/waves.py b/wecopttool/waves.py
index 7603f9b7..97ddc1a2 100644
--- a/wecopttool/waves.py
+++ b/wecopttool/waves.py
@@ -54,6 +54,7 @@ def elevation_fd(
f1: float,
nfreq: int,
directions: Union[float, ArrayLike],
+ nrealizations: int,
amplitudes: Optional[ArrayLike] = None,
phases: Optional[ArrayLike] = None,
attr: Optional[Mapping] = None,
@@ -75,6 +76,8 @@ def elevation_fd(
i.e., :python:`freq = [0, f1, 2*f1, ..., nfreq*f1]`.
directions
Wave directions in degrees. 1D array.
+ nrealizations
+ Number of wave phase realizations.
amplitudes:
Wave elevation amplitude in meters.
phases:
@@ -88,24 +91,33 @@ def elevation_fd(
"""
directions = np.atleast_1d(degrees_to_radians(directions, sort=False))
ndirections = len(directions)
+ realization = range(nrealizations)
freq = frequency(f1, nfreq, False)
omega = freq*2*np.pi
- dims = ('omega', 'wave_direction')
+ 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'}
+ real_attr = {'long_name': 'Phase realization', 'units': ''}
coords = {'omega': (dims[0], omega, omega_attr),
'freq': (dims[0], freq, freq_attr),
- 'wave_direction': (dims[1], directions, dir_attr)}
+ 'wave_direction': (dims[1], directions, dir_attr),
+ 'realization': (dims[2], realization, real_attr)}
if amplitudes is None:
- amplitudes = np.zeros([nfreq, ndirections])
+ 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],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)
@@ -149,7 +161,7 @@ def regular_wave(
# attributes & index
omega = freq*2*np.pi
- tmp_waves = elevation_fd(f1, nfreq, direction)
+ tmp_waves = elevation_fd(f1, nfreq, direction, 1)
iomega = tmp_waves.sel(omega=omega, method='nearest').omega.values
ifreq = iomega/(2*np.pi)
@@ -175,8 +187,8 @@ def regular_wave(
rphase = degrees_to_radians(phase)
# wave elevation
- tmp = np.zeros([nfreq, 1])
- waves = elevation_fd(f1, nfreq, direction, tmp, tmp, attrs)
+ tmp = np.zeros([nfreq, 1, 1])
+ waves = elevation_fd(f1, nfreq, direction, 1, tmp, tmp, attrs)
waves.loc[{'omega': iomega}] = amplitude * np.exp(1j*rphase)
return waves
@@ -184,6 +196,7 @@ def regular_wave(
def long_crested_wave(
efth: DataArray,
+ nrealizations: int,
direction: Optional[float] = 0.0,
seed: Optional[float] = None,
) -> DataArray:
@@ -203,6 +216,9 @@ def long_crested_wave(
efth
Omnidirection wave spectrum in units of m^2/Hz, in the format
used by :py:class:`wavespectra.SpecArray`.
+ nrealizations
+ Number of wave phase realizations to be created for the
+ long-crested wave.
direction
Direction (in degrees) of the long-crested wave.
seed
@@ -221,10 +237,11 @@ def long_crested_wave(
'Direction (degrees)': direction,
}
- return elevation_fd(f1, nfreq, direction, amplitudes, None, attr, seed)
+ return elevation_fd(f1, nfreq, direction, nrealizations, amplitudes, None, attr, seed)
def irregular_wave(efth: DataArray,
+ nrealizations: int,
seed: Optional[float] = None,) -> DataArray:
"""Create a complex frequency-domain wave elevation from a spectrum.
@@ -243,6 +260,9 @@ def irregular_wave(efth: DataArray,
efth
Wave spectrum in units of m^2/Hz/deg, in the format used by
:py:class:`wavespectra.SpecArray`.
+ nrealizations
+ Number of wave phase realizations to be created for the
+ irregular wave.
seed
Seed for random number generator. Used for reproducibility.
Generally should not be used except for testing.
@@ -258,11 +278,11 @@ def irregular_wave(efth: DataArray,
attr = {'Wave type': 'Irregular'}
- return elevation_fd(f1, nfreq, directions, amplitudes, None, attr, seed)
+ return elevation_fd(f1, nfreq, directions, nrealizations, amplitudes, None, attr, seed)
def random_phase(
- shape: Optional[Union[Iterable[int], int]] = None,
+ shape: Optional[Union[Iterable[int], int, int]] = None,
seed: Optional[float] = None,
) -> Union[float , ndarray]:
"""Generate random phases in range [-π, π) radians.