%This test is a test from the ISMP-HOM Intercomparison project %Pattyn and Payne 2006 printingflag=false; L_list={5000,10000,20000,40000,80000,160000}; results={}; minvx=[]; maxvx=[]; for i=1:length(L_list), L=L_list{i}; nx=20; %numberof nodes in x direction ny=20; md=model(); md=squaremesh(md,L,L,nx,ny); md=setmask(md,'',''); %ice sheet test % %Find elements at the corner and extract model % posnodes=find((md.mesh.x==0 | md.mesh.x==max(md.mesh.x)) & (md.mesh.y==0 | md.mesh.y==max(md.mesh.y))); % [a,b]=find(ismember(md.mesh.elements,posnodes)); % elements=ones(md.mesh.numberofelements,1); % elements(a)=0; % md=modelextract(md,elements); md=parameterize(md,'../Par/ISMIPA.par'); md=extrude(md,10,1); md=setflowequation(md,'stokes','all'); %Create dirichlet on the bed only pos=find(md.mesh.vertexonbed); md.diagnostic.spcvx(pos)=0; md.diagnostic.spcvy(pos)=0; md.diagnostic.spcvz(pos)=0; %Create MPCs to have periodic boundary conditions %posx=find(md.mesh.x==0); %posx2=find(md.mesh.x==max(md.mesh.x)); %posx=find(md.mesh.x==0 & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed); %posx2=find(md.mesh.x==max(md.mesh.x) & md.mesh.y~=0 & md.mesh.y~=max(md.mesh.y) & ~md.mesh.vertexonbed); %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 %posy2=find(md.mesh.y==max(md.mesh.y) & md.mesh.x~=0 & md.mesh.x~=max(md.mesh.x) & ~md.mesh.vertexonbed); %md.diagnostic.vertex_pairing=[posx,posx2;posy,posy2]; %Compute the diagnostic md.diagnostic.abstol=NaN; md.diagnostic.reltol=NaN; md.diagnostic.restol=1; md.cluster=generic('name',oshostname(),'np',8); md=solve(md,DiagnosticSolutionEnum()); %Plot the results and save them vx=(md.results.DiagnosticSolution.Vx); vy=(md.results.DiagnosticSolution.Vy); vz=(md.results.DiagnosticSolution.Vz); pressure=(md.results.DiagnosticSolution.Pressure); results{i}=md.results.DiagnosticSolution; minvx(i)=min(vx(end-md.mesh.numberofvertices2d+1:end)); maxvx(i)=max(vx(end-md.mesh.numberofvertices2d+1:end)); %Now plot vx, vy, vz and vx on a cross section plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',2) if printingflag, set(gcf,'Color','w') printmodel(['ismipastokesvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesvx' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',3) if printingflag, set(gcf,'Color','w') printmodel(['ismipastokesvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesvy' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km','figure',4) if printingflag, set(gcf,'Color','w') printmodel(['ismipastokesvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesvz' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end if(L==5000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','') elseif(L==10000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','') elseif(L==20000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','') elseif(L==40000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','') elseif(L==80000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','') elseif(L==160000), plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,... 'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','') end if printingflag, set(gcf,'Color','w') printmodel(['ismipastokesvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesvxsec' num2str(L) '.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end end %Now plot the min and max values of vx for each size of the square plot([5 10 20 40 80 160],minvx);ylim([0 18]) if printingflag, set(gcf,'Color','w') printmodel('ismipastokesminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesminvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end plot([5 10 20 40 80 160],maxvx);ylim([0 120]) if printingflag, set(gcf,'Color','w') printmodel('ismipastokesmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off'); system(['mv ismipastokesmaxvx.png ' ISSM_DIR '/website/doc_pdf/validation/Images/ISMIP/TestA ']); end %Fields and tolerances to track changes field_names ={... 'Vx5km','Vy5km','Vz5km',... 'Vx10km','Vy10km','Vz10km',... 'Vx20km','Vy20km','Vz20km',... 'Vx40km','Vy40km','Vz40km',... 'Vx80km','Vy80km','Vz80km',... 'Vx160km','Vy160km','Vz160km' } field_tolerances={... 1e-12,1e-12,1e-12,... 1e-12,1e-12,1e-12,... 1e-12,1e-11,1e-12,... 1e-12,1e-11,1e-12,... 1e-12,1e-11,1e-12,... 1e-12,1e-11,1e-12,... }; field_values={}; for i=1:6, result=results{i}; field_values={field_values{:},... (result.Vx),... (result.Vy),... (result.Vz),... }; end