Skip to content
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

Jump to bottom

Add multiple phase realizations to waves #297

Merged
merged 13 commits into from
Dec 7, 2023
Merged

Add multiple phase realizations to waves #297

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
29f8077
Add realization to waves variable and adjust tutorials to handle it
jtgrasb Nov 15, 2023
161f367
Update tutorial_3_LUPA.ipynb
jtgrasb Nov 15, 2023
5437203
Update wave and core functionality
jtgrasb Dec 6, 2023
114e019
Update tests for multiple phases
jtgrasb Dec 6, 2023
c8e8c60
Remove demonstration from pioneer
jtgrasb Dec 6, 2023
1647fd6
Update tests
jtgrasb Dec 6, 2023
aee9a7a
Merge branch 'main' into multiple_phases
jtgrasb Dec 6, 2023
45558a8
Test
jtgrasb Dec 6, 2023
551704b
Fix tutorials
jtgrasb Dec 7, 2023
141b61e
Reformat core.py
jtgrasb Dec 7, 2023
7e7edcf
Update AH tutorial
jtgrasb Dec 7, 2023
4f3b514
Update LUPA and Pioneer nrealizations and explanation
jtgrasb Dec 7, 2023
c2fcaa3
Update based on review
jtgrasb Dec 7, 2023
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
20 changes: 10 additions & 10 deletions examples/tutorial_1_WaveBot.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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')"
]
},
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -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)"
]
},
{
Expand Down Expand Up @@ -784,7 +784,7 @@
" scale_x_opt=scale_x_opt,\n",
" scale_obj=scale_obj)\n",
"\n",
" return res.fun"
" return res[0].fun"
]
},
{
Expand Down
16 changes: 8 additions & 8 deletions examples/tutorial_2_AquaHarmonics.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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')"
]
},
{
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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",
Expand Down Expand Up @@ -784,18 +784,18 @@
" 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",
" )\n",
"\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"
]
},
{
Expand Down
50 changes: 30 additions & 20 deletions examples/tutorial_3_LUPA.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 90<sup>th</sup> percentile, maximum percent annual energy, maximum occurrence, and minimum 10<sup>th</sup> 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} $"
]
},
{
Expand All @@ -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",
Expand All @@ -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'])"
Expand All @@ -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",
Expand Down Expand Up @@ -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')"
]
},
{
Expand All @@ -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))"
]
},
Expand All @@ -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)"
]
},
{
Expand Down Expand Up @@ -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"
]
},
{
Expand Down
49 changes: 30 additions & 19 deletions examples/tutorial_4_Pioneer.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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} $"
]
},
{
Expand Down Expand Up @@ -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)"
]
},
{
Expand All @@ -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()"
]
Expand Down Expand Up @@ -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')"
]
},
{
Expand All @@ -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"
]
Expand Down Expand Up @@ -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')"
]
},
{
Expand All @@ -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)"
]
Expand Down
Loading
Loading