source: issm/oecreview/Archive/24684-25833/ISSM-25805-25806.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: 18.8 KB
RevLine 
[25834]1Index: ../trunk-jpl/src/m/classes/frictionschoof.py
2===================================================================
3--- ../trunk-jpl/src/m/classes/frictionschoof.py (nonexistent)
4+++ ../trunk-jpl/src/m/classes/frictionschoof.py (revision 25806)
5@@ -0,0 +1,79 @@
6+import numpy as np
7+
8+from checkfield import checkfield
9+from fielddisplay import fielddisplay
10+from project3d import project3d
11+from structtoobj import structtoobj
12+from WriteData import WriteData
13+
14+
15+class frictionschoof(object):
16+ """FRICTIONSCHOOF class definition
17+
18+ Usage:
19+ friction = frictionschoof()
20+ """
21+
22+ def __init__(self, *args): # {{{
23+ self.C = np.nan
24+ self.Cmax = np.nan
25+ self.m = np.nan
26+ self.effective_pressure_limit = 0
27+
28+ nargs = len(args)
29+ if nargs == 0:
30+ self.setdefaultparameters()
31+ elif nargs == 1:
32+ self = structtoobj(self, args[0])
33+ else:
34+ raise Exception('constructor not supported')
35+ #}}}
36+
37+ def __repr__(self): # {{{
38+ # See Brondex et al. 2019 https://www.the-cryosphere.net/13/177/2019/
39+ s = 'Schoof sliding law parameters:\n'
40+ s += ' Schoof\'s sliding law reads:\n'
41+ s += ' C^2 |u_b|^(m-1) \n'
42+ s += ' tau_b = - _____________________________ u_b \n'
43+ s += ' (1+(C^2/(Cmax N))^1/m |u_b| )^m \n'
44+ s += '\n'
45+ s += "{}\n".format(fielddisplay(self, 'C', 'friction coefficient [SI]'))
46+ s += "{}\n".format(fielddisplay(self, 'Cmax', 'Iken\'s bound (typically between 0.17 and 0.84) [SI]'))
47+ s += "{}\n".format(fielddisplay(self, 'm', 'm exponent (generally taken as m = 1/n = 1/3)'))
48+ s += "{}\n".format(fielddisplay(self, 'effective_pressure_limit', 'fNeff do not allow to fall below a certain limit: effective_pressure_limit*rho_ice*g*thickness (default 0)'))
49+ return s
50+ #}}}
51+
52+ def setdefaultparameters(self): # {{{
53+ self.effective_pressure_limit = 0
54+ return self
55+ #}}}
56+
57+ def extrude(self, md): # {{{
58+ self.C = project3d(md, 'vector', self.C, 'type', 'node')
59+ self.Cmax = project3d(md, 'vector', self.Cmax, 'type', 'node')
60+ self.m = project3d(md, 'vector', self.m, 'type', 'element')
61+ return self
62+ #}}}
63+
64+ def checkconsistency(self, md, solution, analyses): # {{{
65+ # Early return
66+ if 'StressbalanceAnalysis' not in analyses and 'ThermalAnalysis' not in analyses:
67+ return md
68+ md = checkfield(md, 'fieldname', 'friction.C', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>',0.)
69+ md = checkfield(md, 'fieldname', 'friction.Cmax', 'timeseries', 1, 'NaN', 1, 'Inf', 1, '>', 0.)
70+ md = checkfield(md, 'fieldname', 'friction.m', 'NaN', 1, 'Inf', 1, '>', 0., 'size', [md.mesh.numberofelements, 1])
71+ md = checkfield(md, 'fieldname', 'friction.effective_pressure_limit', 'numel', [1], '>=', 0)
72+ return md
73+ # }}}
74+
75+ def marshall(self, prefix, md, fid): # {{{
76+ yts = md.constants.yts
77+
78+ WriteData(fid, prefix, 'name', 'md.friction.law', 'data', 11, 'format', 'Integer')
79+ WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'C', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
80+ WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'Cmax', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
81+ WriteData(fid, prefix, 'class', 'friction', 'object', self, 'fieldname', 'm', 'format', 'DoubleMat', 'mattype', 2)
82+ WriteData(fid, prefix, 'object', self, 'class', 'friction', 'fieldname', 'effective_pressure_limit', 'format', 'Double')
83+ # }}}
84+
85Index: ../trunk-jpl/src/m/classes/frictionshakti.py
86===================================================================
87--- ../trunk-jpl/src/m/classes/frictionshakti.py (revision 25805)
88+++ ../trunk-jpl/src/m/classes/frictionshakti.py (revision 25806)
89@@ -27,7 +27,7 @@
90 def __repr__(self): # {{{
91 s = 'Basal shear stress parameters: Sigma_b = coefficient^2 * Neff * u_b\n'
92 s += '(effective stress Neff = rho_ice * g * thickness + rho_water * g * (head - b))\n'
93- s += "{}\n".format(fielddisplay(self, "coefficient", "friction coefficient [SI]"))
94+ s += '{}\n'.format(fielddisplay(self, 'coefficient', 'friction coefficient [SI]'))
95 return s
96 #}}}
97
98Index: ../trunk-jpl/src/m/classes/calving.py
99===================================================================
100--- ../trunk-jpl/src/m/classes/calving.py (revision 25805)
101+++ ../trunk-jpl/src/m/classes/calving.py (revision 25806)
102@@ -1,31 +1,28 @@
103+import numpy as np
104+
105+from checkfield import checkfield
106 from fielddisplay import fielddisplay
107 from project3d import project3d
108-from checkfield import checkfield
109 from WriteData import WriteData
110
111
112 class calving(object):
113- """
114- CALVING class definition
115+ """CALVING class definition
116
117- Usage:
118- calving = calving()
119+ Usage:
120+ calving = calving()
121 """
122
123 def __init__(self): # {{{
124-
125- self.calvingrate = float('NaN')
126-
127- #set defaults
128+ self.calvingrate = np.nan
129 self.setdefaultparameters()
130
131 #}}}
132
133 def __repr__(self): # {{{
134- string = ' Calving parameters:'
135- string = "%s\n%s" % (string, fielddisplay(self, 'calvingrate', 'calving rate at given location [m / a]'))
136-
137- return string
138+ s = ' Calving parameters:'
139+ s += '{}\n'.format(fielddisplay(self, 'calvingrate', 'calving rate at given location [m / a]'))
140+ return s
141 #}}}
142
143 def extrude(self, md): # {{{
144@@ -39,7 +36,7 @@
145
146 def checkconsistency(self, md, solution, analyses): # {{{
147 #Early return
148- if (solution != 'TransientSolution') or (not md.transient.ismovingfront):
149+ if solution != 'TransientSolution' or not md.transient.ismovingfront:
150 return md
151
152 md = checkfield(md, 'fieldname', 'calving.calvingrate', '>=', 0, 'timeseries', 1, 'NaN', 1, 'Inf', 1)
153Index: ../trunk-jpl/src/m/classes/inversion.py
154===================================================================
155--- ../trunk-jpl/src/m/classes/inversion.py (revision 25805)
156+++ ../trunk-jpl/src/m/classes/inversion.py (revision 25806)
157@@ -53,9 +53,9 @@
158 s += '{}\n'.format(fielddisplay(self, 'step_threshold', 'decrease threshold for misfit, default is 30%'))
159 s += '{}\n'.format(fielddisplay(self, 'min_parameters', 'absolute minimum acceptable value of the inversed parameter on each vertex'))
160 s += '{}\n'.format(fielddisplay(self, 'max_parameters', 'absolute maximum acceptable value of the inversed parameter on each vertex'))
161- s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m / yr]'))
162- s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m / yr]'))
163- s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m / yr]'))
164+ s += '{}\n'.format(fielddisplay(self, 'vx_obs', 'observed velocity x component [m/yr]'))
165+ s += '{}\n'.format(fielddisplay(self, 'vy_obs', 'observed velocity y component [m/yr]'))
166+ s += '{}\n'.format(fielddisplay(self, 'vel_obs', 'observed velocity magnitude [m/yr]'))
167 s += '{}\n'.format(fielddisplay(self, 'thickness_obs', 'observed thickness [m]'))
168 s += '{}\n'.format(fielddisplay(self, 'surface_obs', 'observed surface elevation [m]'))
169 s += '{}\n'.format('Available cost functions:')
170@@ -159,7 +159,7 @@
171 if not self.iscontrol:
172 return
173 WriteData(fid, prefix, 'object', self, 'fieldname', 'nsteps', 'format', 'Integer')
174- WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter_per_step', 'format', 'DoubleMat', 'mattype', 3)
175+ WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter_per_step', 'format', 'IntMat', 'mattype', 3)
176 WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_functions_coefficients', 'format', 'DoubleMat', 'mattype', 1)
177 WriteData(fid, prefix, 'object', self, 'fieldname', 'gradient_scaling', 'format', 'DoubleMat', 'mattype', 3)
178 WriteData(fid, prefix, 'object', self, 'fieldname', 'cost_function_threshold', 'format', 'Double')
179Index: ../trunk-jpl/src/m/classes/solidearthsettings.py
180===================================================================
181--- ../trunk-jpl/src/m/classes/solidearthsettings.py (revision 25805)
182+++ ../trunk-jpl/src/m/classes/solidearthsettings.py (revision 25806)
183@@ -21,7 +21,7 @@
184 self.rotation = 0
185 self.ocean_area_scaling = 0
186 self.runfrequency = 1 # How many time steps we skip before we run grd_core
187- self.computesealevelchange = 0 # Will grd_core compute sea level?
188+ self.computesealevelchange = 1 # Will grd_core compute sea level?
189 self.isgrd = 1 # Will GRD patterns be computed?
190 self.degacc = 0 # Degree increment for resolution of Green tables
191 self.horiz = 0 # Compute horizontal displacement?
192@@ -65,7 +65,7 @@
193 self.rotation = 1
194 self.ocean_area_scaling = 0
195 self.isgrd = 1
196- self.computesealevelchange = 0
197+ self.computesealevelchange = 1
198
199 # Numerical discretization accuracy
200 self.degacc = .01
201@@ -73,11 +73,11 @@
202 # How many time steps we skip before we run solidearthsettings solver during transient
203 self.runfrequency = 1
204
205+ # Fractional contribution
206+ self.glfraction = 1
207+
208 # Horizontal displacement? (not on by default)
209 self.horiz = 0
210-
211- # Fractional contribution
212- self.glfraction = 1
213 #}}}
214
215 def checkconsistency(self, md, solution, analyses): # {{{
216Index: ../trunk-jpl/src/m/inversions/supportedcostfunctions.py
217===================================================================
218--- ../trunk-jpl/src/m/inversions/supportedcostfunctions.py (revision 25805)
219+++ ../trunk-jpl/src/m/inversions/supportedcostfunctions.py (revision 25806)
220@@ -1,2 +1,2 @@
221 def supportedcostfunctions():
222- return [101, 102, 103, 104, 105, 201, 501, 502, 503, 504, 505]
223+ return list(range(101, 105 + 1)) + [201] + list(range(501, 508 + 1)) + [510] + list(range(601, 604 + 1))
224Index: ../trunk-jpl/src/m/inversions/supportedcontrols.py
225===================================================================
226--- ../trunk-jpl/src/m/inversions/supportedcontrols.py (revision 25805)
227+++ ../trunk-jpl/src/m/inversions/supportedcontrols.py (revision 25806)
228@@ -1,2 +1,15 @@
229 def supportedcontrols():
230- return ['BalancethicknessThickeningRate', 'FrictionCoefficient', 'FrictionAs', 'MaterialsRheologyBbar', 'DamageDbar', 'Vx', 'Vy']
231+ return [
232+ 'BalancethicknessThickeningRate',
233+ 'FrictionCoefficient',
234+ 'FrictionC',
235+ 'FrictionAs',
236+ 'MaterialsRheologyBbar',
237+ 'DamageDbar',
238+ 'Vx',
239+ 'Vy',
240+ 'Thickness',
241+ 'BalancethicknessOmega',
242+ 'BalancethicknessApparentMassbalance',
243+ 'MaterialsRheologyB'
244+ ]
245Index: ../trunk-jpl/src/m/qmu/postqmu.m
246===================================================================
247--- ../trunk-jpl/src/m/qmu/postqmu.m (revision 25805)
248+++ ../trunk-jpl/src/m/qmu/postqmu.m (revision 25806)
249@@ -15,7 +15,7 @@
250 disp(sprintf('%s',fline));
251 fline=fgetl(fide);
252 end
253- warning(['Dakota returned error in ''' qmuerrfile ' file. qmu directory retained.'])
254+ warning(['Dakota returned error in ''' qmuerrfile ''' file. QMU directory retained.'])
255 end
256 status=fclose(fide);
257 end
258Index: ../trunk-jpl/test/NightlyRun/test356.py
259===================================================================
260--- ../trunk-jpl/test/NightlyRun/test356.py (nonexistent)
261+++ ../trunk-jpl/test/NightlyRun/test356.py (revision 25806)
262@@ -0,0 +1,56 @@
263+#Test Name: TransientFrictionSchoof
264+import numpy as np
265+
266+from frictionschoof import frictionschoof
267+from MatlabFuncs import oshostname
268+from model import *
269+from parameterize import parameterize
270+from setflowequation import setflowequation
271+from setmask import setmask
272+from solve import solve
273+from transient import transient
274+from triangle import triangle
275+
276+
277+md = triangle(model(), '../Exp/Square.exp', 200000.)
278+md = setmask(md, '', '')
279+md = parameterize(md, '../Par/SquareSheetConstrained.py')
280+md = setflowequation(md, 'SSA', 'all')
281+
282+# Use Schoof's law
283+Cmax = 0.8
284+md.friction = frictionschoof()
285+md.friction.m = 1.0 / 3.0 * np.ones((md.mesh.numberofelements, 1))
286+md.friction.Cmax = Cmax * np.ones((md.mesh.numberofvertices, 1))
287+md.friction.C = 200 * np.ones((md.mesh.numberofvertices, 1))
288+
289+# Control parameters
290+md.inversion.iscontrol = 1
291+md.inversion.control_parameters = ['FrictionC']
292+md.inversion.min_parameters = 1. * np.ones((md.mesh.numberofvertices, 1))
293+md.inversion.max_parameters = 10000. * np.ones((md.mesh.numberofvertices, 1))
294+md.inversion.nsteps = 2
295+md.inversion.cost_functions = [102, 501]
296+md.inversion.cost_functions_coefficients = np.ones((md.mesh.numberofvertices, 2))
297+md.inversion.cost_functions_coefficients[:, 1] = 2e-7
298+md.inversion.gradient_scaling = 3. * np.ones((md.inversion.nsteps, 1))
299+md.inversion.maxiter_per_step = 2 * np.ones((md.inversion.nsteps, 1))
300+md.inversion.step_threshold = 0.3 * np.ones((md.inversion.nsteps, 1))
301+md.inversion.vx_obs = md.initialization.vx
302+md.inversion.vy_obs= md.initialization.vy
303+
304+md.cluster = generic('name', oshostname(), 'np', 3)
305+md = solve(md, 'Stressbalance')
306+
307+#Fields and tolerances to track changes
308+field_names = ['Gradient', 'Misfits', 'FrictionC', 'Pressure', 'Vel', 'Vx', 'Vy']
309+field_tolerances = [1e-12, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
310+field_values = [
311+ md.results.StressbalanceSolution.Gradient1,
312+ md.results.StressbalanceSolution.J,
313+ md.results.StressbalanceSolution.FrictionC,
314+ md.results.StressbalanceSolution.Pressure,
315+ md.results.StressbalanceSolution.Vel,
316+ md.results.StressbalanceSolution.Vx,
317+ md.results.StressbalanceSolution.Vy
318+]
319Index: ../trunk-jpl/test/NightlyRun/test481.py
320===================================================================
321--- ../trunk-jpl/test/NightlyRun/test481.py (nonexistent)
322+++ ../trunk-jpl/test/NightlyRun/test481.py (revision 25806)
323@@ -0,0 +1,60 @@
324+#Test Name: TransientFrictionSchoof
325+import numpy as np
326+
327+from frictionschoof import frictionschoof
328+from MatlabFuncs import oshostname
329+from model import *
330+from parameterize import parameterize
331+from setflowequation import setflowequation
332+from setmask import setmask
333+from solve import solve
334+from transient import transient
335+from triangle import triangle
336+
337+md = triangle(model(), '../Exp/Square.exp', 150000.)
338+md = setmask(md, '../Exp/SquareShelf.exp', '')
339+md = parameterize(md, '../Par/SquareSheetShelf.py')
340+md.extrude(4, 1)
341+md = setflowequation(md, 'HO', 'all')
342+md.transient.isthermal = 0
343+md.friction = frictionschoof(md.friction)
344+md.friction.C = pow(20.e4, 0.5) * np.ones((md.mesh.numberofvertices, 1))
345+md.friction.Cmax = 0.5 * np.ones((md.mesh.numberofvertices, 1))
346+md.friction.m = 1./3.* np.ones((md.mesh.numberofelements, 1))
347+md.cluster = generic('name', oshostname(), 'np', 3)
348+md = solve(md, 'Transient')
349+
350+#Fields and tolerances to track changes
351+field_names = [
352+ 'Vx1', 'Vy1', 'Vel1', 'Pressure1', 'Bed1', 'Surface1', 'Thickness1',
353+ 'Vx2', 'Vy2', 'Vel2', 'Pressure2', 'Bed2', 'Surface2', 'Thickness2',
354+ 'Vx3', 'Vy3', 'Vel3', 'Pressure3', 'Bed3', 'Surface3', 'Thickness3'
355+]
356+field_tolerances = [
357+ 2e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09,
358+ 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09,
359+ 2e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09, 1e-09
360+]
361+field_values = [
362+ md.results.TransientSolution[0].Vx,
363+ md.results.TransientSolution[0].Vy,
364+ md.results.TransientSolution[0].Vel,
365+ md.results.TransientSolution[0].Pressure,
366+ md.results.TransientSolution[0].Base,
367+ md.results.TransientSolution[0].Surface,
368+ md.results.TransientSolution[0].Thickness,
369+ md.results.TransientSolution[1].Vx,
370+ md.results.TransientSolution[1].Vy,
371+ md.results.TransientSolution[1].Vel,
372+ md.results.TransientSolution[1].Pressure,
373+ md.results.TransientSolution[1].Base,
374+ md.results.TransientSolution[1].Surface,
375+ md.results.TransientSolution[1].Thickness,
376+ md.results.TransientSolution[2].Vx,
377+ md.results.TransientSolution[2].Vy,
378+ md.results.TransientSolution[2].Vel,
379+ md.results.TransientSolution[2].Pressure,
380+ md.results.TransientSolution[2].Base,
381+ md.results.TransientSolution[2].Surface,
382+ md.results.TransientSolution[2].Thickness
383+]
384Index: ../trunk-jpl/test/NightlyRun/test350.py
385===================================================================
386--- ../trunk-jpl/test/NightlyRun/test350.py (revision 25805)
387+++ ../trunk-jpl/test/NightlyRun/test350.py (revision 25806)
388@@ -1,15 +1,16 @@
389 #Test Name: SquareSheetHydrologyShakti
390-from operator import itemgetter
391 import numpy as np
392+
393+from frictionshakti import frictionshakti
394+from MatlabFuncs import oshostname
395 from model import *
396-from socket import gethostname
397-from triangle import triangle
398-from setmask import setmask
399+from operator import itemgetter
400 from parameterize import parameterize
401 from setflowequation import setflowequation
402+from setmask import setmask
403 from solve import solve
404-from frictionshakti import frictionshakti
405 from transient import transient
406+from triangle import triangle
407
408 md = triangle(model(), '../Exp/Square.exp', 50000.)
409 md.mesh.x = md.mesh.x / 1000
410Index: ../trunk-jpl/test/NightlyRun/runme.m
411===================================================================
412--- ../trunk-jpl/test/NightlyRun/runme.m (revision 25805)
413+++ ../trunk-jpl/test/NightlyRun/runme.m (revision 25806)
414@@ -106,7 +106,7 @@
415 % }}}
416 %Process Ids according to benchmarks{{{
417 if strcmpi(benchmark,'nightly'),
418- test_ids=intersect(test_ids,[1:999]);
419+ test_ids=intersect(test_ids,union([1:999]));
420 elseif strcmpi(benchmark,'validation'),
421 test_ids=intersect(test_ids,[1001:1999]);
422 elseif strcmpi(benchmark,'ismip'),
423Index: ../trunk-jpl/test/NightlyRun/test444.py
424===================================================================
425--- ../trunk-jpl/test/NightlyRun/test444.py (revision 25805)
426+++ ../trunk-jpl/test/NightlyRun/test444.py (revision 25806)
427@@ -6,16 +6,18 @@
428 #
429
430 import numpy as np
431+
432+from ContourToMesh import *
433+from dmeth_params_set import *
434+from MatlabFuncs import oshostname
435 from model import *
436-from triangle import *
437-from setmask import *
438 from parameterize import *
439+from partitioner import *
440+from regionaloutput import *
441 from setflowequation import *
442+from setmask import *
443 from solve import *
444-from partitioner import *
445-from dmeth_params_set import *
446-from ContourToMesh import *
447-from regionaloutput import *
448+from triangle import *
449
450 #model not consistent: equality thickness = surface-base violated
451
Note: See TracBrowser for help on using the repository browser.