source: issm/trunk/test/NightlyRun/test1501.py@ 14067

Last change on this file since 14067 was 14067, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 14066
File size: 7.8 KB
Line 
1import numpy
2import sys
3from model import *
4from triangle import *
5from setmask import *
6from parameterize import *
7from setflowequation import *
8from EnumDefinitions import *
9from solve import *
10from MatlabFuncs import *
11
12printingflag = False
13
14md=triangle(model(),'../Exp/Square.exp',350000.)
15md=setmask(md,'all','')
16md=parameterize(md,'../Par/SquareShelf.py')
17md=setflowequation(md,'macayeal','all')
18md.cluster=generic('name',oshostname(),'np',3)
19md.transient.isthermal=False
20
21md.timestepping.time_step=1.
22md.settings.output_frequency=1
23md.timestepping.final_time=2000.
24
25#Solve for thinning rate -> -1 * surface mass balance
26smb= 2.*numpy.ones((md.mesh.numberofvertices,1))
27md.surfaceforcings.mass_balance= smb
28md.basalforcings.melting_rate= smb
29
30md=solve(md,PrognosticSolutionEnum())
31
32for i in xrange(1,11):
33 md=solve(md,PrognosticSolutionEnum())
34 md.surfaceforcings.mass_balance= md.surfaceforcings.mass_balance - ((md.results.PrognosticSolution.Thickness)-md.geometry.thickness)
35
36#Set up transient
37smb = md.surfaceforcings.mass_balance
38
39#tooth= [ [ones(400,1)*(smb') - 10.]' [ones(400,1)*(smb')]' ];
40tooth=numpy.hstack((numpy.tile(smb-10.,(1,400)),numpy.tile(smb,(1,400))))
41#smb=[ [ones(399,1)*(smb')]' smb tooth tooth];
42smb=numpy.hstack((numpy.tile(smb,(1,399)),smb,tooth,tooth))
43
44#md.surfaceforcings.mass_balance= smb;
45#md.surfaceforcings.mass_balance(end+1,:)=[1.:2000.];
46md.surfaceforcings.mass_balance=numpy.vstack((smb,numpy.arange(1,2001)))
47
48md=solve(md,TransientSolutionEnum())
49
50#Fields and tolerances to track changes
51field_names=['Vx1','Vy1','Vel1','Pressure1','Bed1','Surface1','Thickness1','SurfaceforcingsMassBalance1', \
52 'Vx2','Vy2','Vel2','Pressure2','Bed2','Surface2','Thickness2','SurfaceforcingsMassBalance2', \
53 'Vx3','Vy3','Vel3','Pressure3','Bed3','Surface3','Thickness3','SurfaceforcingsMassBalance3', \
54 'Vx4','Vy4','Vel4','Pressure4','Bed4','Surface4','Thickness4','SurfaceforcingsMassBalance4', \
55 'Vx5','Vy5','Vel5','Pressure5','Bed5','Surface5','Thickness5','SurfaceforcingsMassBalance5']
56field_tolerances=[1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
57 1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
58 1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
59 1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10,\
60 1e-09,1e-09,1e-09,1e-10,1e-10,1e-10,1e-10,1e-10]
61field_values=[\
62 md.results.TransientSolution[400-1].Vx,\
63 md.results.TransientSolution[400-1].Vy,\
64 md.results.TransientSolution[400-1].Vel,\
65 md.results.TransientSolution[400-1].Pressure,\
66 md.results.TransientSolution[400-1].Bed,\
67 md.results.TransientSolution[400-1].Surface,\
68 md.results.TransientSolution[400-1].Thickness,\
69 md.results.TransientSolution[400-1].SurfaceforcingsMassBalance,\
70 md.results.TransientSolution[800-1].Vx,\
71 md.results.TransientSolution[800-1].Vy,\
72 md.results.TransientSolution[800-1].Vel,\
73 md.results.TransientSolution[800-1].Pressure,\
74 md.results.TransientSolution[800-1].Bed,\
75 md.results.TransientSolution[800-1].Surface,\
76 md.results.TransientSolution[800-1].Thickness,\
77 md.results.TransientSolution[800-1].SurfaceforcingsMassBalance,\
78 md.results.TransientSolution[1200-1].Vx,\
79 md.results.TransientSolution[1200-1].Vy,\
80 md.results.TransientSolution[1200-1].Vel,\
81 md.results.TransientSolution[1200-1].Pressure,\
82 md.results.TransientSolution[1200-1].Bed,\
83 md.results.TransientSolution[1200-1].Surface,\
84 md.results.TransientSolution[1200-1].Thickness,\
85 md.results.TransientSolution[1200-1].SurfaceforcingsMassBalance,\
86 md.results.TransientSolution[1600-1].Vx,\
87 md.results.TransientSolution[1600-1].Vy,\
88 md.results.TransientSolution[1600-1].Vel,\
89 md.results.TransientSolution[1600-1].Pressure,\
90 md.results.TransientSolution[1600-1].Bed,\
91 md.results.TransientSolution[1600-1].Surface,\
92 md.results.TransientSolution[1600-1].Thickness,\
93 md.results.TransientSolution[1600-1].SurfaceforcingsMassBalance,\
94 md.results.TransientSolution[2000-1].Vx,\
95 md.results.TransientSolution[2000-1].Vy,\
96 md.results.TransientSolution[2000-1].Vel,\
97 md.results.TransientSolution[2000-1].Pressure,\
98 md.results.TransientSolution[2000-1].Bed,\
99 md.results.TransientSolution[2000-1].Surface,\
100 md.results.TransientSolution[2000-1].Thickness,\
101 md.results.TransientSolution[2000-1].SurfaceforcingsMassBalance,\
102 ]
103
104if printingflag:
105 pass
106
107 """
108 starttime = 360;
109 endtime = 2000;
110 res = 40;
111 ts = [starttime:res:endtime];
112
113 index = md.mesh.elements;
114 x1=md.mesh.x(index(:,1)); x2=md.mesh.x(index(:,2)); x3=md.mesh.x(index(:,3));
115 y1=md.mesh.y(index(:,1)); y2=md.mesh.y(index(:,2)); y3=md.mesh.y(index(:,3));
116 areas=(0.5*((x2-x1).*(y3-y1)-(y2-y1).*(x3-x1)));
117
118 thickness = [];
119 volume = [];
120 massbal = [];
121 velocity = [];
122 for t=starttime:endtime
123 thickness = [thickness (md.results.TransientSolution(t).Thickness)];
124 volume = [volume mean(md.results.TransientSolution(t).Thickness.value,2).*areas];
125 massbal = [massbal (md.results.TransientSolution(t).SurfaceforcingsMassBalance)];
126 velocity = [velocity (md.results.TransientSolution(t).Vel)];
127 end
128
129 figure('Position', [0 0 860 932])
130
131 options = plotoptions('data','transient_movie','unit','km');
132 options = options.list{1};
133 options = checkplotoptions(md,options);
134
135 %loop over the time steps
136 results=md.results.TransientSolution;
137 count = 1;
138 for i=ts
139
140 subplot(5,9,[28:31 37:40])
141 set(gca,'pos',get(gca,'pos')+[-0.08 -0.08 0.07 0.08])
142 field = 'Thickness';
143
144 %process data
145 [x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
146 [data datatype]=processdata(md,results(i).(field),options);
147
148 titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
149 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
150 options=changefieldvalue(options,'title',titlestring);
151 options=addfielddefault(options,'colorbar',1);
152 options=changefieldvalue(options,'caxis',[0 max(max(thickness))]);
153 applyoptions(md,[],options);
154
155 subplot(5,9,[33:36 42:45])
156 set(gca,'pos',get(gca,'pos')+[-0.00 -0.08 0.07 0.08])
157 field = 'Vel';
158
159 %process data
160 [x y z elements is2d isplanet]=processmesh(md,results(i).(field),options);
161 [data datatype]=processdata(md,results(i).(field),options);
162
163 titlestring=[field ' at time ' num2str(results(i).time/md.constants.yts) ' year'];
164 plot_unit(x,y,z,elements,data,is2d,isplanet,datatype,options)
165 options=changefieldvalue(options,'title',titlestring);
166 options=addfielddefault(options,'colorbar',1);
167 options=changefieldvalue(options,'caxis',[0 max(max(velocity))]);
168 applyoptions(md,[],options);
169
170 subplot(5,4,1:4)
171 cla
172 set(gca,'pos',get(gca,'pos')+[-0.07 0.03 0.12 0.015])
173 plot(starttime:endtime,mean(massbal),'k','LineWidth', 4)
174 hold on
175 ya = ylim;
176 plot([i i], ya, 'r', 'LineWidth',6)
177 ylim(ya); xlim([starttime endtime]);
178 title('Surface Mass Balance','FontSize',14)
179 ylabel('m/year','FontSize',14)
180
181 subplot(5,4,5:8)
182 cla
183 set(gca,'pos',get(gca,'pos')+[-0.07 0.015 0.12 0.015])
184 plot(starttime:endtime,sum(volume)/1000/1000/1000,'LineWidth',4)
185 hold on
186 ya = ylim;
187 plot([i i], ya, 'r', 'LineWidth',6)
188 ylim(ya); xlim([starttime endtime]);
189 title('Ice Volume','FontSize',14)
190 ylabel('km^3','FontSize',14)
191
192 subplot(5,4,9:12)
193 cla
194 set(gca,'pos',get(gca,'pos')+[-0.07 0 0.12 0.015])
195 plot(starttime:endtime,mean(velocity)/1000, 'LineWidth', 4)
196 hold on
197 ya = ylim;
198 plot([i i], ya, 'r', 'LineWidth',6)
199 ylim(ya); xlim([starttime endtime]);
200 title('Mean Velocity','FontSize', 14)
201 ylabel('km/year','FontSize', 14)
202 xlabel('year','FontSize', 14)
203
204 set(gcf,'Renderer','zbuffer','color','white'); %fixes a bug on Mac OS X (not needed in future Matlab version)
205 if i==starttime,
206 %initialize images and frame
207 frame=getframe(gcf);
208 [images,map]=rgb2ind(frame.cdata,256,'nodither');
209 images(1,1,1,length(ts))=0;
210 else
211 frame=getframe(gcf);
212 images(:,:,1,count) = rgb2ind(frame.cdata,map,'nodither');
213 end
214
215 count = count+1;
216
217 end
218
219 filename='transawtooth2d.gif';
220 imwrite(images,map,filename,'DelayTime',1.0,'LoopCount',inf)
221 """
222
Note: See TracBrowser for help on using the repository browser.