Processor.m

%% PreProcessor & PostProcessor for the Ocean Engineering Toolbox (OET)

% Note: Licensing information appended to this file.

% Building and executing models in the OET require the hydrodynamic
% coefficients for all bodies. These parameters can be obtained from a
% diffraction analysis tool such as WAMIT, Nemoh, Capytaine, or Aqwa. To
% convert the different output formats into a standard convention, the
% BEMIO Code from WEC-Sim must be run on MATLAB, and the results saved to
% a Hierarchical Data Format (HDF) file with the '.h5' extension.

% As of early 2024, Modelica does not have a file import protocol for HDF
% formats. Th PreProcessor script extracts the relevant data points from
% the HDF file and saves it to a MATLAB structure that can be imported by
% the OET.

% Once simulations are built and executed in the OET, validation is
% typically performed against WEC-Sim output. For this, the wave elevation
% profile generated by the OET must be supplied to WEC-Sim. Additionally,
% the meshed STL file of all bodies must be supplied for visualization.
% Note that the OET does not require a geometry file at present. Modelica
% can export the dynamic response, forces, and wave profile to a CSV file.

% The PostProcessor script then processes the wave profile into the WEC-Sim
% format. After running WEC-Sim, the script compares the dynamic response
% and forces (in heave) between both tools.

% OET documentation <https://github.com/ajay-menon-iitkgp/OET_Sys-MoDEL>.

% Clear workspace and command line
clear
clc

%% PreProcessor Script

% Export required hydrodynamic coefficients from a BEMIO HDF file to a
% MATLAB structure. This structure is then imported by the OET.

% As of v0.3, the following parameters must be exported:
% 
%   - Total body mass (m33)
%   - Infinite frequency added mass (Ainf33)
%   - Hydrostatic stiffness (Khs)
%   - State space matrices (A, B, C, D)
%   - BEM frequency components (w)
%   - Wave excitation coefficient, real component (FexcRe)
%   - Wave excitation coefficient, imaginary component (FexcIm)
%
% The default filenames are:
%
%   BEMIO output:   'rm3.h5'
%   Structure:      'hydroCoeff.mat'
%
% The hydro file name must match the output file from BEMIO, while the
% structure name must match the import statements in the OET.
%
% The h5read() function reads a HDF file. As the BEMIO output may include
% all degrees of freedom, a matrix transformation removes the non-heave
% modes for data saved as a matrix.

filename = 'F:\...\WEC-Sim\examples\RM3\hydroData\rm3.h5';
% CHECK DIALOG SELECTOR
% CHECK APP

% Hydrostatics
rho = h5read(filename,'/simulation_parameters/rho');                          % Density of water
g = h5read(filename,'/simulation_parameters/g');                              % Acceleration due to gravity
hydroCoeff.m33 = rho * h5read(filename,'/body1/properties/disp_vol');         % Equilibrium mass
hydroCoeff.Ainf33 = rho * h5read(filename,'/body1/hydro_coeffs/added_mass/inf_freq',[3 3],[1 1]);           % Infinite-frequency added mass
hydroCoeff.Khs33 = rho * g * h5read(filename,'/body1/hydro_coeffs/linear_restoring_stiffness',[3 3],[1 1]); % Linear hydrostatic stiffness

% Radiation state-space matrices
hydroCoeff.ss_rad33.A = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/A/all');   % Time-invariant state-space state matrix
hydroCoeff.ss_rad33.B = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/B/all');   % Time-invariant state-space input matrix
hydroCoeff.ss_rad33.C = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/C/all');   % Time-invariant state-space output matrix
hydroCoeff.ss_rad33.D = h5read(filename,'/body1/hydro_coeffs/radiation_damping/state_space/D/all');   % Time-invariant state-space feed-through matrix

% Excitation coefficients
hydroCoeff.w = h5read(filename,'/simulation_parameters/w');                   % Frequency values
hydroCoeff.FexcRe2 = h5read(filename,'/body1/hydro_coeffs/excitation/re');    % Real component of wave excitation force coefficient
hydroCoeff.FexcIm2 = h5read(filename,'/body1/hydro_coeffs/excitation/im');    % Imaginary component of wave excitation force coefficient

% Matrix transformation to remove non-heave modes
hydroCoeff.ss_rad33.A = squeeze(hydroCoeff.ss_rad33.A(:,:,3,3));
hydroCoeff.ss_rad33.B = squeeze(hydroCoeff.ss_rad33.B(:,:,3,3));
hydroCoeff.ss_rad33.C = squeeze(hydroCoeff.ss_rad33.C(:,:,3,3));
hydroCoeff.ss_rad33.D = squeeze(hydroCoeff.ss_rad33.D(3,3));

hydroCoeff.FexcRe(1,:) = hydroCoeff.FexcRe2(:,:,3);
hydroCoeff.FexcIm(1,:) = hydroCoeff.FexcIm2(:,:,3);

% Strip rows and columns with all 0 values
hydroCoeff.ss_rad33.A( ~any(hydroCoeff.ss_rad33.A,2), : ) = [];  % SS matrix A - rows
hydroCoeff.ss_rad33.A( :, ~any(hydroCoeff.ss_rad33.A,1) ) = [];  % SS matrix A - columns
hydroCoeff.ss_rad33.A = hydroCoeff.ss_rad33.A.';

hydroCoeff.ss_rad33.B( ~any(hydroCoeff.ss_rad33.B,2), : ) = [];  % SS matrix B - rows
hydroCoeff.ss_rad33.B( :, ~any(hydroCoeff.ss_rad33.B,1) ) = [];  % SS matrix B - columns
hydroCoeff.ss_rad33.B = hydroCoeff.ss_rad33.B.';

hydroCoeff.ss_rad33.C( ~any(hydroCoeff.ss_rad33.C,2), : ) = [];  % SS matrix C - rows
hydroCoeff.ss_rad33.C( :, ~any(hydroCoeff.ss_rad33.C,1) ) = [];  % SS matrix C - columns
hydroCoeff.ss_rad33.C = rho.*hydroCoeff.ss_rad33.C.';

hydroCoeff.ss_rad33.D( ~any(hydroCoeff.ss_rad33.D,2), : ) = [];  % SS matrix D - rows
hydroCoeff.ss_rad33.D( :, ~any(hydroCoeff.ss_rad33.D,1) ) = [];  % SS matrix D - columns
hydroCoeff.ss_rad33.D = rho.*hydroCoeff.ss_rad33.D;

% Export the MATLAB structure.
save('hydroCoeff.mat')

%% Next steps

% After running this section, build the simulation in the OET using
% the MATLAB structure by specifying the file location in the import
% statements. After simulation, export the OET variables by referring to
% the user documentation.

% Run the PostProcessor script below section-by-section.

%% PostProcessor Script - 1

% Import wave elevation profile into WEC-Sim
M = readmatrix('exportedVariables.csv');
elevationData(:,1) = M(:,1);        % Time
elevationData(:,2) = M(:,2);        % Wave elevation

% Remove duplicate time records from OET output file
[~,uidx] = unique(elevationData(:,1),'stable');
elevationData = elevationData(uidx,:);

% Save the elevation variable to a MATLAB structure. Note, this structure
% must be saved inside the WEC-Sim simulation directory to enable import.
save('elevationData.mat','elevationData')

%% PostProcessor Script - 2

% Preparing WEC-Sim for comparison with the OET. Read the WEC-Sim
% documentation here <https://wec-sim.github.io/WEC-Sim/master/index.html>.
%
% Run WEC-Sim with the elevation profile to obtain validation dataset. The
% following aspects must be kept in mind when preparing the WEC-Sim input:
%
%   1. The wecSimInputFile wave class must be modified to import data. A
%   sample code snippet to perform this is:
%
%                   |10|   %% Wave Information
%                   |11|   waves = waveClass('elevationImport');
%                   |12|   waves.elevationFile = 'elevationData.mat';
%
%   2. When setting up the '\hydroData' directory for WEC-Sim, use the same
%   BEMIO HDF file supplied to the OET.
%
% Navigate to the local directory with the WEC-Sim simulation files.
% Only run this section once the wecSimInputFile.m file, '\hydroData'
% directory, '\geometry' directory, and Simulink model have been prepared.

wecSim

%% PostProcessor Script - 3

% Plot a comparison of various WEC-Sim and OET results. The following plots
% are generated by this section:
%
%   1.      Wave elevation profile (1,2) - COMPARISON NOT APPLICABLE
%   2.1.    Wave excitation force (1,3)
%   2.2.    Wave radiation force (1,4)
%   3.1.    Heave displacement response (1,5)
%   3.2.    Heave velocity response (1,6)
%   3.3.    Heave acceleration response (1,7)
%   4.1.    Absolute difference - excitation force (1,3)
%   4.2.    Absolute difference - radiation force (1,4)
%   4.3.    Absolute difference - heave velocity (1,6)
%
% Note, the WEC-Sim results must be present in the workspace. If not, the
% user can navigate to the simulation '\output' directory and double-click
% the workspace file to add the variables.

% Collect WEC-Sim output variables (heave mode only)
wecTime = output.bodies.time(:);
wecFexc = output.bodies.forceExcitation(:,3);
wecFrad = output.bodies.forceRadiationDamping(:,3);
wecZ = output.bodies.position(:,3);
wecV = output.bodies.velocity(:,3);
wecA = output.bodies.acceleration(:,3);

% Remove duplicate time records from OET output file
[~,uidx] = unique(M(:,1),'stable');
M = M(uidx,:);

figure
plot(wecTime,wecFrad,'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,4),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,4)))
hold off
title('Radiation force comparison');
xlabel('Time (s)');
ylabel('Wave radiation force (N)');
legend;

%% Figure 1 - Wave elevation profile (1,2)
f1 = figure("Name","Wave elevation profile");
plot(elevationData(:,1),elevationData(:,2));
title('Wave elevation profile');
xlabel('Time (s)');
ylabel('Wave elevation (m)');

% Figure 2
% Plot 2.1 - Wave excitation force (1,3)
f2 = figure("Name","Force comparisons");
subplot(2,1,1);
plot(wecTime,wecFexc,'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,3),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,3)))
hold off
title('Wave excitation force comparison');
xlabel('Time (s)');
ylabel('Wave excitation force (N)');
legend;

% Plot 2.2 - Wave radiation force (1,4)
subplot(2,1,2);
plot(wecTime,wecFrad,'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,4),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,4)))
hold off
title('Radiation force comparison');
xlabel('Time (s)');
ylabel('Wave radiation force (N)');
legend;

% Figure 3
% Plot 3.1 - Heave displacement response (1,5)
f3 = figure("Name","Dynamic response comparisons");
subplot(3,1,1);
plot(wecTime,wecZ-body.centerGravity(3),'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,5),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,5)))
hold off
title('Heave displacement comparison');
xlabel('Time (s)');
ylabel('Heave displacement (m)');
legend;

% Plot 3.2 - Heave velocity response (1,6)
subplot(3,1,2);
plot(wecTime,wecV,'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,6),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,6)))
hold off
title('Heave velocity comparison');
xlabel('Time (s)');
ylabel('Heave velocity (m/s)');
legend;

% Plot 3.3 - Heave acceleration response (1,7)
subplot(3,1,3);
plot(wecTime,wecA,'-b','DisplayName','WEC-Sim');
hold on
plot(M(:,1),M(:,7),'--xr','DisplayName','OET','MarkerIndices',1:15:length(M(:,7)))
hold off
title('Heave acceleration comparison');
xlabel('Time (s)');
ylabel('Heave acceleration (m/s ^{2})');
legend;

% Figure 4
% Plot 4.1 - Absolute error in excitation force
f4 = figure("Name","Difference measurements");
subplot(3,1,1);
plot(wecTime,abs(wecFexc - interp1(M(:,1),M(:,3),wecTime)));
title('Excitation force absolute difference');
xlabel('Time (s)');
ylabel('(N)');

% Plot 4.2 - Absolute error in radiation force
subplot(3,1,2);
plot(wecTime,abs(wecFrad - interp1(M(:,1),M(:,4),wecTime)));
title('Radiation force absolute difference');
xlabel('Time (s)');
ylabel('(N)');

% Plot 4.3 - Absolute error in heave velocity
subplot(3,1,3);
plot(wecTime,abs(wecV - interp1(M(:,1),M(:,6),wecTime)));
title('Heave velocity absolute difference');
xlabel('Time (s)');
ylabel('(m/s)');

% Save figures to the WEC-Sim 'output' folder
savefig(f1,'output\Fig1.fig');
saveas(f1,'output\Fig1.png');
savefig(f2,'output\Fig2.fig');
saveas(f2,'output\Fig2.png');
savefig(f3,'output\Fig3.fig');
saveas(f3,'output\Fig3.png');

%% Licensing information

%   Modelica Ocean Engineering Toolbox v0.3
%   Copyright (C) 2024  Ajay Menon, Ali Haider, Kush Bubbar
% 
%   This program is free software: you can redistribute it and/or modify
%   it under the terms of the GNU General Public License as published by
%   the Free Software Foundation, either version 3 of the License, or
%   (at your option) any later version.
% 
%   This program is distributed in the hope that it will be useful,
%   but WITHOUT ANY WARRANTY; without even the implied warranty of
%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%   GNU General Public License for more details.
% 
%   Your copy of the GNU General Public License can be viewed
%   here <https://www.gnu.org/licenses/>.