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

Last change on this file was 24313, checked in by Mathieu Morlighem, 5 years ago

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