showMUD.m 6.92 KB
%% ------------------ EKman Non Linear -------------
clear all

% Mud dimention
nnx= 3585;
nny=65;

dim = [nnx nny];

Lx= 10755;
Lz= 195;
ddx = Lx/nnx;
ddz = Lx/nny;
%% Checking simulation result
home = pwd;

if strcmp(fgetMachin, 'mingan')
    path2Sim = '/share/work/duqk0001/run/2DAh/Ro0.003';
elseif strcmp(fgetMachin, 'mac')
    path2Sim = '/Users/Administrator/Documents/uqar/Leads/result/data/varRoEqAh_v2/Ro0.003';
else
   error('Wrong location in genVreal.m ') 
end

cd(path2Sim)

dynDiag = rdmds('dynDiag', Inf);
U = squeeze(dynDiag(:,1,:,1));
V = squeeze(dynDiag(:,1,:,2));  
W = squeeze(dynDiag(:,1,:,3));

% Reverse upside down U,V,W
%U = flipud(U);
%V = flipud(V);
%W = flipud(W);

[nx,nz] = size(U);

dx=3;

curlY = nan(nx,nz-1);

for ix=1:nx-1
    for iz=1:nz-1
        curlY(ix,iz) = -((U(ix,iz+1) - U(ix,iz)) / dx) - ...
                        ((W(ix+1,iz) - W(ix,iz)) / dx) ;
    end
end
for iz=1:nz-1
    curlY(nx,iz) = -((U(nx,iz+1) - U(nx,iz)) / dx) - ((W(1,iz) - W(nx,iz)) / dx) ;
end

% Stream Function
%{
psiT =  zeros(nx,nz);
psiT(:,2)= -psiT(:,1) + dx*U(:,2);
for i=2:nz
    psiT(:,i+1) = -psiT(:,i-1) + 2*dx*U(:,i);
end
%}

[phiT, psiT] =  streamFunction('showFig',0);
psiT = -psiT ;

cd(home)


figure('unit','centimeters','position',[5,5,40,12])
ax(1) = subplot(1,9,1:3);
pcolor(psiT');
ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('Psi')
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

ax(2) = subplot(1,9,4:6);
pcolor(curlY');
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('Zeta')
set(gca, 'yticklabel',[])
%ax(2).CLim = ax(1).CLim;
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

ax(3) = subplot(1,9,7:9);
pcolor(V')
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('V')
set(gca, 'yticklabel',[])
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

% Set position
print('-dpng','-r150',['result/theoryVZP.png'])


%% Ploting SolverField

psi  = readFortran('psi.bin',dim);
zeta  = readFortran('zeta.bin',dim);
v  = readFortran('v.bin',dim);

figure('unit','centimeters','position',[5,5,40,12])
ax(1) = subplot(1,9,1:3);
pcolor(psi');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('Psi')
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

ax(2) = subplot(1,9,4:6);
pcolor(zeta');
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('Zeta')
set(gca, 'yticklabel',[])
%ax(2).CLim = ax(1).CLim;
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

ax(3) = subplot(1,9,7:9);
pcolor(v')
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('V')
set(gca, 'yticklabel',[])
ax =gca;
colormap(ax,b2r(-max(abs(ax.CLim)),max(abs(ax.CLim))));
colorbar

% Set position
print('-dpng','-r150',['result/showVZP.png'])



%% Compare all Fields

psi  = readFortran('psi.bin',dim);
zeta  = readFortran('zeta.bin',dim);
v  = readFortran('v.bin',dim);

figure('unit','centimeters','position',[5,5,40,12])

%  --------- MIT ---------------
ax(1) = subplot(2,9,1:3);
pcolor(psiT');
ylabel('y')
%xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('Psi')
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar

ax(2) = subplot(2,9,4:6);
pcolor(curlY');
%ylabel('y')
%xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('Zeta')
set(gca, 'yticklabel',[])
%ax(2).CLim = ax(1).CLim;
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar

ax(3) = subplot(2,9,7:9);
pcolor(V')
%ylabel('y')
%xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('V')
set(gca, 'yticklabel',[])
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar

% --------- MUD -------------- 
ax(4) = subplot(2,9,10:12);
pcolor(psi');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('Psi')
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar

ax(5) = subplot(2,9,13:15);
pcolor(zeta');
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('Zeta')
set(gca, 'yticklabel',[])
%ax(2).CLim = ax(1).CLim;
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar

ax(6) = subplot(2,9,16:18);
pcolor(v')
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('V')
set(gca, 'yticklabel',[])
axg =gca;
colormap(axg,b2r(-max(abs(axg.CLim)),max(abs(axg.CLim))));
colorbar


% Set position
print('-dpng','-r150',['result/compareVZP.png'])


%% Compare one field

fieldName = 'psi'

if strcmp(fieldName, 'psi')
    theo = psiT;
elseif strcmp(fieldName,'v')
    theo = V;
elseif strcmp(fieldName,'zeta')
    theo = curlY;
else
    error('wrong field')
end

% Reshape Theo as field
theo = resizem(theo,dim);

% Open result 
field  = readFortran([ fieldName '.bin'],dim);

% Compare with analytique
fDiff= theo-field;


figure('unit','centimeters','position',[5,5,30,12])
ax(1) = subplot(1,9,1:3);
pcolor(theo');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('MIT')
colormap(ax(1),b2r(-max(abs(ax(1).CLim)),max(abs(ax(1).CLim))));
%colorbar

ax(2) = subplot(1,9,4:6);
pcolor(field');
%ylabel('y')
xlabel('x')
shading flat
set(gca,'ydir','reverse')
title('MUD')
set(gca, 'yticklabel',[])
%ax(2).CLim = ax(1).CLim;
colormap(ax(2),b2r(-max(abs(ax(1).CLim)),max(abs(ax(1).CLim))));
%colorbar

ax(3) = subplot(1,9,7:9);
pcolor(fDiff')
%ylabel('y')
xlabel('x')
shading flat
title('Res')
set(gca,'ydir','reverse')
set(gca, 'yticklabel',[])
colormap(ax(3),b2r(-max(abs(ax(3).CLim)),max(abs(ax(3).CLim))));
%colorbar


% Set position
set(ax(1),'position',get(ax(1),'position') - [ 0.05 0 0 0]  )
set(ax(2),'position',get(ax(2),'position') - [ 0.05 0 0 0]  )
set(ax(3),'position',get(ax(3),'position') - [ 0.0 0 0 0]  )


tmpPos = get(ax(3),'Position');

cb1 = colorbar(ax(2),'Position',[tmpPos(1)-0.06  tmpPos(2) 0.02 tmpPos(4)]);
cb2 = colorbar(ax(3),'Position',[tmpPos(1)+tmpPos(3)+0.01 tmpPos(2) 0.02 tmpPos(4)]);

print('-dpng','-r150',['result/compare_' fieldName '.png'])


%% Compare V,U optain by psi

w_psi= -diff(psi,1); %diff(psi,1);
u_psi= diff(psi,2); %diff(psi,2);  

figure('unit','centimeters','position',[5,5,40,12])
ax(1) = subplot(2,2,1);
pcolor(U');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('U - MIT')
colormap(ax(1),b2r(-max(abs(ax(1).CLim)),max(abs(ax(1).CLim))));
colorbar

ax(2) = subplot(2,2,2);
pcolor(u_psi');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('U - MUD')
colormap(ax(2),b2r(-max(abs(ax(2).CLim)),max(abs(ax(2).CLim))));
colorbar

ax(3) = subplot(2,2,3);
pcolor(W');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('W - MIT')
colormap(ax(3),b2r(-max(abs(ax(3).CLim)),max(abs(ax(3).CLim))));
colorbar

ax(4) = subplot(2,2,4);
pcolor(w_psi');
ylabel('y')
xlabel('x')
shading flat
set(gca, 'ydir','reverse')
title('W - MUD')
colormap(ax(4),b2r(-max(abs(ax(4).CLim)),max(abs(ax(4).CLim))));
colorbar



% Set position
print('-dpng','-r150',['result/compareVel.png'])