Draft: Utilities module

examples/tutorial_1_WaveBot.ipynb

@@ -164,6 +164,31 @@
     "bem_data = wot.run_bem(fb, freq)"
    ]
   },
+  {
+   "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 it's hydrodynamic coefficients and inherent inertial properties.\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).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hd = wot.add_linear_friction(bem_data, friction = None) #we're not actually adding friction\n",
+    "hd = wot.check_linear_damping(hd)\n",
+    "\n",
+    "intrinsic_impedance = wot.hydrodynamic_impedance(hd)\n",
+    "fig, axes = wot.utilities.plot_bode_impedance(intrinsic_impedance,'WaveBot Intrinsic Impedance')"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -417,6 +442,36 @@
     "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",
+    "\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 useable output power?\", e.g. *Electrical* power. For this optimization case we used a loss-less 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 powerflow looks fundamentally different when taking the PTO dynamics into account!!\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p_flows = wot.utilities.calculate_power_flows(wec, pto, results[0], waves.sel(realization=0), intrinsic_impedance)\n",
+    "wot.utilities.plot_power_flow(p_flows)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -611,7 +666,36 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The attentive user might have noticed that the amplitude of the position, force and power signals is about half the magnitude of the signals we plotted in the first part of the tutorial. We can see that optimizing for electrical power requires optimal state trajectories with smaller amplitudes. For most WECs the electrical power is the usable form of power, thus the WEC should be designed for electrical power and we can avoid over-designing, which would results from expecting the forces associated with the optimal trajectories for mechanical power maximisation."
+    "The attentive user might have noticed that the amplitude of the position, force and power signals is about half the magnitude of the signals we plotted in the first part of the tutorial. We can see that optimizing for electrical power requires optimal state trajectories with smaller amplitudes. For most WECs the electrical power is the usable form of power, thus the WEC should be designed for electrical power and we can avoid over-designing, which would results from expecting the forces associated with the optimal trajectories for mechanical power maximisation.\n",
+    "\n",
+    "The Sankey power flow diagrams 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)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Another very interesting observation can be drawn from the power flow diagram for the unconstrained case, but optimized for electrical power (below). Although we absorbed 40% more kinetic energy we only converted 13% more *Electrical* power. For the case of the WaveBot increasing the range of applicablt force does not seem to be a good design decision.  "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p_flows_2_nocon = wot.utilities.calculate_power_flows(wec_2_nocon, pto_2, results_2_nocon[0], waves.sel(realization = 0), intrinsic_impedance)\n",
+    "wot.utilities.plot_power_flow(p_flows_2_nocon)"
    ]
   },
   {
@@ -865,7 +949,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.6"
   },
   "vscode": {
    "interpreter": {

examples/tutorial_2_AquaHarmonics.ipynb

@@ -143,23 +143,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]')"
+    "wot.utilities.plot_hydrodynamic_coefficients(bem_data)"
    ]
   },
   {
@@ -915,7 +899,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.6"
   },
   "vscode": {
    "interpreter": {

examples/tutorial_3_LUPA.ipynb

@@ -380,65 +380,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)"
    ]
   },
   {
@@ -1228,6 +1170,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",
@@ -1485,7 +1428,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.4"
+   "version": "3.11.6"
   },
   "vscode": {
    "interpreter": {

examples/tutorial_4_Pioneer.ipynb

@@ -769,17 +769,24 @@
     "Vel = wec_fdom.vel\n",
     "\n",
     "P_max_absorbable = (np.abs(Fex)**2/(8*Rad_res) ).squeeze().sum('omega').item() # after Falnes Eq. 6.10\n",
-    "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)\n",
+    "\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "We also calculate the power dissipated due to flywheel friction and make sure that the absorbed power (calculated as the difference between excited and radiated power) matches the sum of mechanical power captured by the PTO and the power dissipated due to flywheel friction. We use a relative tolerance of 1%."
+    "We also calculate the power dissipated due to flywheel friction and make sure that the absorbed power (calculated as the difference between excited and radiated power) matches the sum of mechanical power captured by the PTO and the power dissipated due to flywheel friction. The utilities function to plot the Sankey diagram lumps the dissipated power due to flywheel friction into the PTO loss."
    ]
   },
   {
@@ -802,7 +809,9 @@
     "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"
+    "power_flows['Electrical (solver)']= avg_elec_power  #solver determins sign\n",
+    "power_flows['Mechanical (solver)']= avg_mech_power - avg_fw_fric_power #solver determins sign, friction is dissipated\n",
+    "power_flows['PTO Loss'] = power_flows['Mechanical (solver)'] -  power_flows['Electrical (solver)'] "
    ]
   },
   {
@@ -811,39 +820,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 = [6,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.2,0.6,0.6,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)"
    ]
   },
   {
@@ -958,13 +935,6 @@
     "    axi.label_outer()\n",
     "    axi.autoscale(axis='x', tight=True)"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
@@ -983,7 +953,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.5"
+   "version": "3.11.6"
   }
  },
  "nbformat": 4,