Commit 04db7baa authored by Jérémy Baudry's avatar Jérémy Baudry

mr proper

parent 43a95ea8
This diff is collapsed.
#!/bin/sh
#PBS -q default
#PBS -N WIM2
#PBS -j oe
#PBS -l walltime=00:15:00
#PBS -l nodes=1:ppn=1
#PBS -l nice=19
module rm gcc/4.9.2
#module load openmpi/1.8.3_new-gcc-4.9.2
cd $PBS_O_WORKDIR
./WIM > std.out
rm -rf WIM2.o*
......@@ -37,7 +37,7 @@ disp =0
nbin =100
dx =5000
Cfl =1
name_sim ='temp2'
name_sim ='test'
root = 'output/'
/
......
/share/work/bauj0001/output
\ No newline at end of file
!
! _______________________
! | |
! | WIM PARAMETERS |
! |_______________________|
!______________________________________________________________________________
! WAVES PARAMETERS:
!
! Tm -> Peak period [s]
! Hs -> Significant wave height [m]
! disp -> Allowing wave dispersion
! 0: Wave dispersion is not allowed,
! group speed is the same at all spectrum
! frequency (cg=max[cg(w)])
! 1: Wave dispersion is allowed
!------------------------------------------------------------------------------
&waves_parameters
Tm =8
Hs =1
disp =0
/
!______________________________________________________________________________
! MODEL PARAMETERS:
!
! nbin -> Number of grid bin
! dx -> Spatial resolution [m]
! Cfl -> Courant–Friedrichs–Lewy condition (0<Cfl<1)
! Only in the case where disp=0. The CFL condition
! is needed to calculate the time step.
! name_sim -> name of the output file
! root -> destination folder for the output file
!------------------------------------------------------------------------------
&model_parameter
nbin =100
dx =5000
Cfl =1
name_sim ='test'
root = 'output/'
/
!______________________________________________________________________________
! SPECTRUM PARAMETERS:
!
! init_spec -> method to build the wave spectrum
! 2: Swell
! 1: JONSWAP spectrum
! 0: Bretschneider spectrum
! nfreq -> number of frequency bin
! Tmin -> Minimum period [s]
! Tmax -> Maximum period [s]
! alpha_s -> parameter for jonswap spectrum (init_spec=1)
! beta_s -> parameter for jonswap spectrum (init_spec=1)
! gamma_s -> parameter for jonswap spectrum (init_spec=1)
! swell_T -> swell period (init_spec=2)
! swell_Hs -> swell significant height (init_spec=2)
!------------------------------------------------------------------------------
&spectrum_parameters
init_spec =0
nfreq =800
Tmin =2.5
Tmax =20
alpha_s =0.2044
beta_s =1.2500
gamma_s =3.3
swell_T =19
swell_Hs=0.09
/
!______________________________________________________________________________
! ICE PARAMETERS
!
! X_ice -> Distance of the ice edge [m]
! c_cice -> Ice concentration
! ice_thick -> method for compute the ice thickness
! 0: constant thickness
! 1: thickness is a function of distance
! from ice edge
! hice -> Ice thickness (if ice_thick=0) [m]
! hmax -> Maximum ice thickness (if ice_thick=1) [m]
! Xh -> distance where h=hmax/2 (if ice_thicl=1) [m]
!------------------------------------------------------------------------------
&ice_parameters
X_ice =50000
cice =1
ice_thick =0
hice =3
hmax =3
Xh =200000
strain_crit =3e-5
P_c =0.55
/
!________________________________________________________________________________
! FSD PARAMETERS
! FSD_sheme -> method for compute <D>
! 0: dumont et al (2011)
! 1: power law
!
! minfloe -> minimum size floe to build the floe size categories [m]
! maxfloe -> maximum size floe to build the floe size categories [m]
! nbcat -> number of floe size categories
!--------------------------------------------------------------------------------
&fsd_parameters
FSD_scheme =1
minfloe =5
maxfloe =400
nbcat =60
/
!_________________________________________________________________________________
! IDT PARAMETERS
!IDT_scheme -> compute the ice thickness distribution
0: no distribution
1: distribution (rayleigh)
!mu_IDT -> parameter for the distribution
!mincat_h -> minimum ice thickness category
!maxcat_h -> maximum ice thickness category
!nbcat_h -> number of ice thickness categories
!---------------------------------------------------------------------------------
&idt_parameters
IDT_scheme =1
mu_IDT =0.5
mincat_h =0.1
maxcat_h =10
nbcat_h =50
/
!________________________________________________________________________________
function FSTD
evalin('base','clear all')
evalin('base','close all')
[filename,filepath]= uigetfile('*.*','All files');
path=[filepath,'/',filename];
x=ncread(path,'x_axis');
t=ncread(path,'time');
om=ncread(path,'omega');
spectre=ncread(path,'Spectrum');
Dave=ncread(path,'Dave');
FSD=ncread(path,'Floe size distribution');
Fsize=ncread(path,'floe size');
thick=ncread(path,'Ice thickness');
conc=ncread(path,'Ice concentration');
Hs=ncread(path,'significant height');
IDT=ncread(path,'Ice Thickness Distribution');
hcat=ncread(path,'thickness categories');
dx=(x(2)-x(1))*1000;
assignin('base','x',x)
assignin('base','FSD',FSD)
assignin('base','IDT',IDT)
assignin('base','Fsize',Fsize)
assignin('base','hcat',hcat)
assignin('base','dx',dx)
fig=figure('Name','Floe size and Thickness distribution','position',[200 200 1200 1000],'Numbertitle','off');
set(fig,'color',[0.2 0.2 0.2])
cmap=coolcolor;
assignin('base','cmap',cmap)
plot1=subplot(2,2,3:4);
[hAx,hLine1,hLine2] =plotyy(x,thick,x,Dave);
set(hAx,'color','none')
set(hAx(1),'ycolor','b')
set(hAx(2),'ycolor','g')
set(hLine1,'linewidth',2,'color','b')
set(hLine2,'linewidth',2,'color','g')
h2=gca;
assignin('base','h2',h2)
axesPosition=get(h2,'Position');
ylabel(hAx(2),'Mean Floe size [m]','interpreter','latex');
ylabel(hAx(1),'Mean Ice thickness [m]','interpreter','latex');
xlabel('Distance [km]','interpreter','latex');
hold on
xLimit = [min(x) max(x)];
yWidth = 0.05;
xOffset = -yWidth*diff(xLimit)/axesPosition(3);
h1 = axes('Position',axesPosition+yWidth.*[-1 0 1 0],...
'Color','none','XColor','k','YColor','r',...
'XLim',xLimit+[xOffset 0],...
'XTick',[],'XTickLabel',[],'NextPlot','add');
plot(h1,x,Hs(find(Hs>0)),'r','linewidth',2);
assignin('base','h1',h1)
ylabel('Significant height [m]','interpreter','latex');
hold on
curs=datacursormode(fig);
set(curs,'enable','on');
set(curs,'UpdateFcn', @plotother)
function txt = plotother(~,event)
x=evalin('base','x');
dx=evalin('base','dx');
h1=evalin('base','h1');
IDT=evalin('base','IDT');
FSD=evalin('base','FSD');
Fsize=evalin('base','Fsize');
hcat=evalin('base','hcat');
cmap=evalin('base','cmap');
pos = get(event,'Position');
txt = {['',num2str(pos(1))]};
[~,posx]=min(abs(x-pos(1)));
if (gca==h1)
ise = evalin( 'base', 'exist(''linex'') == 1' );
if ise
evalin('base', 'set(linex,''visible'',''off'')')
end
FSD1=reshape(FSD(posx,posx,:,:),length(Fsize),length(hcat));
IDT1=IDT(posx,:);
for i=1:length(hcat)
FSD1(:,i)=FSD1(:,i)*IDT1(i);
nfloe(:,i)=((FSD1(:,i)*(dx)^2)./Fsize.^2) ;
end
nfloe=nfloe/sum(sum(nfloe));
subplot(2,2,1)
title('FSTD')
colormap(cmap)
pcolor(hcat,Fsize,nfloe);
shading interp
c=colorbar;
xlabel('Floe thickness [m]','interpreter','latex')
ylabel('Floe size','interpreter','latex')
ylabel(c,'normalized number of floes','interpreter','latex')
subplot(2,2,2)
title('FSD')
FSD2(posx,:)=sum(FSD1,2);
nfloe2=((sum(FSD1,2)*dx^2)./Fsize.^2)/sum((sum(FSD1,2)*dx^2)./Fsize.^2);
plot(Fsize,nfloe2,'linewidth',2,'color','m')
set(gca,'color','none')
xlabel('Floe size [m]','interpreter','latex')
ylabel('normalized number of floes','interpreter','latex')
grid on
get(h1,'ylim')
linex=plot(h1,[pos(1) pos(1)],get(h1,'ylim'),'linewidth',1.5,'color',[0.2 0.2 0.6],'linestyle','--');
assignin('base','linex',linex)
end
end
end
clear all;
close all;
fig=figure;
% fig.Units='inches';
% fig.Position=[6 6 9 3.2];
% fig.OuterPosition=[6 6 9 3.2];
%
% fig.PaperUnits = 'inches';
% fig.PaperPosition = [0 0 9 3.2];
% fig.PaperPositionMode = 'manual';
% fig.PaperOrientation='landscape';
% fig.PaperSize=[9 3.2];
%
for n=3
titre1=['10^{-6}';'10^{-4}';'10^{-3}';'10^{-2}'];
titre2=cellstr(titre1);
dos1=['1e6 ';'1e4 ';'test ';'1e2 '];
dos2=cellstr(dos1);
dossier=[num2str(n),'mu1';num2str(n),'mu2';num2str(n),'mu3';dos2(n);num2str(n),'mu4'];
dossier=cellstr(dossier);
diffus=['1e-2_uniforme'];
mix=[1e-6 1e-4 1e-3 1e-2];
mu=[0.5 1 1.5 2 2.5]/86400;
col=['mo';'ro';'bo';'ko'];
col=cellstr(col);
for i=1:5
dos=char(dossier(i));
number(:,n)=sqrt(mu/mix(n));
lagbiomass=load(['~/projets/correction_memo/output/',dos,'/weight_',diffus,'.dat']);
eulbiomass=load(['~/projets/correction_memo/output/',dos,'/eulerian_',diffus,'.dat']);
biolag_tot=sum(lagbiomass,2)*0.5;
bioeul_tot=sum(eulbiomass,2)*0.5;
std=sqrt((biolag_tot-bioeul_tot).^2./bioeul_tot.^2)*100;
std(1:3)=0;
maxstd(i,n)=max(std);
maxbio(i,n)=max(biolag_tot);
plot(std)
hold on
end
end
clear all;
close all;
dossier=['temp1'];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
x=ncread(path,'x_axis');
t=ncread(path,'time');
om=ncread(path,'omega');
spectre=ncread(path,'Spectrum');
Dave=ncread(path,'Dave');
Dmax=ncread(path,'Dmax');
FSD=ncread(path,'Floe size distribution');
Fsize=ncread(path,'floe size');
thick=ncread(path,'Ice thickness');
conc=ncread(path,'Ice concentration');
Hs=ncread(path,'significant height');
IDT=ncread(path,'Ice Thickness Distribution');
hcat=ncread(path,'thickness categories');
fig=figure(1);
u = fig.Units;
fig.Units = 'inches';
set(fig,'outerposition',[3 3 12 10])
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
w = waitforbuttonpress;
%figure
%filename = 'joint_distribution_anim.gif';
for j=1:length(t)
FSD1=reshape(FSD(j,j,:,:),length(Fsize),length(hcat));
IDT1=IDT(j,:);
for i=1:length(hcat)
FSD1(:,i)=FSD1(:,i)*IDT1(i);
nfloe(:,i)=((FSD1(:,i)*(5000)^2)./Fsize.^2) ;
end
nfloe=nfloe/sum(sum(nfloe));
FSD2(j,:)=sum(FSD1,2);
subplot(2,2,1)
hh=surf(hcat,Fsize,nfloe);
shading interp
%
zoom(1.5)
caxis([0 0.06])
colormap(cmap)
axis([0.1 10 0 400 0 0.1])
view(108,60)
xlabel('Ice thickness [m]')
ylabel('Floe size [m]')
zlabel('area fraction')
subplot(2,2,2)
nfloe2=((FSD2(j,:)*5000^2)./Fsize.^2')/sum((FSD2(j,:)*5000^2)./Fsize.^2');
plot((Fsize),nfloe2)
axis([0 400 0 1])
xlabel('Floe size [m]')
ylabel('normalized number of floe')
subplot(2,2,3:4)
plot(x,thick)
hold on
plot(x(j),thick(j),'ro','markerfacecolor','r')
hold off
xlabel('distance [km]')
ylabel('Average ice thickness [m]')
%pause(0.02)
w = waitforbuttonpress;
% drawnow
% frame = getframe(1);
% im = frame2im(frame);
% [A,map] = rgb2ind(im,256);
% if j == 1;
% imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
% else
% imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
% end
end
clear all;
close all;
x=ncread('simulation1.nc','x_axis');
t=ncread('simulation1.nc','time');
om=ncread('simulation1.nc','omega');
spectre=ncread('simulation1.nc','Spectrum');
Dave=ncread('simulation1.nc','Dave');
Dmax=ncread('simulation1.nc','Dmax');
FSD=ncread('simulation1.nc','Floe size distribution');
Fsize=ncread('simulation1.nc','floe size');
thick=ncread('simulation1.nc','Ice thickness');
conc=ncread('simulation1.nc','Ice concentration');
Hs=ncread('simulation1.nc','significant height');
IDT=ncread('simulation1.nc','Ice Thickness Distribution');
hcat=ncread('simulation1.nc','thickness categories');
f=om/(2*pi);
E=reshape(spectre(end,end,:),length(om),1);
Ei=reshape(spectre(30,:,:),length(x),length(om));
E1=reshape(spectre(1,1,:),length(om),1);
for i=1:length(f)
thick2D(:,i)=thick;
end
thick2D(:,1)=0;
thick2D(:,end)=0;
thick2D(end,:)=0;
fig=figure(1);
u = fig.Units;
fig.Units = 'inches';
set(fig,'outerposition',[3 3 12 7])
cm=colormap('jet');
cm1(1,:)=[1 1 1];
for i=2:10
cm1(i,:)=[1-0.08*i 1-0.06*i 1];
end
cmap(1:10,:)=cm1;
cmap(11:59,:)=cm(16:end,:);
filename = 'spectrum_attenuation_anim.gif';
% w = waitforbuttonpress;
for i=1:length(t)
figure(1)
h2=surf(f,x,thick2D*0.05);
view(-56,22)
hold on
sp=reshape(spectre(i,:,:),length(x),length(om));
h=surf(f,x,sp);
view(-56,22)
colormap(jet)
% caxis([-0.1 0.2])
shading interp
set(h2,'Facealpha',0.3,'facecolor',[0.9 0.9 0.9])
ll=light;
set(ll,'position',[0.25 -20 0.4])
lighting phong
axis([min(f) max(f) min(x) max(x) 0 max(E1)])
xlabel('Frequency [s^{-1}]')
zlabel('Energy')
hold off
ylabel('x [km]')
% pause(0.001)
%w = waitforbuttonpress;
drawnow
frame = getframe(1);
im = frame2im(frame);
[A,map] = rgb2ind(im,256);
if i == 1;
imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',0);
else
imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0);
end
end
%
%
% for i=1:size(Hs,1)
%
% hs2(i)=Hs(i,i);
% end
%
% figure
% subplot(4,4,1:4)
% [ax,p1,p2] = plotyy(x,thick,x,conc,'plot','plot');
%
% ylabel(ax(1),'Ice thickness [m]')
% ylabel(ax(2),'Ice concentration')
% hold on
% line1=line([55 55], ylim);
% line2=line([150 150], ylim);
% line3=line([250 250], ylim);
% line4=line([400 400], ylim);
% set(line1,'color','k','linestyle','--')
% set(line2,'color','k','linestyle','--')
% set(line3,'color','k','linestyle','--')
% set(line4,'color','k','linestyle','--')
%
%
% subplot(4,4,5:8)
%
% plot(x,hs2)
%
% ylabel('Hs [m]')
%
%
% subplot(4,4,9)
% fsd1=reshape(FSD(:,find(x==55),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==55)),'-')
% xlabel('Floe size [m]')
% ylabel('area fraction')
% subplot(4,4,10)
% fsd1=reshape(FSD(:,find(x==150),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==150)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,11)
% fsd1=reshape(FSD(:,find(x==250),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==250)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,12)
% fsd1=reshape(FSD(:,find(x==400),:),size(FSD,1),length(Fsize));
% fsd1=fsd1';
% plot(Fsize,fsd1(:,find(x==400)),'-')
% xlabel('Floe size [m]')
%
% subplot(4,4,13:16)
% plot(x,Dave)
% xlabel('Distance [km]')
% ylabel('<D> [m]')
% saveas(gcf,'new_model.png')
%
clear all;
close all;
fig=figure;
fig.Units='inches';
fig.Position=[7.0312 6.6875 4.0312 3.3021];
fig.OuterPosition=[7.0312 6.6875 4.0312 4.3021];
fig.PaperUnits = 'inches';
fig.PaperPosition = [0 0 4.0312 3.3021];
fig.PaperPositionMode = 'manual';
fig.PaperOrientation='landscape';
fig.PaperSize=[4.0312 3.3021];
strain=[0:1e-6:3e-4];
prob=exp(-2*3e-5^2./strain.^2);
plot(strain(find(prob<0.37)),prob(find(prob<0.37)),'k--','linewidth',1.5)
hold on
plot(strain(find(prob>0.37)),prob(find(prob>0.37)),'k-','linewidth',1.5)
a=xlabel('$E_s$');
b=ylabel('$P_{\varepsilon}=P(E_w>\varepsilon_c)$');
set(a,'interpreter','latex','fontsize',11)
set(b,'interpreter','latex','fontsize',11)
scrit=3e-5;
probcrit=0.37;
plot(xlim,[probcrit probcrit])
plot([scrit scrit],ylim)
c=text(1.5e-4,0.4,'$P_c$');
d=text(3.5e-5,0.65,'$\varepsilon_c$');
set(c,'interpreter','latex','fontsize',11)
set(d,'interpreter','latex','rotation',90,'fontsize',11)
saveas(fig,'Pe_Es','pdf')
clear all;
close all;
dossier=['fig1_1'];
path=['/home/bauj0001/projets/WIM/output/',dossier,'/',dossier,'.nc'];
Fsize=ncread(path,'floe size');
strain1=load('/home/bauj0001/projets/WIM/output/strain1.dat');
beta1=load('/home/bauj0001/projets/WIM/output/beta1.dat');
Q1=load('/home/bauj0001/projets/WIM/output/Q1.dat');
strain2=load('/home/bauj0001/projets/WIM/output/strain2.dat');
beta2=load('/home/bauj0001/projets/WIM/output/beta2.dat');
Q2=load('/home/bauj0001/projets/WIM/output/Q2.dat');
strain3=load('/home/bauj0001/projets/WIM/output/strain3.dat');
beta3=load('/home/bauj0001/projets/WIM/output/beta3.dat');
Q3=load('/home/bauj0001/projets/WIM/output/Q3.dat');
Lmin=pi*0.5*((5e6*2^3)/(3*10*(1-0.3^2)))^0.25;
Lmin=Lmin*2;
wavelength=Fsize*2;