[5049] | 1 | %This test is a test from the ISMP-HOM Intercomparison project
|
---|
| 2 | %Pattyn and Payne 2006
|
---|
[5606] | 3 | printingflag=false;
|
---|
[5049] | 4 |
|
---|
[5055] | 5 | L_list={5000,10000,20000,40000,80000,160000};
|
---|
[5049] | 6 | results={};
|
---|
[5606] | 7 | minvx=[];
|
---|
| 8 | maxvx=[];
|
---|
[5049] | 9 |
|
---|
| 10 | for 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] | 106 | end
|
---|
[5098] | 107 |
|
---|
[5606] | 108 | %Now plot the min and max values of vx for each size of the square
|
---|
| 109 | plot([5 10 20 40 80 160],minvx);ylim([0 18])
|
---|
| 110 | if 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] | 114 | end
|
---|
| 115 | plot([5 10 20 40 80 160],maxvx);ylim([0 120])
|
---|
| 116 | if 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] | 120 | end
|
---|
[5098] | 121 | %Fields and tolerances to track changes
|
---|
| 122 | field_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 | }
|
---|
| 130 | field_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 | };
|
---|
| 138 | field_values={};
|
---|
| 139 | for i=1:6,
|
---|
| 140 | result=results{i};
|
---|
| 141 | field_values={field_values{:},...
|
---|
[10976] | 142 | (result.Vx),...
|
---|
| 143 | (result.Vy),...
|
---|
| 144 | (result.Vz),...
|
---|
[5098] | 145 | };
|
---|
| 146 | end
|
---|