Wave elevation generation

[ckberinger]

Hello,

I am trying to make a Pierson Moskowitz wave elevation 20 minute time series with a peak period of 16s and a significant wave height of 1.5 meters. The code below is what I have but I know it has several issues.
What is the frequency input supposed to look like in the "pierson_moskowitz_spectrum" function? Should it be linearly spaced? Include a zero frequency? Right now it is based on having 32 frequencies between 2 and 30 seconds which is based on how model data reports frequencies. The "surface_elevation" code requires a 0 frequency for the ifft, do you suggest padding with a zero? Example code for this function would be helpful as well!

swell_Tp = 16; %s
swell_Hs = 1.5; %m
swell_S = pierson_moskowitz_spectrum(linspace(1/30,1/2,32),swell_Tp,swell_Hs);

t = 0:dt:20*60; % seconds
swell_eta = surface_elevation(swell_S,t,"seed",1);

Thank you,
Courtney

[simmsa]

Hi @ckberinger,

Note, this post was edited by @simmsa on 4/23/2024 to use the correct units, Spectral Density [m^2/Hz], on the y axis of the spectrum Visualizations

Thank you for reaching out! I am able to replicate this issue. I think I understand the core issue and have found a solution which is detailed below. If you have any questions or need further assistance, don't hesitate to reach out.

Issue Description

Generating a wave spectrum works with a frequency index array starting above zero, but generating a surface elevation using that spectrum fails with Python Error: AssertionError: ifft method must have zero frequency defined.

Solution Summary

Create a frequency index array that starts at zero. Define the frequency index delta based on the desired surface elevation duration.

Issue

@ckberinger, here is what I am seeing with your code:

Generating the wave spectra

swell_Tp = 16; %s
swell_Hs = 1.5; %m
swell_S = pierson_moskowitz_spectrum(linspace(1/30,1/2,32),swell_Tp,swell_Hs);
plot(swell_S.frequency, swell_S.spectrum);
xlabel("Frequency [Hz]");
ylabel("Wave Height [m]")

Yields a valid wave height spectrum:

But generating the surface elevation yields an error:

dt = 0.001;
t = 0:dt:1*60; % seconds
swell_eta = surface_elevation(swell_S,t,"seed",1)

Error using [resource>surface_elevation](matlab:matlab.internal.language.introspective.errorDocCallback('resource>surface_elevation')) Python Error: AssertionError: ifft method must have zero frequency defined

Solution

Define a frequency range that includes zero and use a frequency index derived from the desired surface elevation seconds.

1. Generate the Pierson-Moskowitz spectrum (I used one minute surface elevation duration to create a plot that is easier to inspect. At the end I create a plot that uses the full 20 minutes):

• Note: this spectrum appears to be very "peaky", this is due to the usage of the ifft method and the relatively short time scale. In the full 20 minute surface elevation generation the spectra is not as "peaky"

% Time Definition
surface_elevation_duration_minutes = 1;
surface_elevation_duration_seconds = surface_elevation_duration_minutes * 60;
surface_elevation_delta_time_seconds = 0.001;
surface_elevation_times = 0:surface_elevation_delta_time_seconds:surface_elevation_duration_seconds;

% Frequency Definition
wave_spectrum_frequency_bin_size = 1 / surface_elevation_duration_seconds;
wave_spectrum_frequencies = 0:wave_spectrum_frequency_bin_size:1;

% Resource Definition
swell_Tp = 16; %s
swell_Hs = 1.5; %m

% Spectrum Generation
swell_S = pierson_moskowitz_spectrum(wave_spectrum_frequencies,swell_Tp,swell_Hs);

% Spectrum Plot - Wave Height vs. Frequency
plot(swell_S.frequency, swell_S.spectrum);
xlabel("Frequency [Hz]");
ylabel("Wave Height [m]");
title(sprintf("Pierson Moskowitz Wave Height Spectrum [%.1fs, %.2fm]", swell_Tp, swell_Hs))

Which yields:

2. Generate a valid surface elevation:

% Surface Elevation Generation
swell_eta = surface_elevation(swell_S,surface_elevation_times ,"seed",1);

% Surface Elevation Plot - Wave Height vs. Time
plot(swell_eta.elevation);
xlabel("Time [s]");
ylabel("Wave Height [m]")
xtickformat('%,.0f');
title(sprintf("Pierson Moskowitz %ds Surface Elevation [%.1fs, %.2fm]", surface_elevation_duration_seconds, swell_Tp, swell_Hs));

Which yields the following plot and appears to match expectations for one minute Pierson Moskowitz surface elevations:

3. Generate the full 20 minute surface elevation:

% Time Definition
surface_elevation_duration_minutes = 20;
surface_elevation_duration_seconds = surface_elevation_duration_minutes * 60;
surface_elevation_delta_time_seconds = 0.001;
surface_elevation_times = 0:surface_elevation_delta_time_seconds:surface_elevation_duration_seconds;

% Frequency Definition
wave_spectrum_frequency_bin_size = 1 / surface_elevation_duration_seconds;
wave_spectrum_frequencies = 0:wave_spectrum_frequency_bin_size:0.5;

% Resource Definition
swell_Tp = 16; %s
swell_Hs = 1.5; %m

% Spectrum Generation
swell_S = pierson_moskowitz_spectrum(wave_spectrum_frequencies,swell_Tp,swell_Hs);

% Spectrum Plot - Wave Height vs. Frequency
plot(swell_S.frequency, swell_S.spectrum);
xlabel("Frequency [Hz]");
ylabel("Wave Height [m]");

% Surface Elevation Generation
swell_eta = surface_elevation(swell_S,surface_elevation_times ,"seed",1);

% Surface Elevation Plot - Wave Height vs. Time
plot(swell_eta.time, swell_eta.elevation);
xlabel("Time [s]");
ylabel("Wave Height [m]")
xtickformat('%,.0f');
title(sprintf("Pierson Moskowitz %ds Surface Elevation [%.1fs, %.2fm]", surface_elevation_duration_seconds, swell_Tp, swell_Hs));

Which yields the spectrum:

And the surface elevation:

Which should be the desired result.

One note, for completeness MHKiT-MATLAB needs to add an argument to surface_elevation that allows the user to select the sum_of_sines method for surface elevation generation. This would allow generation of a surface elevation from a spectrum without 0 as a frequency index. We are addressing this in #110 and we plan to add a PR to include this functionality as soon as possible.

Also creating some code examples is a great suggestion! Do you or your team have any suggestions they would like to see?

[ckberinger]

Thank you for working through this and showing a working version! What is the justification for the time step of 0.001 seconds? I had a similarly good-looking result for the surface elevation 20 minute time series plot when I used 1000 frequencies, but it seems not based on physics. This surface elevation plot is very sensitive to the number of frequencies when including frequencies down to 0. For example, should only 32 frequencies (the amount collected from model data/buoys) produce a reasonable time series? It does not though. When using the ifft method, should one perhaps use non-linearly spaced frequencies as (models and buoys do)?
I guess I am just wondering why the ifft method is so sensitive to the number of frequencies and how do you justify without just doing a sanity check on the time series? What if the eventual goal is to input this time series into WEC-Sim or an experiment or a different numerical model? I just want to make sure it is grounded in physics. Thanks for thinking through these questions!

I may just end up using the sum_of_sines with 32 frequencies between 1/30 and 1/2 and see if that produces a reasonable result! Thanks for the suggestion.

[simmsa]

@ckberinger, thank you for the clarification, hopefully I am not leading you in the wrong direction. My justification for using a 0.001 second timestamp is my own perspective. I typically work with time-series data that is collected at high frequency so that is why I chose to use a 1000 Hz sample rate. You can certainly drop that number down to something more reasonable.

In #126 we are adding the method option to the surface elevation function, which allows for choosing between the ifft and sum_of_sines method. The iftt method is the default. The sum_of_sines method allows for more flexibility in the frequency index. Can you try downloading the new toolbox in #126 from here (upper right corner two buttons to the right of raw)? You should be able to choose between the two methods.

Here is the latest code using the sum of sines method:

swell_Tp = 16; %s
swell_Hs = 1.5; %m
swell_S = pierson_moskowitz_spectrum(linspace(1/30,1/2,32),swell_Tp,swell_Hs);
swell_method = "sum_of_sines"

% Spectrum Plot - Wave Height vs. Frequency
plot(swell_S.frequency, swell_S.spectrum);
xlabel("Frequency [Hz]");
ylabel("Wave Height [m]");
title(sprintf("Pierson Moskowitz Wave Height Spectrum [%.1fs, %.2fm]", swell_Tp, swell_Hs))

dt = 0.1;
t = 0:dt:20*60; % seconds
swell_eta = surface_elevation(swell_S,t,"seed",1, "method", "sum_of_sines");

% Surface Elevation Plot - Wave Height vs. Time
plot(swell_eta.time, swell_eta.elevation);
xlabel("Time [s]");
ylabel("Wave Height [m]")

Spectrum:

Surface Elevation:

Does this align better with your use case?

If we need more clarification on the physics we may need to reach out to the developers who originally developed the ifft method. Here is the pull request where this feature was originally implemented. The details of the ifft method will be better explained by them.

Let us know if this is helpful!

Note: this post was edited by @simmsa on 4/23/2024 to use the correct units, Spectral Density [m^2/Hz], on the y axis of the spectrum visualizations

[simmsa]

@ckberinger, one more follow up. We should be handling your use case more gracefully. I submitted Issue 308 in the MHKiT-Python repo to address this issue further.

Please let us know if these fixes are a suitable solution for your use case.

[ckberinger]

Hello,
I downloaded and installed the new toolbox. But it did not update my surface elevation code.
So I got this error:

Error using surface_elevation (line 1)
Invalid argument name 'method'. Name must be 'seed', 'frequency_bins', or 'phases'.

Error in WEC_Workshop2024 (line 76)
swell_eta = surface_elevation(swell_S,t,"seed",1, "method", "sum_of_sines");

I tried restarting MATLAB and my computer. Should I disable the old mhkit toolbox? I tried doing that but I still got the same error.

Oregon State University, TEAMER, NREL, and Sandia are hosting a workshop this coming Monday and I would like the students to be able to use this sine method if it could be in the master before then.

[simmsa]

@ckberinger, this could be due to issues with your MATLAB path sourcing a previously installed add-on of MHKiT-MATLAB. Can you verify that the output of running path in the Command Window does not have two versions of the mhkit addon in the path. On my machine this manifests as mhkit and mhkit(2) (screenshot below) showing up in the path. MATLAB uses the first surface_elevation function in your path and ignores anything after. If you have this condition or something similar navigate to the MATLAB default Add-On Installation folder on your machine, delete any mhkit folders, and then reinstall the add-on. The MATLAB Add-Ons documentation has more info on the default Add-On installation folders.

Let us know if this is not the issue.

For merging into master, I have a review out to both @rpauly18 and @hivanov-nrel on our team. I will reach out to both of them directly to see if they have time to review this before Monday. We appreciate your patience!

[simmsa]

@ckberinger, #126 was merged today. Please let us know if this is an adequate fix of the surface elevation function for your use case.

[ckberinger]

This works great, thank you for adding!

[simmsa]

Awesome! Thanks you for bringing this to our attention and working to get this fixed. Let us know if you find any other issues or have suggestions for improvements.