diff --git a/.github/workflows/pr.yml b/.github/workflows/pr.yml index 5b57220f0..509960977 100644 --- a/.github/workflows/pr.yml +++ b/.github/workflows/pr.yml @@ -101,4 +101,4 @@ jobs: run: | git fetch --force --prune --tags git fetch --force --prune origin '+refs/heads/*:refs/remotes/origin/*' - python3 docs/build_docs.py --build production + python3 docs/build_docs.py --build debug diff --git a/docs/source/conf.py b/docs/source/conf.py index 506577fbb..62d1f2b33 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -191,6 +191,11 @@ def _copy_examples() -> None: ignore=_all_but_nc, dirs_exist_ok=True, ) + shutil.copytree( + os.path.join(project_root, 'examples/M4E_inputs'), + os.path.join(examples_dst, 'M4E_inputs'), + dirs_exist_ok=True, + ) def _generate_theory_animations() -> None: diff --git a/docs/source/tutorials.rst b/docs/source/tutorials.rst index ca25f1204..5106b7f7b 100644 --- a/docs/source/tutorials.rst +++ b/docs/source/tutorials.rst @@ -13,13 +13,17 @@ The third part builds on the second to solve a design optimization problem, usin The second and third tutorials use the `AquaHarmonics`_ and `LUPA`_ WECs, respectively, and provide more robust optimization examples illustrating more complex cases of the WEC body, power take-off, applied forces, and constraints. These tutorials each have two parts, which are similar in structure to the second and third parts of Tutorial 1. - - :doc:`_examples/tutorial_2_AquaHarmonics`: Two part example with a realistic power take-off system and co-design study. You can find an extended version + - :doc:`_examples/tutorial_2_AquaHarmonics`: Two part example with a realistic power take-off system and co-design study. You can find an extended version - :doc:`_examples/tutorial_3_LUPA`: Two part example with multiple bodies, multiple degrees of freedom, irregular waves, a mooring system, and realistic constraints. The fourth tutorial uses the `Pioneer WEC` model, which includes a unique pitch resonator PTO system. This tutorial illustrates how to use WecOptTool to implement and optimize control strategies for less common PTO archetypes. - :doc:`_examples/tutorial_4_Pioneer`: Example with custom PTO physics and modeling both hydrodynamic and non-hydrodynamic degrees of freedom. +The fifth tutorial uses the `FOSWEC` model, which includes two flaps pitching relative to a floating platform. This tutorial illustrates how to use Multibody for Everybody to generate reduced-coordinate matrices and forces for WecOptTool. + + - :doc:`_examples/tutorial_5_FOSWEC`: Example with multiple coupled bodies, code-to-code verification using a static solve, a multi-dof PTO representation, and a study of top-mounted ballast effects. + .. toctree:: :maxdepth: 3 :hidden: @@ -28,7 +32,7 @@ The fourth tutorial uses the `Pioneer WEC` model, which includes a unique pitch _examples/tutorial_2_AquaHarmonics _examples/tutorial_3_LUPA _examples/tutorial_4_Pioneer - + _examples/tutorial_5_FOSWEC Simulating WEC Dynamics without optimization -------------------------------------------- diff --git a/examples/M4E_inputs/FOSWEC_M4E_inputs.py b/examples/M4E_inputs/FOSWEC_M4E_inputs.py new file mode 100644 index 000000000..e94470c1f --- /dev/null +++ b/examples/M4E_inputs/FOSWEC_M4E_inputs.py @@ -0,0 +1,28 @@ +"""FOSWEC Multibody for Everybody input file. +""" + +import numpy as np +from multibody import JointSystem, normalize_prismatic + +# Set FOSWEC parameters +platform_cg = [0, 0, -0.8] +flap_center_distance_apart = 1.44 +flap_draft = 0.59 +cg_height_above_hinge = 0.17 +flap1_cg = [-flap_center_distance_apart/2, 0, -flap_draft+cg_height_above_hinge] # from water surface/origin +flap2_cg = [flap_center_distance_apart/2, 0, -flap_draft+cg_height_above_hinge] # from water surface/origin +flap_hinge_depth = 0.59 + +# Define joints +joints = [[0, 1],[1, 2],[1, 3]] # Joint connectivity: [parent, child] +types = ['F', 'R', 'R'] # Joint types: 'R' for revolute, 'P' for prismatic, 'F' for floating +parent_cg_to_joint = [[0, platform_cg[2]],[flap1_cg[0],-platform_cg[2]-flap_hinge_depth],[flap2_cg[0],-platform_cg[2]-flap_hinge_depth]] # Vectors from parent's center-of-gravity (CG) to the joint location. +joint_to_child_cg = [[np.nan, np.nan],[0,flap_hinge_depth+flap1_cg[2]],[0,flap_hinge_depth+flap2_cg[2]]] # Vectors from the joint to the child's CG. +prismatic_direction = [[np.nan, np.nan],[np.nan, np.nan],[np.nan, np.nan]] # For prismatic joints, the direction vector; for others, [nan, nan] is used. +prismatic_direction = normalize_prismatic(prismatic_direction) + +# Create the JointSystem using the from_data class method and define initial conditions. +joint_system = JointSystem.from_data(joints, types, parent_cg_to_joint, joint_to_child_cg, prismatic_direction) # DO NOT MODIFY +_, _, _, NDOF, _ = joint_system.coordinate_finder() # DO NOT MODIFY +ic = np.zeros(2*sum(NDOF)) # Multiplied by 2 because is position and velocity +ic[2] = 0 * np.pi / 180 # Initial position of the first joint \ No newline at end of file diff --git a/examples/data/FOSWEC_WS_data.nc b/examples/data/FOSWEC_WS_data.nc new file mode 100644 index 000000000..2398ad0b8 Binary files /dev/null and b/examples/data/FOSWEC_WS_data.nc differ diff --git a/examples/tutorial_5_FOSWEC.ipynb b/examples/tutorial_5_FOSWEC.ipynb new file mode 100644 index 000000000..a48bbea25 --- /dev/null +++ b/examples/tutorial_5_FOSWEC.ipynb @@ -0,0 +1,1025 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4b3149a8", + "metadata": {}, + "source": [ + "# Tutorial 5 - FOSWEC\n", + "This tutorial models the floating oscillating surge wave energy converter (FOSWEC) and examines adding top-mounted ballast to the FOSWEC flaps. \n", + "The FOSWEC has been tested in [multiple experimental campaigns at the Oregon State University (OSU) O.H. Hinsdale Wave Research Laboratory (HWRL) Directional Wave Basin (DWB)](https://mhkdr.openei.org/submissions/350) and a [WEC-Sim numerical model](https://github.com/WEC-Sim/FOSWEC2) has been previously developed and validated against experimental data.\n", + "\n", + "The FOSWEC is a coupled multibody system consisting of two flaps pitching relative to a floating platform. \n", + "Efficiently modeling this system in WecOptTool requires a reduced-coordinate formulation. \n", + "In this tutorial, that reduction is performed using [Multibody for Everybody](https://github.com/Project-SEA-Stack/Python_Multibody_for_Everybody) (M4E), a package developed by the National Laboratory of the Rockies (NLR).\n", + "\n", + "![FOSWEC image](https://live.staticflickr.com/65535/55515103812_77d48575c1_z.jpg)\n", + "\n", + "While previous WecOptTool tutorials focus on a single degree of freedom (DOF) or a small set of coupled DOFs, many WEC archetypes have more complex and coupled dynamics that would otherwise require manually defining kinematic constraints or reduced-coordinate equations of motion. \n", + "M4E automates this process by transforming the hydrodynamic, hydrostatic, kinematic, and external force terms from global to joint coordinates.\n", + "For the FOSWEC, this reduces the planar model from 9 global coordinates to 5 joint coordinates, preserving important platform-flap coupling while reducing the number of optimization state variables significantly.\n", + "M4E eliminates the need for manually deriving reduced-coordinate equations, which can be time-consuming and error prone, and greatly improves the computational efficiency, since the WecOptTool optimization cost roughly scales with the square of the number of optimization variables.\n", + "\n", + "Using Multibody for Everybody with WecOptTool requires the following main steps:\n", + "1. Create an [M4E input file](M4E_inputs/FOSWEC_M4E_inputs.py). This defines the body properties, joint locations, and joint types and needs to be imported into the current notebook.\n", + "2. After running BEM for all rigid body degrees of freedom, call `wot.utilities.setup_from_M4E()` to assemble the reduced impedance matrix, reduced excitation coefficients, and coordinate transformation matrix.\n", + "3. Convert PTO kinematics and external force matrices to reduced coordinates using `wot.utilities.reduce_PTO_kinematics_M4E()` and `wot.utilities.reduce_damping_stiffness_M4E()`.\n", + "4. Set up the WEC object using the `from_impedance` method and solve the optimization problem. \n", + "\n", + "This tutorial consists of three parts. \n", + "First, the model is set up and converted to reduced coordinates using M4E.\n", + "Next, the unoptimized response is compared to an equivalent WEC-Sim model for code-to-code verification. \n", + "Then, the impact of adding top-mounted ballast to the flaps is analyzed.\n", + "\n", + "1. [Model setup](#1.-Model-setup)\n", + "2. [Model verification](#2.-Model-verification)\n", + "3. [Top-mounted ballast study](#3.-Top-mounted-ballast-study)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4620c0a", + "metadata": {}, + "outputs": [], + "source": [ + "import capytaine as cpy\n", + "from capytaine.io.meshio import load_from_meshio\n", + "import numpy as np\n", + "import jax.numpy as jnp\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.lines import Line2D\n", + "import xarray as xr\n", + "import gmsh, pygmsh\n", + "\n", + "import wecopttool as wot\n", + "\n", + "## set colorblind-friendly colormap for plots\n", + "plt.style.use('tableau-colorblind10')" + ] + }, + { + "cell_type": "markdown", + "id": "8cbd8440", + "metadata": {}, + "source": [ + "## 1. Model setup" + ] + }, + { + "cell_type": "markdown", + "id": "1cba365f", + "metadata": {}, + "source": [ + "### 1.1 Wave Conditions\n", + "For the initial code-to-code verification, three difference regular wave amplitudes are tested.\n", + "The first wave has a very small amplitude to ensure the response remains within the linear range of operation. \n", + "The intermediate wave amplitude matches a condition previously tested in the wave tank, while the largest wave height is included to further evaluate the effect of larger-amplitude motion on the linearized M4E formulation. \n", + "\n", + "Since we are using only regular waves and linear proportional control gains in this tutorial, the system is entirely linear.\n", + "Thus, only 2 frequencies are needed, enabling a computationally efficient optimization.\n", + "For irregular waves or any control introducing nonlinearities (e.g., unstructured), additional frequencies would be necessary to properly model the dynamics." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d2efcd37", + "metadata": {}, + "outputs": [], + "source": [ + "# Frequency object\n", + "wavefreq = 1/2.61 # Hz\n", + "f1 = wavefreq\n", + "nfreq = 2\n", + "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n", + "\n", + "# Define wave cases\n", + "wave_heights = [0.001, 0.136, 0.272] # m\n", + "wave_amplitudes = [height / 2 for height in wave_heights]\n", + "phase = 0 # degrees\n", + "wavedir = 0 # degrees\n", + "waves = [\n", + " wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n", + " for amplitude in wave_amplitudes\n", + "]" + ] + }, + { + "cell_type": "markdown", + "id": "d6bec81d", + "metadata": {}, + "source": [ + "### 1.2 WEC Geometry\n", + "\n", + "The FOSWEC consists of two flaps pitching relative to a floating platform. \n", + "The geometry is defined below including each of the three meshes, their hydrostatic properties, and the degrees of freedom (DOFs) for solving the boundary element method (BEM) problem.\n", + "\n", + "Usually, only the relevant DOFs are defined for each floating body in Capytaine.\n", + "However, since M4E is handling the transformation from global to reduced coordinates, we add all 6 rigid body DOFs defined about the center of gravity using Capytaine." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a6e18218", + "metadata": {}, + "outputs": [], + "source": [ + "# flap dimensions and properties\n", + "flap_thickness_bottom = 0.04\n", + "flap_thickness_top = 0.1\n", + "flap_center_distance_apart = 1.44 # distance between two flap centers in the direction of wave propagation\n", + "flap_width = 0.76\n", + "flap_height = 0.58\n", + "flap_draft = 0.59 # Fully submerged\n", + "cg_height_above_hinge = 0.17\n", + "flap1_cg = [-flap_center_distance_apart / 2, 0, -flap_draft + cg_height_above_hinge] # from water surface/origin\n", + "flap2_cg = [flap_center_distance_apart / 2, 0, -flap_draft + cg_height_above_hinge] # from water surface/origin\n", + "flap_mass = 23.1\n", + "flap_inertia = 1.19\n", + "\n", + "# platform dimensions and properties; consists of a frame, a DAQ box, and four columns\n", + "platform_frame_length = 1.44\n", + "platform_frame_thickness = 0.05\n", + "platform_frame_width = 1.06\n", + "platform_frame_top_depth = 0.53 + 0.1\n", + "platform_cg = [0, 0, -0.8]\n", + "platform_mass = 189.8\n", + "platform_inertia = 30\n", + "cutout_width = platform_frame_width - 2 * platform_frame_thickness\n", + "cutout_length = platform_frame_length - 2 * platform_frame_thickness\n", + "\n", + "full_width = 1.63\n", + "columns_radius = 0.1\n", + "columns_draft = 0.53 + 0.1 + 0.4\n", + "columns_xy = [flap_center_distance_apart / 2, full_width / 2 - columns_radius]\n", + "\n", + "daq_length = 0.72\n", + "daq_width = 0.6\n", + "daq_height = 0.2\n", + "daq_top_depth = 0.53 + 0.1 + 0.05\n", + "\n", + "# Define the meshes\n", + "with pygmsh.occ.Geometry() as geom:\n", + " gmsh.option.setNumber(\"Mesh.MeshSizeFactor\", 0.3)\n", + " platform = geom.add_box(\n", + " [-platform_frame_length / 2, -platform_frame_width / 2, -platform_frame_top_depth - platform_frame_thickness],\n", + " [platform_frame_length, platform_frame_width, platform_frame_thickness],\n", + " )\n", + " cutout = geom.add_box(\n", + " [-cutout_length / 2, -cutout_width / 2, -platform_frame_top_depth - platform_frame_thickness],\n", + " [cutout_length, cutout_width, platform_frame_thickness],\n", + " )\n", + " cyl1 = geom.add_cylinder(\n", + " [-columns_xy[0], -columns_xy[1], -columns_draft],\n", + " [0, 0, columns_draft + 0.01],\n", + " columns_radius,\n", + " )\n", + " cyl2 = geom.add_cylinder(\n", + " [-columns_xy[0], columns_xy[1], -columns_draft],\n", + " [0, 0, columns_draft + 0.01],\n", + " columns_radius,\n", + " )\n", + " cyl3 = geom.add_cylinder(\n", + " [columns_xy[0], -columns_xy[1], -columns_draft],\n", + " [0, 0, columns_draft + 0.01],\n", + " columns_radius,\n", + " )\n", + " cyl4 = geom.add_cylinder(\n", + " [columns_xy[0], columns_xy[1], -columns_draft],\n", + " [0, 0, columns_draft + 0.01],\n", + " columns_radius,\n", + " )\n", + " daq = geom.add_box(\n", + " [-daq_length / 2, -daq_width / 2, -daq_top_depth - daq_height],\n", + " [daq_length, daq_width, daq_height],\n", + " )\n", + " platform_frame = geom.boolean_difference(platform, cutout)\n", + " geom.boolean_union([platform_frame, daq, cyl1, cyl2, cyl3, cyl4])\n", + " platform_mesh = geom.generate_mesh()\n", + "\n", + "with pygmsh.geo.Geometry() as geom:\n", + " flap1 = geom.add_polygon(\n", + " [\n", + " [-flap_center_distance_apart / 2 - flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [-flap_center_distance_apart / 2 + flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [-flap_center_distance_apart / 2 + flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],\n", + " [-flap_center_distance_apart / 2 - flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],\n", + " ],\n", + " mesh_size=0.1,\n", + " )\n", + " geom.extrude(flap1, [0, flap_width, 0])\n", + " flap1_mesh = geom.generate_mesh()\n", + "\n", + "with pygmsh.geo.Geometry() as geom:\n", + " flap2 = geom.add_polygon(\n", + " [\n", + " [flap_center_distance_apart / 2 - flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [flap_center_distance_apart / 2 + flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [flap_center_distance_apart / 2 + flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],\n", + " [flap_center_distance_apart / 2 - flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],\n", + " ],\n", + " mesh_size=0.1,\n", + " )\n", + " geom.extrude(flap2, [0, flap_width, 0])\n", + " flap2_mesh = geom.generate_mesh()\n", + "\n", + "# define functions for creating inertia matrix and floating body\n", + "def make_inertia_matrix(body, mass, inertia):\n", + " return xr.DataArray(\n", + " data=np.diag([mass, mass, mass, inertia, inertia, inertia]),\n", + " dims=[\"influenced_dof\", \"radiating_dof\"],\n", + " coords={\n", + " \"influenced_dof\": list(body.dofs),\n", + " \"radiating_dof\": list(body.dofs),\n", + " },\n", + " name=\"inertia_matrix\",\n", + " )\n", + "\n", + "def make_floating_body(mesh, name, center_of_mass, mass, inertia):\n", + " body = cpy.FloatingBody(mesh=mesh, name=name, center_of_mass=center_of_mass)\n", + " body.add_all_rigid_body_dofs()\n", + " body.rotation_center = body.center_of_mass\n", + " body.inertia_matrix = make_inertia_matrix(body, mass, inertia)\n", + " body.hydrostatic_stiffness = body.immersed_part().compute_hydrostatic_stiffness()\n", + " return body\n", + "\n", + "# Create floating bodies for platform and flaps and combine\n", + "platform = make_floating_body(platform_mesh, \"platform\", platform_cg, platform_mass, platform_inertia)\n", + "flap1 = make_floating_body(flap1_mesh, \"flap1\", flap1_cg, flap_mass, flap_inertia)\n", + "flap2 = make_floating_body(flap2_mesh, \"flap2\", flap2_cg, flap_mass, flap_inertia)\n", + "\n", + "foswec_fb = platform + flap1 + flap2" + ] + }, + { + "cell_type": "markdown", + "id": "e7d371bc", + "metadata": {}, + "source": [ + "### 1.3 Hydrodynamics\n", + "\n", + "Next, we run the BEM using [Capytaine](https://github.com/capytaine/capytaine) to obtain the hydrodynamic coefficients for the FOSWEC.\n", + "Capytaine returns the hydrodynamic matrices and coefficients for the full 18 rigid body DOFs." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c86c8cf", + "metadata": {}, + "outputs": [], + "source": [ + "# run BEM calculations\n", + "bem_data = wot.run_bem(foswec_fb, freq, depth=2)" + ] + }, + { + "cell_type": "markdown", + "id": "9d4091f8", + "metadata": {}, + "source": [ + "### 1.4 Multibody for Everybody\n", + "\n", + "After generating BEM coefficients for all rigid body DOFs, Multibody for Everybody is called to reduce from global coordinates to joint coordinates.\n", + "M4E is only set up for planar DOFs, so it ignores the sway, roll, and yaw coefficients for each body and reduces the remaining coefficients (surge, heave, and pitch) to joint coordinates.\n", + "Thus, the FOSWEC is reduced from 9 global planar coordinates to 5 joint coordinates.\n", + "`wot.utilities.setup_from_M4E()` calls M4E, inputting the M4E input file, BEM coefficients, floating bodies, and mass proerties, and sets up the reduced-coordinate impedance matrix, excitation coefficients, and hydrostatic stiffness matrix to be used in setting up the WEC object. \n", + "The utilities function also returns `R0`, which is the transformation matrix used for coordinate conversion, and the number of M4E DOFs per body, `M4E_ndof`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "101bb65e", + "metadata": {}, + "outputs": [], + "source": [ + "from M4E_inputs import FOSWEC_M4E_inputs\n", + "\n", + "# create multibody system for wavebot using M4E\n", + "m0 = [platform_mass, flap_mass, flap_mass]\n", + "J0 = [platform_inertia, flap_inertia, flap_inertia]\n", + "\n", + "# return the reduced coordinate matrices\n", + "impedance_reduced, excitation_reduced, hydrostatic_stiffness_reduced, R0, M4E_ndof = wot.utilities.setup_from_M4E(FOSWEC_M4E_inputs, bem_data, [platform, flap1, flap2], m0, J0)" + ] + }, + { + "cell_type": "markdown", + "id": "5f8b7361", + "metadata": {}, + "source": [ + "### 1.5 Additional forces\n", + "\n", + "Additional forces acting on the FOSWEC include the mooring stiffness and damping on the platform and linear damping on the flaps.\n", + "The damping and stiffness matrices are defined in the global planar DOFs and transformed to reduced coordinates using `wot.utilities.reduce_damping_stiffness_M4E()` and the coordinate transformation matrix, `R0`.\n", + "Then, the force functions are defined normally using the reduced-coordinate matrices.\n", + "For the code-to-code verification, no PTO is included here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fa3da84", + "metadata": {}, + "outputs": [], + "source": [ + "# set mooring stiffness and damping matrices and reduce to joint coordinates\n", + "mooring_stiffness_full = np.diag([8e3, 2e5, 2e5, 0, 0, 0, 0, 0, 0])\n", + "mooring_stiffness = wot.utilities.reduce_damping_stiffness_M4E(mooring_stiffness_full, R0)\n", + "\n", + "mooring_damping_full = np.diag([8e2, 1e4, 1e4, 0, 0, 0, 0, 0, 0])\n", + "mooring_damping = wot.utilities.reduce_damping_stiffness_M4E(mooring_damping_full, R0)\n", + "\n", + "def f_mooring(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pos = wec.vec_to_dofmat(x_wec)\n", + " vel = jnp.dot(wec.derivative_mat,pos)\n", + " time_matrix = wec.time_mat_nsubsteps(nsubsteps)\n", + " mooring = -pos @ mooring_stiffness - vel @ mooring_damping\n", + " mooring_force = jnp.dot(time_matrix,mooring)\n", + " return mooring_force\n", + "\n", + "# set linear damping matrix (flap linear damping) and reduce to joint coordinates\n", + "flap_linear_damping_full = np.diag([0, 0, 0, 0, 0, 10, 0, 0, 10])\n", + "flap_linear_damping = wot.utilities.reduce_damping_stiffness_M4E(flap_linear_damping_full, R0)\n", + "\n", + "def f_linear_damping(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pos = wec.vec_to_dofmat(x_wec)\n", + " vel = jnp.dot(wec.derivative_mat,pos)\n", + " time_matrix = wec.time_mat_nsubsteps(nsubsteps)\n", + " damping = - vel @ flap_linear_damping\n", + " linear_damping_force = jnp.dot(time_matrix, damping)\n", + " return linear_damping_force\n", + "\n", + "# define added forces\n", + "f_add = {'f_mooring': f_mooring, 'f_linear_damping': f_linear_damping} " + ] + }, + { + "cell_type": "markdown", + "id": "3db9762e", + "metadata": {}, + "source": [ + "### 1.6 WEC object\n", + "\n", + "The WEC object is defined using the `from_impedance` method with the reduced-coordinate matrices and coefficients." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2496fb1a", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the WEC object using the from_impedance method\n", + "wec = wot.WEC.from_impedance(\n", + " waves[0].freq.values,\n", + " impedance=impedance_reduced,\n", + " exc_coeff=excitation_reduced,\n", + " hydrostatic_stiffness=hydrostatic_stiffness_reduced,\n", + " constraints=None,\n", + " f_add=f_add,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "fb6407aa", + "metadata": {}, + "source": [ + "### 1.7 Objective function\n", + "\n", + "For the code-to-code verification, we are not applying any PTO/control forces. \n", + "The objective function is therefore set to zero with no optimization variables." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c75fd950", + "metadata": {}, + "outputs": [], + "source": [ + "# Objective function - none for now\n", + "obj_fun = lambda wec, x_wec, x_opt, waves : 0\n", + "nstate_opt = 0" + ] + }, + { + "cell_type": "markdown", + "id": "b74e4bb2", + "metadata": {}, + "source": [ + "## 2. Model verification" + ] + }, + { + "cell_type": "markdown", + "id": "9fc94441", + "metadata": {}, + "source": [ + "### 2.1 Solve\n", + "\n", + "With the objective function set to zero, WecOptTool solves for a feasible periodic response by enforcing the dynamic residual constraint, effectively solving the multibody dynamics problem.\n", + "The problem is solved for each verification wave condition.\n", + "For post-processing, `wot.utilities.post_process_M4E()` formats the results and converts them back to the global planar coordinates." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "289a5ea5", + "metadata": {}, + "outputs": [], + "source": [ + "# Solve\n", + "scale_x_wec = 1e1\n", + "scale_x_opt = 1e-3\n", + "scale_obj = 1e-2\n", + "\n", + "results_all = []\n", + "wec_fdom_full_all = []\n", + "wec_tdom_full_all = []\n", + "\n", + "nsubsteps = 10\n", + "\n", + "for wave in waves:\n", + " # Solve\n", + " results = wec.solve(\n", + " wave,\n", + " obj_fun,\n", + " nstate_opt,\n", + " scale_x_wec=scale_x_wec,\n", + " scale_x_opt=scale_x_opt,\n", + " scale_obj=scale_obj,\n", + " )\n", + "\n", + " # Post-process results\n", + " wec_fdom_full, wec_tdom_full = wot.utilities.post_process_M4E(\n", + " wec, results, wave, nsubsteps, R0\n", + " )\n", + "\n", + " # Store outputs\n", + " results_all.append(results)\n", + " wec_fdom_full_all.append(wec_fdom_full)\n", + " wec_tdom_full_all.append(wec_tdom_full)" + ] + }, + { + "cell_type": "markdown", + "id": "cdaa03e3", + "metadata": {}, + "source": [ + "### 2.2 Code-to-code comparison\n", + "\n", + "To verify the WecOptTool+M4E model, the results are compared to an equivalent WEC-Sim model.\n", + "The small waves ensure the system remains entirely within the linear range of operation since M4E linearizes the dynamics while WEC-Sim models the nonlinear dynamics.\n", + "Within the small wave conditions, a near perfect match is achieved between the WecOptTool+M4E and WEC-Sim models. \n", + "As the wave height is increased, the difference between the WecOptTool and WEC-Sim results increases. \n", + "Still, the relatively low error across the tested wave conditions verifies the WecOptTool+M4E model and supporting its use in co-design and optimization studies." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1fb37e6f", + "metadata": {}, + "outputs": [], + "source": [ + "# Load combined WEC-Sim results\n", + "ws_data = xr.open_dataset(\"data/FOSWEC_WS_data.nc\")\n", + "\n", + "for i in range(ws_data.sizes[\"case\"]):\n", + " ws_output = ws_data.isel(case=i)\n", + " wave_height = ws_output.wave_height.values\n", + " wec_tdom_full = wec_tdom_full_all[i]\n", + "\n", + " fig, axs = plt.subplots(2, 1, sharex=True, figsize=(6, 7),gridspec_kw={\"hspace\": 0.05})\n", + "\n", + " ws_time_period = ws_output.time\n", + " wec_time_period = wec_tdom_full[\"time\"]\n", + "\n", + " # Platform response\n", + " axs[0].plot(wec_time_period,wec_tdom_full.sel(realization=0).pos[0, :],\"C0-\",linewidth=1,label=\"surge, WecOptTool\",)\n", + " axs[0].plot(ws_time_period,ws_output.platform_response[:, 0],\"C0--\",linewidth=2,label=\"surge, WEC-Sim\",)\n", + " axs[0].plot(wec_time_period,wec_tdom_full.sel(realization=0).pos[1, :],\"C1-\",linewidth=1,label=\"heave, WecOptTool\",)\n", + " axs[0].plot(ws_time_period,ws_output.platform_response[:, 2] - platform_cg[2],\"C1--\",linewidth=2,label=\"heave, WEC-Sim\",)\n", + " axs[0].plot(wec_time_period,wec_tdom_full.sel(realization=0).pos[2, :],\"C2-\",linewidth=1,label=\"pitch, WecOptTool\",)\n", + " axs[0].plot(ws_time_period,ws_output.platform_response[:, 4],\"C2--\",linewidth=2,label=\"pitch, WEC-Sim\",)\n", + " axs[0].set_ylabel(\"Platform response [m, rad]\")\n", + " axs[0].set_title(f\"H = {wave_height:.3f} m\")\n", + "\n", + " # Flap response\n", + " axs[1].plot(wec_time_period,wec_tdom_full.sel(realization=0).pos[5, :],color=\"C0\",linewidth=1,label=\"WecOptTool - flap 1\",)\n", + " axs[1].plot(wec_time_period,wec_tdom_full.sel(realization=0).pos[8, :],color=\"C1\",linewidth=1,label=\"WecOptTool - flap 2\",)\n", + " axs[1].plot(ws_time_period,ws_output.pto1_response[:, 4],\"--\",color=\"C0\",linewidth=2,label=\"WEC-Sim - flap 1\",)\n", + " axs[1].plot(ws_time_period,ws_output.pto2_response[:, 4],\"--\",color=\"C1\",linewidth=2,label=\"WEC-Sim - flap 2\",)\n", + " axs[1].set_xlabel(\"Time [s]\")\n", + " axs[1].set_ylabel(\"Flap pitch [rad]\")\n", + " axs[1].legend()\n", + "\n", + " plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "3f505714", + "metadata": {}, + "source": [ + "## 3. Top-mounted ballast study" + ] + }, + { + "cell_type": "markdown", + "id": "f334fb0a", + "metadata": {}, + "source": [ + "### 3.1 Top-mounted ballast impact on hydrostatic stiffness\n", + "\n", + "Adding ballast to the top of the flaps moves their center of gravity upward (increase center of gravity height).\n", + "As shown below, this reduces the distance between the center of gravity and center of buoyancy, creating a destabilizing effect that reduces the hydrostatic stiffness.\n", + "The reduced hydrostatic stiffness may increase the WEC response and improve power capture, particularly if the flap rotation constraints are not already active.\n", + "\n", + "\"flap_cg_gs\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8151a1b7", + "metadata": {}, + "outputs": [], + "source": [ + "# plot mass on top effects\n", + "top_mounted_ballast_vec = np.linspace(0, 5, 11)\n", + "\n", + "# Define mesh for single centered flap\n", + "with pygmsh.geo.Geometry() as geom:\n", + " flap = geom.add_polygon(\n", + " [[-flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [flap_thickness_bottom / 2, -flap_width / 2, -flap_draft],\n", + " [flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],\n", + " [-flap_thickness_top / 2, -flap_width / 2, flap_height - flap_draft],],mesh_size=0.1,)\n", + " geom.extrude(flap, [0, flap_width, 0])\n", + " flap_mesh = geom.generate_mesh()\n", + "\n", + "# determine new center of gravity height and hydrostatic stiffness for different top-mounted ballast values\n", + "cg_height_vec = []\n", + "hydrostatic_stiffness_vec = []\n", + "for top_mounted_ballast in top_mounted_ballast_vec:\n", + " # Define FloatingBody\n", + " cg_height_above_hinge_new = (cg_height_above_hinge*flap_mass + top_mounted_ballast*flap_height)/(flap_mass + top_mounted_ballast)\n", + " flap_cg = [0, 0, -flap_draft + cg_height_above_hinge_new]\n", + " flap = cpy.FloatingBody(mesh=flap_mesh, name=\"flap\", center_of_mass=flap_cg)\n", + "\n", + " # Define rotation dof\n", + " flap.rotation_center = (0, 0, -flap_draft)\n", + " flap.add_rotation_dof(name='Pitch')\n", + "\n", + " # Set FloatingBody inertia matrix\n", + " flap.mass = flap_mass\n", + " flap_inertia_hinge = flap_inertia + flap_mass*cg_height_above_hinge**2\n", + " rigid_inertia_matrix_xr = xr.DataArray(data=np.asarray((np.diag([flap_inertia_hinge]))),\n", + " dims=['influenced_dof', 'radiating_dof'],\n", + " coords={'influenced_dof': list(flap.dofs),\n", + " 'radiating_dof': list(flap.dofs)},\n", + " name=\"inertia_matrix\")\n", + " flap.inertia_matrix = rigid_inertia_matrix_xr\n", + "\n", + " # record cg height and hydrostatic stiffness\n", + " flap.hydrostatic_stiffness = flap.immersed_part().compute_hydrostatic_stiffness()\n", + " cg_height_vec.append(cg_height_above_hinge_new)\n", + " hydrostatic_stiffness_vec.append(flap.hydrostatic_stiffness)\n", + "\n", + "fig, axs = plt.subplots(2, 1, sharex=True, figsize=(6, 5))\n", + "\n", + "axs[0].plot(top_mounted_ballast_vec, cg_height_vec)\n", + "axs[0].set_ylabel(\"Center of gravity height [m]\")\n", + "axs[0].grid(True)\n", + "\n", + "axs[1].plot(top_mounted_ballast_vec, np.squeeze(hydrostatic_stiffness_vec))\n", + "axs[1].set_xlabel(\"Top-mounted ballast [kg]\")\n", + "axs[1].set_ylabel(\"Hydrostatic stiffness [N·m/rad]\")\n", + "axs[1].grid(True)\n" + ] + }, + { + "cell_type": "markdown", + "id": "506ef00c", + "metadata": {}, + "source": [ + "### 3.2 Top-mounted ballast study loop\n", + "\n", + "The top-mounted ballast study loops through ballast mass, wave period, and wave amplitude.\n", + "Capytaine is rerun for each ballast mass and wave period to generate the hydrodynamic coefficients using the updated mass properties and required frequencies.\n", + "The PTO, additional forces, and constraints are defined for each case, and the optimization is run to determine the optimal controller gains and added stiffness coefficients.\n", + "\n", + "*Note that a few details differ from the corresponding paper, including the case discretization and controller formulation: this tutorial uses a proportional controller, while the paper uses an unstructured controller.*\n", + "\n", + "#### 3.2.1 PTO\n", + "\n", + "Since no PTO was included in the verification setup, the PTO dynamics, kinematics, and constraints are now included here.\n", + "Although there are two PTO units in reality, they can be represented together by one kinematics matrix and one impedance matrix.\n", + "The PTOs are connected between the flap and the platform and resist the relative pitch motion between them. \n", + "\n", + "The PTO kinematics matrix maps the global planar coordinates to the two PTO coordinates: flap 1 pitch relative to the platform and flap 2 pitch relative to the platform. \n", + "This kinematics matrix is reduced to joint coordinates using `wot.utilities.reduce_PTO_kinematics_M4E()`.\n", + "\n", + "The PTO dynamics are defined using the two-port impedance model based on the drivetrain and generator parameters.\n", + "Because there are two PTOs, the combined PTO impedance matrix is 4x4: two mechanical ports and two electrical ports.\n", + "In this formulation, rows/columns 1 and 3 correspond to the mechanical and electrical ports of the first PTO, while rows/columns 2 and 4 correspond to the mechanical and electrical ports of the second PTO.\n", + "\n", + "#### 3.2.2 Optimal drivetrain stiffness\n", + "\n", + "Along with optimizing the PTO controller gains, an optimized stiffness coefficient is added to each flap to represent an added torsional spring. \n", + "The added stiffness coefficients give the optimizer an additional design variable that can compensate for changes in hydrostatic stiffness while still satisfying the maximum flap rotation constraints.\n", + "\n", + "#### 3.2.3 Controller\n", + "\n", + "A proportional controller is used here for efficient optimization, providing linear damping control. There are two proportional gains: one for each PTO. \n", + "\n", + "#### 3.2.4 Constraints\n", + "\n", + "The constraints include rotational constraints for each flap and peak and rms generator torque limits.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfdc3584", + "metadata": {}, + "outputs": [], + "source": [ + "top_mounted_ballast_vec = np.linspace(0, 2, 3)\n", + "wave_periods = [1, 3, 6]\n", + "wave_amplitudes = [0.05, 0.15, 0.25]\n", + "\n", + "case_outputs = []\n", + "\n", + "for top_mounted_ballast in top_mounted_ballast_vec:\n", + "\n", + " # redefine floating bodies with top-mounted ballast\n", + " platform = make_floating_body(platform_mesh, \"platform\", platform_cg, platform_mass, platform_inertia)\n", + " total_mass = flap_mass + top_mounted_ballast\n", + " cg_height_above_hinge_new = (flap_mass * cg_height_above_hinge + top_mounted_ballast * flap_height) / total_mass\n", + " flap1_cg = [-flap_center_distance_apart / 2, 0, -flap_draft + cg_height_above_hinge_new]\n", + " flap1 = make_floating_body(flap1_mesh, \"flap1\", flap1_cg, flap_mass, flap_inertia)\n", + " flap2_cg = [flap_center_distance_apart / 2, 0, -flap_draft + cg_height_above_hinge_new]\n", + " flap2 = make_floating_body(flap2_mesh, \"flap2\", flap2_cg, flap_mass, flap_inertia)\n", + " foswec_fb = platform + flap1 + flap2\n", + " \n", + " for wave_period in wave_periods:\n", + "\n", + " # define frequency object\n", + " wavefreq = 1/wave_period # Hz\n", + " f1 = wavefreq\n", + " nfreq = 2\n", + " freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n", + "\n", + " # run BEM calculations for wave period\n", + " bem_data = wot.run_bem(foswec_fb, freq, depth=2)\n", + "\n", + " # create multibody system for wavebot using M4E\n", + " m0 = [platform_mass, flap_mass, flap_mass]\n", + " J0 = [platform_inertia, flap_inertia, flap_inertia]\n", + "\n", + " # return the reduced coordinate matrices\n", + " impedance_reduced, excitation_reduced, hydrostatic_stiffness_reduced, R0, M4E_ndof = wot.utilities.setup_from_M4E(FOSWEC_M4E_inputs, bem_data, [platform, flap1, flap2], m0, J0)\n", + "\n", + " for wave_amplitude in wave_amplitudes:\n", + "\n", + " # define wave case\n", + " amplitude = wave_amplitude\n", + " phase = 0 # degrees\n", + " wavedir = 0 # degrees\n", + " waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n", + "\n", + " ## PTO impedance definition\n", + " omega = bem_data.omega.values\n", + " gear_ratio = 3.75\n", + " torque_constant = 1.021\n", + " winding_resistance = 1.028\n", + " winding_inductance = 0.0\n", + " drivetrain_inertia = 0.05\n", + " drivetrain_friction_aft = 3.27 / 3.75**2\n", + " drivetrain_friction_bow = 2.7 / 3.75**2\n", + " drivetrain_stiffness = 0.0\n", + "\n", + " drivetrain_impedance_aft = (\n", + " 1j * omega * drivetrain_inertia\n", + " + drivetrain_friction_aft\n", + " + 1 / (1j * omega) * drivetrain_stiffness\n", + " )\n", + " drivetrain_impedance_bow = (\n", + " 1j * omega * drivetrain_inertia\n", + " + drivetrain_friction_bow\n", + " + 1 / (1j * omega) * drivetrain_stiffness\n", + " )\n", + "\n", + " winding_impedance = winding_resistance + 1j * omega * winding_inductance\n", + "\n", + " pto_impedance_11_aft = -1 * gear_ratio**2 * drivetrain_impedance_aft\n", + " pto_impedance_11_bow = -1 * gear_ratio**2 * drivetrain_impedance_bow\n", + " off_diag = np.sqrt(3.0 / 2.0) * torque_constant * gear_ratio\n", + " pto_impedance_12 = -1 * (off_diag + 0j) * np.ones(omega.shape)\n", + " pto_impedance_21 = -1 * (off_diag + 0j) * np.ones(omega.shape)\n", + " pto_impedance_22 = winding_impedance\n", + "\n", + " pto_impedance = np.zeros((4, 4, len(omega)), dtype=complex)\n", + " pto_impedance[0, 0, :] = pto_impedance_11_aft\n", + " pto_impedance[1, 1, :] = pto_impedance_11_bow\n", + " pto_impedance[0, 2, :] = pto_impedance_12\n", + " pto_impedance[1, 3, :] = pto_impedance_12\n", + " pto_impedance[2, 0, :] = pto_impedance_21\n", + " pto_impedance[3, 1, :] = pto_impedance_21\n", + " pto_impedance[2, 2, :] = pto_impedance_22\n", + " pto_impedance[3, 3, :] = pto_impedance_22\n", + "\n", + " # PTO settings\n", + " pto_ndof = 2\n", + " controller = wot.controllers.pid_controller(ndof_pto=pto_ndof, proportional=True, integral=False, derivative=False)\n", + " loss = None\n", + " names = [\"flap 1 PTO\",\"flap 2 PTO\"]\n", + "\n", + " # PTO kinematics\n", + " pto_kinematics_full = np.array(\n", + " [[0, 0, -1, 0, 0, 1, 0, 0, 0],\n", + " [0, 0, -1, 0, 0, 0, 0, 0, 1],]\n", + " )\n", + " pto_kinematics_reduced = wot.utilities.reduce_PTO_kinematics_M4E(pto_kinematics_full, R0)\n", + "\n", + " pto = wot.pto.PTO(pto_ndof, pto_kinematics_reduced, controller, pto_impedance, loss, names)\n", + "\n", + " mooring_stiffness = wot.utilities.reduce_damping_stiffness_M4E(mooring_stiffness_full, R0)\n", + " mooring_damping = wot.utilities.reduce_damping_stiffness_M4E(mooring_damping_full, R0)\n", + " flap_linear_damping = wot.utilities.reduce_damping_stiffness_M4E(flap_linear_damping_full, R0)\n", + "\n", + " def f_mooring(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pos = wec.vec_to_dofmat(x_wec)\n", + " vel = jnp.dot(wec.derivative_mat, pos)\n", + " time_matrix = wec.time_mat_nsubsteps(nsubsteps)\n", + " force = -pos @ mooring_stiffness - vel @ mooring_damping\n", + " return jnp.dot(time_matrix, force)\n", + " def f_linear_damping(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pos = wec.vec_to_dofmat(x_wec)\n", + " vel = jnp.dot(wec.derivative_mat, pos)\n", + " time_matrix = wec.time_mat_nsubsteps(nsubsteps)\n", + " force = -vel @ flap_linear_damping\n", + " return jnp.dot(time_matrix, force)\n", + " def f_drivetrain_stiffness_flap1(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pto_pos = pto.position(wec, x_wec, x_opt, wave, nsubsteps)\n", + " flap1_stiffness = x_opt[-2]\n", + " spring_force_pto = jnp.zeros_like(pto_pos)\n", + " spring_force_pto = spring_force_pto.at[:, 0].set(-flap1_stiffness * pto_pos[:, 0])\n", + " spring_force_wec = jnp.dot(spring_force_pto, pto_kinematics_reduced)\n", + " return spring_force_wec\n", + " def f_drivetrain_stiffness_flap2(wec, x_wec, x_opt, wave, nsubsteps=1):\n", + " pto_pos = pto.position(wec, x_wec, x_opt, wave, nsubsteps)\n", + " flap2_stiffness = x_opt[-1]\n", + " spring_force_pto = jnp.zeros_like(pto_pos)\n", + " spring_force_pto = spring_force_pto.at[:, 1].set(-flap2_stiffness * pto_pos[:, 1])\n", + " spring_force_wec = jnp.dot(spring_force_pto, pto_kinematics_reduced)\n", + " return spring_force_wec\n", + "\n", + " f_add = {\n", + " \"f_mooring\": f_mooring,\n", + " \"f_linear_damping\": f_linear_damping,\n", + " \"f_drivetrain_stiffness_flap1\": f_drivetrain_stiffness_flap1,\n", + " \"f_drivetrain_stiffness_flap2\": f_drivetrain_stiffness_flap2,\n", + " }\n", + "\n", + " # Constraints\n", + " max_flap_pitch = np.deg2rad(30.0)\n", + " max_pto_torque = 40.0 \n", + " max_rms_pto_torque = 20.0\n", + "\n", + " def const_flap1_rotation(wec, x_wec, x_opt, wave, nsubsteps=nsubsteps):\n", + " pto_pos = pto.position(wec, x_wec, x_opt, wave, nsubsteps)\n", + " flap1_pos = pto_pos[:, 0]\n", + " return max_flap_pitch - jnp.abs(flap1_pos.flatten())\n", + " def const_flap2_rotation(wec, x_wec, x_opt, wave, nsubsteps=nsubsteps):\n", + " pto_pos = pto.position(wec, x_wec, x_opt, wave, nsubsteps)\n", + " flap2_pos = pto_pos[:, 1]\n", + " return max_flap_pitch - jnp.abs(flap2_pos.flatten())\n", + " def const_motor_torque(wec, x_wec, x_opt, wave):\n", + " pto_force = pto.force(wec, x_wec, x_opt, wave, nsubsteps)\n", + " motor_torque = pto_force / gear_ratio\n", + " return max_pto_torque - jnp.abs(motor_torque.flatten())\n", + " def const_motor_rms_torque(wec, x_wec, x_opt, wave):\n", + " pto_force = pto.force(wec, x_wec, x_opt, wave, nsubsteps)\n", + " motor_torque = pto_force / gear_ratio\n", + " rms_by_flap = jnp.sqrt(jnp.mean(motor_torque**2, axis=0) + 1e-12)\n", + " return max_rms_pto_torque - rms_by_flap\n", + "\n", + " constraints = [\n", + " {\"type\": \"ineq\", \"fun\": const_flap1_rotation},\n", + " {\"type\": \"ineq\", \"fun\": const_flap2_rotation},\n", + " {\"type\": \"ineq\", \"fun\": const_motor_torque},\n", + " {\"type\": \"ineq\", \"fun\": const_motor_rms_torque},\n", + " ]\n", + "\n", + " # define optimization states and bounds\n", + " n_force_states = 1 * 2 # 1 P gain * 2 PTO DOFs\n", + " nstate_opt = n_force_states + 2 # 2 optimal stiffnesses for the two flaps\n", + " bounds_opt = (((-1e8, 1e8),) * n_force_states + ((0.0, 1e4),) + ((0.0, 1e4),))\n", + "\n", + " # create WEC object\n", + " wec = wot.WEC.from_impedance(\n", + " waves.freq.values,\n", + " impedance=impedance_reduced,\n", + " exc_coeff=excitation_reduced,\n", + " hydrostatic_stiffness=hydrostatic_stiffness_reduced,\n", + " constraints=constraints,\n", + " f_add=f_add,\n", + " )\n", + "\n", + " obj_fun = pto.average_power\n", + "\n", + " wave_coeffs = waves.sel(realization=0, wave_direction=0).values\n", + " x_wec_mat = np.zeros((2 * len(wave_coeffs), 5))\n", + " x_wec_mat[0::2, :] = np.real(wave_coeffs)[:, None]\n", + " x_wec_mat[1::2, :] = np.imag(wave_coeffs)[:, None]\n", + " x_wec_0 = x_wec_mat.flatten() # Set x_wec_0 based on wave coefficients for initial guess of WEC states\n", + " x_opt_0 = [0, 0, 250, 250] # Set last 2 elements of x_opt_0 to 250 for initial guess of optimal stiffnesses\n", + "\n", + " results = wec.solve(\n", + " waves,\n", + " obj_fun,\n", + " nstate_opt,\n", + " scale_x_wec=1e1,\n", + " scale_x_opt=1e-2,\n", + " scale_obj=1.0,\n", + " x_wec_0=x_wec_0,\n", + " x_opt_0=x_opt_0,\n", + " bounds_opt=bounds_opt,\n", + " optim_options={\"maxiter\": 500},\n", + " )\n", + "\n", + " nsubsteps_post = 10\n", + " wec_fdom_full, wec_tdom_full = wot.utilities.post_process_M4E(wec, results, waves, nsubsteps_post, R0)\n", + " pto_fdom, pto_tdom = pto.post_process(wec, results, waves, nsubsteps=nsubsteps_post)\n", + "\n", + " avg_power = np.mean(pto_tdom[\"power\"][0, 1, :, 0].values + pto_tdom[\"power\"][0, 1, :, 1].values)\n", + "\n", + " x_full = np.asarray(results[0].x)\n", + " nstate_wec = int(wec.nstate_wec)\n", + " x_opt_sol = np.array(x_full[nstate_wec:], copy=True)\n", + " opt_stiffness_flap1 = x_opt_sol[-2]\n", + " opt_stiffness_flap2 = x_opt_sol[-1]\n", + "\n", + " pos = np.squeeze(pto_tdom[\"pos\"].sel(realization=0).values)\n", + " max_abs_pos = np.max(np.abs(pos), axis=0)\n", + " flap1_peak_rotation_utilization = max_abs_pos[0] / max_flap_pitch\n", + " flap2_peak_rotation_utilization = max_abs_pos[1] / max_flap_pitch\n", + "\n", + " case_outputs.append(\n", + " {\n", + " \"top_mounted_ballast\": top_mounted_ballast,\n", + " \"wave_period\": wave_period,\n", + " \"wave_amplitude\": wave_amplitude,\n", + " \"electrical_power\": -avg_power, # positive generated power\n", + " \"flap1_peak_rotation_utilization\": flap1_peak_rotation_utilization,\n", + " \"flap2_peak_rotation_utilization\": flap2_peak_rotation_utilization,\n", + " \"opt_stiffness_flap1\": opt_stiffness_flap1,\n", + " \"opt_stiffness_flap2\": opt_stiffness_flap2,\n", + " \"results\": results,\n", + " \"wec_fdom_full\": wec_fdom_full,\n", + " \"wec_tdom_full\": wec_tdom_full,\n", + " \"pto_fdom\": pto_fdom,\n", + " \"pto_tdom\": pto_tdom,\n", + " }\n", + " )" + ] + }, + { + "cell_type": "markdown", + "id": "6a04e28a", + "metadata": {}, + "source": [ + "### 3.3 Top-mounted ballast study results\n", + "\n", + "The results of the top-mounted ballast study are summarized in the plots below. \n", + "Adding ballast can improve average power by reducing the flap hydrostatic stiffness and increasing flap motion. \n", + "However, the benefits are limited by the peak flap pitch contraints.\n", + "As a result, top-mounted ballast is most beneficial in smaller-amplitude waves, while larger-amplitude waves may already reach the rotation constraint and therefore see little additional benefit from destabilization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d56981a", + "metadata": {}, + "outputs": [], + "source": [ + "fig, axs = plt.subplots(\n", + " 4,\n", + " len(wave_periods),\n", + " sharex=\"col\",\n", + " sharey=\"row\",\n", + " figsize=(2.6 * len(wave_periods), 9),\n", + " gridspec_kw={\"hspace\": 0.15},\n", + " constrained_layout=True,\n", + " squeeze=False,\n", + ")\n", + "\n", + "colors = [\"#9ecae1\", \"#4292c6\", \"#08306b\"]\n", + "markers = [\"o\", \"s\", \"^\"]\n", + "\n", + "for j, period in enumerate(wave_periods):\n", + " for amplitude, color, marker in zip(wave_amplitudes, colors, markers):\n", + "\n", + " case_list = sorted(\n", + " [\n", + " c for c in case_outputs\n", + " if np.isclose(c[\"wave_period\"], period)\n", + " and np.isclose(c[\"wave_amplitude\"], amplitude)\n", + " ],\n", + " key=lambda c: c[\"top_mounted_ballast\"],\n", + " )\n", + "\n", + " if len(case_list) == 0:\n", + " continue\n", + "\n", + " ballast = np.array([c[\"top_mounted_ballast\"] for c in case_list])\n", + " mean_elec_power = np.array([c[\"electrical_power\"] for c in case_list])\n", + " flap1_util = np.array([c[\"flap1_peak_rotation_utilization\"] for c in case_list])\n", + " flap2_util = np.array([c[\"flap2_peak_rotation_utilization\"] for c in case_list])\n", + " stiffness_flap1 = np.array([c[\"opt_stiffness_flap1\"] for c in case_list])\n", + " stiffness_flap2 = np.array([c[\"opt_stiffness_flap2\"] for c in case_list])\n", + "\n", + " baseline_power = mean_elec_power[np.where(np.isclose(ballast, 0))[0][0]]\n", + " power_improvement = 100 * (mean_elec_power - baseline_power) / baseline_power\n", + "\n", + " axs[0, j].plot(ballast, mean_elec_power, marker=marker, color=color, alpha=0.85)\n", + " axs[1, j].plot(ballast, power_improvement, marker=marker, color=color, alpha=0.85)\n", + "\n", + " axs[2, j].plot(ballast, flap1_util, marker=marker, linestyle=\"-\", color=color, alpha=0.85)\n", + " axs[2, j].plot(ballast, flap2_util, marker=marker, linestyle=\"--\", linewidth=2, color=color, markerfacecolor=\"none\", alpha=0.85,)\n", + "\n", + " axs[3, j].plot(ballast, stiffness_flap1, marker=marker, linestyle=\"-\", color=color, alpha=0.85)\n", + " axs[3, j].plot(ballast,stiffness_flap2, marker=marker, linestyle=\"--\", linewidth=2, color=color, markerfacecolor=\"none\", alpha=0.85,)\n", + "\n", + " axs[0, j].set_title(f\"T = {period:g} s\")\n", + " axs[3, j].set_xlabel(\"Top-mounted ballast [kg]\")\n", + "\n", + "for ax in axs.flatten():\n", + " ax.grid(True, alpha=0.3)\n", + "\n", + "axs[0, 0].set_ylabel(\"Electrical power [W]\")\n", + "axs[1, 0].set_ylabel(\"Power improvement [%]\")\n", + "axs[2, 0].set_ylabel(\"Peak rotation utilization [-]\")\n", + "axs[3, 0].set_ylabel(\"Optimal stiffness [N·m/rad]\")\n", + "\n", + "amplitude_handles = [\n", + " Line2D([0], [0], color=color, marker=marker, linestyle=\"-\", label=f\"A = {amp:.2f} m\")\n", + " for amp, color, marker in zip(wave_amplitudes, colors, markers)\n", + "]\n", + "\n", + "axs[0, -1].legend(handles=amplitude_handles, fontsize=8, framealpha=0.9, loc=\"upper right\")\n", + "axs[1, -1].legend(handles=amplitude_handles, fontsize=8, framealpha=0.9, loc=\"upper left\")\n", + "\n", + "flap_handles = [\n", + " Line2D([0], [0], color=\"k\", linestyle=\"-\", label=\"Flap 1\"),\n", + " Line2D([0], [0], color=\"k\", linestyle=\"--\", label=\"Flap 2\"),\n", + "]\n", + "\n", + "amp_legend = axs[2, -1].legend(handles=amplitude_handles, fontsize=8, framealpha=0.9, loc=\"lower left\")\n", + "axs[2, -1].add_artist(amp_legend)\n", + "axs[2, -1].legend(handles=flap_handles, fontsize=8, framealpha=0.9, loc=\"lower right\")\n", + "\n", + "amp_legend = axs[3, -1].legend(handles=amplitude_handles, fontsize=8, framealpha=0.9, loc=\"upper left\")\n", + "axs[3, -1].add_artist(amp_legend)\n", + "axs[3, -1].legend(handles=flap_handles, fontsize=8, framealpha=0.9, loc=\"upper right\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "wot_dev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.8" + }, + "vscode": { + "interpreter": { + "hash": "a3e13d9eb6391ec8c830b5b864d7e2cac256aef786c5e95ba02dc5258710976f" + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/pyproject.toml b/pyproject.toml index be12f5832..ce0674efb 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -27,6 +27,7 @@ dependencies = [ "joblib", "wavespectra>=4.0", "netcdf4", + "multibody4everybody @ git+https://github.com/Project-SEA-Stack/Python_Multibody_for_Everybody.git@v2.0.0" ] [project.optional-dependencies] diff --git a/tests/test_utilities.py b/tests/test_utilities.py index 4b7847e11..e6d1cfe06 100644 --- a/tests/test_utilities.py +++ b/tests/test_utilities.py @@ -9,12 +9,7 @@ from pytest import approx import capytaine as cpy - - -# test function in the utilities.py - - - +# test functions in utilities.py @pytest.fixture(scope="module") def power_flows(): """Dictionary of power flows.""" @@ -34,7 +29,6 @@ def f1(): """Fundamental frequency [Hz].""" return 0.1 - @pytest.fixture(scope="module") def nfreq(): """Number of frequencies in frequency vector.""" @@ -125,7 +119,6 @@ def fb(): fb.add_translation_dof(name="Heave") return fb - @pytest.fixture(scope="module") def wb_bem(f1, nfreq, fb): """Boundary elemement model (Capytaine) results""" @@ -140,7 +133,16 @@ def wb_hydro_impedance(wb_bem): Zi = wot.hydrodynamic_impedance(hd) return Zi - +@pytest.fixture(scope="module") +def R0(): + """Transformation matrix from reduced to full coordinates.""" + return np.array( + [ + [1.0, 0.0], + [0.0, 1.0], + [1.0, 1.0], + ] + ) def test_plot_hydrodynamic_coefficients(bem_data,ndof): @@ -166,7 +168,6 @@ def test_plot_bode_impedance(intrinsic_impedance, ndof): 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) @@ -232,3 +233,98 @@ def test_linear_solve(wb_bem, regular_wave): power, _, _, _ = wot.utilities.linear_solve(wb_bem, pto_impedance, regular_wave.isel(realization=0), np.eye(1)) assert power == approx(-29.2, abs=0.05) + +def test_create_dataarray(): + omega = np.array([1.0, 2.0, 3.0]) + directions = np.array([0.0]) + dof_names = ["DOF_1"] + + impedance = np.array([1 + 1j, 2 + 2j, 3 + 3j]) + exc_coeff = np.array([4 + 1j, 5 + 2j, 6 + 3j]) + + exc_da, zi_da = wot.utilities.create_dataarray( + impedance, exc_coeff, omega, directions, dof_names + ) + + assert exc_da.dims == ("omega", "wave_direction", "influenced_dof") + assert zi_da.dims == ("omega", "radiating_dof", "influenced_dof") + + assert exc_da.shape == (3, 1, 1) + assert zi_da.shape == (3, 1, 1) + + assert exc_da.omega.values == approx(omega) + assert zi_da.omega.values == approx(omega) + assert np.squeeze(exc_da.values) == approx(exc_coeff) + assert np.squeeze(zi_da.values) == approx(impedance) + +def test_reduce_PTO_kinematics_M4E(R0): + kinematics_mat = np.array([[1.0, 2.0, 3.0]]) + + result = wot.utilities.reduce_PTO_kinematics_M4E(kinematics_mat, R0) + expected = kinematics_mat @ R0 + + assert result == approx(expected) + +def test_reduce_damping_stiffness_M4E(R0): + damping_stiffness_mat = np.diag([1.0, 2.0, 3.0]) + + result = wot.utilities.reduce_damping_stiffness_M4E(damping_stiffness_mat, R0) + expected = R0.T @ damping_stiffness_mat @ R0 + + assert result == approx(expected) + +def test_expand_and_reduce_timeseries_M4E(R0): + + timeseries_reduced = np.array( + [ + [1.0, 2.0], + [3.0, 4.0], + ] + ) + + timeseries_full = wot.utilities.expand_timeseries_M4E(timeseries_reduced, R0) + timeseries_reduced_back = wot.utilities.reduce_timeseries_M4E(timeseries_full, R0) + + assert timeseries_reduced_back == approx(timeseries_reduced) + +def test_post_process_M4E(R0): + + coords = { + "realization": [0], + "time": [0.0, 1.0], + "influenced_dof": ["red_0", "red_1"], + "type": ["test"], + } + + pos = xr.DataArray( + np.array([[[1.0, 2.0], [3.0, 4.0]]]), + dims=("realization", "influenced_dof", "time"), + coords={k: coords[k] for k in ["realization", "influenced_dof", "time"]}, + name="pos", + ) + + force = xr.DataArray( + np.ones((1, 2, 1, 2)), + dims=("realization", "influenced_dof", "type", "time"), + coords=coords, + name="force", + ) + + ds = xr.Dataset({"pos": pos, "force": force}, coords=coords) + + class DummyWEC: + def post_process(self, *args, **kwargs): + return ds, ds + + wec_fdom_full, wec_tdom_full = wot.utilities.post_process_M4E( + DummyWEC(), None, None, 1, R0 + ) + + assert wec_tdom_full.sizes["influenced_dof"] == 3 + assert wec_tdom_full.sizes["reduced_influenced_dof"] == 2 + assert "force_full" in wec_tdom_full + + expected_pos_t0 = np.array([1.0, 3.0, 4.0]) + assert wec_tdom_full["pos"].isel(realization=0, time=0).values == approx(expected_pos_t0) + + assert wec_tdom_full["force"].dims == ("realization","reduced_influenced_dof","type","time") \ No newline at end of file diff --git a/wecopttool/utilities.py b/wecopttool/utilities.py index 152cb216e..6b49fd2f1 100644 --- a/wecopttool/utilities.py +++ b/wecopttool/utilities.py @@ -22,7 +22,7 @@ import numpy as np from numpy.linalg import inv from numpy.typing import ArrayLike -from xarray import DataArray, concat +from xarray import Dataset, DataArray, concat import matplotlib.pyplot as plt from matplotlib.figure import Figure from matplotlib.axes import Axes @@ -32,6 +32,11 @@ from wecopttool.core import hydrodynamic_impedance, frequency_parameters from wecopttool.core import fd_to_td, time +# add function from Multibody for Everbody +from multibody import MbdSystem +from multibody.linearization.linearization_main import LinearizationManager +from multibody.linearization.hydro_linear_mckf import HydroLinearMCKF + # logger _log = logging.getLogger(__name__) @@ -106,9 +111,9 @@ def plot_hydrodynamic_coefficients(bem_data, 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(r'$\omega$', fontsize=10) - ax_rd[i, j].set_xlabel(r'$\omega$', fontsize=10) - ax_ex[j, 0].set_xlabel(r'$\omega$', fontsize=10) + 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('') @@ -175,12 +180,12 @@ def plot_bode_impedance(impedance: DataArray, 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(r'Frequency [Hz]', fontsize=10) + 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(r'Phase. [deg]', fontsize=10) + axes[2*i+1, j].set_ylabel(f'Phase. [deg]', fontsize=10) else: axes[i, j].set_ylabel('') if i == 0: @@ -252,9 +257,9 @@ def calculate_power_flows(wec, 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), dtype=np.float64),#eq 6.68 - 'Radiated': -1*np.sum(np.real(P_r), dtype=np.float64), - 'Actual Excitation': -1*np.sum(np.real(P_e), dtype=np.float64), + '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, } @@ -274,22 +279,19 @@ def calculate_power_flows(wec, return power_flows -def plot_power_flow( - power_flows: dict[str, float], +def plot_power_flow(power_flows: dict[str, float], tolerance: Optional[float] = None, -) -> tuple[Figure, Axes]: +)-> tuple(Figure, Axes): """Plot power flow through a WEC as Sankey diagram. Parameters ---------- power_flows - Power flow dictionary produced by + Power flow dictionary produced by for example by :py:func:`wecopttool.utilities.calculate_power_flows`. - - Required keys are ``'Optimal Excitation'``, ``'Radiated'``, - ``'Actual Excitation'``, ``'Electrical (solver)'``, - ``'Mechanical (solver)'``, ``'Absorbed'``, - ``'Unused Potential'``, and ``'PTO Loss'``. + Required keys: 'Optimal Excitation', 'Radiated', 'Actual Excitation', + 'Electrical (solver)', 'Mechanical (solver)', + 'Absorbed', 'Unused Potential', 'PTO Loss' tolerance Tolerance value for sankey diagram. """ @@ -512,3 +514,177 @@ def create_dataarray( Zi = DataArray(Zi, dims=dims_imp, coords=coords_imp, attrs=attrs_imp, name='Intrinsic impedance') return exc_coeff, Zi + +def setup_from_M4E(M4E_input_file, bem_data, fb_list, mass_list, inertia_list): + """Setup WOT from M4E""" + + body_inputs = { + i + 1: { + 'mesh': fb.mesh, + 'name': fb.name, + 'mesh_reference': 'absolute', + 'inertia_diag': [m, m, m, J, J, J], + } + for i, (fb, m, J) in enumerate(zip(fb_list, mass_list, inertia_list)) + } + + # Define the multibody system + MBDsys = MbdSystem.from_example(M4E_input_file) + mainNumVars = M4E_input_file.ic.copy() + + # Linearization + q0 = (MBDsys.ic - M4E_input_file.ic)[:len(MBDsys.Q)] + + # Instantiate the linearized MBD system + LinManager = LinearizationManager(MBDsys, q0, mainNumVars, mass_list, inertia_list, print_sym_matrices=False) + + # Instanciate and register hydrodynamics adapter + hydroAdapter = HydroLinearMCKF( + MBDsys, + (mainNumVars, mass_list, inertia_list), + body_inputs = body_inputs, + data = bem_data, + equilibrium_pos=q0 + ) + + LinManager.register(hydroAdapter) + + # Instantiate the total linearized system with attributes M,C,K,F + omega = bem_data.omega.values + FD_system = LinManager.assemble_frequency_domain(omega) + + # Extract linearized matrices + M = FD_system.M + C = FD_system.C + K = FD_system.K + F = FD_system.Fhat + + # Impedance matrix in joint coordinates + impedance_reduced = M * (1j * omega[:, None, None]) + C - 1j * K[None, :, :] * (1/omega)[:, None, None] + hydrostatic_stiffness_reduced = K + + # Convert excitation force to xarray + excitation_reduced = DataArray(F.transpose(0,2,1), coords=hydroAdapter.coords_exc, attrs={}) + + # transformation variables + R0 = LinManager.R0 + + return impedance_reduced, excitation_reduced, hydrostatic_stiffness_reduced, R0, MBDsys.NDOF + +def reduce_PTO_kinematics_M4E(kinematics_mat, R0): + """Convert PTO kinematics from full dofs to reduced dofs using the transformation matrix R0 from M4E linearization""" + + kinematics_reduced = kinematics_mat @ R0 + return kinematics_reduced + +def reduce_damping_stiffness_M4E(damping_stiffness_mat, R0): + """Convert damping or stiffness matrix from full dofs to reduced dofs using the transformation matrix R0 from M4E linearization""" + + damping_stiffness_reduced = R0.T @ damping_stiffness_mat @ R0 + return damping_stiffness_reduced + +def reduce_timeseries_M4E(timeseries_full, R0): + """Convert full timeseries to reduced DOFs using R0""" + + timeseries_reduced = np.linalg.solve(R0.T @ R0, R0.T @ timeseries_full.T).T + return timeseries_reduced + +def expand_timeseries_M4E(timeseries_reduced, R0): + """Convert reduced-coordinate time series to full DOFs using R0.""" + + timeseries_full = timeseries_reduced @ R0.T + return timeseries_full + +def _full_dof_labels(n_full): + """Create generic full-DOF labels.""" + return [f"DOF_{i}" for i in range(n_full)] + +def _expand_reduced_da(da, R0, dof_dim="influenced_dof", new_labels=None): + + n_full, n_red = R0.shape + if da.sizes[dof_dim] != n_red: + raise ValueError( + f"DataArray '{da.name}' has {da.sizes[dof_dim]} reduced DOFs, " + f"but R0 has shape {R0.shape}." + ) + + # Move reduced dof to last axis for matrix multiply + other_dims = [d for d in da.dims if d != dof_dim] + da_last = da.transpose(*other_dims, dof_dim) + + # (..., n_red) @ (n_red, n_full) -> (..., n_full) + data_full = expand_timeseries_M4E(da_last.data, R0) + + coords = {d: da_last.coords[d] for d in other_dims} + coords[dof_dim] = new_labels if new_labels is not None else _full_dof_labels(n_full) + + da_full = DataArray( + data=data_full, + dims=other_dims + [dof_dim], + coords=coords, + attrs=da.attrs, + name=da.name, + ) + + # Restore original dim ordering + return da_full.transpose(*da.dims) + +def _rename_reduced_dof_dim(da, old_dim="influenced_dof", new_dim="reduced_influenced_dof"): + if old_dim not in da.dims: + return da + return da.rename({old_dim: new_dim}) + +# Post-process and convert reduced coordinate results back into full coordinates +def post_process_M4E(wec, results, waves, nsubsteps, R0): + + wec_fdom, wec_tdom = wec.post_process(wec, results, waves, nsubsteps=nsubsteps) + + n_full, _ = R0.shape + full_labels = _full_dof_labels(n_full) + + def transform_dataset(ds): + data_vars = {} + + for name, da in ds.data_vars.items(): + if name in ["pos", "vel", "acc"] and "influenced_dof" in da.dims: + data_vars[name] = _expand_reduced_da( + da, R0, dof_dim="influenced_dof", new_labels=full_labels + ) + + elif name == "force" and "influenced_dof" in da.dims: + force_reduced = _rename_reduced_dof_dim( + da, old_dim="influenced_dof", new_dim="reduced_influenced_dof" + ) + data_vars["force"] = force_reduced + + force_full = _expand_reduced_da( + da, R0, dof_dim="influenced_dof", new_labels=full_labels + ) + force_full.name = "force_full" + force_full.attrs = dict(force_full.attrs) + force_full.attrs["note"] = ( + "Reconstructed from reduced generalized forces using R0. " + "This is not guaranteed to be the unique physical full-DOF force." + ) + data_vars["force_full"] = force_full + + else: + data_vars[name] = da + + # Only keep original coords that do NOT use influenced_dof + coords = {} + for cname, c in ds.coords.items(): + if "influenced_dof" not in c.dims: + coords[cname] = c + + # Add both coordinate sets explicitly + coords["influenced_dof"] = full_labels + if "influenced_dof" in ds.coords: + coords["reduced_influenced_dof"] = ds.coords["influenced_dof"].values + + return Dataset(data_vars=data_vars, coords=coords, attrs=ds.attrs) + + wec_fdom_full = transform_dataset(wec_fdom) + wec_tdom_full = transform_dataset(wec_tdom) + + return wec_fdom_full, wec_tdom_full \ No newline at end of file