source: issm/trunk/test/NightlyRun/test1101.py@ 21729

Last change on this file since 21729 was 21729, checked in by Mathieu Morlighem, 8 years ago

merged trunk-jpl and trunk for revision 21727
File size: 5.4 KB
Line 
1#Test Name: ISMIPAHO
2import numpy as np
3import shutil
4from model import *
5from socket import gethostname
6from squaremesh import *
7from setmask import *
8from parameterize import *
9from setflowequation import *
10from solve import *
11
12"""
13This test is a test from the ISMP-HOM Intercomparison project.
14Pattyn and Payne 2006
15"""
16
17printingflag=False
18
19#L_list=[80000.]
20L_list=[80000.]
21results=[]
22minvx=[]
23maxvx=[]
24
25for L in L_list:
26 nx=20 #numberof nodes in x direction
27 ny=20
28 md=model()
29 md=squaremesh(md,L,L,nx,ny)
30 md=setmask(md,'','') #ice sheet test
31 md=parameterize(md,'../Par/ISMIPA.py')
32 md.extrude(9,1.)
33
34 md=setflowequation(md,'HO','all')
35
36 #Create dirichlet on the bed only
37 md.stressbalance.spcvx=np.nan*np.ones((md.mesh.numberofvertices))
38 md.stressbalance.spcvy=np.nan*np.ones((md.mesh.numberofvertices))
39 md.stressbalance.spcvz=np.nan*np.ones((md.mesh.numberofvertices))
40
41 pos=np.where(md.mesh.vertexonbase)
42 md.stressbalance.spcvx[pos]=0.
43 md.stressbalance.spcvy[pos]=0.
44
45 #Create MPCs to have periodic boundary conditions
46 posx=np.where(md.mesh.x==0.)[0]
47 posx2=np.where(md.mesh.x==np.max(md.mesh.x))[0]
48
49 posy=np.where(np.logical_and.reduce((md.mesh.y==0.,md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x))))[0] #Don't take the same nodes two times
50 posy2=np.where(np.logical_and.reduce((md.mesh.y==np.max(md.mesh.y),md.mesh.x!=0.,md.mesh.x!=np.max(md.mesh.x))))[0]
51
52 md.stressbalance.vertex_pairing=np.vstack((np.vstack((posx+1,posx2+1)).T,np.vstack((posy+1,posy2+1)).T))
53
54 #Compute the stressbalance
55 md.cluster=generic('name',gethostname(),'np',8)
56 md=solve(md,'Stressbalance')
57
58 #Plot the results and save them
59 vx=md.results.StressbalanceSolution.Vx
60 vy=md.results.StressbalanceSolution.Vy
61 vz=md.results.StressbalanceSolution.Vz
62 results.append(md.results.StressbalanceSolution)
63 minvx.append(np.min(vx[-md.mesh.numberofvertices2d:]))
64 maxvx.append(np.max(vx[-md.mesh.numberofvertices2d:]))
65
66 #Now plot vx, vy, vz and vx on a cross section
67# plotmodel(md,'data',vx,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
68 if printingflag:
69 pass
70# set(gcf,'Color','w')
71# printmodel(['ismipaHOvx' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
72# shutil.move("ismipaHOvx%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
73# plotmodel(md,'data',vy,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
74 if printingflag:
75 pass
76# set(gcf,'Color','w')
77# printmodel(['ismipaHOvy' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
78# shutil.move("ismipaHOvy%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
79# plotmodel(md,'data',vz,'layer#all',md.mesh.numberoflayers,'xlim',[0 L/10^3],'ylim',[0 L/10^3],'unit','km')
80 if printingflag:
81 pass
82# set(gcf,'Color','w')
83# printmodel(['ismipaHOvz' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
84# shutil.move("ismipaHOvz%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
85
86 if (L==5000.):
87 pass
88# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP5000.exp','layer',md.mesh.numberoflayers,...
89# 'resolution',[10 10],'ylim',[10 18],'xlim',[0 5000],'title','','xlabel','')
90 elif (L==10000.):
91 pass
92# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP10000.exp','layer',md.mesh.numberoflayers,...
93# 'resolution',[10 10],'ylim',[10 30],'xlim',[0 10000],'title','','xlabel','')
94 elif (L==20000.):
95 pass
96# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP20000.exp','layer',md.mesh.numberoflayers,...
97# 'resolution',[10 10],'ylim',[0 50],'xlim',[0 20000],'title','','xlabel','')
98 elif (L==40000.):
99 pass
100# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP40000.exp','layer',md.mesh.numberoflayers,...
101# 'resolution',[10 10],'ylim',[0 80],'xlim',[0 40000],'title','','xlabel','')
102 elif (L==80000.):
103 pass
104# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP80000.exp','layer',md.mesh.numberoflayers,...
105# 'resolution',[10 10],'ylim',[0 100],'xlim',[0 80000],'title','','xlabel','')
106 elif (L==160000.):
107 pass
108# plotmodel(md,'data',vx,'sectionvalue','../Exp/ISMIP160000.exp','layer',md.mesh.numberoflayers,...
109# 'resolution',[10 10],'ylim',[0 120],'xlim',[0 160000],'title','','xlabel','')
110 if printingflag:
111 pass
112# set(gcf,'Color','w')
113# printmodel(['ismipaHOvxsec' num2str(L)],'png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
114# shutil.move("ismipaHOvxsec%d.png" % L,ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
115
116#Now plot the min and max values of vx for each size of the square
117#plot([5 10 20 40 80 160],minvx)ylim([0 18])xlim([0 160])
118if printingflag:
119 pass
120# set(gcf,'Color','w')
121# printmodel('ismipaHOminvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
122# shutil.move('ismipaHOminvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
123#plot([5 10 20 40 80 160],maxvx)ylim([0 120])xlim([0 160])
124if printingflag:
125 pass
126# set(gcf,'Color','w')
127# printmodel('ismipaHOmaxvx','png','margin','on','marginsize',25,'frame','off','resolution',1.5,'hardcopy','off')
128# shutil.move('ismipaHOmaxvx.png',ISSM_DIR+'/website/doc_pdf/validation/Images/ISMIP/TestA')
129
130#Fields and tolerances to track changes
131field_names =['Vx80km','Vy80km','Vz80km']
132field_tolerances=[1e-08,1e-08,1e-08]
133field_values=[]
134for result in results:
135 field_values=field_values+[result.Vx,result.Vy,result.Vz]
Note: See TracBrowser for help on using the repository browser.