% Pseudo-periodic convective heat transport in channel
%
% This model simulates convective heat transfer in a channel filled with 
% water. To reduce memory requirements, the model is solved repeatedly on a
% pseudo-periodic section of the channel. Each solution corresponds to a 
% different section, and before each solution step the temperature at the 
% outlet boundary from the previous solution is mapped to the inlet 
% boundary.

% Use CTRL+SHIFT+ENTER keys to run the script in sequence.

%%
% Load the model object pseudoperiodic_llmatlab.mph that contains the 
% initial settings.
model = mphload('pseudoperiodicity_llmatlab');

%%
% Run a loop over 6 iterations.
for i = 1:6  
    % Run the study.
    model.study('std1').run;
    %%
    % At the first iteration, modify the model settings.
    if i==1
        % Add an interpolation function node that read an external data 
        % file(pseudoperiodic_data.txt) to interpolate the temperature 
        % function of space coordinates (x- and y-axis). The function name 
        % is inletTemp.
        int1 = model.func.create('int1', 'Interpolation');
        int1.model('mod1');
        int1.set('source', 'file');
        %%
        % Define the file name. Here the file is stored in the MATLAB
        % temporary directory (tempdir).
        filename = fullfile(tempdir,'pseudoperiodic_data.txt');
        int1.set('filename', filename);
        int1.set('nargs', '2');
        model.func('int1').setIndex('funcs', 'inletTemp', 0, 0);
        %%
        % Set the temperature boundary condition with the interpolation
        % function inletTemp(y,z).
        temp1 = model.physics('ht').feature('temp1');
        temp1.set('T0',1,'inletTemp(y,z)');
        %%
        % Create a 1D plot group defined with a line graph along edge 
        % number 5. The x axis set using an expression.
        pg1 = model.result.create('pg1', 'PlotGroup1D');
        lngr1 = pg1.feature.create('lngr1', 'LineGraph');
        lngr1.selection.set(5);
        lngr1.set('xdata', 'expr');
        % Create a 3D plot group that display the temperature distribution
        % on the geometry surface. The color range is set manually between
        % the value 283.15[K] and 306.5[K].
        pg2 = model.result.create('pg2', 3);
        surf1 = pg2.feature.create('surf1', 'Surface');
        surf1.set('expr', {'T'});
        surf1.set('rangecoloractive', 'on');
        surf1.set('rangecolormin', '283.15');
        surf1.set('rangecolormax', '306.5');
    % End of the IF statement.
    end
    
    %% 
    % Set the 1D x-axis expression based on the iteration number in order
    % to represent of the position in the real channel and not the position
    % in the reduced model.
    str = sprintf('x+10*R*(%d-1)',i);
    lngr1.set('xdataexpr', str);
    % Plot the 3D surface plot in a multiple plot figure.
    figure(1)
    subplot(3,2,i)
    mphplot(model,'pg2')
    % Plot the 1D plot on a different figure.
    figure(2)
    mphplot(model,'pg1')
    hold on
    %%
    % Extract the the Temperature and y- and z-coordinates on the x=1
    % plane.
    x = 1;
    y = -1e-1:2e-2:1e-1;
    z = -1e-1:2e-2:1e-1;
    [xx,yy,zz] = meshgrid(x,y,z);
    coord = [xx(:),yy(:),zz(:)]';
    [y0,z0,T] = mphinterp(model,{'y','z','T'},'coord',coord);
    %% 
    % Open the file pseudoperiodic_data.txt for writting and enter the
    % evaluated data with the format '%%y z T'. Exclude in the file the NaN
    % value that result of evaluating the temperature outside the domain.
    fid = fopen(filename,'wt');
    fprintf(fid,'%%y z T\n');
    for j = 1:length(T)
        if ~isnan(y0(j))||~isnan(z0(j))
           fprintf(fid,'%f %f %f\n',y0(j),z0(j),T(j));
        end
    end
    fclose(fid);
    %% 
    % End of the iteration.
    disp(sprintf('End of iteration No.%d',i));
end

%% 
% The following commands show how to use exact solution at the 
% computational node instead of using interpolated data. Make sure to run
% the previous sequence before.
%%
% Initialize the temperature boundary condition.
temp1.set('T0',1,'T0');
% Change the line color.
lngr1.set('linecolor', 'red');
%% 
% Run the model in loop for 6 iterations. 
for i = 1:6
    %%
    % Run the study.
    model.study('std1').run;
    %% 
    % At the first iteration, set the temperature condition with the 
    % function inletTemp(y,z).
    if i==1
       temp1.set('T0',1,'inletTemp(y,z)');
    end
    %%
    % Extract the temperature data at the computational node on boundary 6.
    data = mpheval(model,'T','edim',2,'selection',6,'refine',2);
    T = data.d1;
    y = data.p(2,:);
    z = data.p(3,:);
    %%
    % Write the data in the file pseudoperiodic_data.txt.
    fid = fopen(filename,'w');
    fprintf(fid,'%%y z T\n');
    for j = 1:length(T)
       fprintf(fid,'%f %f %f\n',y(j),z(j),T(j));
    end
    fclose(fid);
    %%
    % Display the 1D plot in the second figure. 
    figure(2)
    str = sprintf('x+10*R*(%d-1)',i);
    lngr1.set('xdataexpr', str);
    mphplot(model,'pg1')
    %%
    % End of the iteration.
    disp(sprintf('End of iteration No.%d',i));
end