showMUD.m~ 5.05 KB
%% ------------------ EKman Non Linear -------------

% Mud dimention
nnx= 97;
nny=6;

%% Checking simulation result
home = pwd;
path2Sim = '/share/work/duqk0001/run/2DSmag/Ro0.003';
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,:) = ((U(nx,iz+1) - U(nx,iz)) / dx) - ((W(1,iz) - W(nx,iz)) / dx) ;
end

%Psi_z = u
%Psi_x = -w
%cx  = nan(nx,nz);
%cz  = nan(nx,nz);

%for ix=1:nx
%   cz = cumsimp(U(ix,:)) 
%end

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


cd(home)

V2Write = [V; V(1,:)];
fid=fopen('Vini.bin','w','b'); fwrite(fid,V2Write,'real*4'); fclose(fid);


figure('unit','centimeters','position',[5,5,40,12])
ax(1) = subplot(1,9,1:3);
pcolor(phiT');
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


dim = [nnx nny];
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

fieldName = 'v'

if strcmp(fieldName, 'psi')
    theo = phiT;
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(2).CLim)),max(abs(ax(2).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);
u_psi=  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(1).CLim)),max(abs(ax(1).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(1).CLim)),max(abs(ax(1).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(1).CLim)),max(abs(ax(1).CLim))));
colorbar



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