source: issm/trunk-jpl/test/NightlyRun/test1205.py@ 21056

Last change on this file since 21056 was 21056, checked in by Mathieu Morlighem, 9 years ago

CHG: do not request solution in solution string

  • Property svn:executable set to *
File size: 3.6 KB
Line 
1#Test Name: EISMINTRoundIceSheetStaticSIA
2import numpy
3from model import *
4from roundmesh import *
5from setmask import *
6from parameterize import *
7from setflowequation import *
8from EnumDefinitions import *
9from solve import *
10from MatlabFuncs import *
11
12"""
13The aim of this program is to compare a model with an analytical solution given in SSA EISMINT : Lessons in Ice-Sheet Modeling.
14"""
15
16printingflag=False
17
18numlayers=10
19resolution=30000.
20
21#To begin with the numerical model
22md=model()
23md=roundmesh(md,750000.,resolution)
24md=setmask(md,'','') #We can not test iceshelves nor ice rises with this analytical solution
25md=parameterize(md,'../Par/RoundSheetStaticEISMINT.py')
26
27#Calculation of the analytical 2d velocity field
28constant=0.3
29vx_obs=constant/2.*md.mesh.x*(md.geometry.thickness)**-1
30vy_obs=constant/2.*md.mesh.y*(md.geometry.thickness)**-1
31vel_obs=numpy.sqrt((md.inversion.vx_obs)**2+(md.inversion.vy_obs)**2)
32
33#We extrude the model to have a 3d model
34md.extrude(numlayers,1.);
35md=setflowequation(md,'SIA','all')
36
37#Spc the nodes on the bed
38pos=numpy.nonzero(md.mesh.vertexonbase)
39md.stressbalance.spcvx[pos]=0.
40md.stressbalance.spcvy[pos]=0.
41md.stressbalance.spcvz[pos]=0.
42
43#Now we can solve the problem
44md.cluster=generic('name',oshostname(),'np',8)
45md=solve(md,'Stressbalance')
46
47#Calculate the depth averaged velocity field (2d):
48vx=md.results.StressbalanceSolution.Vx
49vy=md.results.StressbalanceSolution.Vy
50vel=numpy.zeros((md.mesh.numberofvertices2d,1))
51
52for 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
95if 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
102field_names =[ \
103 'Vx','Vy','Vel', \
104]
105field_tolerances=[ \
106 1e-13,1e-13,1e-13, \
107]
108field_values=[ \
109 vx,vy,vel, \
110]
Note: See TracBrowser for help on using the repository browser.