source: issm/trunk-jpl/test/NightlyRun/test1102.m@ 13663

Last change on this file since 13663 was 13663, checked in by Mathieu Morlighem, 12 years ago

CHG: cosmetics on matlab file (model -> model(); Enum -> Enum())

File size: 5.9 KB
RevLine 
[5049]1%This test is a test from the ISMP-HOM Intercomparison project
2%Pattyn and Payne 2006
[5606]3printingflag=false;
[5049]4
[5055]5L_list={5000,10000,20000,40000,80000,160000};
[5049]6results={};
[5606]7minvx=[];
8maxvx=[];
[5049]9
10for i=1:length(L_list),
11 L=L_list{i};
12 nx=20; %numberof nodes in x direction
13 ny=20;
[13663]14 md=model();
[5049]15 md=squaremesh(md,L,L,nx,ny);
[9641]16 md=setmask(md,'',''); %ice sheet test
[5606]17
18% %Find elements at the corner and extract model
[9734]19% posnodes=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y)));
[9733]20% [a,b]=find(ismember(md.mesh.elements,posnodes));
[9725]21% elements=ones(md.mesh.numberofelements,1);
[5606]22% elements(a)=0;
23% md=modelextract(md,elements);
24
[5049]25 md=parameterize(md,'../Par/ISMIPA.par');
[5606]26 md=extrude(md,10,1);
[9664]27 md=setflowequation(md,'stokes','all');
[5049]28
29 %Create dirichlet on the bed only
[9729]30 pos=find(md.mesh.vertexonbed);
[9679]31 md.diagnostic.spcvx(pos)=0;
32 md.diagnostic.spcvy(pos)=0;
33 md.diagnostic.spcvz(pos)=0;
[5049]34
35 %Create MPCs to have periodic boundary conditions
[9734]36 %posx=find(md.mesh.x==0);
37 %posx2=find(md.mesh.x==max(md.mesh.x));
38 %posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
39 %posx2=find(md.mesh.x==max(md.mesh.x) & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed);
[5049]40
[9734]41 %posy=find(md.mesh.y==0 & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %Don't take the same nodes two times
42 %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed);
[5049]43
[9679]44 %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2];
[5049]45
46 %Compute the diagnostic
[9763]47 md.diagnostic.abstol=NaN;
[9679]48 md.diagnostic.reltol=NaN;
49 md.diagnostic.restol=1;
[8630]50 md.cluster=generic('name',oshostname(),'np',8);
[13663]51 md=solve(md,DiagnosticSolutionEnum());
[5049]52
53 %Plot the results and save them
[10976]54 vx=(md.results.DiagnosticSolution.Vx);
55 vy=(md.results.DiagnosticSolution.Vy);
56 vz=(md.results.DiagnosticSolution.Vz);
57 pressure=(md.results.DiagnosticSolution.Pressure);
[5049]58 results{i}=md.results.DiagnosticSolution;
[9725]59 minvx(i)=min(vx(end-md.mesh.numberofvertices2d+1:end));
60 maxvx(i)=max(vx(end-md.mesh.numberofvertices2d+1:end));
[5049]61
[5606]62 %Now plot vx, vy, vz and vx on a cross section
[9725]63 plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2)
[5606]64 if printingflag,
65 set(gcf,'Color','w')
66 printmodel(['ismipastokesvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]67 system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]68 end
[9725]69 plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3)
[5606]70 if printingflag,
71 set(gcf,'Color','w')
72 printmodel(['ismipastokesvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]73 system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]74 end
[9725]75 plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4)
[5606]76 if printingflag,
77 set(gcf,'Color','w')
78 printmodel(['ismipastokesvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]79 system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]80 end
[5049]81
[5606]82 if(L==5000),
[9725]83 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
[5606]84 'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
85 elseif(L==10000),
[9725]86 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
[5606]87 'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
88 elseif(L==20000),
[9725]89 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
[5606]90 'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
91 elseif(L==40000),
[9725]92 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
[5606]93 'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
94 elseif(L==80000),
[9725]95 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
[5606]96 'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
97 elseif(L==160000),
[9725]98 plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
[5606]99 'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
100 end
101 if printingflag,
102 set(gcf,'Color','w')
103 printmodel(['ismipastokesvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]104 system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]105 end
[5049]106end
[5098]107
[5606]108%Now plot the min and max values of vx for each size of the square
109plot([5 10 20 40 80 160],minvx);ylim([0 18])
110if printingflag,
111 set(gcf,'Color','w')
112 printmodel('ismipastokesminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]113 system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]114end
115plot([5 10 20 40 80 160],maxvx);ylim([0 120])
116if printingflag,
117 set(gcf,'Color','w')
118 printmodel('ismipastokesmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off');
[6088]119 system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']);
[5606]120end
[5098]121%Fields and tolerances to track changes
122field_names ={...
123 'Vx5km','Vy5km','Vz5km',...
124 'Vx10km','Vy10km','Vz10km',...
125 'Vx20km','Vy20km','Vz20km',...
126 'Vx40km','Vy40km','Vz40km',...
127 'Vx80km','Vy80km','Vz80km',...
128 'Vx160km','Vy160km','Vz160km'
129}
130field_tolerances={...
[8228]131 1e-12,1e-12,1e-12,...
132 1e-12,1e-12,1e-12,...
[8384]133 1e-12,1e-11,1e-12,...
134 1e-12,1e-11,1e-12,...
135 1e-12,1e-11,1e-12,...
[9877]136 1e-12,1e-11,1e-12,...
[5098]137};
138field_values={};
139for i=1:6,
140 result=results{i};
141 field_values={field_values{:},...
[10976]142 (result.Vx),...
143 (result.Vy),...
144 (result.Vz),...
[5098]145 };
146end
Note: See TracBrowser for help on using the repository browser.