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
sandialabs · Dec 7, 2023 · 598e875 · 598e875
1 parent 9a36dca
commit 598e875
Show file tree

Hide file tree

Showing 9 changed files with 265 additions and 192 deletions.
diff --git 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
@@ -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,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"
    ]
   },
   {

diff --git 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 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} $"
    ]
   },
   {
@@ -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
@@ -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)"
    ]