BUG: error creating wave if nfreq=1

[rebeccamccabe]

Describe the bug
Attribute error when creating waves with nfreq=1

To Reproduce
Steps to reproduce the behavior:

w = np.array([1])
f1 = w/(2*np.pi)
nfreq = 1
wavefreq = w/(2*np.pi)
amplitude = 1
phase = 30
wavedir = 0

waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)

Expected behavior
Runs without error.

Observed behavior
WecOptTool 2.6.0: runs without error
newest commit: attribute error

Traceback (most recent call last):

  File ~\Anaconda3\envs\wot-pooch\Lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\users\rgm222\documents\github\sea-lab\wec-decider\force_saturation_sweep.py:253
    try_different_nfreqs()

  File c:\users\rgm222\documents\github\sea-lab\wec-decider\force_saturation_sweep.py:240 in try_different_nfreqs
    X[idx] = -inner_function(zeta_u, w_u_star, f_max_Fp, m, w, F_h, amplitude=1, nfreq=nfreqs[idx], plot_on=True)

  File c:\users\rgm222\documents\github\sea-lab\wec-decider\force_saturation_sweep.py:112 in inner_function
    waves = wot.waves.regular_wave(f1, nfreq, wavefreq, amplitude, phase, wavedir)

  File C:\Users/rgm222/Documents/Github/SEA-Lab/WecOptTool\wecopttool\waves.py:191 in regular_wave
    waves = elevation_fd(f1, nfreq, direction, 1, tmp, tmp, attrs)

  File C:\Users/rgm222/Documents/Github/SEA-Lab/WecOptTool\wecopttool\waves.py:120 in elevation_fd
    assert phases.shape == (nfreq, ndirections, nrealizations)

AttributeError: 'float' object has no attribute 'shape'

System:

• OS: windows 10
• Python version: 3.11.7
• WecOptTool version: main, ahead of current release 598e875

Additional information
Seems like this was introduced when the assert statement was added in #297
The root cause seems to be that degrees_to_radians outputs a float instead of a np array when the input size is (1,1,1) as it is here.
Relevant source code in waves.py:

phases = degrees_to_radians(phases, False)
assert phases.shape == (nfreq, ndirections, nrealizations)

[rebeccamccabe]

I forgot to mention that my current workaround is just to reshape in elevation_fd, but probably the long term solution is to fix degrees_to_radians.

phases_float = degrees_to_radians(phases, False)
phases = np.reshape(phases_float,phases.shape)
assert phases.shape == (nfreq, ndirections, nrealizations)

[cmichelenstrofer]

I think the reshape should be done for the input to degrees_to_radians. It makes sense for this function to return a float if the input was a float.

We should add a test that uses a single frequency (not just for waves, for a whole problem). This used to work and then broke without us realizing it.

[jtgrasb]

I think the reshape should be done for the input to degrees_to_radians.

By this do you mean the assert statement should be done before calling degrees_to_radians?

I'll start working on a test case as well.

[rebeccamccabe]

It makes sense for this function to return a float if the input was a float.

I think the issue was, it returns a float when the input was not a float, it was an ndarray of size (1,1,1).

[jtgrasb]

After looking into it more, I'm not sure that setting nfreq = 1 will give a reasonable answer. When nfreq = 1, the derivative matrix is [[0, 0] [0, 0]], which means the velocity (and therefore power) is always equal to zero.

[cmichelenstrofer]

For optimizing power there is no optimal solution (always zero) but I guess you could have a different objective function that depends on position or acceleration instead. I think we should fix the size issue, but not do anything else since all our examples use power for the objective.

[rebeccamccabe]

Could you explain why velocity and power are always zero? If nfreq=1 corresponds to a regular wave where everything is linear and there is only one frequency present in all signals, it would seem to me that velocity and power could be well-defined signals at that frequency, the optimum would exist, and could be verified analytically since it's all linear.

[cmichelenstrofer]

With only one frequency that frequency is both the fundamental and Nyquist frequency. The Nyquist frequency is only sampled at the peaks/troughs (see image below for n=4). In WecOptTool we represent the position in x_wec. Because position and velocity are 90 degrees out of phase, the peaks/troughs of position correspond to the zeros of velocity. By the same logic, they correspond to the peaks/troughs of acceleration.

So in WecOptTool, the velocity component for the largest (Nyquist) frequency is always zero. With only one frequency, velocity is zero.

(ignore the top n=10... with 8 points this example would only have 4 frequencies... the image is using n=10 just as an example of aliasing for frequency components above the Nyquist frequency)

[rebeccamccabe]

Oh I see, I think this was also partially the source of my confusion in #307. I originally assumed the conversion from freq domain to time domain would work like this:

• take the continuous-time inverse fourier transform of the frequency harmonics, so $a + bi$ at $w$ becomes $a \cos(wt) + b \sin(wt)$
• this produces a sum of nfreqs sin/cosines in continuous time
• sample that at a bunch of points in time (which can be independent of the frequency resolution) to get the time domain signals

But it seems like instead of a discrete sampling of the continuous inverse fourier transform, the PS method uses the discrete inverse fourier transform directly (np.fft.irfft), which is not equivalent, and has the property you described above where the derivative of the highest frequency appears zero. Is that a correct understanding now?

If I want to simulate (for the purpose of validating against analytical) a linear system where I know there will be nonzero harmonic content at only one frequency f, is it correct that the freq vector must consist of f and 2f, not just f, so that there are sample points at the right places in time to pick up the derivative? And if I want to simulate a system with a constraint, where I expect the solution to have odd harmonics of $f, 3f, 5f, ...$, then I need the freq vector to include their doubles so $f, 2f, 3f, 5f, 6f, 10f, ...$?

Is there a benefit of using the ifft, which appears to require frequency doubles to get accurate derivatives, compared to the way I misunderstood it to work, which would not?

My experience in signals has been almost all continuous time, so sorry if I'm misunderstanding some nuance of discrete time.

[cmichelenstrofer]

You would still need to sample the same number of points, else the problem would be over/under determined. So the number of frequency determines the number of collocation points. The only difference would be that you would not be limited to equally spaced points in time (cannot use the DFT). I am sure there are theoretical reasons why equally-spaced frequency and time steps & DFT are the best choice for some predetermined finite number of samples... but I don't have a good explanation.

The derivative components of the highest frequency are always zero, so your last frequency needs to be larger than the largest frequency you want to use for the derivative (not just the main signal: position). I don't believe you need to include double the frequency of each harmonic you are interested in.

[rebeccamccabe]

Ok, I think that makes sense now.

My reasoning for the factor of two thing was this:

• We want to resolve the signal (sin) and its derivative (cos) for a signal with freq $f$ and period $T=1/f$.
• To resolve the signal, we need time points at its peak and trough, at $T/4$ and $3T/4$
• To resolve the derivative, we need time points at its peak and trough, at $0$ and $T/2$
• To resolve both the signal and its derivative, we need points at $0, T/4, T/2, 3T/4$ which requires a sample rate of $4f$
• The nyquist frequency is $2f$ to obtain a sample rate of $4f$
• Therefore the freq vector must contain $2f$ to accurately resolve a signal with freq $f$ and its derivative
  Does that reasoning seem valid?

[cmichelenstrofer]

I don't think that's right. You don't need points at the troughs/peaks for every frequency (see n=3 in image above). And for all frequencies (except the largest) you have more than just two points per period. The "issue" with the Nyquist frequency is that there are only two points per period and they fall at the peak/trough, so the derivative is zero at all the collocation points for that frequency. This is the true value, and is not an issue, except if you try to use only one frequency and have an objective that depends on the derivative. In general you just need your final frequency to be larger than the largest frequency you want to use for the derivative.