Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve execution time for surface elevation computation #229

Closed
mbruggs opened this issue Mar 15, 2023 · 3 comments
Closed

Improve execution time for surface elevation computation #229

mbruggs opened this issue Mar 15, 2023 · 3 comments
Labels
enhancement New feature or request wave module

Comments

@mbruggs
Copy link
Contributor

mbruggs commented Mar 15, 2023

This is primarily placeholder, I'm happy to contribute this improvement when I get a moment. The code below needs generalising to properly use the Dataframe that is passed in.

Describe the bug:

The computation of the surface elevation from a spectrum uses a 'sum of sines' approach which is prohibitively slow. During some experiments, trying calculate a 3hr sea state at 20Hz resolution crashed the python process on my machine, a 1hr sea state took 11s.

Alternative approach

Alternatively, we can use an approach using an ifft. For the same calculation that previously took 11s, the ifft version took 0.007s.

The code below (loosely adapted to match current MHKit structure) gave the same results as the previous surface_elevation method using the sum of sines.

def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None)
    duration_s = time_index[-1]
    dt = time_index[1] - time_index[0]

    # see notes section in
    # https://numpy.org/doc/stable/reference/generated/numpy.fft.irfft.html
    if time.size % 2 != 0:
        time = np.append(time, [duration_s + dt])

    f = S.index
    df = f[1] - f[0]
    if phases is None:
        np.random.seed(seed)
        phases = 2 * np.pi * np.random.rand(f.size)

    # Convert from energy to amplitude spectrum
    amp = np.sqrt(2 * S.values.squeeze() * df)
    amp[0] = 0  # 0 frequency gives NaN value in spectrum

    amp_re = amp * np.cos(phases)
    amp_im = amp * np.sin(phases)
    amp_complex = amp_re + 1j * amp_im

    eta = np.fft.irfft(0.5 * amp_complex * time.size, time.size)
@ssolson ssolson added enhancement New feature or request wave module labels Mar 17, 2023
@ssolson
Copy link
Contributor

ssolson commented Mar 17, 2023

That is an incredible time improvement and we should use this assuming there are no other tradeoffs. @cmichelenstrofer could you comment on @mbruggs proposed approach?

@cmichelenstrofer
Copy link
Contributor

@ssolson The IFFT is the better approach for this. The only limitation is that the frequencies have to be evenly spaced (but I believe the current function enforces this anyways). In some applications (rare in practice tbh) you might want uneven spacing, such as the "equal energy binning". So maybe there is still an advantage to keeping a separate function that uses the sum of sines? @ryancoe thoughts?

@ryancoe
Copy link
Contributor

ryancoe commented Mar 20, 2023

I'd say it'd be relatively easy to check that frequency spacing and then select either the ifft or summation of sinusoids depending on whether the spacing is equal.

mbruggs added a commit to mbruggs/MHKiT-Python that referenced this issue Jul 13, 2023
Previously the surface elevation was only computed using the 'sum of sines' method. This has been found to be prohibitively slow when computing long duration surface elevation traces.

This commit introduces the capability to compute the surface elevation using the more computationally efficient IFFT. The IFFT routine is used by default if no frequency bins are provided and the input frequency vector is equally spaced.

For example a 1hr sea state realisation at 20Hz took 11s to compute using the 'sum of sines' and 0.007s using the IFFT.

Fixes MHKiT-Software#229
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request wave module
Projects
None yet
Development

No branches or pull requests

4 participants