% Parametric sweep script example using parallel processing (see section 5.3 of User Manual - Part II)
% Plot efficiency of the motor versus Tooth Width and Slot Depth
cd ../
initMXPscriptingAPI
file = getpmpath('SimFiles\example_innerrotor\example_innerrotor.mxp');   % initial mxp-project file
[motorProps,settingsMagnetostatic] = openMXP(file);    % open mxp-project
% motorProps - structure with motor model parameters (geometry, winding, mesh)
% settingsMagnetostatic - magnetostatic FEA settings

%% Simulation settings

ToothWidth = [3.5 4 4.5 5 5.5 6];        % sweep through a range of Tooth Width values
SlotDepth = [15 17 19 21 23 25];         % sweep through a range of Slot Depth values

% Creating a vector with all combinations for ToothWidth and SlotDepth
k = 1;
combinationSweep = zeros(length(ToothWidth)*length(SlotDepth),2);

for i = 1: length(ToothWidth)
    for j = 1:length(SlotDepth)
        combinationSweep(k,:) = [SlotDepth(j) ToothWidth(i)];
        k = k + 1;
    end
end

targetTorque = 100;  % N*m - all designs to be compared with the same output torque 100 N*m and speed 4500 RPM
initCurrent = 200;   % Arms - initial supply current value
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Save each step solution',0);
settingsMagnetostatic=setParampm(settingsMagnetostatic,'Compute cogging torque',0);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Convergence tolerance',0.01);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Current waveform','Sinusoidal');
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Advance angle',0);
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Simulation time','One electrical period');
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Mechanical speed',4500);           % RPM
settingsMagnetostatic = setParampm(settingsMagnetostatic,'Number of points',1);   % define single rotor position simulation setup

motorProps = setParampm(motorProps,'Wire size method','Fill factor');
motorProps = setParampm(motorProps,'Coil fill factor',0.5);            % keep the same coil fill factor for all the designs
motorProps = setParampm(motorProps,'End winding inductance','Auto');   % end winding inductance is automatically calculated
motorProps = setParampm(motorProps,'Phase resistance','Auto');         % phase resistance is automatically calculated

% initialize with NaN values
Efficiency = NaN*ones(length(ToothWidth),length(SlotDepth));

%% Maximum number of processes working in parallel
maxNumberOfProcesses = feature('numcores');                    % equal to the number of physical cpu cores
% maxNumberOfProcesses = 7;                                      % or can be defined manually

%% assemble all designs first
nDesigns = length(ToothWidth)*length(SlotDepth);

% Create the cell arrays to store data for all designs
cell_Geometry_ready = cell(nDesigns,1);
cell_Windings_ready = cell(nDesigns,1);
cell_Mesh_ready = cell(nDesigns,1);
cell_Materials_ready = cell(nDesigns,1);
cell_SubdomainProperty_ready = cell(nDesigns,1);
cell_Error_ready = cell(nDesigns,1);
cell_motorProps_output_ready = cell(nDesigns,1);
cell_Weight_ready = cell(nDesigns,1);
numberOfProcesses = 0; % current number of processes
i = 1;  % design number
initAssembleMXP_par();  % init design assembling in parallel
tic
% Assembling the designs using parallel processing
% keep number of processes equal to maxNumberOfProcesses
while 1
    if numberOfProcesses<maxNumberOfProcesses
        fprintf('assemble design #%i', i);
        % set up geometry parameters
        motorProps = setParampm(motorProps,'Slot Depth',combinationSweep(i,1));
        motorProps = setParampm(motorProps,'Tooth Width',combinationSweep(i,2));
        
        % assemble designs in parallel - add one process
        runAssembleMXP_par(file, motorProps);
        numberOfProcesses = numberOfProcesses+1;
        i = i+1;
    end
    if numberOfProcesses==maxNumberOfProcesses
        % wait when one of processes finish assembling the design
        while 1
            % check if any design(s) are already assembled
            [cell_Geometry, cell_Windings, cell_Mesh, cell_Materials, cell_SubdomainProperty,...
                cell_Error, cell_motorProps_output, cell_Weight, numberOfProcesses] = ...
                getAssembleMXP_par();   % get assembled design(s)
            if numberOfProcesses<maxNumberOfProcesses
                break;
            end
            pause(1); % wait for 1 second
        end
        % put assembled design(s)
        for n = 1:length(cell_motorProps_output)
            if ~isempty(cell_motorProps_output{n})
                cell_Geometry_ready{n} = cell_Geometry{n};
                cell_Windings_ready{n} = cell_Windings{n};
                cell_Mesh_ready{n} = cell_Mesh{n};
                cell_Materials_ready{n} = cell_Materials{n};
                cell_SubdomainProperty_ready{n} = cell_SubdomainProperty{n};
                cell_Error_ready{n} = cell_Error{n};
                cell_motorProps_output_ready{n} = cell_motorProps_output{n};
                cell_Weight_ready{n} = cell_Weight{n};
            end
        end
    end
    if i>nDesigns
        break
    end
end

% get all (the rest) designs
numberOfProcesses0 = numberOfProcesses;
while numberOfProcesses>0
    % check if any design(s) are finished assembling
    [cell_Geometry, cell_Windings, cell_Mesh, cell_Materials, cell_SubdomainProperty,...
        cell_Error, cell_motorProps_output, cell_Weight, numberOfProcesses] = ...
        getAssembleMXP_par();   % get assembled design(s)
    if numberOfProcesses~=numberOfProcesses0  % if new design(s) are finished assembling (number of processes changed)
        % put assembled design(s)
        for n = 1:length(cell_motorProps_output)
            if ~isempty(cell_motorProps_output{n})
                cell_Geometry_ready{n} = cell_Geometry{n};
                cell_Windings_ready{n} = cell_Windings{n};
                cell_Mesh_ready{n} = cell_Mesh{n};
                cell_Materials_ready{n} = cell_Materials{n};
                cell_SubdomainProperty_ready{n} = cell_SubdomainProperty{n};
                cell_Error_ready{n} = cell_Error{n};
                cell_motorProps_output_ready{n} = cell_motorProps_output{n};
                cell_Weight_ready{n} = cell_Weight{n};
            end
        end
        numberOfProcesses0 = numberOfProcesses;
    end
    pause(1); % wait for 1 second
end
%% Check designs
% check if any design failed assembling and try again
for i = 1:nDesigns
    if isempty(cell_Geometry_ready{i})
        motorProps = setParampm(motorProps,'Slot Depth',combinationSweep(i,1));
        motorProps = setParampm(motorProps,'Tooth Width',combinationSweep(i,2));
        
        [Geometry, Windings, Mesh, Materials, SubdomainProperty, Error, motorProps_output, Weight] = assembleMXP(file,motorProps);
        if isempty(Geometry)
            cell_Error_ready{i} = 'Unknown error';
        else
            cell_Geometry_ready{i} = Geometry;
            cell_Windings_ready{i} = Windings;
            cell_Mesh_ready{i} = Mesh;
            cell_Materials_ready{i} = Materials;
            cell_SubdomainProperty_ready{i} = SubdomainProperty;
            cell_Error_ready{i} = Error;
            cell_motorProps_output_ready{i} = motorProps_output;
            cell_Weight_ready{i} = Weight;
        end
    end
end


%% Adjust current for targetTorque
arr_Current = initCurrent*ones(nDesigns,1);
arr_Torque = zeros(nDesigns,1);
enableRun = ones(nDesigns,1);
enableRun=logical(enableRun);
for n = 1:nDesigns
    if ~isempty(cell_Error_ready{n})
        enableRun(n) = 0;           % do not run simulation if there is an error while assembling
    end
end

% Determine current required to produce target torque for each design
cell_Results_ready = cell(nDesigns,1);
cell_Timesteppingdata_ready = cell(nDesigns,1);
cell_Private_ready = cell(nDesigns,1);
iter = 1;

while 1
    initMStimestepping_par(); % init running magnetostatic FE simulations in parallel
    numberOfProcesses = 0; % current number of processes
    i = 1;
    % keep number of processes equal to maxNumberOfProcesses
    while 1
        if numberOfProcesses<maxNumberOfProcesses
            Geometry = cell_Geometry_ready{i};
            Windings = cell_Windings_ready{i};
            Mesh = cell_Mesh_ready{i};
            Materials = cell_Materials_ready{i};
            SubdomainProperty = cell_SubdomainProperty_ready{i};
            Error = cell_Error_ready{i};
            if ~isempty(cell_Private_ready{i})
                Private = cell_Private_ready{i};
            else
                Private = [];
            end
            if  enableRun(i)
                settingsMagnetostatic=setParampm(settingsMagnetostatic,'RMS supply current',arr_Current(i));
                fprintf('run simulation #%i', i);
                % run magnetostatic FE simulation in parallel - add one process
                runMStimestepping_par(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty,[]);
                numberOfProcesses = numberOfProcesses+1;
            end
            i = i+1;
        end
        if numberOfProcesses == maxNumberOfProcesses
            % wait when one of processes finish calculation
            while 1
                % check if any results are ready
                [cell_Results, cell_Timesteppingdata, ~, cell_Private, numberOfProcesses] = getMStimestepping_par();   % get magnetostatic analysis results
                if numberOfProcesses<maxNumberOfProcesses
                    break;
                end
                pause(1); % wait for 1 second
            end
            % put results
            k = 1;
            for n = 1:length(cell_Results)
                while enableRun(k) == 0
                    k = k+1;
                end
                if ~isempty(cell_Results{n})
                    cell_Results_ready{k} = cell_Results{n};
                    cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
                    cell_Private_ready{k} = cell_Private{n};
                end
                k = k+1;
            end
        end
        if i>nDesigns
            break
        end
    end
    
    % get all (the rest) results
    numberOfProcesses0 = numberOfProcesses;
    while numberOfProcesses>0
        % check if any results are ready
        [cell_Results, cell_Timesteppingdata, ~, cell_Private, numberOfProcesses] = getMStimestepping_par();
        if numberOfProcesses ~= numberOfProcesses0  % if new results are ready (number of processes changed)
            % put results
            k = 1;
            for n = 1:length(cell_Results)
                while enableRun(k) == 0
                    k = k+1;
                end
                if ~isempty(cell_Results{n})
                    cell_Results_ready{k} = cell_Results{n};
                    cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
                    cell_Private_ready{k} = cell_Private{n};
                end
                k = k+1;
            end
            numberOfProcesses0 = numberOfProcesses;
        end
        pause(1); % wait for 1 second
    end

    % Adjust the current to reach target torque or stop calculating if
    %the targetTorque is already reached
    for n = 1:length(cell_Results_ready)
        if enableRun(n)
            if isempty(cell_Results_ready{n})
                error(['Simulation results for design #' num2str(n) ' are empty (1)']);
            end
            arr_Torque(n) = cell_Results_ready{n}.torque;
            if 100*abs((arr_Torque(n)-targetTorque)/targetTorque)<0.5          % 0.5% torque error tolerance
                enableRun(n) = 0;
            end
%             if isempty(cell_Private_output{n})
%                 cell_Private_input{n} = cell_Private{n};
%             end
        end
    end
    if any(enableRun)
        arr_Current(enableRun) = (targetTorque./arr_Torque(enableRun)).*arr_Current(enableRun);
    else
        break
    end
    iter = iter+1;
    if iter>10
        break     % interrupt if more than 10 iterations 
    end
end


%% Run multi-step magnetostatic FE simulations
nPolePairs = getParampm(motorProps,'Number of pole pairs');
Ns = getParampm(motorProps,'Number of slots');
numberOfRotorPositions = setNumberOfPoints('Low accuracy',nPolePairs,Ns);
settingsMagnetostatic = setParampm(settingsMagnetostatic, 'Number of points', numberOfRotorPositions);

enableRun = ones(nDesigns,1);
for n = 1:nDesigns
    if ~isempty(cell_Error_ready{n})
        enableRun(n) = 0;
    end
end

initMStimestepping_par(); % init running magnetostatic FE simulations in parallel
numberOfProcesses = 0; % current number of processes
i = 1;

while 1
    if numberOfProcesses<maxNumberOfProcesses
        Geometry = cell_Geometry_ready{i};
        Windings = cell_Windings_ready{i};
        Mesh = cell_Mesh_ready{i};
        Materials = cell_Materials_ready{i};
        SubdomainProperty = cell_SubdomainProperty_ready{i};
        Error = cell_Error_ready{i};
        Private = cell_Private_ready{i};
        if  enableRun(i)
            settingsMagnetostatic=setParampm(settingsMagnetostatic,'RMS supply current',arr_Current(i));
            fprintf('run simulation #%i', i);
            % run magnetostatic FE simulation in parallel - add one process
            runMStimestepping_par(Private,Geometry,Mesh,Windings,settingsMagnetostatic,Materials,SubdomainProperty);
            numberOfProcesses = numberOfProcesses+1;
        end
        i = i+1;
    end
    if numberOfProcesses == maxNumberOfProcesses
        % wait when one of processes finish calculation
        while 1
            % check if any results are ready
            [cell_Results, cell_Timesteppingdata, ~, ~, numberOfProcesses] = getMStimestepping_par();   % get magnetostatic analysis results
            if numberOfProcesses<maxNumberOfProcesses
                break;
            end
            pause(1); % wait for 1 second
        end
        % put results
        k = 1;
        for n = 1:length(cell_Results)
            while enableRun(k) == 0
                k = k+1;
            end
            if ~isempty(cell_Results{n})
                cell_Results_ready{k} = cell_Results{n};
                cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
            end
            k = k+1;
        end
    end
    if i>nDesigns
        break
    end
end

% get all (the rest) results
numberOfProcesses0 = numberOfProcesses;
while numberOfProcesses>0
    % check if any results are ready
    [cell_Results, cell_Timesteppingdata, ~, ~, numberOfProcesses] = getMStimestepping_par();
    if numberOfProcesses ~= numberOfProcesses0  % if new results are ready (number of processes changed)
        % put results
        k = 1;
        for n = 1:length(cell_Results)
            while enableRun(k) == 0
                k = k+1;
            end
            if ~isempty(cell_Results{n})
                cell_Results_ready{k} = cell_Results{n};
                cell_Timesteppingdata_ready{k} = cell_Timesteppingdata{n};
            end
            k = k+1;
        end
        numberOfProcesses0 = numberOfProcesses;
    end
    pause(1); % wait for 1 second
end
toc

%% Save and plot results
Results = cell_Results_ready;
k = 1;
for i = 1: length(ToothWidth)
    for j = 1:length(SlotDepth)
        if enableRun(k)
            if isempty(Results{k})
                error(['Simulation results for design #' num2str(n) ' are empty (2)'])
            end
            Efficiency(i,j) = Results{k}.Efficiency;
            k = k + 1;
        end
    end
end

save('ExParSw_par.mat')
figure('Name','Efficiency vs. Tooth Width and Slot Depth'); surf(ToothWidth,SlotDepth,Efficiency); xlabel('Tooth Width'); ylabel('Slot Depth'); zlabel('Efficiency');

