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