Merge pull request #179 from hivanov-nrel/loads_bug_fix

Loads mler bug fix

Author: simmsa

examples/extreme_response_MLER_example.html

@@ -46,7 +46,7 @@
 % would like to find the wave profile that will generate a heave response of 1
 % meter for our WEC device.
 
-mler_data = mler_coefficients(RAO', js, 1);
+mler_data = mler_coefficients(RAO, js, 1);
 plot(mler_data.frequency, mler_data.conditioned_spectrum); xlabel('Frequency [Hz]'); ylabel('Conditoned wave spectrum [m^2-s]')
 plot(mler_data.frequency, mler_data.phase); xlabel('Frequency [Hz]'); ylabel('Phase [rad]')
 
@@ -82,11 +82,6 @@
 
 peakHeightDesired = Hs/2 * 1.9;
 
-% Transpose MLER output
-mler_data.conditioned_spectrum = mler_data.conditioned_spectrum';
-mler_data.phase = mler_data.phase';
-mler_data.frequency = mler_data.frequency';
-
 mler_norm = mler_wave_amp_normalize(peakHeightDesired, mler_data, sim, k.values);
 
 %%

examples/extreme_response_MLER_example.m

@@ -7,18 +7,20 @@
 %
 %     Parameters
 %     ----------
-%         RAO : array
+%         RAO : array (N x 1)
 %             Response amplitude operator [-]
 %         wave_spectrum: struct
-%             Struct with wave spectral density [m^2/Hz] and frequency [Hz]
-%         response_desired: int or float
-%             Latitude longitude pairs at which to extract data.
+%             wave_spectrum.spectrum - Spectral density [m^2/Hz] (N x 1)
+%             wave_spectrum.frequency - Frequency [Hz] (N x 1)
+%         response_desired: scalar
+%             Desired response amplitude
 %
 %      Returns
 %      -------
 %         mler : struct
-%             containing conditioned wave spectral amplitude
-%             coefficient [m^2-s], and phase [rad] indexed by frequency [Hz].
+%             mler.conditioned_spectrum - Conditioned wave spectral amplitude [m^2-s] (N x 1)
+%             mler.phase - Phase [rad] (N x 1)
+%             mler.frequency - Frequency [Hz] (N x 1)
 %
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
@@ -32,16 +34,21 @@
     error('ERROR: response_desired must be an int or double')
 end
 
+% validate all inputs have same length
+N = length(RAO);
+if length(wave_spectrum.frequency) ~= N || length(wave_spectrum.spectrum) ~= N
+    error('MHKiT:loads:mler_coefficients: RAO, frequency, and spectrum must have same length');
+end
+
 % convert from Hz to rad/s
-freq = wave_spectrum.frequency * (2*pi);
-freq_hz = wave_spectrum.frequency;
-wave_spectrum = wave_spectrum.spectrum / (2*pi);
-dw = (2*pi - 0) / (length(freq)-1);
+freq_rad = wave_spectrum.frequency * (2*pi);
+wave_spectrum_rad = wave_spectrum.spectrum / (2*pi);
+dw = (2*pi - 0) / (N-1);
 
 % response spectrum
-R.spectrum = abs(RAO).^2 .* (2*wave_spectrum);
+R.spectrum = abs(RAO).^2 .* (2*wave_spectrum_rad);
 R.type = 'response';
-R.frequency = freq;
+R.frequency = freq_rad;
 
 % spectral moment calculations
 m0 = frequency_moment(R, 0);
@@ -50,16 +57,16 @@
 wBar = m1/m0;
 
 % calculate coefficient_a from Quon2016 Eqn.8
-coeff_a_rn = abs(RAO) .* sqrt(2*wave_spectrum*dw) .* ((m2 - freq*m1) + wBar*(freq*m0 - m1)) ./ (m0*m2 - m1^2);
+coeff_a_rn = abs(RAO) .* sqrt(2*dw.*wave_spectrum_rad) .* ((m2 - freq_rad.*m1) + wBar.*(freq_rad.*m0 - m1)) ./ (m0*m2 - m1^2);
 % phase delay should be positive number
 phase = -unwrap(angle(RAO));
 % for negative values of Amp, add pi phase shift, flip sign
 phase(coeff_a_rn < 0) = phase(coeff_a_rn < 0) - pi;
 coeff_a_rn(coeff_a_rn < 0) = coeff_a_rn(coeff_a_rn < 0) * -1;
 
 % calculate conditioned spectrum [m^2-s/rad]
-S = wave_spectrum .* coeff_a_rn.^2 .* response_desired^2;
-S(isnan(S)) = 0; % replace nans with zero
+conditioned_spectrum = wave_spectrum_rad .* coeff_a_rn.^2 .* response_desired^2;
+conditioned_spectrum(isnan(conditioned_spectrum)) = 0; % replace nans with zero
 % if the response amplitude we ask for is negative, we will add
 % a pi phase shift to the phase information.  This is because
 % the sign of response_desired is lost in the squaring above.
@@ -71,8 +78,8 @@
 end
 
 % outputs
-mler.conditioned_spectrum = S;
+mler.conditioned_spectrum = conditioned_spectrum;
 mler.phase = phase;
-mler.frequency = freq_hz;
+mler.frequency = wave_spectrum.frequency;
 
 end

examples/extreme_response_MLER_example.mlx

@@ -3,25 +3,18 @@
     methods (Test)
 
         function test_mler_coefficients(testCase)
-            % `mler_coefficients` Line 54, 55
-            % phase delay should be positive number
-            % phase = -unwrap(angle(RAO));
-            % Per @simmsa, this does not perform this conversion correctly for positive RAO values
-            assumeFail(testCase, "Per @simmsa, ask @hivanov about setting RAO to negative values per the above comment");
-
             % create inputs and load validation data
             fpath = '../../examples/data/loads/mler.csv';
             validation = readtable(fpath);
             wave_freq = linspace(0,1,500);
-            js = jonswap_spectrum(wave_freq,15.1,9);
+            js = pierson_moskowitz_spectrum(wave_freq,15.1,9);
             response_desired = 1;
             RAO = validation.RAO;
-            RAO = RAO';
             % execute function
             mler = mler_coefficients(RAO, js, response_desired);
             % assertions
-            assertEqual(testCase, mler.conditioned_spectrum, validation.Res_Spec', 'RelTol',0.005)
-            assertEqual(testCase, mler.phase, validation.phase', 'RelTol',0.001)
+            assertEqual(testCase, mler.conditioned_spectrum, validation.Res_Spec, 'RelTol',0.005)
+            assertEqual(testCase, mler.phase, validation.phase, 'RelTol',0.001)
         end
 
         function test_mler_simulation(testCase)