Utils dev - replaces old utilities PR (sandialabs#343)

* run CI on PRs against dev branch * coppied fundamental utility files * import utilities module * added utilities funtions to tut1 * added bem plot from utils * added bem plot from utils * updated sankey plot * updated check_radiation_damping * cleared outputs * corrected bug * changed Zi to hydro_impedance to be consistent with our variables name python convention * PR review edits * add grid to plots * removed draft functions in utilities.py * typo * Fixed one more typo I found while reviewing Daniel's changes --------- Co-authored-by: Ryan Coe <rcoe@sandia.gov> Co-authored-by: Michael Devin <michaelcdevin@outlook.com>
dtgaebe · May 2, 2024 · b8506d8 · b8506d8
1 parent ea039c5
commit b8506d8
Show file tree

Hide file tree

Showing 8 changed files with 689 additions and 120 deletions.
diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb
@@ -267,6 +267,33 @@
     "# wot.write_netcdf('bem_data.nc', bem_data) # saves BEM data to file"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next, we use the utilities class `wecopttool.utilities`, which has functions that help you analyze and design WECs, but are not needed for the core function of WecOptTool.\n",
+    "For example, we calculate the WEC's intrinsic impedance with its hydrodynamic coefficients and inherent inertial properties. We make use of `wecopttool.add_linear_friction()` to convert the `bem_data` into a dataset which contains all zero friction data, because `wecopttool.check_radiation_damping()` and `wecopttool.hydrodynamic_impedance()` excpect a data variable called 'friction'.\n",
+    "\n",
+    "The intrinsic impedance tells us how a hydrodynamic body responds to different excitation frequencies. For oscillating systems we are intersted in both, the amplitude of the resulting velocity as well as the phase between velocity and excitation force. Bode plots are useful tool to visualize the frequency response function.\n",
+    "\n",
+    "The natural frequency reveals itself as a trough in the intrinsic impedance for restoring degrees of freedom (heave, roll, pitch)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hd = wot.add_linear_friction(bem_data, friction = None) \n",
+    "#we're not actually adding friction, but need the datavariables in hd \n",
+    "hd = wot.check_radiation_damping(hd)\n",
+    "\n",
+    "intrinsic_impedance = wot.hydrodynamic_impedance(hd)\n",
+    "fig, axes = wot.utilities.plot_bode_impedance(intrinsic_impedance,\n",
+    "                                              'WaveBot Intrinsic Impedance')"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -546,6 +573,46 @@
     "pto_tdom"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Lastly, we will visualize the average power at different stages and how the power flows through the system.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p_flows = wot.utilities.calculate_power_flows(wec,\n",
+    "                                              pto, \n",
+    "                                              results[0], \n",
+    "                                              waves.sel(realization=0), \n",
+    "                                              intrinsic_impedance)\n",
+    "wot.utilities.plot_power_flow(p_flows)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "On the very left of the power flow Sankey diagram we start with the *Optimal Excitation*, which is only determined by the incident wave and the hydrodynamic damping and friction. \n",
+    "In order for the actual *Excited* power to equal the *Optimal Excitation* the absorbed *Mechanical* power would need to equal the *Radiated* power (power that is radiated away from the WEC). In this case the *Unused Potential* would be zero.\n",
+    "In other words, you can never convert more than 50% of the incident wave energy into kinetic energy. For more information on this we refer to Johannes Falnes Book - Ocean Waves and Oscillating System, specifically Chapter 6.\n",
+    "\n",
+    "However, the optimal 50% absorption is an overly optimistic goal considering that real world systems are likely constrained in their oscillation velocity amplitude and phase, due to limitations in stroke and applicable force.\n",
+    "\n",
+    "\n",
+    "The PTO force constraint used in this optimization also stopped us from absorbing the maximal potential wave energy. It is notable that we absorbed approximately 3/4 of the maximal absorbable power (*Mechanical* / (1/2 *Optimal Excitation*)), with relatively little *Radiated* power (about 1/2 of the absorbed power). To also absorb the last 1/4 of the potential wave power, we would need to increase the *Radiated* power three times!!\n",
+    "\n",
+    "\n",
+    "A more important question than \"How to achieve optimal absorption?\" is \"How do we optimize usable output power?\", i.e. *Electrical* power. For this optimization case we used a lossless PTO that has no impedance in itself. Therefore, the *Electrical* power equals the absorbed *Mechanical* power.\n",
+    "We'll show in the following sections that this a poor assumption and that the power flow looks fundamentally different when taking the PTO dynamics into account!!\n"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -752,7 +819,19 @@
     "\n",
     "The PTO force trajectory for optimizing mechanical power is saturated at the maximum for longer compared to the electrical power.\n",
     "This could inform the WEC designer optimizing for mechanical power to consider larger components that would not be utilized at their limit as frequently.\n",
-    "However, the electrical power (_not_ the mechanical power) is the usable form of power, thus designing the WEC for optimal electrical power does not indicate a need for larger components and prevents this over-design."
+    "However, the electrical power (_not_ the mechanical power) is the usable form of power, thus designing the WEC for optimal electrical power does not indicate a need for larger components and prevents this over-design.\n",
+    "\n",
+    "The Sankey power flow diagram confirm this observation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p_flows_2 = wot.utilities.calculate_power_flows(wec_2, pto_2, results_2[0], waves.sel(realization = 0), intrinsic_impedance)\n",
+    "wot.utilities.plot_power_flow(p_flows_2)"
    ]
   },
   {
@@ -1003,7 +1082,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   },
   "vscode": {
    "interpreter": {

diff --git a/examples/tutorial_2_AquaHarmonics.ipynb b/examples/tutorial_2_AquaHarmonics.ipynb
@@ -170,23 +170,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, axes = plt.subplots(3,1)\n",
-    "bem_data['added_mass'].plot(ax = axes[0])\n",
-    "bem_data['radiation_damping'].plot(ax = axes[1])\n",
-    "axes[2].plot(bem_data['omega'],np.abs(np.squeeze(bem_data['diffraction_force'].values)), color = 'orange')\n",
-    "axes[2].set_ylabel('abs(diffraction_force)', color = 'orange')\n",
-    "axes[2].tick_params(axis ='y', labelcolor = 'orange')\n",
-    "ax2r = axes[2].twinx()\n",
-    "ax2r.plot(bem_data['omega'], np.abs(np.squeeze(bem_data['Froude_Krylov_force'].values)), color = 'blue')\n",
-    "ax2r.set_ylabel('abs(FK_force)', color = 'blue')\n",
-    "ax2r.tick_params(axis ='y', labelcolor = 'blue')\n",
-    "\n",
-    "for axi in axes:\n",
-    "    axi.set_title('')\n",
-    "    axi.label_outer()\n",
-    "    axi.grid()\n",
-    "    \n",
-    "axes[-1].set_xlabel('Frequency [rad/s]')"
+    "[(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)] = wot.utilities.plot_hydrodynamic_coefficients(bem_data)"
    ]
   },
   {
@@ -916,7 +900,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   },
   "vscode": {
    "interpreter": {

diff --git a/examples/tutorial_3_LUPA.ipynb b/examples/tutorial_3_LUPA.ipynb
@@ -472,65 +472,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot coefficients\n",
-    "radiating_dofs = bem_data.radiating_dof.values\n",
-    "influenced_dofs = bem_data.influenced_dof.values\n",
-    "\n",
-    "# plots\n",
-    "fig_am, ax_am = plt.subplots(len(radiating_dofs), len(influenced_dofs),\n",
-    "                             tight_layout=True, sharex=True)\n",
-    "fig_rd, ax_rd = plt.subplots(len(radiating_dofs), len(influenced_dofs),\n",
-    "                             tight_layout=True, sharex=True)\n",
-    "fig_ex, ax_ex = plt.subplots(len(influenced_dofs), 1,\n",
-    "                             tight_layout=True, sharex=True, figsize=(4, 6))\n",
-    "\n",
-    "# plot titles\n",
-    "fig_am.suptitle('Added Mass Coefficients', fontweight='bold')\n",
-    "fig_rd.suptitle('Radiation Damping Coefficients', fontweight='bold')\n",
-    "fig_ex.suptitle('Wave Excitation Coefficients', fontweight='bold')\n",
-    "\n",
-    "# subplotting across 4DOF\n",
-    "sp_idx = 0\n",
-    "for i, rdof in enumerate(radiating_dofs):\n",
-    "    for j, idof in enumerate(influenced_dofs):\n",
-    "        sp_idx += 1\n",
-    "        if i == 0:\n",
-    "            np.abs(bem_data.diffraction_force.sel(influenced_dof=idof)).plot(\n",
-    "                ax=ax_ex[j], linestyle='dashed', label='Diffraction force')\n",
-    "            np.abs(bem_data.Froude_Krylov_force.sel(influenced_dof=idof)).plot(\n",
-    "                ax=ax_ex[j], linestyle='dashdot', label='Froude-Krylov force')\n",
-    "            ex_handles, ex_labels = ax_ex[j].get_legend_handles_labels()\n",
-    "            ax_ex[j].set_title(f'{idof}')\n",
-    "            ax_ex[j].set_xlabel('')\n",
-    "            ax_ex[j].set_ylabel('')\n",
-    "        if j <= i:\n",
-    "            bem_data.added_mass.sel(\n",
-    "                radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_am[i, j])\n",
-    "            bem_data.radiation_damping.sel(\n",
-    "                radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_rd[i, j])\n",
-    "            if i == len(radiating_dofs)-1:\n",
-    "                ax_am[i, j].set_xlabel(f'$\\omega$', fontsize=10)\n",
-    "                ax_rd[i, j].set_xlabel(f'$\\omega$', fontsize=10)\n",
-    "                ax_ex[j].set_xlabel(f'$\\omega$', fontsize=10)\n",
-    "            else:\n",
-    "                ax_am[i, j].set_xlabel('')\n",
-    "                ax_rd[i, j].set_xlabel('')\n",
-    "            if j == 0:\n",
-    "                ax_am[i, j].set_ylabel(f'{rdof}', fontsize=10)\n",
-    "                ax_rd[i, j].set_ylabel(f'{rdof}', fontsize=10)\n",
-    "            else:\n",
-    "                ax_am[i, j].set_ylabel('')\n",
-    "                ax_rd[i, j].set_ylabel('')\n",
-    "            if j == i:\n",
-    "                ax_am[i, j].set_title(f'{idof}', fontsize=10)\n",
-    "                ax_rd[i, j].set_title(f'{idof}', fontsize=10)\n",
-    "            else:\n",
-    "                ax_am[i, j].set_title('')\n",
-    "                ax_rd[i, j].set_title('')\n",
-    "        else:\n",
-    "            fig_am.delaxes(ax_am[i, j])\n",
-    "            fig_rd.delaxes(ax_rd[i, j])\n",
-    "fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False)"
+    "wot.utilities.plot_hydrodynamic_coefficients(bem_data)"
    ]
   },
   {
@@ -1241,6 +1183,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "influenced_dofs = bem_data.influenced_dof.values\n",
     "## Forces\n",
     "fig_res2, ax_res2 = plt.subplots(len(influenced_dofs), 1, sharex=True, figsize=(8, 10))\n",
     "for j, idof in enumerate(influenced_dofs):\n",
@@ -1499,7 +1442,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   },
   "vscode": {
    "interpreter": {

diff --git a/examples/tutorial_4_Pioneer.ipynb b/examples/tutorial_4_Pioneer.ipynb
@@ -838,7 +838,14 @@
     "P_opt_excitation = 2*P_max_absorbable # after Falnes Eq. 6.10\n",
     "P_radiated = ((1/2)*(Rad_res * np.abs(Vel)**2 ).squeeze().sum('omega').item()) # after Falnes Eq. 6.4\n",
     "P_excited= (1/4)*(Fex*np.conjugate(Vel) + np.conjugate(Fex)*Vel ).squeeze().sum('omega').item() # after Falnes Eq. 6.3\n",
-    "P_absorbed = P_excited - P_radiated # after Falnes Eq. 6.2 absorbed by WEC-PTO (difference between actual excitation power and actual radiated power needs to be absorbed by PTO)"
+    "\n",
+    "power_flows = {'Optimal Excitation' : -2* (np.real(P_max_absorbable)),    \n",
+    "            'Radiated': -1*(np.real(P_radiated)), \n",
+    "            'Actual Excitation': -1*(np.real(P_excited)), \n",
+    "}\n",
+    "power_flows['Unused Potential'] =  power_flows['Optimal Excitation'] - power_flows['Actual Excitation']\n",
+    "power_flows['Absorbed'] =  power_flows['Actual Excitation'] - power_flows['Radiated']  # after Falnes Eq. 6.2 \n",
+    "#absorbed by WEC-PTO (difference between actual excitation power and actual radiated power needs to be absorbed by PTO)"
    ]
   },
   {
@@ -868,9 +875,18 @@
     "    return vel_td * force_td\n",
     "\n",
     "fw_fric_power = fw_friction_power(wec, x_wec, x_opt, waves_regular, nsubsteps)\n",
-    "avg_fw_fric_power = np.mean(fw_fric_power)\n",
-    "\n",
-    "assert(np.isclose(P_absorbed, -1*(avg_mech_power -avg_fw_fric_power), rtol=0.01)) # assert that solver solution matches"
+    "avg_fw_fric_power = np.mean(fw_fric_power)\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "power_flows['Electrical (solver)']= avg_elec_power  #solver determines sign\n",
+    "power_flows['Mechanical (solver)']= avg_mech_power - avg_fw_fric_power #solver determines sign, friction is dissipated\n",
+    "power_flows['PTO Loss'] = power_flows['Mechanical (solver)'] -  power_flows['Electrical (solver)'] "
    ]
   },
   {
@@ -879,39 +895,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "from matplotlib.sankey import Sankey\n",
-    "\n",
-    "P_PTO_loss = avg_mech_power - avg_elec_power  \n",
-    "P_unused = P_opt_excitation - P_excited # Difference between the theoretical optimum excitation, if the WEC velocity would be in resonance with the excitation force\n",
-    "\n",
-    "Power_flows = [P_opt_excitation, P_PTO_loss, -1*avg_fw_fric_power, -1*P_radiated, -1*P_unused, avg_elec_power, ]\n",
-    "\n",
-    "fig = plt.figure(figsize = [8,4])\n",
-    "ax = fig.add_subplot(1, 1, 1,)\n",
-    "sankey = Sankey(ax=ax, \n",
-    "                scale=0.5/P_max_absorbable,\n",
-    "                offset= 0,\n",
-    "                format = '%.2f W',shoulder = 0.02)\n",
-    "\n",
-    "sankey.add(flows=Power_flows, \n",
-    "           labels = ['Optimal Excitation \\n $ 2 \\\\frac{\\left | F_e \\\\right |^2}{8Re(Z_i)} = 2*P_{mech}^{max}$ ', \n",
-    "                      'PTO-Loss \\n $ P_{mech} - P_{elec}$', \n",
-    "                      'Flywheel friction',\n",
-    "                      'Radiated \\n $ \\\\frac{1}{2} Re(Z_i) \\left | U \\\\right |^2  $ ', \n",
-    "                      'Unused Potential \\n(neither absorbed nor radiated)', \n",
-    "                      'Electrical'], \n",
-    "           orientations=[0, -1, -1,-1, -1, 0], # arrow directions\n",
-    "           pathlengths = [0.4,0.4,0.8,1.1,0.7,0.5,],\n",
-    "           trunklength = 1.5,\n",
-    "           edgecolor = '#2b8cbe',\n",
-    "           facecolor = '#2b8cbe' )\n",
-    "\n",
-    "diagrams = sankey.finish()\n",
-    "for diagram in diagrams:\n",
-    "    for text in diagram.texts:\n",
-    "        text.set_fontsize(10)\n",
-    "plt.axis(\"off\") \n",
-    "plt.show()"
+    "wot.utilities.plot_power_flow(power_flows)"
    ]
   },
   {
@@ -1178,7 +1162,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.6"
+   "version": "3.11.7"
   }
  },
  "nbformat": 4,