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