From b8506d84ed41b9f02281deae3d0b72ce5b08a68e Mon Sep 17 00:00:00 2001 From: Daniel Gaebele <86246113+dtgaebe@users.noreply.github.com> Date: Thu, 2 May 2024 09:43:56 -0700 Subject: [PATCH] Utils dev - replaces old utilities PR (#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 Co-authored-by: Michael Devin --- examples/tutorial_1_WaveBot.ipynb | 83 +++++- examples/tutorial_2_AquaHarmonics.ipynb | 20 +- examples/tutorial_3_LUPA.ipynb | 63 +--- examples/tutorial_4_Pioneer.ipynb | 60 ++-- tests/test_utilities.py | 209 ++++++++++++++ wecopttool/__init__.py | 1 + wecopttool/core.py | 4 +- wecopttool/utilities.py | 369 ++++++++++++++++++++++++ 8 files changed, 689 insertions(+), 120 deletions(-) create mode 100644 tests/test_utilities.py create mode 100644 wecopttool/utilities.py diff --git a/examples/tutorial_1_WaveBot.ipynb b/examples/tutorial_1_WaveBot.ipynb index 25d7119c..feefdcc5 100644 --- 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 index 1eedc1f4..a32f7e10 100644 --- 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 index 1fc6c3cf..ae1fd2ac 100644 --- 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 index 6a27f931..5404931d 100644 --- 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, diff --git a/tests/test_utilities.py b/tests/test_utilities.py new file mode 100644 index 00000000..ee28e36f --- /dev/null +++ b/tests/test_utilities.py @@ -0,0 +1,209 @@ +""" Unit tests for functions in the :python:`utilities.py` module. +""" + +import pytest +import numpy as np +import xarray as xr +from matplotlib.pyplot import Figure, Axes +import wecopttool as wot +from pytest import approx +import capytaine as cpy + + + +# test function in the utilities.py + + + +@pytest.fixture(scope="module") +def power_flows(): + """Dictionary of power flows.""" + pflows = {'Optimal Excitation': -100, + 'Radiated': -20, + 'Actual Excitation': -70, + 'Electrical (solver)': -40, + 'Mechanical (solver)': -50, + 'Absorbed': -50, + 'Unused Potential': -30, + 'PTO Loss': -10 + } + return pflows + +@pytest.fixture(scope="module") +def f1(): + """Fundamental frequency [Hz].""" + return 0.1 + + +@pytest.fixture(scope="module") +def nfreq(): + """Number of frequencies in frequency vector.""" + return 5 + +@pytest.fixture(scope="module") +def ndof(): + """Number of degrees of freedom.""" + return 2 + +@pytest.fixture(scope="module") +def ndir(): + """Number of wave directions.""" + return 3 + +@pytest.fixture(scope='module') +def bem_data(f1, nfreq, ndof, ndir): + """Synthetic BEM data.""" + # TODO - start using single BEM solution across entire test suite + coords = { + 'omega': [2*np.pi*(ifreq+1)*f1 for ifreq in range(nfreq)], + 'influenced_dof': [f'DOF_{idof+1}' for idof in range(ndof)], + 'radiating_dof': [f'DOF_{idof+1}' for idof in range(ndof)], + 'wave_direction': [2*np.pi/ndir*idir for idir in range(ndir)], + } + radiation_dims = ['omega', 'radiating_dof', 'influenced_dof'] + excitation_dims = ['omega', 'influenced_dof', 'wave_direction'] + hydrostatics_dims = ['radiating_dof', 'influenced_dof'] + + added_mass = np.ones([nfreq, ndof, ndof]) + radiation_damping = np.ones([nfreq, ndof, ndof]) + diffraction_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j + Froude_Krylov_force = np.ones([nfreq, ndof, ndir], dtype=complex) + 1j + inertia_matrix = np.ones([ndof, ndof]) + hydrostatic_stiffness = np.ones([ndof, ndof]) + + data_vars = { + 'added_mass': (radiation_dims, added_mass), + 'radiation_damping': (radiation_dims, radiation_damping), + 'diffraction_force': (excitation_dims, diffraction_force), + 'Froude_Krylov_force': (excitation_dims, Froude_Krylov_force), + 'inertia_matrix': (hydrostatics_dims, inertia_matrix), + 'hydrostatic_stiffness': (hydrostatics_dims, hydrostatic_stiffness) + } + return xr.Dataset(data_vars=data_vars, coords=coords) + +@pytest.fixture(scope='module') +def intrinsic_impedance(bem_data): + bem_data = wot.add_linear_friction(bem_data) + intrinsic_impedance = wot.hydrodynamic_impedance(bem_data) + return intrinsic_impedance + +@pytest.fixture(scope='module') +def pi_controller_pto(): + """Basic PTO: proportional-integral (PI) controller, 1DOF, mechanical + power.""" + ndof = 1 + pto = wot.pto.PTO(ndof=ndof, kinematics=np.eye(ndof), + controller=wot.pto.controller_pi, + names=["PI controller PTO"]) + return pto + +@pytest.fixture(scope='module') +def regular_wave(f1, nfreq): + """Single frequency wave""" + wfreq = 0.3 + wamp = 0.0625 + wphase = 0 + wdir = 0 + waves = wot.waves.regular_wave(f1, nfreq, wfreq, wamp, wphase, wdir) + return waves + +@pytest.fixture(scope="module") +def fb(): + """Capytaine FloatingBody object""" + try: + import wecopttool.geom as geom + except ImportError: + pytest.skip( + 'Skipping integration tests due to missing optional geometry ' + + 'dependencies. Run `pip install wecopttool[geometry]` to run ' + + 'these tests.' + ) + mesh_size_factor = 0.5 + wb = geom.WaveBot() + mesh = wb.mesh(mesh_size_factor) + fb = cpy.FloatingBody.from_meshio(mesh, name="WaveBot") + fb.add_translation_dof(name="Heave") + return fb + + +@pytest.fixture(scope="module") +def wb_bem(f1, nfreq, fb): + """Boundary elemement model (Capytaine) results""" + freq = wot.frequency(f1, nfreq, False) + return wot.run_bem(fb, freq) + +@pytest.fixture(scope='class') +def wb_hydro_impedance(wb_bem): + """Intrinsic hydrodynamic impedance""" + hd = wot.add_linear_friction(wb_bem) + hd = wot.check_radiation_damping(hd) + Zi = wot.hydrodynamic_impedance(hd) + return Zi + + + + +def test_plot_hydrodynamic_coefficients(bem_data,ndof): + bem_figure_list = wot.utilities.plot_hydrodynamic_coefficients(bem_data) + correct_len = ndof*(ndof+1)/2 #using only the subdiagonal elements + #added mass + fig_am = bem_figure_list[0][0] + assert correct_len == len(fig_am.axes) + assert isinstance(fig_am,Figure) + #radiation damping + fig_rd = bem_figure_list[1][0] + assert correct_len == len(fig_rd.axes) + assert isinstance(fig_rd,Figure) + #radiation damping + fig_ex = bem_figure_list[2][0] + assert ndof == len(fig_ex.axes) + assert isinstance(fig_ex,Figure) + +def test_plot_bode_impedance(intrinsic_impedance, ndof): + fig_Zi, axes_Zi = wot.utilities.plot_bode_impedance(intrinsic_impedance) + + assert 2*ndof*ndof == len(fig_Zi.axes) + assert isinstance(fig_Zi,Figure) + assert all([isinstance(ax, Axes) for ax in np.reshape(axes_Zi,-1)]) + + +def test_plot_power_flow(power_flows): + fig_sankey, ax_sankey = wot.utilities.plot_power_flow(power_flows) + + assert isinstance(fig_sankey, Figure) + assert isinstance(ax_sankey, Axes) + +def test_calculate_power_flow(wb_bem, + regular_wave, + pi_controller_pto, + wb_hydro_impedance): + """PI controller matches optimal for any regular wave, + thus we check if the radiated power is equal the absorber power + and if the Optimal excitation is equal the actual excitation""" + + f_add = {"PTO": pi_controller_pto.force_on_wec} + wec = wot.WEC.from_bem(wb_bem, f_add=f_add) + + res = wec.solve(waves=regular_wave, + obj_fun=pi_controller_pto.average_power, + nstate_opt=2, + x_wec_0=1e-1*np.ones(wec.nstate_wec), + x_opt_0=[-1e3, 1e4], + scale_x_wec=1e2, + scale_x_opt=1e-3, + scale_obj=1e-2, + optim_options={'maxiter': 50}, + bounds_opt=((-1e4, 0), (0, 2e4),) + ) + + pflows = wot.utilities.calculate_power_flows(wec, + pi_controller_pto, + res[0], + regular_wave.sel(realization = 0), + wb_hydro_impedance) + + assert pflows['Absorbed'] == approx(pflows['Radiated'], rel=1e-4) + assert pflows['Optimal Excitation'] == approx(pflows['Actual Excitation'], rel=1e-4) + + + diff --git a/wecopttool/__init__.py b/wecopttool/__init__.py index 9b3d6b66..3cd72035 100644 --- a/wecopttool/__init__.py +++ b/wecopttool/__init__.py @@ -32,6 +32,7 @@ from wecopttool.core import * from wecopttool import waves from wecopttool import pto +from wecopttool import utilities try: from wecopttool import geom diff --git a/wecopttool/core.py b/wecopttool/core.py index 2db40b24..5ddfb9ff 100644 --- a/wecopttool/core.py +++ b/wecopttool/core.py @@ -2288,11 +2288,11 @@ def hydrodynamic_impedance(hydro_data: Dataset) -> Dataset: :py:func:`wecopttool.add_linear_friction`. """ - Zi = (hydro_data['inertia_matrix'] \ + hydro_impedance = (hydro_data['inertia_matrix'] \ + hydro_data['added_mass'])*1j*hydro_data['omega'] \ + hydro_data['radiation_damping'] + hydro_data['friction'] \ + hydro_data['hydrostatic_stiffness']/1j/hydro_data['omega'] - return Zi.transpose('omega', 'radiating_dof', 'influenced_dof') + return hydro_impedance.transpose('omega', 'radiating_dof', 'influenced_dof') def atleast_2d(array: ArrayLike) -> ndarray: diff --git a/wecopttool/utilities.py b/wecopttool/utilities.py new file mode 100644 index 00000000..ff013d5c --- /dev/null +++ b/wecopttool/utilities.py @@ -0,0 +1,369 @@ +"""Functions that are useful for WEC analysis and design. +""" + + +from __future__ import annotations + + +__all__ = [ + "plot_hydrodynamic_coefficients", + "plot_bode_impedance", + "calculate_power_flows", + "plot_power_flow", +] + + +from typing import Optional, Union +import logging +from pathlib import Path + +import autograd.numpy as np +from autograd.numpy import ndarray + +from xarray import DataArray +from numpy.typing import ArrayLike +# from autograd.numpy import ndarray +from xarray import DataArray, concat +import matplotlib.pyplot as plt +from matplotlib.figure import Figure +from matplotlib.axes import Axes + +from matplotlib.sankey import Sankey + + + +# logger +_log = logging.getLogger(__name__) + + + + +def plot_hydrodynamic_coefficients(bem_data, + wave_dir: Optional[float] = 0.0 + )-> list(tuple(Figure, Axes)): + """Plots hydrodynamic coefficients (added mass, radiation damping, + and wave excitation) based on BEM data. + + + Parameters + ---------- + bem_data + Linear hydrodynamic coefficients obtained using the boundary + element method (BEM) code Capytaine, with sign convention + corrected. + wave_dir + Wave direction(s) to plot the + + """ + + bem_data = bem_data.sel(wave_direction = wave_dir, method='nearest') + radiating_dofs = bem_data.radiating_dof.values + influenced_dofs = bem_data.influenced_dof.values + + # plots + fig_am, ax_am = plt.subplots( + len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs),3*len(influenced_dofs)), + squeeze=False + ) + fig_rd, ax_rd = plt.subplots( + len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs), 3*len(influenced_dofs)), + squeeze=False + ) + fig_ex, ax_ex = plt.subplots( + len(influenced_dofs), + 1, + tight_layout=True, + sharex=True, + figsize=(3, 3*len(radiating_dofs)), + squeeze=False + ) + [ax.grid(True) for axs in (ax_am, ax_rd, ax_ex) for ax in axs.flatten()] + # plot titles + fig_am.suptitle('Added Mass Coefficients', fontweight='bold') + fig_rd.suptitle('Radiation Damping Coefficients', fontweight='bold') + fig_ex.suptitle('Wave Excitation Coefficients', fontweight='bold') + + sp_idx = 0 + for i, rdof in enumerate(radiating_dofs): + for j, idof in enumerate(influenced_dofs): + sp_idx += 1 + if i == 0: + np.abs(bem_data.diffraction_force.sel(influenced_dof=idof)).plot( + ax=ax_ex[j,0], linestyle='dashed', label='Diffraction') + np.abs(bem_data.Froude_Krylov_force.sel(influenced_dof=idof)).plot( + ax=ax_ex[j,0], linestyle='dashdot', label='Froude-Krylov') + ex_handles, ex_labels = ax_ex[j,0].get_legend_handles_labels() + ax_ex[j,0].set_title(f'{idof}') + ax_ex[j,0].set_xlabel('') + ax_ex[j,0].set_ylabel('') + if j <= i: + bem_data.added_mass.sel( + radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_am[i, j]) + bem_data.radiation_damping.sel( + radiating_dof=rdof, influenced_dof=idof).plot(ax=ax_rd[i, j]) + if i == len(radiating_dofs)-1: + ax_am[i, j].set_xlabel(f'$\omega$', fontsize=10) + ax_rd[i, j].set_xlabel(f'$\omega$', fontsize=10) + ax_ex[j, 0].set_xlabel(f'$\omega$', fontsize=10) + else: + ax_am[i, j].set_xlabel('') + ax_rd[i, j].set_xlabel('') + if j == 0: + ax_am[i, j].set_ylabel(f'{rdof}', fontsize=10) + ax_rd[i, j].set_ylabel(f'{rdof}', fontsize=10) + else: + ax_am[i, j].set_ylabel('') + ax_rd[i, j].set_ylabel('') + if j == i: + ax_am[i, j].set_title(f'{idof}', fontsize=10) + ax_rd[i, j].set_title(f'{idof}', fontsize=10) + else: + ax_am[i, j].set_title('') + ax_rd[i, j].set_title('') + else: + fig_am.delaxes(ax_am[i, j]) + fig_rd.delaxes(ax_rd[i, j]) + fig_ex.legend(ex_handles, ex_labels, loc=(0.08, 0), ncol=2, frameon=False) + return [(fig_am,ax_am), (fig_rd,ax_rd), (fig_ex,ax_ex)] + +def plot_bode_impedance(impedance: DataArray, + title: Optional[str]= '', + #plot_natural_freq: Optional[bool] = False, +)-> tuple(Figure, Axes): + """Plot Bode graph from wecoptool impedance data array. + + Parameters + ---------- + impedance: DataArray + Complex impedance matrix produced by for example by + :py:func:`wecopttool.hydrodynamic_impedance`. + Dimensions: omega, radiating_dofs, influenced_dofs + title: String + Title string to be displayed in the plot. + + """ + radiating_dofs = impedance.radiating_dof.values + influenced_dofs = impedance.influenced_dof.values + mag = 20.0 * np.log10(np.abs(impedance)) + phase = np.rad2deg(np.unwrap(np.angle(impedance))) + freq = impedance.omega.values/2/np.pi + fig, axes = plt.subplots( + 2*len(radiating_dofs), + len(influenced_dofs), + tight_layout=True, + sharex=True, + figsize=(3*len(radiating_dofs), 3*len(influenced_dofs)), + squeeze=False + ) + fig.suptitle(title + ' Bode Plots', fontweight='bold') + + sp_idx = 0 + for i, rdof in enumerate(radiating_dofs): + for j, idof in enumerate(influenced_dofs): + sp_idx += 1 + axes[2*i, j].semilogx(freq, mag[:, i, j]) # Bode magnitude plot + axes[2*i+1, j].semilogx(freq, phase[:, i, j]) # Bode phase plot + axes[2*i, j].grid(True, which = 'both') + axes[2*i+1, j].grid(True, which = 'both') + if i == len(radiating_dofs)-1: + axes[2*i+1, j].set_xlabel(f'Frequency (Hz)', fontsize=10) + else: + axes[i, j].set_xlabel('') + if j == 0: + axes[2*i, j].set_ylabel(f'{rdof} \n Mag. (dB)', fontsize=10) + axes[2*i+1, j].set_ylabel(f'Phase. (deg)', fontsize=10) + else: + axes[i, j].set_ylabel('') + if i == 0: + axes[i, j].set_title(f'{idof}', fontsize=10) + else: + axes[i, j].set_title('') + return fig, axes + + + + +def calculate_power_flows(wec, + pto, + results, + waves, + intrinsic_impedance)-> dict[str, float]: + """Calculate power flows into a :py:class:`wecopttool.WEC` + and through a :py:class:`wecopttool.pto.PTO` based on the results + of :py:meth:`wecopttool.WEC.solve` for a single wave realization. + + Parameters + ---------- + wec + WEC object of :py:class:`wecopttool.WEC` + pto + PTO object of :py:class:`wecopttool.pto.PTO` + results + Results produced by :py:func:`scipy.optimize.minimize` for a single wave + realization. + waves + :py:class:`xarray.Dataset` with the structure and elements + shown by :py:mod:`wecopttool.waves`. + intrinsic_impedance: DataArray + Complex intrinsic impedance matrix produced by + :py:func:`wecopttool.hydrodynamic_impedance`. + Dimensions: omega, radiating_dofs, influenced_dofs + + + """ + wec_fdom, _ = wec.post_process(results, waves) + x_wec, x_opt = wec.decompose_state(results.x) + + #power quntities from solver + P_mech = pto.mechanical_average_power(wec, x_wec, x_opt, waves) + P_elec = pto.average_power(wec, x_wec, x_opt, waves) + + #compute analytical power flows + Fex_FD = wec_fdom.force.sel(type=['Froude_Krylov', 'diffraction']).sum('type') + Rad_res = np.real(intrinsic_impedance.squeeze()) + Vel_FD = wec_fdom.vel + + P_max, P_e, P_r = [], [], [] + + #This solution requires radiation resistance matrix Rad_res to be invertible + # TODO In the future we might want to add an entirely unconstrained solve + # for optimized mechanical power + + for om in Rad_res.omega.values: + #use frequency vector from intrinsic impedance (no zero freq) + #Eq. 6.69 + #Dofs are row vector, which is transposed in standard convention + Fe_FD_t = np.atleast_2d(Fex_FD.sel(omega = om)) + Fe_FD = np.transpose(Fe_FD_t) + R_inv = np.linalg.inv(np.atleast_2d(Rad_res.sel(omega= om))) + P_max.append((1/8)*(Fe_FD_t@R_inv)@np.conj(Fe_FD)) + #Eq.6.57 + U_FD_t = np.atleast_2d(Vel_FD.sel(omega = om)) + U_FD = np.transpose(U_FD_t) + R = np.atleast_2d(Rad_res.sel(omega= om)) + P_r.append((1/2)*(U_FD_t@R)@np.conj(U_FD)) + #Eq. 6.56 (replaced pinv(Fe)*U with U'*conj(Fe) + # as suggested in subsequent paragraph) + P_e.append((1/4)*(Fe_FD_t@np.conj(U_FD) + U_FD_t@np.conj(Fe_FD))) + + power_flows = { + 'Optimal Excitation' : -2* np.sum(np.real(P_max)),#eq 6.68 + 'Radiated': -1*np.sum(np.real(P_r)), + 'Actual Excitation': -1*np.sum(np.real(P_e)), + 'Electrical (solver)': P_elec, + 'Mechanical (solver)': P_mech, + } + + power_flows['Absorbed'] = ( + power_flows['Actual Excitation'] + - power_flows['Radiated'] + ) + power_flows['Unused Potential'] = ( + power_flows['Optimal Excitation'] + - power_flows['Actual Excitation'] + ) + power_flows['PTO Loss'] = ( + power_flows['Mechanical (solver)'] + - power_flows['Electrical (solver)'] + ) + return power_flows + + +def plot_power_flow(power_flows: dict[str, float])-> tuple(Figure, Axes): + """Plot power flow through a WEC as Sankey diagram. + + Parameters + ---------- + power_flows: dictionary + Power flow dictionary produced by for example by + :py:func:`wecopttool.utilities.calculate_power_flows`. + Required keys: 'Optimal Excitation', 'Radiated', 'Actual Excitation' + 'Electrical (solver)', 'Mechanical (solver)', + 'Absorbed', 'Unused Potential', 'PTO Loss' + + + """ + # fig = plt.figure(figsize = [8,4]) + # ax = fig.add_subplot(1, 1, 1,) + fig, ax = plt.subplots(1, 1, + tight_layout=True, + figsize=(8, 4), + ) + + # plt.viridis() + sankey = Sankey(ax=ax, + scale= -1/power_flows['Optimal Excitation'], + offset= 0, + format = '%.1f', + shoulder = 0.02, + tolerance=-1e-03*power_flows['Optimal Excitation'], + unit = 'W' + ) + + sankey.add(flows=[-1*power_flows['Optimal Excitation'], + power_flows['Unused Potential'], + power_flows['Actual Excitation']], + labels = ['Optimal Excitation', + 'Unused Potential ', + 'Excited'], + orientations=[0, -1, -0],#arrow directions, + pathlengths = [0.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.253935, 0.265254, 0.529983, 1.0) #viridis(0.2) + ) + + sankey.add(flows=[ + -1*(power_flows['Absorbed'] + power_flows['Radiated']), + power_flows['Radiated'], + power_flows['Absorbed'], + ], + labels = ['Excited', + 'Radiated', + ''], + prior= (0), + connect=(2,0), + orientations=[0, -1, -0],#arrow directions, + pathlengths = [0.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.127568, 0.566949, 0.550556, 1.0) #viridis(0.5) + ) + + sankey.add(flows=[-1*(power_flows['Mechanical (solver)']), + power_flows['PTO Loss'], + power_flows['Electrical (solver)'], + ], + labels = ['Mechanical', + 'PTO-Loss' , + 'Electrical'], + prior= (1), + connect=(2,0), + orientations=[0, -1, -0],#arrow directions, + pathlengths = [.2,0.3,0.2], + trunklength = 1.0, + edgecolor = 'None', + facecolor = (0.741388, 0.873449, 0.149561, 1.0) #viridis(0.9) + ) + + + diagrams = sankey.finish() + for diagram in diagrams: + for text in diagram.texts: + text.set_fontsize(10) + #remove text label from last entries + for diagram in diagrams[0:2]: + diagram.texts[2].set_text('') + + plt.axis("off") + # plt.show() + + return fig, ax