Phase Realizations Demo and Docs

docs/source/theory.rst

@@ -146,6 +146,16 @@ This relationship can be derived from conservation of energy in both frames, and
     f_w = K^T f_p \\
     :label: conservation_energy
 
+Phase Realizations
+^^^^^^^^^^^^^^^^^^
+Irregular waves are defined in WecOptTool as a spectrum of complex frequency-domain wave elevations. The phase of each of the elevation elements impacts the time-domain result. Thus, the standard calculation of the objective function (average power) may change across a range of phase realizations. The amount of variance in power depends on multiple factors such as the optimization problem and the frequency array. When creating an irregular wave using :py:meth:`wecopttool.waves.long_crested_wave` or :py:meth:`wecopttool.waves.irregular_wave`, :code:`nrealizations` can be used to select the number of phase realizations to be used. By default, random realizations will be used to create the selected number of wave elevation spectra. The :py:meth:`wecopttool.WEC.solve` function will automatically iterate through and solve the optimization problem for each realization, and the overall result can be found by averaging the value of the objective function across all realizations. A general recommendation is to use sufficient random phase realizations such that the total simulation time sums to 20 minutes. 
+
+.. math::
+    t_{total} = \frac{nrealizations}{f1}
+    :label: total_time
+
+The selection of the number of realizations is further detailed in :doc:`_examples/tutorial_4_Pioneer`.
+
 Troubleshooting
 ---------------
 If your simulation is not behaving as expected, consider some of the general troubleshooting steps below:

examples/tutorial_1_WaveBot.ipynb

@@ -115,12 +115,43 @@
     "_ = wb.plot_cross_section(show=True)  # specific to WaveBot"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Waves\n",
+    "In this case we will use a regular wave with a frequency of 0.3 Hz and an amplitude of 0.0625 m. \n",
+    "\n",
+    "In regular waves, the WEC system response will occur at the wave frequency and corresponding odd harmonics (3rd, 5th, etc.) of the wave frequency (due to nonlinearities such as PTO force constraints). Thus, we can set the fundamental frequency, $f_1$, (also equal to the array spacing) equal to the wave frequency. We analyzed the system response in the frequency domain and determined that 10 frequencies was sufficient to capture the relevant dynamics including the nonlinear effects. Using the lowest number of frequencies possible (while still maintaining accuracy) is very important to minimize degrees of freedom and computation time for the optimization solver.\n",
+    "\n",
+    "The wave environment must be specified as a 2-dimensional `xarray.DataArray` containing the complex amplitude (m). \n",
+    "The two coordinates are the radial frequency ``omega`` (rad/s)  and the direction ``wave_direction`` (rad). \n",
+    "The `wecopttool.waves` submodule contains functions for creating this `xarray.DataArray` for different types of wave environments. We will use the `wecopttool.waves.regular_wave` function to create the regular wave. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "amplitude = 0.0625  \n",
+    "wavefreq = 0.3\n",
+    "phase = 30\n",
+    "wavedir = 0\n",
+    "\n",
+    "f1 = wavefreq\n",
+    "nfreq = 10\n",
+    "\n",
+    "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "#### Frequency and mesh check\n",
-    "We will analyze 50 frequencies with a spacing of 0.05 Hz. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n",
+    "We already defined the fundamental frequency equal to the wave frequency and the number of frequencies to include all energetic harmonics. These frequencies will be used for the Fourier representation of both the wave and the desired PTO force in the pseudo-spectral problem. See the Theory section of the Documentation for more details on the pseudo-spectral problem formulation.\n",
     "\n",
     "The `fb.minimal_computable_wavelength` parameter checks the mesh to determine the minimum wavelength that can be reliably computed using Capytaine. This warning is ignored here because the BEM results have been validated, but can be used as a guide for mesh refinement to ensure accurate BEM results."
    ]
@@ -131,8 +162,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "f1 = 0.05\n",
-    "nfreq = 50\n",
     "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
     "\n",
     "min_computable_wavelength = fb.minimal_computable_wavelength\n",
@@ -252,32 +281,6 @@
     "Note: We might receive a warning regarding negative linear damping values. Per default, WecOptTool ensures that the BEM data does not contain non-negative damping values. If you would like to correct the BEM solution manually to a minimum damping value you can specify `min_damping`. "
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Waves\n",
-    "The wave environment must be specified as a 2-dimensional `xarray.DataArray` containing the complex amplitude (m). \n",
-    "The two coordinates are the radial frequency ``omega`` (rad/s)  and the direction ``wave_direction`` (rad). \n",
-    "The `wecopttool.waves` submodule contains functions for creating this `xarray.DataArray` for different types of wave environments. \n",
-    "\n",
-    "In this case we will use a regular wave with a frequency of 0.3 Hz and an amplitude of 0.0625 m. \n",
-    "We will use the `wecopttool.waves.regular_wave` function. "
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "amplitude = 0.0625  \n",
-    "wavefreq = 0.3\n",
-    "phase = 30\n",
-    "wavedir = 0\n",
-    "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -820,8 +823,7 @@
    "metadata": {},
    "source": [
     "### Results\n",
-    "From a quick plot of the results, we see that the power absorption (where negative power is power absorbed by the device) generally improves for smaller values of `r2`.\n",
-    "It is also clear that when the WEC is cylindrical (where `r2=0.88`), power absorption is reduced."
+    "From a quick plot of the results, we see that the power absorption (where negative power is power absorbed by the device) generally improves for larger values of `r2`."
    ]
   },
   {

examples/tutorial_2_AquaHarmonics.ipynb

@@ -27,9 +27,6 @@
     "import matplotlib.pyplot as plt\n",
     "from matplotlib import cm\n",
     "from scipy.optimize import brute\n",
-    "import logging\n",
-    "logging.basicConfig()\n",
-    "logging.getLogger().setLevel(logging.DEBUG)\n",
     "\n",
     "import wecopttool as wot"
    ]
@@ -106,14 +103,41 @@
     "fb.mass = np.atleast_2d(5e3) # [kg]"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Waves\n",
+    "A regular wave will allow us to get a good initial understanding of the optimal control trajectory.\n",
+    "Note that we need to choose an appropriate fundamental frequency, $f_1$, and number of frequencies, nfreq, to ensure that \n",
+    "the wave frequency and harmonic components are within the frequency array we use to calculate the hydrodynamic data."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "amplitude = 0.5\n",
+    "wavefreq =  0.24/2 \n",
+    "phase = 30\n",
+    "wavedir = 0\n",
+    "\n",
+    "f1 = wavefreq # [Hz]\n",
+    "nfreq = 10\n",
+    "\n",
+    "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
     "#### Hydrodynamics\n",
     "\n",
     "Next we will run the boundary element model (BEM) [Capytaine](https://github.com/capytaine/capytaine) to obtain the hydrodynamic coefficients for the hull.\n",
-    "Before running the BEM solver, we must specify a set of frequencies at which to perform the calculations.\n",
+    "Before running the BEM solver, we must create a set of frequencies at which to perform the calculations.\n",
     "For `WecOptTool`, these must be a regularly spaced frequency array."
    ]
   },
@@ -123,8 +147,6 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "f1 = 0.06 # [Hz]\n",
-    "nfreq = 20\n",
     "freq = wot.frequency(f1, nfreq, False) # False -> no zero frequency\n",
     "\n",
     "bem_data = wot.run_bem(fb, freq, rho=rho, g=g)"
@@ -489,28 +511,6 @@
     ")"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Waves\n",
-    "A regular wave will allow us to get a good initial understanding of the optimal control trajectory.\n",
-    "Note that we'll want to choose a wave frequency that is within the frequency array we used to calculate the hydrodynamic data."
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "amplitude = 0.5\n",
-    "wavefreq =  0.24/2 \n",
-    "phase = 30\n",
-    "wavedir = 0\n",
-    "waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},

examples/tutorial_3_LUPA.ipynb

@@ -343,8 +343,18 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### BEM\n",
-    "With the LUPA geometry and physical properties fully defined, we can now run Capytaine to calculate the hydrodynamic and hydrostatic coefficients of the device, as done in previous tutorials. Capytaine can handle generalized modes and will calculate the coefficients for our 4 degrees of freedom. Because we are now using irregular waves, we need significantly more frequencies to capture the entire wave spectrum. The BEM coefficients have been pre-calculated and are saved in a file. To re-run the BEM, which takes about 1 hour, simply move or delete the existing `data/bem.nc` file."
+    "### Waves\n",
+    "Oregon State University has defined two sets of wave testing conditions for the LUPA: one corresponding to the PacWave South site and a scaling factor of 25, and one for the PacWave North site and a scaling factor of 20. \n",
+    "For each site/scale they provide four wave conditions to test at the Oregon State Large Wave Flume (LWF): the maximum 90<sup>th</sup> percentile, maximum percent annual energy, maximum occurrence, and minimum 10<sup>th</sup> percentile. \n",
+    "\n",
+    "In this tutorial we will use the PacWave South conditions scaled to the LWF and will design for the maximum occurrence wave. \n",
+    "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum. \n",
+    "\n",
+    "The irregular wave is created with multiple phase realizations. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes. For this tutorial, the number of realizations has been set to 2 to reduce the total run-time. \n",
+    "\n",
+    "$t_{total} = \\frac{nrealizations}{f1} $\n",
+    "\n",
+    "Because we are now using irregular waves, we need significantly more frequencies to capture the entire wave spectrum and WEC response. "
    ]
   },
   {
@@ -353,9 +363,57 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# compute hydrodynamic coefficients\n",
+    "waves = {}\n",
+    "\n",
     "f1 = 0.02\n",
     "nfreq = 50\n",
+    "\n",
+    "# regular (for testing/setup)\n",
+    "amplitude = 0.1\n",
+    "wavefreq = 0.4\n",
+    "phase = 0\n",
+    "wavedir = 0\n",
+    "waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n",
+    "\n",
+    "nrealizations = 2\n",
+    "\n",
+    "# irregular wave cases from OSU\n",
+    "wave_cases = {\n",
+    "    'south_max_90': {'Hs': 0.21, 'Tp': 3.09}, \n",
+    "    'south_max_annual': {'Hs': 0.13, 'Tp': 2.35},\n",
+    "    'south_max_occurrence': {'Hs': 0.07, 'Tp': 1.90},\n",
+    "    'south_min_10': {'Hs': 0.04, 'Tp': 1.48},  \n",
+    "    'north_max_90': {'Hs': 0.25, 'Tp': 3.46}, \n",
+    "    'north_max_annual': {'Hs': 0.16, 'Tp': 2.63},\n",
+    "    'north_max_occurrence': {'Hs': 0.09, 'Tp': 2.13},\n",
+    "    'north_min_10': {'Hs': 0.05, 'Tp': 1.68}, \n",
+    "}\n",
+    "\n",
+    "def irregular_wave(hs, tp):\n",
+    "    fp = 1/tp\n",
+    "    spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
+    "    efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
+    "    return wot.waves.long_crested_wave(efth,nrealizations=nrealizations)\n",
+    "\n",
+    "for case, data in wave_cases.items():\n",
+    "    waves[case] = irregular_wave(data['Hs'], data['Tp'])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### BEM\n",
+    "With the LUPA geometry and physical properties fully defined, we can now run Capytaine to calculate the hydrodynamic and hydrostatic coefficients of the device, as done in previous tutorials. Capytaine can handle generalized modes and will calculate the coefficients for our 4 degrees of freedom. The BEM coefficients have been pre-calculated and are saved in a file. To re-run the BEM, which takes about 1 hour, simply move or delete the existing `data/bem.nc` file."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# compute hydrodynamic coefficients\n",
     "freq = wot.frequency(f1, nfreq, False)\n",
     "\n",
     "# read BEM data file if it exists\n",
@@ -947,65 +1005,6 @@
     ")"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "### Waves\n",
-    "Oregon State University has defined two sets of wave testing conditions for the LUPA: one corresponding to the PacWave South site and a scaling factor of 25, and one for the PacWave North site and a scaling factor of 20. \n",
-    "For each site/scale they provide four wave conditions to test at the Oregon State Large Wave Flume (LWF): the maximum 90<sup>th</sup> percentile, maximum percent annual energy, maximum occurrence, and minimum 10<sup>th</sup> percentile. \n",
-    "\n",
-    "In this tutorial we will use the PacWave South conditions scaled to the LWF and will design for the maximum occurrence wave. \n",
-    "The wave conditions are specified in terms of significant wave height and peak period. Waves are mostly fully developed at the PacWave site, so we will use a Pierson-Moskowitz wave spectrum. \n",
-    "\n",
-    "The irregular wave is created with multiple phase realizations. For this tutorial, the number of realizations has been overwritten to 2 to reduce the total run-time. The solver will be run once for each wave phase realization. Each of these phase realizations leads to a slightly different result for optimal average power. Thus, for irregular wave conditions, it is recommended to include multiple phase realizations. The number of phase realizations required is dependent on the desired accuracy of the result, but it is generally recommended to include enough realizations for the total simulation time to equal 20 minutes.\n",
-    "\n",
-    "$t_{total} = \\frac{nrealizations}{f1} $"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "waves = {}\n",
-    "\n",
-    "# regular (for testing/setup)\n",
-    "amplitude = 0.1\n",
-    "wavefreq = 0.4\n",
-    "phase = 0\n",
-    "wavedir = 0\n",
-    "waves['regular'] = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)\n",
-    "\n",
-    "# number of realizations to reach 20 minutes of total simulation time\n",
-    "minutes_needed = 20\n",
-    "nrealizations = minutes_needed*60*f1\n",
-    "print(f'Number of realizations for a 20 minute total simulation time: {nrealizations}')\n",
-    "nrealizations = 2 # overwrite nrealizations to reduce run-time\n",
-    "\n",
-    "# irregular wave cases from OSU\n",
-    "wave_cases = {\n",
-    "    'south_max_90': {'Hs': 0.21, 'Tp': 3.09}, \n",
-    "    'south_max_annual': {'Hs': 0.13, 'Tp': 2.35},\n",
-    "    'south_max_occurrence': {'Hs': 0.07, 'Tp': 1.90},\n",
-    "    'south_min_10': {'Hs': 0.04, 'Tp': 1.48},  \n",
-    "    'north_max_90': {'Hs': 0.25, 'Tp': 3.46}, \n",
-    "    'north_max_annual': {'Hs': 0.16, 'Tp': 2.63},\n",
-    "    'north_max_occurrence': {'Hs': 0.09, 'Tp': 2.13},\n",
-    "    'north_min_10': {'Hs': 0.05, 'Tp': 1.68}, \n",
-    "}\n",
-    "\n",
-    "def irregular_wave(hs, tp):\n",
-    "    fp = 1/tp\n",
-    "    spectrum = lambda f: wot.waves.pierson_moskowitz_spectrum(f, fp, hs)\n",
-    "    efth = wot.waves.omnidirectional_spectrum(f1, nfreq, spectrum, \"Pierson-Moskowitz\")\n",
-    "    return wot.waves.long_crested_wave(efth,nrealizations)\n",
-    "\n",
-    "for case, data in wave_cases.items():\n",
-    "    waves[case] = irregular_wave(data['Hs'], data['Tp'])"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},