source: issm/oecreview/Archive/24684-25833/ISSM-25171-25172.diff@ 25834

Last change on this file since 25834 was 25834, checked in by Mathieu Morlighem, 4 years ago

CHG: added 24684-25833
File size: 11.8 KB
RevLine 
[25834]1Index: ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py
2===================================================================
3--- ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py (revision 25171)
4+++ ../trunk-jpl/test/Par/GiaIvinsBenchmarksAB.py (revision 25172)
5@@ -14,29 +14,29 @@
6 rad = 600000.
7 nv = md.mesh.numberofvertices
8 if (np.isnan(md.geometry.thickness)):
9- md.geometry.thickness = np.zeros((md.mesh.numberofvertices, ))
10+ md.geometry.thickness = np.zeros((md.mesh.numberofvertices, 1))
11 for i in range(nv):
12 dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2)
13 if (dist <= rad):
14 md.geometry.thickness[i] = 2000.0
15 else:
16- md.geometry.thickness[i] = 1.0 # non - zero thickness
17+ md.geometry.thickness[i] = 1.0 # non-zero thickness
18
19 md.geometry.thickness = md.geometry.thickness.reshape(-1, 1)
20-md.geometry.base = np.zeros((md.mesh.numberofvertices, ))
21+md.geometry.base = np.zeros((md.mesh.numberofvertices, 1))
22 md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1, 1) #would otherwise create a 91x91 matrix
23
24-#Ice density used for benchmarking, not 917 kg / m^3
25-md.materials.rho_ice = 1000 #kg m^3
26+#Ice density used for benchmarking, not 917 kg/m^3
27+md.materials.rho_ice = 1000 #kg/m^3
28
29-#GIA parameters specific to Experiments A and B
30+#GIA parameters specific to Experiments A and B
31 md.gia=giaivins();
32-md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, )) #in Pa.s
33-md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, )) #in km
34+md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, 1)) #in Pa.s
35+md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, 1)) #in km
36 md.materials.lithosphere_shear_modulus = 6.7 * 1e10 #in Pa
37-md.materials.lithosphere_density = 3.36 #in g / cm^3
38+md.materials.lithosphere_density = 3.36 #in g/cm^3
39 md.materials.mantle_shear_modulus = 1.45 * 1e11 #in Pa
40-md.materials.mantle_density = 3.38 #in g / cm^3
41+md.materials.mantle_density = 3.38 #in g/cm^3
42
43 #Initial velocity
44 x = archread('../Data/SquareSheetConstrained.arch', 'x')
45@@ -52,19 +52,19 @@
46 x = None
47 y = None
48 index = None
49-md.initialization.vz = np.zeros((md.mesh.numberofvertices, ))
50-md.initialization.pressure = np.zeros((md.mesh.numberofvertices, ))
51+md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1))
52+md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1))
53
54 #Materials
55-md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, ))
56+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, 1))
57 md.materials.rheology_B = paterson(md.initialization.temperature)
58-md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, ))
59+md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, 1))
60
61 #Friction
62-md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, ))
63+md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, 1))
64 md.friction.coefficient[np.where(md.mask.ocean_levelset < 0.)] = 0.
65-md.friction.p = np.ones((md.mesh.numberofelements, ))
66-md.friction.q = np.ones((md.mesh.numberofelements, ))
67+md.friction.p = np.ones((md.mesh.numberofelements, 1))
68+md.friction.q = np.ones((md.mesh.numberofelements, 1))
69
70 #Numerical parameters
71 md.groundingline.migration = 'None'
72@@ -75,7 +75,7 @@
73 md.stressbalance.restol = 0.05
74 md.steadystate.reltol = 0.05
75 md.stressbalance.reltol = 0.05
76-md.stressbalance.abstol = float('NaN')
77+md.stressbalance.abstol = np.nan
78 md.timestepping.time_step = 1.
79 md.timestepping.final_time = 3.
80
81Index: ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py
82===================================================================
83--- ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py (revision 25171)
84+++ ../trunk-jpl/test/Par/GiaIvinsBenchmarksCD.py (revision 25172)
85@@ -1,18 +1,22 @@
86 #Geometry specific to Experiments C and D
87
88+import inspect
89 import os.path
90-import inspect
91+
92 import numpy as np
93+
94 from arch import *
95+from giaivins import giaivins
96 from InterpFromMeshToMesh2d import *
97 from paterson import *
98 from verbose import *
99 from SetIceSheetBC import *
100
101+
102 rad = 800000
103 nv = md.mesh.numberofvertices
104 if (np.isnan(md.geometry.thickness)):
105- md.geometry.thickness = np.zeros((md.mesh.numberofvertices, ))
106+ md.geometry.thickness = np.zeros((md.mesh.numberofvertices, 1))
107 for i in range(nv):
108 dist = np.sqrt(md.mesh.x[i]**2 + md.mesh.y[i]**2)
109 if (dist <= rad):
110@@ -21,15 +25,16 @@
111 md.geometry.thickness[i] = 1.0 # non - zero thickness
112
113 md.geometry.thickness = md.geometry.thickness.reshape(-1, 1)
114-md.geometry.base = np.zeros((md.mesh.numberofvertices, ))
115+md.geometry.base = np.zeros((md.mesh.numberofvertices, 1))
116 md.geometry.surface = md.geometry.thickness + md.geometry.base.reshape(-1, 1) #would otherwise create a 91x91 matrix
117
118 #Ice density used for benchmarking, not 917 kg / m^3
119-md.materials.rho_ice = 1000 #kg m^3
120+md.materials.rho_ice = 1000 #kg m^3
121
122-#GIA parameters specific to Experiments A and B
123-md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, )) #in Pa.s
124-md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, )) #in km
125+#GIA parameters specific to Experiments A and B
126+md.gia=giaivins();
127+md.gia.mantle_viscosity = 1e21 * np.ones((md.mesh.numberofvertices, 1)) #in Pa.s
128+md.gia.lithosphere_thickness = 100 * np.ones((md.mesh.numberofvertices, 1)) #in km
129 md.materials.lithosphere_shear_modulus = 6.7e10 #in Pa
130 md.materials.lithosphere_density = 3.32 #in g / cm^3
131 md.materials.mantle_shear_modulus = 1.45e11 #in Pa
132@@ -49,19 +54,19 @@
133 x = None
134 y = None
135 index = None
136-md.initialization.vz = np.zeros((md.mesh.numberofvertices, ))
137-md.initialization.pressure = np.zeros((md.mesh.numberofvertices, ))
138+md.initialization.vz = np.zeros((md.mesh.numberofvertices, 1))
139+md.initialization.pressure = np.zeros((md.mesh.numberofvertices, 1))
140
141 #Materials
142-md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, ))
143+md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, 1))
144 md.materials.rheology_B = paterson(md.initialization.temperature)
145-md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, ))
146+md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, 1))
147
148 #Friction
149-md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, ))
150+md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices, 1))
151 md.friction.coefficient[np.where(md.mask.ocean_levelset < 0.)] = 0.
152-md.friction.p = np.ones((md.mesh.numberofelements, ))
153-md.friction.q = np.ones((md.mesh.numberofelements, ))
154+md.friction.p = np.ones((md.mesh.numberofelements, 1))
155+md.friction.q = np.ones((md.mesh.numberofelements, 1))
156
157 #Numerical parameters
158 md.groundingline.migration = 'None'
159@@ -72,7 +77,7 @@
160 md.stressbalance.restol = 0.05
161 md.steadystate.reltol = 0.05
162 md.stressbalance.reltol = 0.05
163-md.stressbalance.abstol = float('NaN')
164+md.stressbalance.abstol = np.nan
165 md.timestepping.time_step = 1.
166 md.timestepping.final_time = 3.
167
168Index: ../trunk-jpl/test/NightlyRun/test2071.py
169===================================================================
170--- ../trunk-jpl/test/NightlyRun/test2071.py (revision 25171)
171+++ ../trunk-jpl/test/NightlyRun/test2071.py (revision 25172)
172@@ -16,11 +16,11 @@
173 md = parameterize(md, '../Par/GiaIvinsBenchmarksCD.py')
174
175 #indicate what you want to compute
176-md.gia.cross_section_shape = 1 # for square-edged x - section
177+md.gia.cross_section_shape = 1 # for square-edged x-section
178
179 #define loading history
180 md.timestepping.start_time = 0.3 # for t \approx 0 kyr : to get eleastic response!
181-md.timestepping.final_time = 2500000 # 2, 500 kyr
182+md.timestepping.final_time = 2500000 # 2,500 kyr
183 md.geometry.thickness = np.array([np.append(md.geometry.thickness * 0.0, 0.0),
184 np.append(md.geometry.thickness / 2.0, 0.1),
185 np.append(md.geometry.thickness, 0.2),
186Index: ../trunk-jpl/src/m/classes/model.m
187===================================================================
188--- ../trunk-jpl/src/m/classes/model.m (revision 25171)
189+++ ../trunk-jpl/src/m/classes/model.m (revision 25172)
190@@ -1263,7 +1263,6 @@
191 md.calving = calving();
192 md.frontalforcings = frontalforcings();
193 md.gia = giamme();
194- md.esa = esa();
195 md.love = fourierlove();
196 md.esa = esa();
197 md.autodiff = autodiff();
198Index: ../trunk-jpl/src/m/classes/model.py
199===================================================================
200--- ../trunk-jpl/src/m/classes/model.py (revision 25171)
201+++ ../trunk-jpl/src/m/classes/model.py (revision 25172)
202@@ -86,29 +86,27 @@
203 '''
204 self.mesh = mesh2d()
205 self.mask = mask()
206+ self.constants = constants()
207 self.geometry = geometry()
208- self.constants = constants()
209+ self.initialization = initialization()
210 self.smb = SMBforcing()
211 self.basalforcings = basalforcings()
212+ self.friction = friction()
213+ self.rifts = rifts()
214+ self.solidearth = solidearth()
215+ self.dsl = dsl()
216+ self.timestepping = timestepping()
217+ self.groundingline = groundingline()
218 self.materials = matice()
219 self.damage = damage()
220- self.friction = friction()
221 self.flowequation = flowequation()
222- self.timestepping = timestepping()
223- self.initialization = initialization()
224- self.rifts = rifts()
225- self.dsl = dsl()
226- self.solidearth = solidearth()
227-
228 self.debug = debug()
229 self.verbose = verbose()
230 self.settings = issmsettings()
231 self.toolkits = toolkits()
232 self.cluster = generic()
233-
234 self.balancethickness = balancethickness()
235 self.stressbalance = stressbalance()
236- self.groundingline = groundingline()
237 self.hydrology = hydrologyshreve()
238 self.masstransport = masstransport()
239 self.thermal = thermal()
240@@ -124,10 +122,9 @@
241 self.inversion = inversion()
242 self.qmu = qmu()
243 self.amr = amr()
244-
245+ self.radaroverlay = radaroverlay()
246 self.results = results()
247 self.outputdefinition = outputdefinition()
248- self.radaroverlay = radaroverlay()
249 self.miscellaneous = miscellaneous()
250 self.private = private()
251 #}}}
252@@ -136,19 +133,20 @@
253 # ordered list of properties since vars(self) is random
254 return ['mesh',
255 'mask',
256+ 'constants',
257 'geometry',
258- 'constants',
259+ 'initialization',
260 'smb',
261 'basalforcings',
262+ 'friction',
263+ 'rifts',
264+ 'solidearth',
265+ 'dsl',
266+ 'timestepping',
267+ 'groundingline',
268 'materials',
269 'damage',
270- 'friction',
271 'flowequation',
272- 'timestepping',
273- 'initialization',
274- 'rifts',
275- 'dsl',
276- 'solidearth',
277 'debug',
278 'verbose',
279 'settings',
280@@ -156,7 +154,6 @@
281 'cluster',
282 'balancethickness',
283 'stressbalance',
284- 'groundingline',
285 'hydrology',
286 'masstransport',
287 'thermal',
288@@ -165,16 +162,16 @@
289 'levelset',
290 'calving',
291 'frontalforcings',
292+ 'gia',
293 'love',
294- 'gia',
295 'esa',
296 'autodiff',
297 'inversion',
298 'qmu',
299 'amr',
300+ 'radaroverlay',
301+ 'results',
302 'outputdefinition',
303- 'results',
304- 'radaroverlay',
305 'miscellaneous',
306 'private']
307 # }}}
Note: See TracBrowser for help on using the repository browser.