Compute surface elevation using IFFT (#250)

* Prevent nan values for zero-frequency component in spectrum

Previously providing a zero frequency when creating a spectrum (JONSWAP or Pierson Moskowitz) produced a NaN value. This is due the f^-5 and f^-4 terms in the calculation.

The zero frequency is a valid input and is important when calculating the surface elevation from the spectrum. We know, however, that the zero frequency should __always__ have zero amplitude as the surface elevation has a mean of zero.

This change ensures that if a zero frequency is provided, the amplitude is set to zero.

* Compute surface elevation using IFFT

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 #229

* add docstrings to surface_elevation

* update surface_elevation calls in examples

---------

Co-authored-by: akeeste <akeeste@sandia.gov>

Author: mbruggs

.gitignore

@@ -32,3 +32,6 @@ $ cat .gitignore
 # Exemptions
 !**/examples/data/wave/*.mat
 !**/tests/data/wave/*.mat
+
+# Files created during tests
+mhkit/tests/wave/plots/

examples/extreme_response_contour_example.ipynb

@@ -32,12 +32,11 @@ class TestResourceSpectrum(unittest.TestCase):
 
     @classmethod
     def setUpClass(self):
-        omega = np.arange(0.1,3.5,0.01)
-        self.f = omega/(2*np.pi)
+        Trep = 600
+        df = 1 / Trep
+        self.f = np.arange(0, 1, df)
         self.Hs = 2.5
         self.Tp = 8
-        df = self.f[1] - self.f[0]
-        Trep = 1/df
         self.t = np.arange(0, Trep, 0.05)
 
     @classmethod
@@ -55,6 +54,17 @@ def test_pierson_moskowitz_spectrum(self):
         self.assertLess(errorHm0, 0.01)
         self.assertLess(errorTp0, 0.01)
 
+    def test_pierson_moskowitz_spectrum_zero_freq(self):
+        df = 0.1
+        f_zero = np.arange(0, 1, df)
+        f_nonzero = np.arange(df, 1, df)
+
+        S_zero = wave.resource.pierson_moskowitz_spectrum(f_zero, self.Tp, self.Hs)
+        S_nonzero = wave.resource.pierson_moskowitz_spectrum(f_nonzero, self.Tp, self.Hs)
+
+        self.assertEqual(S_zero.values.squeeze()[0], 0.0)
+        self.assertGreater(S_nonzero.values.squeeze()[0], 0.0)
+
     def test_jonswap_spectrum(self):
         S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
         Hm0 = wave.resource.significant_wave_height(S).iloc[0,0]
@@ -66,6 +76,17 @@ def test_jonswap_spectrum(self):
         self.assertLess(errorHm0, 0.01)
         self.assertLess(errorTp0, 0.01)
 
+    def test_jonswap_spectrum_zero_freq(self):
+        df = 0.1
+        f_zero = np.arange(0, 1, df)
+        f_nonzero = np.arange(df, 1, df)
+
+        S_zero = wave.resource.jonswap_spectrum(f_zero, self.Tp, self.Hs)
+        S_nonzero = wave.resource.jonswap_spectrum(f_nonzero, self.Tp, self.Hs)
+
+        self.assertEqual(S_zero.values.squeeze()[0], 0.0)
+        self.assertGreater(S_nonzero.values.squeeze()[0], 0.0)
+
     def test_surface_elevation_phases_np_and_pd(self):
         S0 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs)
         S1 = wave.resource.jonswap_spectrum(self.f,self.Tp,self.Hs*1.1)
@@ -129,6 +150,14 @@ def test_surface_elevation_rmse(self):
 
         self.assertLess(rmse_sum, 0.02)
 
+    def test_ifft_sum_of_sines(self):
+        S = wave.resource.jonswap_spectrum(self.f, self.Tp, self.Hs)
+
+        eta_ifft = wave.resource.surface_elevation(S, self.t, seed=1, method='ifft')
+        eta_sos = wave.resource.surface_elevation(S, self.t, seed=1, method='sum_of_sines')
+
+        assert_allclose(eta_ifft, eta_sos)     
+
     def test_plot_spectrum(self):
         filename = abspath(join(plotdir, 'wave_plot_spectrum.png'))
         if isfile(filename):

examples/extreme_response_full_sea_state_example.ipynb

@@ -90,7 +90,17 @@ def pierson_moskowitz_spectrum(f, Tp, Hs):
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
     A_PM = B_PM*(Hs/2)**2
-    Sf  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
+
+    # Avoid a divide by zero if the 0 frequency is provided
+    # The zero frequency should always have 0 amplitude, otherwise
+    # we end up with a mean offset when computing the surface elevation.
+    Sf = np.zeros(f.size)
+    if f[0] == 0.0:
+        inds = range(1, f.size)
+    else:
+        inds = range(0, f.size)
+    
+    Sf[inds]  = A_PM*f[inds]**(-5)*np.exp(-B_PM*f[inds]**(-4))
 
     col_name = 'Pierson-Moskowitz ('+str(Tp)+'s)'
     S = pd.DataFrame(Sf, index=f, columns=[col_name])
@@ -132,7 +142,17 @@ def jonswap_spectrum(f, Tp, Hs, gamma=None):
     f.sort()
     B_PM = (5/4)*(1/Tp)**4
     A_PM = B_PM*(Hs/2)**2
-    S_f  = A_PM*f**(-5)*np.exp(-B_PM*f**(-4))
+
+    # Avoid a divide by zero if the 0 frequency is provided
+    # The zero frequency should always have 0 amplitude, otherwise
+    # we end up with a mean offset when computing the surface elevation.
+    S_f = np.zeros(f.size)
+    if f[0] == 0.0:
+        inds = range(1, f.size)
+    else:
+        inds = range(0, f.size)
+
+    S_f[inds]  = A_PM*f[inds]**(-5)*np.exp(-B_PM*f[inds]**(-4))
 
     if not gamma:
         TpsqrtHs = Tp/np.sqrt(Hs);
@@ -162,7 +182,7 @@ def jonswap_spectrum(f, Tp, Hs, gamma=None):
     return S
 
 ### Metrics
-def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None):
+def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None, method='ifft'):
     """
     Calculates wave elevation time-series from spectrum
 
@@ -175,9 +195,18 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
         for example, time = np.arange(0,100,0.01)
     seed: int (optional)
         Random seed
+    frequency_bins: numpy array or pandas DataFrame (optional)
+        Bin widths for frequency of S. Required for unevenly sized bins
     phases: numpy array or pandas DataFrame (optional)
         Explicit phases for frequency components (overrides seed)
         for example, phases = np.random.rand(len(S)) * 2 * np.pi
+    method: str (optional)
+        Method used to calculate the surface elevation. 'ifft'
+        (Inverse Fast Fourier Transform) used by default if the
+        given frequency_bins==None.
+        'sum_of_sines' explicitly sums each frequency component
+        and used by default if frequency_bins are provided.
+        The 'ifft' method is significantly faster.
 
     Returns
     ---------
@@ -194,13 +223,21 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
             "frequency_bins must be of type None, np.ndarray, or pd,DataFrame")
     assert isinstance(phases, (type(None), np.ndarray, pd.DataFrame)), (
             'phases must be of type None, np.ndarray, or pd,DataFrame')
+    assert isinstance(method, str)
 
     if frequency_bins is not None:
         assert frequency_bins.squeeze().shape == (S.squeeze().shape[0],),(
             'shape of frequency_bins must match shape of S')
     if phases is not None:
         assert phases.squeeze().shape == S.squeeze().shape,(
             'shape of phases must match shape of S')
+        
+    if method is not None:
+        assert method == 'ifft' or method == 'sum_of_sines',(
+            f"unknown method {method}, options are 'ifft' or 'sum_of_sines'")
+        
+    if method == 'ifft':
+        assert S.index.values[0] == 0, ('ifft method must have zero frequency defined')        
 
     f = pd.Series(S.index)
     f.index = f
@@ -209,10 +246,12 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
         assert np.allclose(f.diff()[1:], delta_f)
     elif isinstance(frequency_bins, np.ndarray):
         delta_f = pd.Series(frequency_bins, index=S.index)
+        method = 'sum_of_sines'
     elif isinstance(frequency_bins, pd.DataFrame):
         assert len(frequency_bins.columns) == 1, ('frequency_bins must only'
                 'contain 1 column')
         delta_f = frequency_bins.squeeze()
+        method = 'sum_of_sines'
 
     if phases is None:
         np.random.seed(seed)
@@ -231,17 +270,28 @@ def surface_elevation(S, time_index, seed=None, frequency_bins=None, phases=None
     A = A.multiply(delta_f, axis=0)
     A = np.sqrt(A)
 
-    # Product of omega and time
-    B = np.outer(time_index, omega)
-    B = B.reshape((len(time_index), len(omega)))
-    B = pd.DataFrame(B, index=time_index, columns=omega.index)
-
-    # wave elevation
-    eta = pd.DataFrame(columns=S.columns, index=time_index)
-    for mcol in eta.columns:
-        C = np.cos(B+phase[mcol])
-        C = pd.DataFrame(C, index=time_index, columns=omega.index)
-        eta[mcol] = (C*A[mcol]).sum(axis=1)
+    if method == 'ifft':
+        A_cmplx = A * (np.cos(phase) + 1j*np.sin(phase))
+
+        def func(v):
+            eta = np.fft.irfft(0.5 * v.values.squeeze() * time_index.size, time_index.size)
+            return pd.Series(data=eta, index=time_index)
+        
+        eta = A_cmplx.apply(func)
+
+    elif method == 'sum_of_sines':
+        # Product of omega and time
+        B = np.outer(time_index, omega)
+        B = B.reshape((len(time_index), len(omega)))
+        B = pd.DataFrame(B, index=time_index, columns=omega.index)
+
+        # wave elevation
+        eta = pd.DataFrame(columns=S.columns, index=time_index)
+        for mcol in eta.columns:
+            C = np.cos(B+phase[mcol])
+            C = pd.DataFrame(C, index=time_index, columns=omega.index)
+            eta[mcol] = (C*A[mcol]).sum(axis=1)
+    
     return eta