source: issm/trunk-jpl/test/NightlyRun/test1207.m@ 10976

Last change on this file since 10976 was 10976, checked in by Mathieu Morlighem, 13 years ago

removed all PatchToVec not needed anymore as default is md.settings.results_as_patches = 0

  • Property svn:executable set to *
File size: 3.3 KB
Line 
1%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
2printingflag=false;
3
4numlayers=10;
5resolution=30000;
6
7%To begin with the numerical model
8md=model;
9md=roundmesh(md,750000,resolution);
10md=setmask(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution
11md=parameterize(md,'../Par/RoundSheetStaticEISMINT.par');
12
13%Calculation of the analytical 2d velocity field
14constant=0.3;
15vx_obs=constant/2*md.mesh.x.*(md.geometry.thickness).^-1;
16vy_obs=constant/2*md.mesh.y.*(md.geometry.thickness).^-1;
17vel_obs=(sqrt((md.inversion.vx_obs).^2+(md.inversion.vy_obs).^2));
18
19%We extrude the model to have a 3d model
20md=extrude(md,numlayers,1);
21md=setflowequation(md,'stokes','all');
22
23%Spc the nodes on the bed
24pos=find(md.mesh.vertexonbed);
25md.diagnostic.spcvx(pos)=0;
26md.diagnostic.spcvy(pos)=0;
27md.diagnostic.spcvz(pos)=0;
28
29%Now we can solve the problem
30md.cluster=generic('name',oshostname(),'np',8);
31md=solve(md,DiagnosticSolutionEnum);
32
33%Calculate the depth averaged velocity field (2d):
34vx=(md.results.DiagnosticSolution.Vx);
35vy=(md.results.DiagnosticSolution.Vy);
36vel=zeros(md.mesh.numberofvertices2d,1);
37
38node_vel=0;
39for i=1:md.mesh.numberofvertices2d
40 for j=1:(md.mesh.numberoflayers-1)
41 node_vel=node_vel+1/(2*(md.mesh.numberoflayers-1))*(sqrt(vx(i+j*md.mesh.numberofvertices2d,1).^2+...
42 vy(i+j*md.mesh.numberofvertices2d,1).^2)+...
43 sqrt(vx(i+(j-1)*md.mesh.numberofvertices2d,1).^2+vy(i+(j-1)*md.mesh.numberofvertices2d,1).^2));
44 end
45 vel(i,1)=node_vel;
46 node_vel=0;
47end
48
49%Plot of the velocity from the exact and calculated solutions
50figure(1)
51subplot(2,2,1)
52p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',...
53vel,'FaceColor','interp','EdgeColor','none');
54title('Modelled velocity','FontSize',14,'FontWeight','bold')
55colorbar;
56caxis([0 200]);
57
58subplot(2,2,2)
59p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',...
60vel_obs,'FaceColor','interp','EdgeColor','none');
61title('Analytical velocity','FontSize',14,'FontWeight','bold')
62colorbar;
63caxis([0 200]);
64
65subplot(2,2,3)
66hold on;
67plot(sqrt((md.mesh.x(1:md.mesh.numberofvertices2d)).^2+(md.mesh.y(1:md.mesh.numberofvertices2d)).^2),vel,'r.');
68plot(sqrt((md.mesh.x2d).^2+(md.mesh.y2d).^2),vel_obs,'b.');
69title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold');
70xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold');
71ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold');
72legend('calculated velocity','exact velocity');
73axis([0 750000 0 200]);
74hold off;
75
76subplot(2,2,4)
77p=patch('Faces',md.mesh.elements2d,'Vertices',[md.mesh.x2d md.mesh.y2d],'FaceVertexCData',...
78abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none');
79title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
80colorbar;
81caxis([0 100]);
82
83if printingflag,
84 set(gcf,'Color','w')
85 printmodel('stokesstatic','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
86 system(['mv stokesstatic.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT/IceSheet']);
87end
88
89%Fields and tolerances to track changes
90field_names ={ ...
91 'Vx','Vy','Vel', ...
92};
93field_tolerances={...
94 1e-12,1e-12,1e-12, ...
95};
96field_values={
97 vx,vy,vel, ...
98};
Note: See TracBrowser for help on using the repository browser.