source: issm/oecreview/Archive/25834-26739/ISSM-26300-26301.diff

Last change on this file was 26740, checked in by Mathieu Morlighem, 3 years ago

CHG: added 25834-26739
File size: 167.9 KB
RevLine 
[26740]1Index: ../trunk-jpl/src/m/boundaryconditions/getlovenumbers.py
2===================================================================
3--- ../trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 26300)
4+++ ../trunk-jpl/src/m/boundaryconditions/getlovenumbers.py (revision 26301)
5@@ -4,23 +4,27 @@
6
7
8 def getlovenumbers(*args): #{{{
9- '''
10- GETLOVENUMBERS - provide love numbers retrieved from:
11+ """GETLOVENUMBERS - provide love numbers retrieved from:
12 http://www.srosat.com/iag-jsg/loveNb.php in a chosen reference frame
13
14- Usage:
15- series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
16+ Usage:
17+ series = love_numbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', 1001)
18
19- - type = one of 'loadingverticaldisplacement',
20- 'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
21- 'tidalverticaldisplacement', 'tidalgravitationalpotential',
22- 'tidalhorizontaldisplacement'
23- - reference_frame = one of 'CM' (default) and 'CF'
24- - maxdeg = default 1000
25+ - type = one of 'loadingverticaldisplacement',
26+ 'loadinggravitationalpotential', 'loadinghorizontaldisplacement',
27+ 'tidalverticaldisplacement', 'tidalgravitationalpotential',
28+ 'tidalhorizontaldisplacement'
29+ - reference_frame = one of 'CM' (default) and 'CF'
30+ - maxdeg = default 1000
31
32- Example:
33- h = love_number
34- '''
35+ Example:
36+ h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
37+ k = getlovenumbers('type', 'loadinggravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
38+ l = getlovenumbers('type', 'loadinghorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
39+ th = getlovenumbers('type', 'tidalverticaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
40+ tk = getlovenumbers('type', 'tidalgravitationalpotential', 'referenceframe', 'CM', 'maxdeg', maxdeg)
41+ tl = getlovenumbers('type', 'tidalhorizontaldisplacement', 'referenceframe', 'CM', 'maxdeg', maxdeg)
42+ """
43
44 TYPES = [
45 'loadingverticaldisplacement',
46@@ -31,12 +35,15 @@
47 'tidalhorizontaldisplacement'
48 ]
49
50- #recover options
51+ # Recover options
52 options = pairoptions(*args)
53 type = options.getfieldvalue('type')
54 frame = options.getfieldvalue('referenceframe', 'CM')
55 maxdeg = options.getfieldvalue('maxdeg', 1000)
56
57+ if (maxdeg > 10000):
58+ raise Exception('PREM love numbers computed only for deg < 10,000. Request lower maxdeg')
59+
60 if type not in TYPES:
61 raise Exception('type should be one of {}'.format(TYPES))
62
63@@ -10044,10 +10051,10 @@
64 [-6.27342778, -0.00030945, 0.00018956, 0.00031100, -0.00000116, 0.00000000]
65 ])
66
67- #cut off series at maxdeg
68+ # Cut off series at maxdeg
69 love_numbers = np.delete(love_numbers, range(maxdeg + 1, len(love_numbers)), axis=0)
70
71- #retrieve right type
72+ # Retrieve right type
73 if type == 'loadingverticaldisplacement':
74 series = love_numbers[:, 0]
75 elif type == 'loadinggravitationalpotential':
76@@ -10061,7 +10068,7 @@
77 elif type == 'tidalhorizontaldisplacement':
78 series = love_numbers[:, 5]
79 else:
80- raise Exception("love_numbers error message: unknown type: {}".format(type))
81+ raise Exception('love_numbers error message: unknown type: {}'.format(type))
82
83 #choose degree 1 term for CF reference system
84 if frame == 'CM':
85@@ -10074,7 +10081,7 @@
86 elif type == 'loadinghorizontaldisplacement':
87 series[1] = 0.134
88 else:
89- raise Exception("love_numbers error message: unknown reference frame: {}".format(frame))
90+ raise Exception('love_numbers error message: unknown reference frame: {}'.format(frame))
91
92 return series
93 #}}}
94Index: ../trunk-jpl/src/m/classes/amr.m
95===================================================================
96--- ../trunk-jpl/src/m/classes/amr.m (revision 26300)
97+++ ../trunk-jpl/src/m/classes/amr.m (revision 26301)
98@@ -83,7 +83,7 @@
99 %control of element lengths
100 self.gradation=1.5;
101
102- %other criterias
103+ %other criteria
104 self.groundingline_resolution=500.;
105 self.groundingline_distance=0.;
106 self.icefront_resolution=500.;
107Index: ../trunk-jpl/src/m/classes/amr.py
108===================================================================
109--- ../trunk-jpl/src/m/classes/amr.py (revision 26300)
110+++ ../trunk-jpl/src/m/classes/amr.py (revision 26301)
111@@ -4,84 +4,93 @@
112
113
114 class amr(object):
115- """
116- AMR Class definition
117+ """AMR Class definition
118
119 Usage:
120 amr = amr()
121 """
122
123- def __init__(self): # {{{
124- self.hmin = 0.
125- self.hmax = 0.
126+ def __init__(self): #{{{
127+ self.hmin = 0
128+ self.hmax = 0
129 self.fieldname = ''
130- self.err = 0.
131- self.keepmetric = 0.
132- self.gradation = 0.
133- self.groundingline_resolution = 0.
134- self.groundingline_distance = 0.
135- self.icefront_resolution = 0.
136- self.icefront_distance = 0.
137- self.thicknesserror_resolution = 0.
138- self.thicknesserror_threshold = 0.
139- self.thicknesserror_groupthreshold = 0.
140- self.thicknesserror_maximum = 0.
141- self.deviatoricerror_resolution = 0.
142- self.deviatoricerror_threshold = 0.
143- self.deviatoricerror_groupthreshold = 0.
144- self.deviatoricerror_maximum = 0.
145- self.restart = 0.
146- #set defaults
147+ self.err = 0
148+ self.keepmetric = 0
149+ self.gradation = 0
150+ self.groundingline_resolution = 0
151+ self.groundingline_distance = 0
152+ self.icefront_resolution = 0
153+ self.icefront_distance = 0
154+ self.thicknesserror_resolution = 0
155+ self.thicknesserror_threshold = 0
156+ self.thicknesserror_groupthreshold = 0
157+ self.thicknesserror_maximum = 0
158+ self.deviatoricerror_resolution = 0
159+ self.deviatoricerror_threshold = 0
160+ self.deviatoricerror_groupthreshold = 0
161+ self.deviatoricerror_maximum = 0
162+ self.restart = 0
163+
164 self.setdefaultparameters()
165 #}}}
166
167- def __repr__(self): # {{{
168- string = " amr parameters:"
169- string = "%s\n%s" % (string, fielddisplay(self, "hmin", "minimum element length"))
170- string = "%s\n%s" % (string, fielddisplay(self, "hmax", "maximum element length"))
171- string = "%s\n%s" % (string, fielddisplay(self, "fieldname", "name of input that will be used to compute the metric (should be an input of FemModel)"))
172- string = "%s\n%s" % (string, fielddisplay(self, "keepmetric", "indicates whether the metric should be kept every remeshing time"))
173- string = "%s\n%s" % (string, fielddisplay(self, "gradation", "maximum ratio between two adjacent edges"))
174- string = "%s\n%s" % (string, fielddisplay(self, "groundingline_resolution", "element length near the grounding line"))
175- string = "%s\n%s" % (string, fielddisplay(self, "groundingline_distance", "distance around the grounding line which elements will be refined"))
176- string = "%s\n%s" % (string, fielddisplay(self, "icefront_resolution", "element length near the ice front"))
177- string = "%s\n%s" % (string, fielddisplay(self, "icefront_distance", "distance around the ice front which elements will be refined"))
178- string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_resolution", "element length when thickness error estimator is used"))
179- string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_threshold", "maximum threshold thickness error permitted"))
180- string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_groupthreshold", "maximum group threshold thickness error permitted"))
181- string = "%s\n%s" % (string, fielddisplay(self, "thicknesserror_maximum", "maximum thickness error permitted"))
182- string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_resolution", "element length when deviatoric stress error estimator is used"))
183- string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_threshold", "maximum threshold deviatoricstress error permitted"))
184- string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_groupthreshold", "maximum group threshold deviatoric stress error permitted"))
185- string = "%s\n%s" % (string, fielddisplay(self, "deviatoricerror_maximum", "maximum deviatoricstress error permitted"))
186- string = "%s\n%s" % (string, fielddisplay(self, "restart", "indicates if ReMesh() will call before first time step"))
187- return string
188+ def __repr__(self): #{{{
189+ s = ' amr parameters:\n'
190+ s += '{}\n'.format(fielddisplay(self, 'hmin', 'minimum element length'))
191+ s += '{}\n'.format(fielddisplay(self, 'hmax', 'maximum element length'))
192+ s += '{}\n'.format(fielddisplay(self, 'fieldname', 'name of input that will be used to compute the metric (should be an input of FemModel)'))
193+ s += '{}\n'.format(fielddisplay(self, 'keepmetric', 'indicates whether the metric should be kept every remeshing time'))
194+ s += '{}\n'.format(fielddisplay(self, 'gradation', 'maximum ratio between two adjacent edges'))
195+ s += '{}\n'.format(fielddisplay(self, 'groundingline_resolution', 'element length near the grounding line'))
196+ s += '{}\n'.format(fielddisplay(self, 'groundingline_distance', 'distance around the grounding line which elements will be refined'))
197+ s += '{}\n'.format(fielddisplay(self, 'icefront_resolution', 'element length near the ice front'))
198+ s += '{}\n'.format(fielddisplay(self, 'icefront_distance', 'distance around the ice front which elements will be refined'))
199+ s += '{}\n'.format(fielddisplay(self, 'thicknesserror_resolution', 'element length when thickness error estimator is used'))
200+ s += '{}\n'.format(fielddisplay(self, 'thicknesserror_threshold', 'maximum threshold thickness error permitted'))
201+ s += '{}\n'.format(fielddisplay(self, 'thicknesserror_groupthreshold', 'maximum group threshold thickness error permitted'))
202+ s += '{}\n'.format(fielddisplay(self, 'thicknesserror_maximum', 'maximum thickness error permitted'))
203+ s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_resolution', 'element length when deviatoric stress error estimator is used'))
204+ s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_threshold', 'maximum threshold deviatoricstress error permitted'))
205+ s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_groupthreshold', 'maximum group threshold deviatoric stress error permitted'))
206+ s += '{}\n'.format(fielddisplay(self, 'deviatoricerror_maximum', 'maximum deviatoricstress error permitted'))
207+ s += '{}\n'.format(fielddisplay(self, 'restart', 'indicates if ReMesh() will call before first time step'))
208+ return s
209 #}}}
210
211- def setdefaultparameters(self): # {{{
212- self.hmin = 100.
213- self.hmax = 100.e3
214+ def setdefaultparameters(self): #{{{
215+ self.hmin = 100
216+ self.hmax = 100e3
217+
218+ # Fields
219 self.fieldname = 'Vel'
220- self.err = 3.
221+ self.err = 3
222+
223+ # Keep metric?
224 self.keepmetric = 1
225+
226+ # Control of element lengths
227 self.gradation = 1.5
228- self.groundingline_resolution = 500.
229+
230+ # Other criteria
231+ self.groundingline_resolution = 500
232 self.groundingline_distance = 0
233- self.icefront_resolution = 500.
234+ self.icefront_resolution = 500
235 self.icefront_distance = 0
236- self.thicknesserror_resolution = 500.
237+ self.thicknesserror_resolution = 500
238 self.thicknesserror_threshold = 0
239 self.thicknesserror_groupthreshold = 0
240 self.thicknesserror_maximum = 0
241- self.deviatoricerror_resolution = 500.
242+ self.deviatoricerror_resolution = 500
243 self.deviatoricerror_threshold = 0
244 self.deviatoricerror_groupthreshold = 0
245 self.deviatoricerror_maximum = 0
246- self.restart = 0.
247+
248+ # Is restart? This calls femmodel->ReMesh() before first time step.
249+ self.restart = 0
250 return self
251 #}}}
252
253- def checkconsistency(self, md, solution, analyses): # {{{
254+ def checkconsistency(self, md, solution, analyses): #{{{
255 md = checkfield(md, 'fieldname', 'amr.hmax', 'numel', [1], '>', 0, 'NaN', 1)
256 md = checkfield(md, 'fieldname', 'amr.hmin', 'numel', [1], '>', 0, '<', self.hmax, 'NaN', 1)
257 md = checkfield(md, 'fieldname', 'amr.keepmetric', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1)
258@@ -100,9 +109,9 @@
259 md = checkfield(md, 'fieldname', 'amr.deviatoricerror_maximum', 'numel', [1], '>=', 0, 'NaN', 1, 'Inf', 1)
260 md = checkfield(md, 'fieldname', 'amr.restart', 'numel', [1], '>=', 0, '<=', 1, 'NaN', 1)
261 return md
262- # }}}
263+ # }}}
264
265- def marshall(self, prefix, md, fid): # {{{
266+ def marshall(self, prefix, md, fid): #{{{
267 WriteData(fid, prefix, 'name', 'md.amr.type', 'data', 1, 'format', 'Integer')
268 WriteData(fid, prefix, 'object', self, 'fieldname', 'hmin', 'format', 'Double')
269 WriteData(fid, prefix, 'object', self, 'fieldname', 'hmax', 'format', 'Double')
270@@ -123,4 +132,4 @@
271 WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_groupthreshold', 'format', 'Double')
272 WriteData(fid, prefix, 'object', self, 'fieldname', 'deviatoricerror_maximum', 'format', 'Double')
273 WriteData(fid, prefix, 'object', self, 'class', 'amr', 'fieldname', 'restart', 'format', 'Integer')
274- # }}}
275+ #}}}
276Index: ../trunk-jpl/src/m/classes/dsl.m
277===================================================================
278--- ../trunk-jpl/src/m/classes/dsl.m (revision 26300)
279+++ ../trunk-jpl/src/m/classes/dsl.m (revision 26301)
280@@ -12,10 +12,6 @@
281
282 end
283 methods
284- function self = extrude(self,md) % {{{
285- self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
286- self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
287- end % }}}
288 function self = dsl(varargin) % {{{
289 switch nargin
290 case 0
291@@ -33,6 +29,14 @@
292 self.sea_water_pressure_at_sea_floor=NaN;
293
294 end % }}}
295+ function disp(self) % {{{
296+
297+ disp(sprintf(' dsl parameters:'));
298+ fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
299+ fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
300+ fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).');
301+
302+ end % }}}
303 function md = checkconsistency(self,md,solution,analyses) % {{{
304
305 %Early return
306@@ -48,28 +52,17 @@
307 end
308
309 end % }}}
310- function disp(self) % {{{
311-
312- disp(sprintf(' dsl parameters:'));
313- fielddisplay(self,'global_average_thermosteric_sea_level','Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).');
314- fielddisplay(self,'sea_surface_height_above_geoid','Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).');
315- fielddisplay(self,'sea_water_pressure_at_sea_floor','Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).');
316-
317- end % }}}
318 function marshall(self,prefix,md,fid) % {{{
319-
320+ yts=md.constants.yts;
321 WriteData(fid,prefix,'name','md.dsl.model','data',1,'format','Integer');
322- WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',md.constants.yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
323- WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify DSL at vertex locations.
324- WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts); %mattype 1 because we specify bottom pressure at vertex locations.
325+ WriteData(fid,prefix,'object',self,'fieldname','global_average_thermosteric_sea_level','format','DoubleMat','mattype',2,'timeseries',1,'timeserieslength',2,'yts',yts); %mattype 2, because we are sending a GMSL value identical everywhere on each element.
326+ WriteData(fid,prefix,'object',self,'fieldname','sea_surface_height_above_geoid','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify DSL at vertex locations.
327+ WriteData(fid,prefix,'object',self,'fieldname','sea_water_pressure_at_sea_floor','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts); %mattype 1 because we specify bottom pressure at vertex locations.
328
329 end % }}}
330- function savemodeljs(self,fid,modelname) % {{{
331-
332- writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
333- writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
334- writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
335-
336+ function self = extrude(self,md) % {{{
337+ self.sea_surface_height_above_geoid=project3d(md,'vector',self.sea_surface_height_above_geoid,'type','node','layer',1);
338+ self.sea_water_pressure_at_sea_floor=project3d(md,'vector',self.sea_water_pressure_at_sea_floor,'type','node','layer',1);
339 end % }}}
340 function self = initialize(self,md) % {{{
341
342@@ -86,6 +79,12 @@
343 disp(' no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero');
344 end
345 end % }}}
346-
347+ function savemodeljs(self,fid,modelname) % {{{
348+
349+ writejs1Darray(fid,[modelname '.dsl.global_average_thermosteric_sea_level'],self.global_average_thermosteric_sea_level);
350+ writejs1Darray(fid,[modelname '.dsl.sea_surface_height_above_geoid'],self.sea_surface_height_above_geoid);
351+ writejs1Darray(fid,[modelname '.dsl.sea_water_pressure_at_sea_floor'],self.sea_water_pressure_at_sea_floor);
352+
353+ end % }}}
354 end
355 end
356Index: ../trunk-jpl/src/m/classes/fourierlove.m
357===================================================================
358--- ../trunk-jpl/src/m/classes/fourierlove.m (revision 26300)
359+++ ../trunk-jpl/src/m/classes/fourierlove.m (revision 26301)
360@@ -5,25 +5,25 @@
361
362 classdef fourierlove
363 properties (SetAccess=public)
364- nfreq = 0;
365- frequencies = 0;
366- sh_nmax = 0;
367- sh_nmin = 0;
368- g0 = 0;
369- r0 = 0;
370- mu0 = 0;
371- Gravitational_Constant = 0;
372- allow_layer_deletion = 0;
373- underflow_tol = 0;
374- integration_steps_per_layer= 0;
375- istemporal = 0;
376- n_temporal_iterations = 0;
377- time = 0;
378- love_kernels = 0;
379- forcing_type = 0;
380- inner_core_boundary = 0;
381- core_mantle_boundary = 0;
382- complex_computation = 0;
383+ nfreq = 0;
384+ frequencies = 0;
385+ sh_nmax = 0;
386+ sh_nmin = 0;
387+ g0 = 0;
388+ r0 = 0;
389+ mu0 = 0;
390+ Gravitational_Constant = 0;
391+ allow_layer_deletion = 0;
392+ underflow_tol = 0;
393+ integration_steps_per_layer = 0;
394+ istemporal = 0;
395+ n_temporal_iterations = 0;
396+ time = 0;
397+ love_kernels = 0;
398+ forcing_type = 0;
399+ inner_core_boundary = 0;
400+ core_mantle_boundary = 0;
401+ complex_computation = 0;
402
403 end
404 methods (Static)
405@@ -33,8 +33,6 @@
406 end% }}}
407 end
408 methods
409- function self = extrude(self,md) % {{{
410- end % }}}
411 function self = fourierlove(varargin) % {{{
412 switch nargin
413 case 0
414@@ -52,7 +50,7 @@
415 % work on matlab script for computing g0 for given Earth's structure.
416 self.g0=9.81; % m/s^2;
417 self.r0=6371*1e3; %m;
418- self.mu0=10^11; % Pa
419+ self.mu0=1e11; % Pa
420 self.Gravitational_Constant=6.67259e-11; % m^3 kg^-1 s^-2
421 self.allow_layer_deletion=1;
422 self.underflow_tol=1e-16; %threshold of deep to surface love number ratio to trigger the deletion of layer
423@@ -67,24 +65,24 @@
424 self.complex_computation=0;
425 end % }}}
426 function disp(self) % {{{
427- fielddisplay(self,'nfreq','number of frequencies sampled (default 1, elastic) [Hz]');
428+ fielddisplay(self,'nfreq','number of frequencies sampled (default: 1, elastic) [Hz]');
429 fielddisplay(self,'frequencies','frequencies sampled (convention defaults to 0 for the elastic case) [Hz]');
430- fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)');
431- fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default 1)');
432- fielddisplay(self,'g0','adimensioning constant for gravity (default 10) [m/s^2]');
433- fielddisplay(self,'r0','adimensioning constant for radius (default 6371*10^3) [m]');
434- fielddisplay(self,'mu0','adimensioning constant for stress (default 10^11) [Pa]');
435- fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default 6.67259e-11 [m^3 kg^-1 s^-2])');
436- fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)');
437- fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default 2.2204460492503131E-016)');
438- fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default 100)');
439- fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
440- fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default 8)');
441- fielddisplay(self,'time','time vector for deformation if istemporal (default 0) [s]');
442- fielddisplay(self,'love_kernels','compute love numbers at depth? (default 0)');
443- fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
444- fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default 1)');
445- fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default 2)');
446+ fielddisplay(self,'sh_nmax','maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)');
447+ fielddisplay(self,'sh_nmin','minimum spherical harmonic degree (default: 1)');
448+ fielddisplay(self,'g0','adimensioning constant for gravity (default: 10) [m/s^2]');
449+ fielddisplay(self,'r0','adimensioning constant for radius (default: 6371*10^3) [m]');
450+ fielddisplay(self,'mu0','adimensioning constant for stress (default: 10^11) [Pa]');
451+ fielddisplay(self,'Gravitational_Constant','Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])');
452+ fielddisplay(self,'allow_layer_deletion','allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)');
453+ fielddisplay(self,'underflow_tol','threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)');
454+ fielddisplay(self,'integration_steps_per_layer','number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)');
455+ fielddisplay(self,'istemporal',{'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'});
456+ fielddisplay(self,'n_temporal_iterations','max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)');
457+ fielddisplay(self,'time','time vector for deformation if istemporal (default: 0) [s]');
458+ fielddisplay(self,'love_kernels','compute love numbers at depth? (default: 0)');
459+ fielddisplay(self,'forcing_type',{'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):','1: Inner core boundary -- Volumic Potential','2: Inner core boundary -- Pressure','3: Inner core boundary -- Loading','4: Inner core boundary -- Tangential traction','5: Core mantle boundary -- Volumic Potential','6: Core mantle boundary -- Pressure','7: Core mantle boundary -- Loading','8: Core mantle boundary -- Tangential traction','9: Surface -- Volumic Potential','10: Surface -- Pressure','11: Surface -- Loading','12: Surface -- Tangential traction '});
460+ fielddisplay(self,'inner_core_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)');
461+ fielddisplay(self,'core_mantle_boundary','interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)');
462
463 end % }}}
464 function md = checkconsistency(self,md,solution,analyses) % {{{
465@@ -111,8 +109,8 @@
466 md = checkfield(md,'fieldname','love.n_temporal_iterations','NaN',1,'Inf',1,'numel',1,'>',0);
467 md = checkfield(md,'fieldname','love.time','NaN',1,'Inf',1,'numel',md.love.nfreq/2/md.love.n_temporal_iterations);
468 end
469- if md.love.sh_nmin<=1 & (md.love.forcing_type==9 || md.love.forcing_type==5 || md.love.forcing_type==1)
470- error('Degree 1 not supported for Volumetric Potential forcing. Use sh_min>=2 for this kind of calculation.')
471+ if md.love.sh_nmin<=1 & (md.love.forcing_type==1 || md.love.forcing_type==5 || md.love.forcing_type==9)
472+ error(['Degree 1 not supported for forcing type ' num2str(md.love.forcing_type) '. Use sh_min>=2 for this kind of calculation.'])
473 end
474
475 %need 'litho' material:
476@@ -129,7 +127,7 @@
477
478 end % }}}
479 function marshall(self,prefix,md,fid) % {{{
480-
481+
482 WriteData(fid,prefix,'object',self,'fieldname','nfreq','format','Integer');
483 WriteData(fid,prefix,'object',self,'fieldname','frequencies','format','DoubleMat','mattype',3);
484 WriteData(fid,prefix,'object',self,'fieldname','sh_nmax','format','Integer');
485@@ -151,6 +149,8 @@
486 WriteData(fid,prefix,'object',self,'fieldname','core_mantle_boundary','format','Integer');
487
488 end % }}}
489+ function self = extrude(self,md) % {{{
490+ end % }}}
491 function savemodeljs(self,fid,modelname) % {{{
492 error('not implemented yet!');
493 end % }}}
494Index: ../trunk-jpl/src/m/classes/fourierlove.py
495===================================================================
496--- ../trunk-jpl/src/m/classes/fourierlove.py (revision 26300)
497+++ ../trunk-jpl/src/m/classes/fourierlove.py (revision 26301)
498@@ -1,3 +1,5 @@
499+import numpy as np
500+
501 from fielddisplay import fielddisplay
502 from checkfield import checkfield
503 from WriteData import WriteData
504@@ -4,12 +6,13 @@
505
506
507 class fourierlove(object):
508- """FOURIERLOVE - Fourier Love Number class definition
509+ """FOURIERLOVE - class definition
510
511 Usage:
512- fourierlove = fourierlove()
513+ md.love = fourierlove()
514 """
515- def __init__(self): # {{{
516+
517+ def __init__(self): #{{{
518 self.nfreq = 0
519 self.frequencies = 0
520 self.sh_nmax = 0
521@@ -17,96 +20,155 @@
522 self.g0 = 0
523 self.r0 = 0
524 self.mu0 = 0
525+ self.Gravitational_Constant = 0
526 self.allow_layer_deletion = 0
527+ self.underflow_tol = 0
528+ self.integration_steps_per_layer = 0
529+ self.istemporal = 0
530+ self.n_temporal_iterations = 0
531+ self.time = 0
532 self.love_kernels = 0
533 self.forcing_type = 0
534+ self.inner_core_boundary = 0
535+ self.core_mantle_boundary = 0
536+ self.complex_computation = 0
537
538 self.setdefaultparameters()
539 #}}}
540
541- def __repr__(self): # {{{
542- # TODO:
543- # - Convert all formatting to calls to <string>.format (see any
544- # already converted <class>.__repr__ method for examples)
545- #
546- string = ' Fourier Love class:'
547- string = "%s\n%s" % (string, fielddisplay(self, 'nfreq', 'number of frequencies sampled (default 1, elastic) [Hz]'))
548- string = "%s\n%s" % (string, fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
549- string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default 256, .35 deg, or 40 km at equator)'))
550- string = "%s\n%s" % (string, fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default 1)'))
551- string = "%s\n%s" % (string, fielddisplay(self, 'g0', 'adimensioning constant for gravity (default 10) [m / s^2]'))
552- string = "%s\n%s" % (string, fielddisplay(self, 'r0', 'adimensioning constant for radius (default 6378e3) [m]'))
553- string = "%s\n%s" % (string, fielddisplay(self, 'mu0', 'adimensioning constant for stress (default 1.0e11) [Pa]'))
554- string = "%s\n%s" % (string, fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default 1)'))
555- string = "%s\n%s" % (string, fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default 0)'))
556- string = "%s\n%s" % (string, fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default 11) :'))
557- string = "%s\n%s" % (string, ' 1: Inner core boundary -- Volumic Potential')
558- string = "%s\n%s" % (string, ' 2: Inner core boundary -- Pressure')
559- string = "%s\n%s" % (string, ' 3: Inner core boundary -- Loading')
560- string = "%s\n%s" % (string, ' 4: Inner core boundary -- Tangential traction')
561- string = "%s\n%s" % (string, ' 5: Core mantle boundary -- Volumic Potential')
562- string = "%s\n%s" % (string, ' 6: Core mantle boundary -- Pressure')
563- string = "%s\n%s" % (string, ' 7: Core mantle boundary -- Loading')
564- string = "%s\n%s" % (string, ' 8: Core mantle boundary -- Tangential traction')
565- string = "%s\n%s" % (string, ' 9: Surface-- Volumic Potential')
566- string = "%s\n%s" % (string, ' 10: Surface-- Pressure')
567- string = "%s\n%s" % (string, ' 11: Surface-- Loading')
568- string = "%s\n%s" % (string, ' 12: Surface-- Tangential traction ')
569+ def __repr__(self): #{{{
570+ s = ' Fourier Love class:\n'
571+ s += '{}\n'.format(fielddisplay(self, 'nfreq', 'number of frequencies sampled (default: 1, elastic) [Hz]'))
572+ s += '{}\n'.format(fielddisplay(self, 'frequencies', 'frequencies sampled (convention defaults to 0 for the elastic case) [Hz]'))
573+ s += '{}\n'.format(fielddisplay(self, 'sh_nmax', 'maximum spherical harmonic degree (default: 256, .35 deg, or 40 km at equator)'))
574+ s += '{}\n'.format(fielddisplay(self, 'sh_nmin', 'minimum spherical harmonic degree (default: 1)'))
575+ s += '{}\n'.format(fielddisplay(self, 'g0', 'adimensioning constant for gravity (default: 10) [m / s^2]'))
576+ s += '{}\n'.format(fielddisplay(self, 'r0', 'adimensioning constant for radius (default: 6378e3) [m]'))
577+ s += '{}\n'.format(fielddisplay(self, 'mu0', 'adimensioning constant for stress (default: 1.0e11) [Pa]'))
578+ s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
579+ s += '{}\n'.format(fielddisplay(self, 'Gravitational_Constant', 'Newtonian constant of gravitation (default: 6.67259e-11 [m^3 kg^-1 s^-2])'))
580+ s += '{}\n'.format(fielddisplay(self, 'allow_layer_deletion', 'allow for migration of the integration boundary with increasing spherical harmonics degree (default: 1)'))
581+ s += '{}\n'.format(fielddisplay(self, 'underflow_tol', 'threshold of deep to surface love number ratio to trigger the deletion of layers (default: 1e-16)'))
582+ s += '{}\n'.format(fielddisplay(self, 'integration_steps_per_layer', 'number of radial steps to propagate the yi system from the bottom to the top of each layer (default: 100)'))
583+ s += '{}\n'.format(fielddisplay(self, 'istemporal', {'1 for time-dependent love numbers, 0 for frequency-dependent or elastic love numbers (default: 0)', 'If 1: use fourierlove function build_frequencies_from_time to meet consistency'}))
584+ s += '{}\n'.format(fielddisplay(self, 'n_temporal_iterations', 'max number of iterations in the inverse Laplace transform. Also the number of spectral samples per time step requested (default: 8)'))
585+ s += '{}\n'.format(fielddisplay(self, 'time', 'time vector for deformation if istemporal (default: 0) [s]'))
586+ s += '{}\n'.format(fielddisplay(self, 'love_kernels', 'compute love numbers at depth? (default: 0)'))
587+ s += '{}\n'.format(fielddisplay(self, 'forcing_type', 'integer indicating the nature and depth of the forcing for the Love number calculation (default: 11):'))
588+ s += '{}\n'.format(' 1: Inner core boundary -- Volumic Potential')
589+ s += '{}\n'.format(' 2: Inner core boundary -- Pressure')
590+ s += '{}\n'.format(' 3: Inner core boundary -- Loading')
591+ s += '{}\n'.format(' 4: Inner core boundary -- Tangential traction')
592+ s += '{}\n'.format(' 5: Core mantle boundary -- Volumic Potential')
593+ s += '{}\n'.format(' 6: Core mantle boundary -- Pressure')
594+ s += '{}\n'.format(' 7: Core mantle boundary -- Loading')
595+ s += '{}\n'.format(' 8: Core mantle boundary -- Tangential traction')
596+ s += '{}\n'.format(' 9: Surface-- Volumic Potential')
597+ s += '{}\n'.format(' 10: Surface-- Pressure')
598+ s += '{}\n'.format(' 11: Surface-- Loading')
599+ s += '{}\n'.format(' 12: Surface-- Tangential traction ')
600+ s += '{}\n'.format(fielddisplay(self, 'inner_core_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 1--4 (default: 1)'))
601+ s += '{}\n'.format(fielddisplay(self, 'core_mantle_boundary', 'interface index in materials.radius locating forcing. Only used for forcing_type 5--8 (default: 2)'))
602
603- return string
604+ return s
605 #}}}
606
607- def extrude(self, md): # {{{
608- return self
609- #}}}
610-
611- def setdefaultparameters(self): # {{{
612- #we setup an elastic love number computation by default.
613+ def setdefaultparameters(self): #{{{
614+ # We setup an elastic love number computation by default
615 self.nfreq = 1
616- self.frequencies = [0] #Hz
617- self.sh_nmax = 256 # .35 degree, 40 km at the equator.
618+ self.frequencies = [0] # Hz
619+ self.sh_nmax = 256 # .35 degree, 40 km at the equator
620 self.sh_nmin = 1
621+ # Work on matlab script for computing g0 for given Earth's structure
622 self.g0 = 9.81 # m/s^2
623- self.r0 = 6371e3 #m
624+ self.r0 = 6371 * 1e3 # m
625 self.mu0 = 1e11 # Pa
626+ self.Gravitational_Constant = 6.67259e-11 # m^3 kg^-1 s^-2
627 self.allow_layer_deletion = 1
628+ self.underflow_tol = 1e-16 # Threshold of deep to surface love number ratio to trigger the deletion of layer
629+ self.integration_steps_per_layer = 100
630+ self.istemporal = 0
631+ self.n_temporal_iterations = 8
632+ self.time = [0] # s
633 self.love_kernels = 0
634- self.forcing_type = 11
635-
636- return self
637+ self.forcing_type = 11 # Surface loading
638+ self.inner_core_boundary = 1
639+ self.core_mantle_boundary = 2
640+ self.complex_computation = 0
641 #}}}
642
643- def checkconsistency(self, md, solution, analyses): # {{{
644+ def checkconsistency(self, md, solution, analyses): #{{{
645 if 'LoveAnalysis' not in analyses:
646 return md
647- md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
648- md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', [md.love.nfreq])
649- md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
650- md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
651- md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
652- md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
653- md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0)
654+
655+ md = checkfield(md, 'fieldname', 'love.nfreq', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
656+ md = checkfield(md, 'fieldname', 'love.frequencies', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq)
657+ md = checkfield(md, 'fieldname', 'love.sh_nmax', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
658+ md = checkfield(md, 'fieldname', 'love.sh_nmin', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
659+ md = checkfield(md, 'fieldname', 'love.g0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
660+ md = checkfield(md, 'fieldname', 'love.r0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
661+ md = checkfield(md, 'fieldname', 'love.mu0', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
662+ md = checkfield(md, 'fieldname', 'love.Gravitational_Constant', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
663 md = checkfield(md, 'fieldname', 'love.allow_layer_deletion', 'values', [0, 1])
664+ md = checkfield(md, 'fieldname', 'love.underflow_tol', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
665+ md = checkfield(md, 'fieldname', 'love.integration_steps_per_layer', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
666 md = checkfield(md, 'fieldname', 'love.love_kernels', 'values', [0, 1])
667- md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', [1], '>', 0, '<=', 12)
668- if md.love.sh_nmin <= 1 and md.love.forcing_type == 9:
669- raise RuntimeError("Degree 1 not supported for Volumetric Potential forcing. Use sh_min >= 2 for this kind of calculation.")
670+ md = checkfield(md, 'fieldname', 'love.forcing_type', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', 12)
671+ md = checkfield(md, 'fieldname', 'love.complex_computation', 'NaN', 1, 'Inf', 1, 'numel', 1, 'values', [0, 1])
672
673- # need 'litho' material
674+ md = checkfield(md, 'fieldname', 'love.istemporal', 'values', [0, 1])
675+ if md.love.istemporal:
676+ md = checkfield(md, 'fieldname', 'love.n_temporal_iterations', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
677+ md = checkfield(md, 'fieldname', 'love.time', 'NaN', 1, 'Inf', 1, 'numel', md.love.nfreq / 2 / md.love.n_temporal_iterations)
678+ if md.love.sh_nmin <= 1 and (md.love.forcing_type == 1 or md.love.forcing_type == 5 or md.love.forcing_type == 9):
679+ raise RuntimeError('Degree 1 not supported for forcing type {}. Use sh_min >= 2 for this kind of calculation.'.format(md.love.forcing_type))
680+
681+ # Need 'litho' material
682 if not hasattr(md.materials, 'materials') or 'litho' not in md.materials.nature:
683 raise RuntimeError('Need a \'litho\' material to run a Fourier Love number analysis')
684+
685+ mat = np.where(np.array(nature) == 'litho')[0]
686+ if md.love.forcing_type <= 4:
687+ md = checkfield(md, 'fieldname', 'love.inner_core_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
688+ elif md.love.forcing_type <= 8:
689+ md = checkfield(md, 'fieldname', 'love.core_mantle_boundary', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0, '<=', md.materials[mat].numlayers)
690+
691 return md
692- # }}}
693+ #}}}
694
695- def marshall(self, prefix, md, fid): # {{{
696+ def marshall(self, prefix, md, fid): #{{{
697 WriteData(fid, prefix, 'object', self, 'fieldname', 'nfreq', 'format', 'Integer')
698- WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype', 3)
699+ WriteData(fid, prefix, 'object', self, 'fieldname', 'frequencies', 'format', 'DoubleMat', 'mattype',3)
700 WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmax', 'format', 'Integer')
701 WriteData(fid, prefix, 'object', self, 'fieldname', 'sh_nmin', 'format', 'Integer')
702 WriteData(fid, prefix, 'object', self, 'fieldname', 'g0', 'format', 'Double')
703 WriteData(fid, prefix, 'object', self, 'fieldname', 'r0', 'format', 'Double')
704 WriteData(fid, prefix, 'object', self, 'fieldname', 'mu0', 'format', 'Double')
705+ WriteData(fid, prefix, 'object', self, 'fieldname', 'Gravitational_Constant', 'format', 'Double')
706 WriteData(fid, prefix, 'object', self, 'fieldname', 'allow_layer_deletion', 'format', 'Boolean')
707+ WriteData(fid, prefix, 'object', self, 'fieldname', 'underflow_tol', 'format', 'Double')
708+ WriteData(fid, prefix, 'object', self, 'fieldname', 'integration_steps_per_layer', 'format', 'Integer')
709+ WriteData(fid, prefix, 'object', self, 'fieldname', 'istemporal', 'format', 'Boolean')
710+ WriteData(fid, prefix, 'object', self, 'fieldname', 'n_temporal_iterations', 'format', 'Integer')
711+ WriteData(fid, prefix, 'object', self, 'fieldname', 'complex_computation', 'format', 'Boolean')
712+ # Note: no need to marshall the time vector, we have frequencies
713 WriteData(fid, prefix, 'object', self, 'fieldname', 'love_kernels', 'format', 'Boolean')
714 WriteData(fid, prefix, 'object', self, 'fieldname', 'forcing_type', 'format', 'Integer')
715- # }}}
716+ WriteData(fid, prefix, 'object', self, 'fieldname', 'inner_core_boundary', 'format', 'Integer')
717+ WriteData(fid, prefix, 'object', self, 'fieldname', 'core_mantle_boundary', 'format', 'Integer')
718+ #}}}
719+
720+ def extrude(self, md): #{{{
721+ return self
722+ #}}}
723+
724+ def build_frequencies_from_time(self): #{{{
725+ if not self.istemporal:
726+ raise RuntimeError('cannot build frequencies for temporal love numbers if love.istemporal==0')
727+ print('Temporal love numbers: Overriding md.love.nfreq and md.love.frequencies')
728+ self.nfreq = len(self.time) * 2 * self.n_temporal_iterations
729+ self.frequencies = np.zeros((self.nfreq, 1))
730+ for i in range(len(self.time)):
731+ for j in range(2 * self.n_temporal_iterations):
732+ self.frequencies[(i - 1) * 2 * self.n_temporal_iterations + j] = j * np.log(2) / self.time[i] / 2 / np.pi
733+ #}}}
734Index: ../trunk-jpl/src/m/classes/geometry.py
735===================================================================
736--- ../trunk-jpl/src/m/classes/geometry.py (revision 26300)
737+++ ../trunk-jpl/src/m/classes/geometry.py (revision 26301)
738@@ -36,7 +36,7 @@
739 #}}}
740
741 def setdefaultparameters(self): #{{{
742- return self
743+ return
744 #}}}
745
746 def checkconsistency(self, md, solution, analyses): #{{{
747Index: ../trunk-jpl/src/m/classes/hydrologyshreve.py
748===================================================================
749--- ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26300)
750+++ ../trunk-jpl/src/m/classes/hydrologyshreve.py (revision 26301)
751@@ -1,3 +1,5 @@
752+import numpy as np
753+
754 from fielddisplay import fielddisplay
755 from checkfield import checkfield
756 from WriteData import WriteData
757@@ -10,43 +12,35 @@
758 hydrologyshreve = hydrologyshreve()
759 """
760
761- def __init__(self): # {{{
762- self.spcwatercolumn = float('NaN')
763+ def __init__(self, *args): #{{{
764+ self.spcwatercolumn = np.nan
765 self.stabilization = 0
766 self.requested_outputs = []
767
768- self.setdefaultparameters()
769+ nargs = len(args)
770+ if nargs == 0:
771+ self.setdefaultparameters()
772+ elif nargs == 1:
773+ self.setdefaultparameters(args)
774+ else:
775+ raise RuntimeError('constructor not supported')
776 #}}}
777
778- def __repr__(self): # {{{
779- # TODO:
780- # - Convert all formatting to calls to <string>.format (see any
781- # already converted <class>.__repr__ method for examples)
782- #
783- string = ' hydrologyshreve solution parameters:'
784- string = "%s\n%s" % (string, fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
785- string = "%s\n%s" % (string, fielddisplay(self, 'stabilization', 'artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.'))
786- string = "%s\n%s" % (string, fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
787- return string
788+ def __repr__(self): #{{{
789+ s = ' hydrologyshreve solution parameters:\n'
790+ s += '{}\n'.format(fielddisplay(self, 'spcwatercolumn', 'water thickness constraints (NaN means no constraint) [m]'))
791+ s += '{}\n'.format(fielddisplay(self, 'stabilization', 'artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.'))
792+ s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested'))
793+ return s
794 #}}}
795
796- def extrude(self, md): # {{{
797- return self
798- #}}}
799-
800- def setdefaultparameters(self): # {{{
801- #Type of stabilization to use 0:nothing 1:artificial_diffusivity
802+ def setdefaultparameters(self): #{{{
803+ # Type of stabilization to use 0:nothing 1:artificial_diffusivity
804 self.stabilization = 1
805 self.requested_outputs = ['default']
806- return self
807 #}}}
808
809- def defaultoutputs(self, md): # {{{
810- list = ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
811- return list
812- #}}}
813-
814- def checkconsistency(self, md, solution, analyses): # {{{
815+ def checkconsistency(self, md, solution, analyses): #{{{
816 #Early return
817 if 'HydrologyShreveAnalysis' not in analyses or (solution == 'TransientSolution' and not md.transient.ishydrology):
818 return md
819@@ -53,20 +47,25 @@
820
821 md = checkfield(md, 'fieldname', 'hydrology.spcwatercolumn', 'Inf', 1, 'timeseries', 1)
822 md = checkfield(md, 'fieldname', 'hydrology.stabilization', '>=', 0)
823- md = checkfield(md, 'fieldname', 'hydrology.requested_outputs', 'stringrow', 1)
824 return md
825 # }}}
826
827- def marshall(self, prefix, md, fid): # {{{
828+ def defaultoutputs(self, md): #{{{
829+ return ['Watercolumn', 'HydrologyWaterVx', 'HydrologyWaterVy']
830+ #}}}
831+
832+ def marshall(self, prefix, md, fid): #{{{
833 WriteData(fid, prefix, 'name', 'md.hydrology.model', 'data', 2, 'format', 'Integer')
834 WriteData(fid, prefix, 'object', self, 'fieldname', 'spcwatercolumn', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
835 WriteData(fid, prefix, 'object', self, 'fieldname', 'stabilization', 'format', 'Double')
836- #process requested outputs
837 outputs = self.requested_outputs
838 indices = [i for i, x in enumerate(outputs) if x == 'default']
839- if len(indices) > 0:
840+ if len(indices):
841 outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
842 outputs = outputscopy
843 WriteData(fid, prefix, 'data', outputs, 'name', 'md.hydrology.requested_outputs', 'format', 'StringArray')
844+ # }}}
845
846- # }}}
847+ def extrude(self, md): #{{{
848+ return self
849+ #}}}
850Index: ../trunk-jpl/src/m/classes/initialization.m
851===================================================================
852--- ../trunk-jpl/src/m/classes/initialization.m (revision 26300)
853+++ ../trunk-jpl/src/m/classes/initialization.m (revision 26301)
854@@ -26,26 +26,6 @@
855 sample = NaN;
856 end
857 methods
858- function self = extrude(self,md) % {{{
859- self.vx=project3d(md,'vector',self.vx,'type','node');
860- self.vy=project3d(md,'vector',self.vy,'type','node');
861- self.vz=project3d(md,'vector',self.vz,'type','node');
862- self.vel=project3d(md,'vector',self.vel,'type','node');
863- self.temperature=project3d(md,'vector',self.temperature,'type','node');
864- self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
865- self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
866- self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
867- self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
868- self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
869- self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
870- self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
871- self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
872- self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
873- self.str=project3d(md,'vector',self.str,'type','node','layer',1);
874-
875- %Lithostatic pressure by default
876- self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
877- end % }}}
878 function self = initialization(varargin) % {{{
879 switch nargin
880 case 0
881@@ -200,6 +180,26 @@
882 WriteData(fid,prefix,'data',enthalpy,'format','DoubleMat','mattype',1,'name','md.initialization.enthalpy');
883 end
884 end % }}}
885+ function self = extrude(self,md) % {{{
886+ self.vx=project3d(md,'vector',self.vx,'type','node');
887+ self.vy=project3d(md,'vector',self.vy,'type','node');
888+ self.vz=project3d(md,'vector',self.vz,'type','node');
889+ self.vel=project3d(md,'vector',self.vel,'type','node');
890+ self.temperature=project3d(md,'vector',self.temperature,'type','node');
891+ self.enthalpy=project3d(md,'vector',self.enthalpy,'type','node');
892+ self.waterfraction=project3d(md,'vector',self.waterfraction,'type','node');
893+ self.watercolumn=project3d(md,'vector',self.watercolumn,'type','node','layer',1);
894+ self.sediment_head=project3d(md,'vector',self.sediment_head,'type','node','layer',1);
895+ self.epl_head=project3d(md,'vector',self.epl_head,'type','node','layer',1);
896+ self.epl_thickness=project3d(md,'vector',self.epl_thickness,'type','node','layer',1);
897+ self.sealevel=project3d(md,'vector',self.sealevel,'type','node','layer',1);
898+ self.bottompressure=project3d(md,'vector',self.bottompressure,'type','node','layer',1);
899+ self.dsl=project3d(md,'vector',self.dsl,'type','node','layer',1);
900+ self.str=project3d(md,'vector',self.str,'type','node','layer',1);
901+
902+ %Lithostatic pressure by default
903+ self.pressure=md.constants.g*md.materials.rho_ice*(md.geometry.surface-md.mesh.z);
904+ end % }}}
905 function savemodeljs(self,fid,modelname) % {{{
906
907 writejs1Darray(fid,[modelname '.initialization.vx'],self.vx);
908Index: ../trunk-jpl/src/m/classes/lovenumbers.m
909===================================================================
910--- ../trunk-jpl/src/m/classes/lovenumbers.m (revision 26300)
911+++ ../trunk-jpl/src/m/classes/lovenumbers.m (revision 26301)
912@@ -15,14 +15,14 @@
913 l = []; %idem
914
915 %tidal love numbers for computing rotational feedback:
916- th = [];
917- tk = [];
918- tl = [];
919- tk2secular = 0; %deg 2 secular number.
920+ th = [];
921+ tk = [];
922+ tl = [];
923+ tk2secular = 0; %deg 2 secular number.
924
925 %time/frequency for visco-elastic love numbers
926 timefreq = [];
927- istime = 1;
928+ istime = 1;
929
930 end
931 methods
932@@ -32,8 +32,24 @@
933 referenceframe=getfieldvalue(options,'referenceframe','CM');
934 self=setdefaultparameters(self,maxdeg,referenceframe);
935 end % }}}
936+ function disp(self) % {{{
937+ disp(sprintf(' lovenumbers parameters:'));
938+
939+ fielddisplay(self,'h','load Love number for radial displacement');
940+ fielddisplay(self,'k','load Love number for gravitational potential perturbation');
941+ fielddisplay(self,'l','load Love number for horizontal displacements');
942+
943+ fielddisplay(self,'th','tidal load Love number (deg 2)');
944+ fielddisplay(self,'tk','tidal load Love number (deg 2)');
945+ fielddisplay(self,'tl','tidal load Love number (deg 2)');
946+ fielddisplay(self,'tk2secular','secular fluid Love number');
947+
948+ fielddisplay(self,'istime','time (default: 1) or frequency love numbers (0)');
949+ fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
950+
951+ end % }}}
952 function self = setdefaultparameters(self,maxdeg,referenceframe) % {{{
953-
954+
955 %initialize love numbers:
956 self.h=getlovenumbers('type','loadingverticaldisplacement','referenceframe',referenceframe,'maxdeg',maxdeg);
957 self.k=getlovenumbers('type','loadinggravitationalpotential','referenceframe',referenceframe,'maxdeg',maxdeg);
958@@ -44,7 +60,7 @@
959
960 %secular fluid love number:
961 self.tk2secular=0.942;
962-
963+
964 %time:
965 self.istime=1; %temporal love numbers by default
966 self.timefreq=0; %elastic case by default.
967@@ -72,7 +88,7 @@
968 if (size(self.h,1)~=size(self.k,1) | size(self.h,1)~=size(self.l,1)),
969 error('lovenumbers error message: love numbers should be provided at the same level of accuracy');
970 end
971-
972+
973 ntf=length(self.timefreq);
974 if( size(self.h,2) ~= ntf | size(self.k,2) ~= ntf | size(self.l,2) ~= ntf | size(self.th,2) ~= ntf | size(self.tk,2) ~= ntf | size(self.tl,2) ~= ntf ),
975 error('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector');
976@@ -82,22 +98,6 @@
977 function list=defaultoutputs(self,md) % {{{
978 list = {};
979 end % }}}
980- function disp(self) % {{{
981- disp(sprintf(' lovenumbers parameters:'));
982-
983- fielddisplay(self,'h','load Love number for radial displacement');
984- fielddisplay(self,'k','load Love number for gravitational potential perturbation');
985- fielddisplay(self,'l','load Love number for horizontal displacements');
986-
987- fielddisplay(self,'th','tidal load Love number (deg 2)');
988- fielddisplay(self,'tk','tidal load Love number (deg 2)');
989- fielddisplay(self,'tl','tidal load Love number (deg 2)');
990- fielddisplay(self,'tk2secular','secular fluid Love number');
991-
992- fielddisplay(self,'istime','time (default=1) or frequency love numbers (0)');
993- fielddisplay(self,'timefreq','time/frequency vector (yr or 1/yr)');
994-
995- end % }}}
996 function marshall(self,prefix,md,fid) % {{{
997
998 WriteData(fid,prefix,'object',self,'fieldname','h','name','md.solidearth.lovenumbers.h','format','DoubleMat','mattype',1);
999Index: ../trunk-jpl/src/m/classes/materials.m
1000===================================================================
1001--- ../trunk-jpl/src/m/classes/materials.m (revision 26300)
1002+++ ../trunk-jpl/src/m/classes/materials.m (revision 26301)
1003@@ -203,7 +203,7 @@
1004 fielddisplay(self,'ebm_alpha','array describing each layer''s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)');
1005 fielddisplay(self,'ebm_delta','array describing each layer''s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)');
1006 fielddisplay(self,'ebm_taul','array describing each layer''s starting period for transient relaxation, only for EBM rheology (numlayers) [s]');
1007- fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]');
1008+ fielddisplay(self,'ebm_tauh','array describing each layer''s array describing each layer''s end period for transient relaxation, only for Burgers rheology (numlayers) [s]');
1009
1010
1011 fielddisplay(self,'rheologymodel','array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)');
1012@@ -255,7 +255,7 @@
1013 if md.materials.rheologymodel(i)==1 & (isnan(md.materials.burgers_viscosity(i) | isnan(md.materials.burgers_mu(i)))),
1014 error('materials checkconsistency error message: Litho burgers_viscosity or burgers_mu has NaN values, inconsistent with rheologymodel choice');
1015 end
1016- if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i)) ),
1017+ if md.materials.rheologymodel(i)==2 & (isnan(md.materials.ebm_alpha(i)) | isnan(md.materials.ebm_delta(i)) | isnan(md.materials.ebm_taul(i)) | isnan(md.materials.ebm_tauh(i))),
1018 error('materials checkconsistency error message: Litho ebm_alpha, ebm_delta, ebm_taul or ebm_tauh has NaN values, inconsistent with rheologymodel choice');
1019 end
1020 end
1021@@ -283,7 +283,7 @@
1022
1023 %1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
1024 WriteData(fid,prefix,'name','md.materials.nature','data',naturetointeger(self.nature),'format','IntMat','mattype',3);
1025- WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER, this can evolve if you have classes.
1026+ WriteData(fid,prefix,'name','md.materials.type','data',5,'format','Integer'); %DANGER: this can evolve if you have classes.
1027 for i=1:length(self.nature),
1028 nat=self.nature{i};
1029 switch nat
1030@@ -325,6 +325,7 @@
1031 earth_density=earth_density + (self.radius(i+1)^3-self.radius(i)^3)*self.density(i);
1032 end
1033 earth_density=earth_density/self.radius(self.numlayers+1)^3;
1034+ disp(earth_density)
1035 self.earth_density=earth_density;
1036 case 'hydro'
1037 WriteData(fid,prefix,'object',self,'class','materials','fieldname','rho_ice','format','Double');
1038@@ -540,7 +541,7 @@
1039 t4 = s(i,4)/((ra^3)*6.);
1040 vs = vs + t1*(r2^3) + t2*(r2^4) + t3*(r2^5) + t4*(r2^6) - ( t1*(r1^3) + t2*(r1^4) + t3*(r1^5) + t4*(r1^6) );
1041
1042- end
1043+ end
1044 ro = ro*3 / (rad(j+1)^3-rad(j)^3);
1045 vp = vp*3 /(rad(j+1)^3-rad(j)^3);
1046 vs = vs*3 / (rad(j+1)^3-rad(j)^3);
1047@@ -556,7 +557,6 @@
1048 end
1049
1050 end % }}}
1051-
1052 end
1053 end
1054
1055@@ -583,5 +583,3 @@
1056 end
1057 end
1058 end % }}}
1059-
1060-
1061Index: ../trunk-jpl/src/m/classes/model.m
1062===================================================================
1063--- ../trunk-jpl/src/m/classes/model.m (revision 26300)
1064+++ ../trunk-jpl/src/m/classes/model.m (revision 26301)
1065@@ -38,10 +38,10 @@
1066 thermal = 0;
1067 steadystate = 0;
1068 transient = 0;
1069- levelset = 0;
1070+ levelset = 0;
1071 calving = 0;
1072 frontalforcings = 0;
1073- love = 0;
1074+ love = 0;
1075 esa = 0;
1076 sampling = 0;
1077
1078@@ -48,7 +48,7 @@
1079 autodiff = 0;
1080 inversion = 0;
1081 qmu = 0;
1082- amr = 0;
1083+ amr = 0;
1084 results = 0;
1085 outputdefinition = 0;
1086 radaroverlay = 0;
1087@@ -206,6 +206,50 @@
1088
1089 end
1090 %}}}
1091+ function disp(self) % {{{
1092+ disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
1093+ disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
1094+ disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
1095+ disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants'));
1096+ disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
1097+ disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
1098+ disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties'));
1099+ disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
1100+ disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
1101+ disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
1102+ disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
1103+ disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
1104+ disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
1105+ disp(sprintf('%19s: %-22s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
1106+ disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
1107+ disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
1108+ disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
1109+ disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties'));
1110+ disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
1111+ disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
1112+ disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
1113+ disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
1114+ disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
1115+ disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
1116+ disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
1117+ disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
1118+ disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
1119+ disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
1120+ disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
1121+ disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
1122+ disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
1123+ disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
1124+ disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
1125+ disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
1126+ disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
1127+ disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
1128+ disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
1129+ disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
1130+ disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
1131+ disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results'));
1132+ disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
1133+ disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
1134+ end % }}}
1135 function md = setdefaultparameters(md,planet) % {{{
1136
1137 %initialize subclasses
1138@@ -1132,7 +1176,7 @@
1139 if md.mesh.average_vertex_connectivity<=25,
1140 md.mesh.average_vertex_connectivity=100;
1141 end
1142- end % }}}
1143+ end % }}}
1144 function md = structtomodel(md,structmd) % {{{
1145
1146 if ~isstruct(structmd) error('input model is not a structure'); end
1147@@ -1583,97 +1627,54 @@
1148 md.mesh.average_vertex_connectivity=100;
1149 end
1150 end % }}}
1151- function disp(self) % {{{
1152- disp(sprintf('%19s: %-22s -- %s','mesh' ,['[1x1 ' class(self.mesh) ']'],'mesh properties'));
1153- disp(sprintf('%19s: %-22s -- %s','mask' ,['[1x1 ' class(self.mask) ']'],'defines grounded and floating elements'));
1154- disp(sprintf('%19s: %-22s -- %s','geometry' ,['[1x1 ' class(self.geometry) ']'],'surface elevation, bedrock topography, ice thickness,...'));
1155- disp(sprintf('%19s: %-22s -- %s','constants' ,['[1x1 ' class(self.constants) ']'],'physical constants'));
1156- disp(sprintf('%19s: %-22s -- %s','smb' ,['[1x1 ' class(self.smb) ']'],'surface mass balance'));
1157- disp(sprintf('%19s: %-22s -- %s','basalforcings' ,['[1x1 ' class(self.basalforcings) ']'],'bed forcings'));
1158- disp(sprintf('%19s: %-22s -- %s','materials' ,['[1x1 ' class(self.materials) ']'],'material properties'));
1159- disp(sprintf('%19s: %-22s -- %s','damage' ,['[1x1 ' class(self.damage) ']'],'parameters for damage evolution solution'));
1160- disp(sprintf('%19s: %-22s -- %s','friction' ,['[1x1 ' class(self.friction) ']'],'basal friction/drag properties'));
1161- disp(sprintf('%19s: %-22s -- %s','flowequation' ,['[1x1 ' class(self.flowequation) ']'],'flow equations'));
1162- disp(sprintf('%19s: %-22s -- %s','timestepping' ,['[1x1 ' class(self.timestepping) ']'],'time stepping for transient models'));
1163- disp(sprintf('%19s: %-22s -- %s','initialization' ,['[1x1 ' class(self.initialization) ']'],'initial guess/state'));
1164- disp(sprintf('%19s: %-22s -- %s','rifts' ,['[1x1 ' class(self.rifts) ']'],'rifts properties'));
1165- disp(sprintf('%19s: %-22s -- %s','solidearth' ,['[1x1 ' class(self.solidearth) ']'],'solidearth inputs and settings'));
1166- disp(sprintf('%19s: %-22s -- %s','dsl' ,['[1x1 ' class(self.dsl) ']'],'dynamic sea-level '));
1167- disp(sprintf('%19s: %-22s -- %s','debug' ,['[1x1 ' class(self.debug) ']'],'debugging tools (valgrind, gprof)'));
1168- disp(sprintf('%19s: %-22s -- %s','verbose' ,['[1x1 ' class(self.verbose) ']'],'verbosity level in solve'));
1169- disp(sprintf('%19s: %-22s -- %s','settings' ,['[1x1 ' class(self.settings) ']'],'settings properties'));
1170- disp(sprintf('%19s: %-22s -- %s','toolkits' ,['[1x1 ' class(self.toolkits) ']'],'PETSc options for each solution'));
1171- disp(sprintf('%19s: %-22s -- %s','cluster' ,['[1x1 ' class(self.cluster) ']'],'cluster parameters (number of CPUs...)'));
1172- disp(sprintf('%19s: %-22s -- %s','balancethickness',['[1x1 ' class(self.balancethickness) ']'],'parameters for balancethickness solution'));
1173- disp(sprintf('%19s: %-22s -- %s','stressbalance' ,['[1x1 ' class(self.stressbalance) ']'],'parameters for stressbalance solution'));
1174- disp(sprintf('%19s: %-22s -- %s','groundingline' ,['[1x1 ' class(self.groundingline) ']'],'parameters for groundingline solution'));
1175- disp(sprintf('%19s: %-22s -- %s','hydrology' ,['[1x1 ' class(self.hydrology) ']'],'parameters for hydrology solution'));
1176- disp(sprintf('%19s: %-22s -- %s','masstransport' ,['[1x1 ' class(self.masstransport) ']'],'parameters for masstransport solution'));
1177- disp(sprintf('%19s: %-22s -- %s','thermal' ,['[1x1 ' class(self.thermal) ']'],'parameters for thermal solution'));
1178- disp(sprintf('%19s: %-22s -- %s','steadystate' ,['[1x1 ' class(self.steadystate) ']'],'parameters for steadystate solution'));
1179- disp(sprintf('%19s: %-22s -- %s','transient' ,['[1x1 ' class(self.transient) ']'],'parHwoameters for transient solution'));
1180- disp(sprintf('%19s: %-22s -- %s','levelset' ,['[1x1 ' class(self.levelset) ']'],'parameters for moving boundaries (level-set method)'));
1181- disp(sprintf('%19s: %-22s -- %s','calving' ,['[1x1 ' class(self.calving) ']'],'parameters for calving'));
1182- disp(sprintf('%19s: %-22s -- %s','frontalforcings' ,['[1x1 ' class(self.frontalforcings) ']'],'parameters for frontalforcings'));
1183- disp(sprintf('%19s: %-22s -- %s','esa' ,['[1x1 ' class(self.esa) ']'],'parameters for elastic adjustment solution'));
1184- disp(sprintf('%19s: %-22s -- %s','love' ,['[1x1 ' class(self.love) ']'],'parameters for love solution'));
1185- disp(sprintf('%19s: %-22s -- %s','sampling' ,['[1x1 ' class(self.sampling) ']'],'parameters for stochastic sampler'));
1186- disp(sprintf('%19s: %-22s -- %s','autodiff' ,['[1x1 ' class(self.autodiff) ']'],'automatic differentiation parameters'));
1187- disp(sprintf('%19s: %-22s -- %s','inversion' ,['[1x1 ' class(self.inversion) ']'],'parameters for inverse methods'));
1188- disp(sprintf('%19s: %-22s -- %s','qmu' ,['[1x1 ' class(self.qmu) ']'],'Dakota properties'));
1189- disp(sprintf('%19s: %-22s -- %s','amr' ,['[1x1 ' class(self.amr) ']'],'adaptive mesh refinement properties'));
1190- disp(sprintf('%19s: %-22s -- %s','outputdefinition',['[1x1 ' class(self.outputdefinition) ']'],'output definition'));
1191- disp(sprintf('%19s: %-22s -- %s','results' ,['[1x1 ' class(self.results) ']'],'model results'));
1192- disp(sprintf('%19s: %-22s -- %s','radaroverlay' ,['[1x1 ' class(self.radaroverlay) ']'],'radar image for plot overlay'));
1193- disp(sprintf('%19s: %-22s -- %s','miscellaneous' ,['[1x1 ' class(self.miscellaneous) ']'],'miscellaneous fields'));
1194- end % }}}
1195 function memory(self) % {{{
1196
1197- disp(sprintf('\nMemory imprint:\n'));
1198+ disp(sprintf('\nMemory imprint:\n'));
1199
1200- fields=properties('model');
1201- mem=0;
1202+ fields=properties('model');
1203+ mem=0;
1204
1205- for i=1:length(fields),
1206- field=self.(fields{i});
1207- s=whos('field');
1208- mem=mem+s.bytes/1e6;
1209- disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
1210+ for i=1:length(fields),
1211+ field=self.(fields{i});
1212+ s=whos('field');
1213+ mem=mem+s.bytes/1e6;
1214+ disp(sprintf('%19s: %6.2f Mb',fields{i},s.bytes/1e6));
1215+ end
1216+ disp(sprintf('%19s--%10s','--------------','--------------'));
1217+ disp(sprintf('%19s: %g Mb','Total',mem));
1218 end
1219- disp(sprintf('%19s--%10s','--------------','--------------'));
1220- disp(sprintf('%19s: %g Mb','Total',mem));
1221- end % }}}
1222+ % }}}
1223 function netcdf(self,filename) % {{{
1224- %NETCDF - save model as netcdf
1225- %
1226- % Usage:
1227- % netcdf(md,filename)
1228- %
1229- % Example:
1230- % netcdf(md,'model.nc');
1231+ %NETCDF - save model as netcdf
1232+ %
1233+ % Usage:
1234+ % netcdf(md,filename)
1235+ %
1236+ % Example:
1237+ % netcdf(md,'model.nc');
1238
1239- disp('Saving model as NetCDF');
1240- %1. Create NetCDF file
1241- ncid=netcdf.create(filename,'CLOBBER');
1242- netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
1243- netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
1244- netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
1245- netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
1246+ disp('Saving model as NetCDF');
1247+ %1. Create NetCDF file
1248+ ncid=netcdf.create(filename,'CLOBBER');
1249+ netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Conventions','CF-1.4');
1250+ netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Title',['ISSM model (' self.miscellaneous.name ')']);
1251+ netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Author',getenv('USER'));
1252+ netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'),'Date',datestr(now));
1253
1254- %Preallocate variable id, needed to write variables in netcdf file
1255- var_id=zeros(1000,1);%preallocate
1256+ %Preallocate variable id, needed to write variables in netcdf file
1257+ var_id=zeros(1000,1);%preallocate
1258
1259- for step=1:2,
1260- counter=0;
1261- [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
1262- if step==1, netcdf.endDef(ncid); end
1263- end
1264+ for step=1:2,
1265+ counter=0;
1266+ [var_id,counter]=structtonc(ncid,'md',self,0,var_id,counter,step);
1267+ if step==1, netcdf.endDef(ncid); end
1268+ end
1269
1270- if counter>1000,
1271- warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
1272- end
1273+ if counter>1000,
1274+ warning(['preallocation of var_id need to be updated from ' num2str(1000) ' to ' num2str(counter)]);
1275+ end
1276
1277- netcdf.close(ncid)
1278+ netcdf.close(ncid)
1279 end % }}}
1280 function xylim(self) % {{{
1281
1282@@ -1681,40 +1682,40 @@
1283 ylim([min(self.mesh.y) max(self.mesh.y)])
1284 end % }}}
1285 function md=upload(md) % {{{
1286- %the goal of this routine is to upload the model onto a server, and to empty it.
1287- %So first, save the model with a unique name and upload the file to the server:
1288- random_part=fix(rand(1)*10000);
1289- id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload'];
1290- eval(['save ' id ' md']);
1291+ %the goal of this routine is to upload the model onto a server, and to empty it.
1292+ %So first, save the model with a unique name and upload the file to the server:
1293+ random_part=fix(rand(1)*10000);
1294+ id=[md.miscellaneous.name '-' regexprep(datestr(now),'[^\w'']','') '-' num2str(random_part) '-' getenv('USER') '-' oshostname() '.upload'];
1295+ eval(['save ' id ' md']);
1296
1297- %Now, upload the file:
1298- issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
1299+ %Now, upload the file:
1300+ issmscpout(md.settings.upload_server,md.settings.upload_path,md.settings.upload_login,md.settings.upload_port,{id},1);
1301
1302- %Now, empty this model of everything except settings, and record name of file we just uploaded!
1303- settings_back=md.settings;
1304- md=model();
1305- md.settings=settings_back;
1306- md.settings.upload_filename=id;
1307+ %Now, empty this model of everything except settings, and record name of file we just uploaded!
1308+ settings_back=md.settings;
1309+ md=model();
1310+ md.settings=settings_back;
1311+ md.settings.upload_filename=id;
1312
1313- %get locally rid of file that was uploaded
1314- eval(['delete ' id]);
1315+ %get locally rid of file that was uploaded
1316+ eval(['delete ' id]);
1317
1318 end % }}}
1319 function md=download(md) % {{{
1320
1321- %the goal of this routine is to download the internals of the current model from a server, because
1322- %this model is empty, except for the settings which tell us where to go and find this model!
1323+ %the goal of this routine is to download the internals of the current model from a server, because
1324+ %this model is empty, except for the settings which tell us where to go and find this model!
1325
1326- %Download the file:
1327- issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
1328+ %Download the file:
1329+ issmscpin(md.settings.upload_server, md.settings.upload_login, md.settings.upload_port, md.settings.upload_path, {md.settings.upload_filename});
1330
1331- name=md.settings.upload_filename;
1332+ name=md.settings.upload_filename;
1333
1334- %Now, load this model:
1335- md=loadmodel(md.settings.upload_filename);
1336+ %Now, load this model:
1337+ md=loadmodel(md.settings.upload_filename);
1338
1339- %get locally rid of file that was downloaded
1340- eval(['delete ' name]);
1341+ %get locally rid of file that was downloaded
1342+ eval(['delete ' name]);
1343
1344 end % }}}
1345 function savemodeljs(md,modelname,websiteroot,varargin) % {{{
1346Index: ../trunk-jpl/src/m/classes/model.py
1347===================================================================
1348--- ../trunk-jpl/src/m/classes/model.py (revision 26300)
1349+++ ../trunk-jpl/src/m/classes/model.py (revision 26301)
1350@@ -75,31 +75,39 @@
1351
1352
1353 class model(object):
1354- #properties
1355+ """MODEL - class definition
1356+
1357+ Usage:
1358+ md = model()
1359+ """
1360+
1361 def __init__(self, *args): #{{{
1362 self.mesh = None
1363 self.mask = None
1364+
1365+ self.geometry = None
1366 self.constants = None
1367- self.geometry = None
1368- self.initialization = None
1369 self.smb = None
1370 self.basalforcings = None
1371+ self.materials = None
1372+ self.damage = None
1373 self.friction = None
1374+ self.flowequation = None
1375+ self.timestepping = None
1376+ self.initialization = None
1377 self.rifts = None
1378+ self.dsl = None
1379 self.solidearth = None
1380- self.dsl = None
1381- self.timestepping = None
1382- self.groundingline = None
1383- self.materials = None
1384- self.damage = None
1385- self.flowequation = None
1386+
1387 self.debug = None
1388 self.verbose = None
1389 self.settings = None
1390 self.toolkits = None
1391 self.cluster = None
1392+
1393 self.balancethickness = None
1394 self.stressbalance = None
1395+ self.groundingline = None
1396 self.hydrology = None
1397 self.masstransport = None
1398 self.thermal = None
1399@@ -111,81 +119,30 @@
1400 self.love = None
1401 self.esa = None
1402 self.sampling = None
1403+
1404 self.autodiff = None
1405 self.inversion = None
1406 self.qmu = None
1407 self.amr = None
1408- self.radaroverlay = None
1409 self.results = None
1410 self.outputdefinition = None
1411+ self.radaroverlay = None
1412 self.miscellaneous = None
1413 self.private = None
1414
1415- nargs = len(args)
1416-
1417- if nargs == 0:
1418+ if len(args) == 0:
1419 self.setdefaultparameters('earth')
1420 else:
1421- self.setdefaultparameters(args[0])
1422 options = pairoptions(*args)
1423 planet = options.getfieldvalue('planet', 'earth')
1424 self.setdefaultparameters(planet)
1425-#}}}
1426+ #}}}
1427
1428- def properties(self): # {{{
1429- # ordered list of properties since vars(self) is random
1430- return ['mesh',
1431- 'mask',
1432- 'constants',
1433- 'geometry',
1434- 'initialization',
1435- 'smb',
1436- 'basalforcings',
1437- 'friction',
1438- 'rifts',
1439- 'solidearth',
1440- 'dsl',
1441- 'timestepping',
1442- 'groundingline',
1443- 'materials',
1444- 'damage',
1445- 'flowequation',
1446- 'debug',
1447- 'verbose',
1448- 'settings',
1449- 'toolkits',
1450- 'cluster',
1451- 'balancethickness',
1452- 'stressbalance',
1453- 'hydrology',
1454- 'masstransport',
1455- 'thermal',
1456- 'steadystate',
1457- 'transient',
1458- 'levelset',
1459- 'calving',
1460- 'frontalforcings',
1461- 'love',
1462- 'esa',
1463- 'sampling',
1464- 'autodiff',
1465- 'inversion',
1466- 'qmu',
1467- 'amr',
1468- 'radaroverlay',
1469- 'results',
1470- 'outputdefinition',
1471- 'miscellaneous',
1472- 'private']
1473- # }}}
1474-
1475 def __repr__(obj): #{{{
1476 # TODO:
1477 # - Convert all formatting to calls to <string>.format (see any
1478 # already converted <class>.__repr__ method for examples)
1479 #
1480-
1481- #print "Here %s the number: %d" % ("is", 37)
1482 s = "%19s: %-22s -- %s" % ("mesh", "[%s %s]" % ("1x1", obj.mesh.__class__.__name__), "mesh properties")
1483 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("mask", "[%s %s]" % ("1x1", obj.mask.__class__.__name__), "defines grounded and floating elements"))
1484 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("geometry", "[%s %s]" % ("1x1", obj.geometry.__class__.__name__), "surface elevation, bedrock topography, ice thickness, ..."))
1485@@ -229,8 +186,55 @@
1486 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("radaroverlay", "[%s %s]" % ("1x1", obj.radaroverlay.__class__.__name__), "radar image for plot overlay"))
1487 s = "%s\n%s" % (s, "%19s: %-22s -- %s" % ("miscellaneous", "[%s %s]" % ("1x1", obj.miscellaneous.__class__.__name__), "miscellaneous fields"))
1488 return s
1489- # }}}
1490+ #}}}
1491
1492+ def properties(self): #{{{
1493+ # ordered list of properties since vars(self) is random
1494+ return ['mesh',
1495+ 'mask',
1496+ 'geometry',
1497+ 'constants',
1498+ 'smb',
1499+ 'basalforcings',
1500+ 'materials',
1501+ 'damage',
1502+ 'friction',
1503+ 'flowequation',
1504+ 'timestepping',
1505+ 'initialization',
1506+ 'rifts',
1507+ 'dsl',
1508+ 'solidearth',
1509+ 'debug',
1510+ 'verbose',
1511+ 'settings',
1512+ 'toolkits',
1513+ 'cluster',
1514+ 'balancethickness',
1515+ 'stressbalance',
1516+ 'groundingline',
1517+ 'hydrology',
1518+ 'masstransport',
1519+ 'thermal',
1520+ 'steadystate',
1521+ 'transient',
1522+ 'levelset',
1523+ 'calving',
1524+ 'frontalforcings',
1525+ 'love',
1526+ 'esa',
1527+ 'sampling',
1528+ 'autodiff',
1529+ 'inversion',
1530+ 'qmu',
1531+ 'amr',
1532+ 'results',
1533+ 'outputdefinition',
1534+ 'radaroverlay',
1535+ 'miscellaneous',
1536+ 'private']
1537+ #}}}
1538+
1539 def setdefaultparameters(self, planet): #{{{
1540 self.mesh = mesh2d()
1541 self.mask = mask()
1542@@ -277,14 +281,14 @@
1543 self.private = private()
1544 #}}}
1545
1546- def checkmessage(self, string): # {{{
1547+ def checkmessage(self, string): #{{{
1548 print("model not consistent: {}".format(string))
1549 self.private.isconsistent = False
1550 return self
1551- # }}}
1552+ #}}}
1553 #@staticmethod
1554
1555- def extract(self, area): # {{{
1556+ def extract(self, area): #{{{
1557 """EXTRACT - extract a model according to an Argus contour or flag list
1558
1559 This routine extracts a submodel from a bigger model with respect to a given contour
1560@@ -561,9 +565,9 @@
1561 md2.mesh.extractedelements = pos_elem + 1
1562
1563 return md2
1564- # }}}
1565+ #}}}
1566
1567- def extrude(md, *args): # {{{
1568+ def extrude(md, *args): #{{{
1569 """EXTRUDE - vertically extrude a 2d mesh
1570
1571 vertically extrude a 2d mesh and create corresponding 3d mesh.
1572@@ -748,7 +752,7 @@
1573 md.mesh.average_vertex_connectivity = 100
1574
1575 return md
1576- # }}}
1577+ #}}}
1578
1579 def collapse(md): #{{{
1580 """COLLAPSE - collapses a 3d mesh into a 2d mesh
1581Index: ../trunk-jpl/src/m/classes/sampling.m
1582===================================================================
1583--- ../trunk-jpl/src/m/classes/sampling.m (revision 26300)
1584+++ ../trunk-jpl/src/m/classes/sampling.m (revision 26301)
1585@@ -5,107 +5,98 @@
1586
1587 classdef sampling
1588 properties (SetAccess=public)
1589- kappa = NaN;
1590- tau = 0;
1591- beta = NaN;
1592- phi = 0;
1593- alpha = 0;
1594- robin = 0;
1595- seed = 0;
1596- requested_outputs = {};
1597- end
1598+ kappa = NaN;
1599+ tau = 0;
1600+ beta = NaN;
1601+ phi = 0;
1602+ alpha = 0;
1603+ robin = 0;
1604+ seed = 0;
1605+ requested_outputs = {};
1606+ end
1607 methods
1608 function self = sampling(varargin) % {{{
1609- switch nargin
1610+ switch nargin
1611 case 0
1612 self=setdefaultparameters(self);
1613 otherwise
1614 error('constructor not supported');
1615- end
1616+ end
1617 end % }}}
1618- function list = defaultoutputs(self,md) % {{{
1619+ function disp(self) % {{{
1620
1621- list = {};
1622+ disp(sprintf(' Sampling parameters:'));
1623
1624- end % }}}
1625+ disp(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
1626+ fielddisplay(self,'kappa','coefficient of the identity operator');
1627+ fielddisplay(self,'tau','scaling coefficient of the solution (default: 1.0)');
1628+ fielddisplay(self,'alpha','exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)');
1629+
1630+ disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
1631+ fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)');
1632+ fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');
1633+
1634+ disp(sprintf('\n %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
1635+ fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
1636+
1637+ disp(sprintf('\n %s','Other parameters of stochastic sampler:'));
1638+ fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)');
1639+ fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
1640+
1641+ end % }}}','
1642 function self = setdefaultparameters(self) % {{{
1643-
1644- %Scaling coefficient
1645- self.tau=1;
1646
1647- %Apply Robin boundary conditions
1648- self.robin=0;
1649-
1650- %Temporal correlation factor
1651- self.phi=0;
1652-
1653- %Exponent in fraction SPDE (default=2, biLaplacian covariance
1654+ %Scaling coefficient
1655+ self.tau=1;
1656+
1657+ %Apply Robin boundary conditions
1658+ self.robin=0;
1659+
1660+ %Temporal correlation factor
1661+ self.phi=0;
1662+
1663+ %Exponent in fraction SPDE (default: 2, biLaplacian covariance
1664 %operator)
1665 self.alpha=2; % Default
1666-
1667- %Seed for pseudorandom number generator (default -1 for random seed)
1668- self.seed=-1;
1669-
1670+
1671+ %Seed for pseudorandom number generator (default: -1, for random seed)
1672+ self.seed=-1;
1673+
1674 %default output
1675 self.requested_outputs={'default'};
1676
1677 end % }}}
1678- function md = setparameters(self,md,lc,sigma) % {{{
1679-
1680- nu = self.alpha-1;
1681- KAPPA = sqrt(8*nu)/lc;
1682- TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
1683- md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
1684- md.sampling.tau = TAU;
1685-
1686+ function list = defaultoutputs(self,md) % {{{
1687+
1688+ list = {};
1689+
1690 end % }}}
1691 function md = checkconsistency(self,md,solution,analyses) % {{{
1692-
1693- if ~ismember('SamplingAnalysis',analyses), return; end
1694-
1695- md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
1696- md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
1697- md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
1698- if(md.sampling.robin)
1699- md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
1700- end
1701- md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
1702- md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
1703- md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
1704- md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
1705
1706- end % }}}
1707- function disp(self) % {{{
1708+ if ~ismember('SamplingAnalysis',analyses), return; end
1709
1710- disp(sprintf(' Sampling parameters:'));
1711+ md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
1712+ md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0);
1713+ md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0 1]);
1714+ if(md.sampling.robin)
1715+ md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices 1],'>',0);
1716+ end
1717+ md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0);
1718+ md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0);
1719+ md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1);
1720+ md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1);
1721
1722- disp(sprintf('\n %s','Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):'));
1723- fielddisplay(self,'kappa','coefficient of the identity operator');
1724- fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)');
1725- fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)');
1726-
1727- disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():'));
1728- fielddisplay(self,'robin','Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default 0)');
1729- fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)');
1730-
1731- disp(sprintf('\n %s','Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):'));
1732- fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)');
1733-
1734- disp(sprintf('\n %s','Other parameters of stochastic sampler:'));
1735- fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)');
1736- fielddisplay(self,'requested_outputs','additional outputs requested (not implemented yet)');
1737-
1738 end % }}}
1739 function marshall(self,prefix,md,fid) % {{{
1740
1741- WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
1742- WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
1743- WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
1744- WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
1745- WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
1746- WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
1747- WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
1748-
1749+ WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1);
1750+ WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double');
1751+ WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1);
1752+ WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double');
1753+ WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer');
1754+ WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean');
1755+ WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer');
1756+
1757 %process requested outputs
1758 outputs = self.requested_outputs;
1759 pos = find(ismember(outputs,'default'));
1760@@ -115,17 +106,26 @@
1761 end
1762 WriteData(fid,prefix,'data',outputs,'name','md.sampling.requested_outputs','format','StringArray');
1763 end % }}}
1764+ function md = setparameters(self,md,lc,sigma) % {{{
1765+
1766+ nu = self.alpha-1;
1767+ KAPPA = sqrt(8*nu)/lc;
1768+ TAU = sqrt(gamma(nu)/(gamma(self.alpha)*(4*pi)*KAPPA^(2*nu)*sigma^2));
1769+ md.sampling.kappa = KAPPA*ones(md.mesh.numberofvertices,1);
1770+ md.sampling.tau = TAU;
1771+
1772+ end % }}}
1773 function savemodeljs(self,fid,modelname) % {{{
1774-
1775- writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
1776- writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
1777- writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
1778- writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
1779- writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
1780- writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
1781- writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
1782+
1783+ writejsdouble(fid,[modelname '.sampling.kappa'],self.kappa);
1784+ writejsdouble(fid,[modelname '.sampling.tau'],self.tau);
1785+ writejsdouble(fid,[modelname '.sampling.beta'],self.beta);
1786+ writejsdouble(fid,[modelname '.sampling.phi'],self.beta);
1787+ writejsdouble(fid,[modelname '.sampling.alpha'],self.alpha);
1788+ writejsdouble(fid,[modelname '.sampling.robin'],self.robin);
1789+ writejsdouble(fid,[modelname '.sampling.seed'],self.seed);
1790 writejscellstring(fid,[modelname '.sampling.requested_outputs'],self.requested_outputs);
1791
1792 end % }}}
1793 end
1794-end
1795\ No newline at end of file
1796+end
1797Index: ../trunk-jpl/src/m/classes/sampling.py
1798===================================================================
1799--- ../trunk-jpl/src/m/classes/sampling.py (revision 26300)
1800+++ ../trunk-jpl/src/m/classes/sampling.py (revision 26301)
1801@@ -1,5 +1,7 @@
1802 import numpy as np
1803
1804+from math import *
1805+
1806 from checkfield import checkfield
1807 from fielddisplay import fielddisplay
1808 from project3d import project3d
1809@@ -13,11 +15,10 @@
1810 sampling = sampling()
1811 """
1812
1813- def __init__(self): # {{{
1814- self.kappa = float('NaN')
1815+ def __init__(self, *args): #{{{
1816+ self.kappa = np.nan
1817 self.tau = 0
1818- self.beta = float('NaN')
1819- self.hydrostatic_adjustment = 0
1820+ self.beta = np.nan
1821 self.phi = 0
1822 self.alpha = 0
1823 self.robin = 0
1824@@ -24,72 +25,87 @@
1825 self.seed = 0
1826 self.requested_outputs = []
1827
1828- # Set defaults
1829- self.setdefaultparameters()
1830-
1831+ if len(args) == 0:
1832+ self.setdefaultparameters()
1833+ else:
1834+ raise RuntimeError('constructor not supported')
1835 #}}}
1836
1837- def __repr__(self): # {{{
1838+ def __repr__(self): #{{{
1839 s = ' Sampling parameters::\n'
1840- s += '{}\n'.format(fielddisplay(self,'kappa','coefficient of the identity operator'))
1841- s += '{}\n'.format(fielddisplay(self,'tau','scaling coefficient of the solution (default 1.0)'))
1842- s += '{}\n'.format(fielddisplay(self,'alpha','exponent in PDE operator, (default 2.0, BiLaplacian covariance operator)'))
1843- s += '{}\n'.format(disp(sprintf('\n %s','Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():')))
1844- s += '{}\n'.format(fielddisplay(self,'beta','Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
1845- s += '{}\n'.format(fielddisplay(self,'phi','Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
1846- s += '{}\n'.format(fielddisplay(self,'seed','Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default -1)'))
1847+ s += ' Parameters of PDE operator (kappa^2 I-Laplacian)^(alpha/2)(tau):\n'
1848+ s += '{}\n'.format(fielddisplay(self, 'kappa', 'coefficient of the identity operator'))
1849+ s += '{}\n'.format(fielddisplay(self, 'tau', 'scaling coefficient of the solution (default: 1.0)'))
1850+ s += '{}\n'.format(fielddisplay(self, 'alpha', 'exponent in PDE operator, (default: 2.0, BiLaplacian covariance operator)'))
1851+
1852+ s += ' Parameters of Robin boundary conditions nabla () \cdot normvec + beta ():\n'
1853+ s += '{}\n'.format(fielddisplay(self, 'robin', 'Apply Robin boundary conditions (1 if applied and 0 for homogenous Neumann boundary conditions) (default: 0)'))
1854+ s += '{}\n'.format(fielddisplay(self, 'beta', 'Coefficient in Robin boundary conditions (to be defined for robin = 1)'))
1855+
1856+ s += ' Parameters for first-order autoregressive process (X_t = phi X_{t-1} + noise) (if transient):\n'
1857+ s += '{}\n'.format(fielddisplay(self, 'phi', 'Temporal correlation factor (|phi|<1 for stationary process, phi = 1 for random walk process) (default 0)'))
1858+
1859+ s += ' Other parameters of stochastic sampler:\n'
1860+ s += '{}\n'.format(fielddisplay(self, 'seed', 'Seed for pseudorandom number generator (given seed if >=0 and random seed if <0) (default: -1)'))
1861 s += '{}\n'.format(fielddisplay(self, 'requested_outputs', 'additional outputs requested (not implemented yet)'))
1862+
1863 return s
1864 #}}}
1865
1866- def defaultoutputs(self, md): # {{{
1867- return []
1868-
1869- #}}}
1870- def setdefaultparameters(self): # {{{
1871+ def setdefaultparameters(self): #{{{
1872 # Scaling coefficient
1873 self.tau = 1
1874+
1875 # Apply Robin boundary conditions
1876 self.robin = 0
1877+
1878 # Temporal correlation factor
1879 self.phi = 0
1880- # Exponent in fraction SPDE (default=2, biLaplacian covariance operator)
1881+
1882+ # Exponent in fraction SPDE (default: 2, biLaplacian covariance operator)
1883 self.alpha = 2
1884- # Seed for pseudorandom number generator (default -1 for random seed)
1885- self.alpha = -1
1886+
1887+ # Seed for pseudorandom number generator (default: -1, for random seed)
1888+ self.seed = -1
1889+
1890 # Default output
1891 self.requested_outputs = ['default']
1892+
1893 return self
1894 #}}}
1895
1896- def checkconsistency(self, md, solution, analyses): # {{{
1897- # Early return
1898+ def defaultoutputs(self, md): #{{{
1899+ return []
1900+ #}}}
1901+
1902+ def checkconsistency(self, md, solution, analyses): #{{{
1903 if ('SamplingAnalysis' not in analyses):
1904 return md
1905
1906- md = checkfield(md,'fieldname','sampling.kappa','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
1907- md = checkfield(md,'fieldname','sampling.tau','NaN',1,'Inf',1,'numel',1,'>',0)
1908- md = checkfield(md,'fieldname','sampling.robin','numel',1,'values',[0, 1])
1909+ md = checkfield(md, 'fieldname', 'sampling.kappa', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
1910+ md = checkfield(md, 'fieldname', 'sampling.tau', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
1911+ md = checkfield(md, 'fieldname', 'sampling.robin', 'numel', 1, 'values', [0, 1])
1912 if md.sampling.robin:
1913- md = checkfield(md,'fieldname','sampling.beta','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices],'>',0)
1914- md = checkfield(md,'fieldname','sampling.phi','NaN',1,'Inf',1,'numel',1,'>=',0)
1915- md = checkfield(md,'fieldname','sampling.alpha','NaN',1,'Inf',1,'numel',1,'>',0)
1916- md = checkfield(md,'fieldname','sampling.seed','NaN',1,'Inf',1,'numel',1)
1917- md = checkfield(md,'fieldname','sampling.requested_outputs','stringrow',1)
1918+ md = checkfield(md, 'fieldname', 'sampling.beta', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices], '>', 0)
1919+ end
1920+ md = checkfield(md, 'fieldname', 'sampling.phi', 'NaN', 1, 'Inf', 1, 'numel', 1, '>=', 0)
1921+ md = checkfield(md, 'fieldname', 'sampling.alpha', 'NaN', 1, 'Inf', 1, 'numel', 1, '>', 0)
1922+ md = checkfield(md, 'fieldname', 'sampling.seed', 'NaN', 1, 'Inf', 1, 'numel', 1)
1923+ md = checkfield(md, 'fieldname', 'sampling.requested_outputs', 'stringrow', 1)
1924
1925 return md
1926- # }}}
1927+ #}}}
1928
1929- def marshall(self, prefix, md, fid): # {{{
1930- WriteData(fid,prefix,'object',self,'fieldname','kappa','format','DoubleMat','mattype',1)
1931- WriteData(fid,prefix,'object',self,'fieldname','tau','format','Double')
1932- WriteData(fid,prefix,'object',self,'fieldname','beta','format','DoubleMat','mattype',1)
1933- WriteData(fid,prefix,'object',self,'fieldname','phi','format','Double')
1934- WriteData(fid,prefix,'object',self,'fieldname','alpha','format','Integer')
1935- WriteData(fid,prefix,'object',self,'fieldname','robin','format','Boolean')
1936- WriteData(fid,prefix,'object',self,'fieldname','seed','format','Integer')
1937+ def marshall(self, prefix, md, fid): #{{{
1938+ WriteData(fid, prefix, 'object', self, 'fieldname', 'kappa', 'format', 'DoubleMat', 'mattype', 1)
1939+ WriteData(fid, prefix, 'object', self, 'fieldname', 'tau', 'format', 'Double')
1940+ WriteData(fid, prefix, 'object', self, 'fieldname', 'beta', 'format', 'DoubleMat', 'mattype', 1)
1941+ WriteData(fid, prefix, 'object', self, 'fieldname', 'phi', 'format', 'Double')
1942+ WriteData(fid, prefix, 'object', self, 'fieldname', 'alpha', 'format', 'Integer')
1943+ WriteData(fid, prefix, 'object', self, 'fieldname', 'robin', 'format', 'Boolean')
1944+ WriteData(fid, prefix, 'object', self, 'fieldname', 'seed', 'format', 'Integer')
1945
1946- #process requested outputs
1947+ # Process requested outputs
1948 outputs = self.requested_outputs
1949 indices = [i for i, x in enumerate(outputs) if x == 'default']
1950 if len(indices) > 0:
1951@@ -96,5 +112,14 @@
1952 outputscopy = outputs[0:max(0, indices[0] - 1)] + self.defaultoutputs(md) + outputs[indices[0] + 1:]
1953 outputs = outputscopy
1954 WriteData(fid, prefix, 'data', outputs, 'name', 'md.sampling.requested_outputs', 'format', 'StringArray')
1955+ #}}}
1956
1957- # }}}
1958+ def setparameters(self, md, lc, sigma): #{{{
1959+ nu = self.alpha - 1
1960+ KAPPA = pow((8 * nu), 0.5) / lc
1961+ TAU = pow((math.gamma(nu) / math.gamma(self.alpha) * (4 * np.pi) * pow(KAPPA, 2 * nu) * pow(sigma, 2)), 0.5)
1962+ md.sampling.kappa = KAPPA * np.ones((md.mesh.numberofvertices, 1))
1963+ md.sampling.tau = TAU
1964+
1965+ return md
1966+ #}}}
1967Index: ../trunk-jpl/src/m/classes/solidearthsettings.m
1968===================================================================
1969--- ../trunk-jpl/src/m/classes/solidearthsettings.m (revision 26300)
1970+++ ../trunk-jpl/src/m/classes/solidearthsettings.m (revision 26301)
1971@@ -91,6 +91,27 @@
1972 self.grdmodel=0;
1973
1974 end % }}}
1975+ function disp(self) % {{{
1976+ disp(sprintf(' solidearth settings:'));
1977+
1978+ fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
1979+ fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
1980+ fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
1981+ fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
1982+ fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
1983+ fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
1984+ fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
1985+ fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
1986+ fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
1987+ fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
1988+ fielddisplay(self,'elastic','enables elastic deformation from surface loading');
1989+ fielddisplay(self,'viscous','enables viscous deformation from surface loading');
1990+ fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
1991+ fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
1992+ fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
1993+ fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
1994+ fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
1995+ end % }}}
1996 function md = checkconsistency(self,md,solution,analyses) % {{{
1997
1998 if ~ismember('SealevelchangeAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.isslc==0),
1999@@ -134,27 +155,6 @@
2000 end
2001
2002 end % }}}
2003- function disp(self) % {{{
2004- disp(sprintf(' solidearth settings:'));
2005-
2006- fielddisplay(self,'reltol','sea level change relative convergence criterion (default, NaN: not applied)');
2007- fielddisplay(self,'abstol','sea level change absolute convergence criterion(default, NaN: not applied)');
2008- fielddisplay(self,'maxiter','maximum number of nonlinear iterations');
2009- fielddisplay(self,'grdocean','does this planet have an ocean, if set to 1: global water mass is conserved in GRD module (default: 1)');
2010- fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area (default: No correction)');
2011- fielddisplay(self,'sealevelloading','enables surface loading from sea-level change (default: 1)');
2012- fielddisplay(self,'isgrd','compute GRD patterns (default: 1)');
2013- fielddisplay(self,'compute_bp_grd','compute GRD patterns for bottom pressure loads (default: 1)');
2014- fielddisplay(self,'runfrequency','how many time steps we skip before we run solidearthsettings solver during transient (default: 1)');
2015- fielddisplay(self,'selfattraction','enables surface mass load to perturb the gravity field');
2016- fielddisplay(self,'elastic','enables elastic deformation from surface loading');
2017- fielddisplay(self,'viscous','enables viscous deformation from surface loading');
2018- fielddisplay(self,'rotation','enables polar motion to feedback on the GRD fields');
2019- fielddisplay(self,'degacc','accuracy (default: .01 deg) for numerical discretization of the Green''s functions');
2020- fielddisplay(self,'timeacc','time accuracy (default: 1 yr) for numerical discretization of the Green''s functions');
2021- fielddisplay(self,'grdmodel','type of deformation model, 0 for no GRD, 1 for spherical GRD model (SESAW model), 2 for half-space planar GRD (visco-elastic model from Ivins)');
2022- fielddisplay(self,'cross_section_shape','1: square-edged (default). 2: elliptical. See iedge in GiaDeflectionCore');
2023- end % }}}
2024 function marshall(self,prefix,md,fid) % {{{
2025 WriteData(fid,prefix,'object',self,'fieldname','reltol','name','md.solidearth.settings.reltol','format','Double');
2026 WriteData(fid,prefix,'object',self,'fieldname','abstol','name','md.solidearth.settings.abstol','format','Double');
2027Index: ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py
2028===================================================================
2029--- ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26300)
2030+++ ../trunk-jpl/src/m/miscellaneous/MatlabFuncs.py (revision 26301)
2031@@ -1,60 +1,79 @@
2032-def oshostname():
2033- import socket
2034+def acosd(X): #{{{
2035+ """ function acosd - Inverse cosine in degrees
2036
2037- return socket.gethostname()
2038+ Usage:
2039+ Y = acosd(X)
2040+ """
2041+ import numpy as np
2042
2043+ return np.degrees(np.arccos(X))
2044+#}}}
2045
2046-def ispc():
2047- import platform
2048+def asind(X): #{{{
2049+ """ function asind - Inverse sine in degrees
2050
2051- if 'Windows' in platform.system():
2052- return True
2053- else:
2054- return False
2055+ Usage:
2056+ Y = asind(X)
2057+ """
2058+ import numpy as np
2059
2060+ return np.degrees(np.arcsin(X))
2061+#}}}
2062
2063-def ismac():
2064- import platform
2065+def atand(X): #{{{
2066+ """ function atand - Inverse tangent in degrees
2067
2068- if 'Darwin' in platform.system():
2069- return True
2070- else:
2071- return False
2072+ Usage:
2073+ Y = atand(X)
2074+ """
2075+ import numpy as np
2076
2077+ return np.degrees(np.arctan(X))
2078+#}}}
2079
2080-def strcmp(s1, s2):
2081
2082- if s1 == s2:
2083- return True
2084- else:
2085- return False
2086+def atan2d(Y, X): #{{{
2087+ """ function atan2d - Four-quadrant inverse tangent in degrees
2088
2089+ Usage:
2090+ D = atan2d(Y, X)
2091+ """
2092+ import numpy as np
2093
2094-def strncmp(s1, s2, n):
2095+ return np.degrees(np.arctan2(Y, X))
2096+#}}}
2097
2098- if s1[0:n] == s2[0:n]:
2099- return True
2100+def det(a): #{{{
2101+ if a.shape == (1, ):
2102+ return a[0]
2103+ elif a.shape == (1, 1):
2104+ return a[0, 0]
2105+ elif a.shape == (2, 2):
2106+ return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
2107 else:
2108- return False
2109+ raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
2110+#}}}
2111
2112+def heaviside(x): #{{{
2113+ import numpy as np
2114
2115-def strcmpi(s1, s2):
2116+ y = np.zeros_like(x)
2117+ y[np.nonzero(x > 0.)] = 1.
2118+ y[np.nonzero(x == 0.)] = 0.5
2119
2120- if s1.lower() == s2.lower():
2121- return True
2122- else:
2123- return False
2124+ return y
2125+#}}}
2126
2127+def ismac(): #{{{
2128+ import platform
2129
2130-def strncmpi(s1, s2, n):
2131-
2132- if s1.lower()[0:n] == s2.lower()[0:n]:
2133+ if 'Darwin' in platform.system():
2134 return True
2135 else:
2136 return False
2137+#}}}
2138
2139-
2140-def ismember(a, s):
2141+def ismember(a, s): #{{{
2142 import numpy as np
2143
2144 if not isinstance(s, (tuple, list, dict, np.ndarray)):
2145@@ -65,32 +84,33 @@
2146
2147 if not isinstance(a, np.ndarray):
2148 b = [item in s for item in a]
2149-
2150 else:
2151 if not isinstance(s, np.ndarray):
2152 b = np.empty_like(a)
2153 for i, item in enumerate(a.flat):
2154 b.flat[i] = item in s
2155-
2156 else:
2157 b = np.in1d(a.flat, s.flat).reshape(a.shape)
2158
2159 return b
2160+#}}}
2161
2162+def ispc(): #{{{
2163+ import platform
2164
2165-def det(a):
2166-
2167- if a.shape == (1, ):
2168- return a[0]
2169- elif a.shape == (1, 1):
2170- return a[0, 0]
2171- elif a.shape == (2, 2):
2172- return a[0, 0] * a[1, 1] - a[0, 1] * a[1, 0]
2173+ if 'Windows' in platform.system():
2174+ return True
2175 else:
2176- raise TypeError("MatlabFunc.det only implemented for shape (2, 2), not for shape %s." % str(a.shape))
2177+ return False
2178+#}}}
2179
2180+def oshostname(): #{{{
2181+ import socket
2182
2183-def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0):
2184+ return socket.gethostname()
2185+#}}}
2186+
2187+def sparse(ivec, jvec, svec, m=0, n=0, nzmax=0): #{{{
2188 import numpy as np
2189
2190 if not m:
2191@@ -104,13 +124,32 @@
2192 a[i - 1, j - 1] += s
2193
2194 return a
2195+#}}}
2196
2197+def strcmp(s1, s2): #{{{
2198+ if s1 == s2:
2199+ return True
2200+ else:
2201+ return False
2202+#}}}
2203
2204-def heaviside(x):
2205- import numpy as np
2206+def strcmpi(s1, s2): #{{{
2207+ if s1.lower() == s2.lower():
2208+ return True
2209+ else:
2210+ return False
2211+#}}}
2212
2213- y = np.zeros_like(x)
2214- y[np.nonzero(x > 0.)] = 1.
2215- y[np.nonzero(x == 0.)] = 0.5
2216+def strncmp(s1, s2, n): #{{{
2217+ if s1[0:n] == s2[0:n]:
2218+ return True
2219+ else:
2220+ return False
2221+#}}}
2222
2223- return y
2224+def strncmpi(s1, s2, n): #{{{
2225+ if s1.lower()[0:n] == s2.lower()[0:n]:
2226+ return True
2227+ else:
2228+ return False
2229+#}}}
2230Index: ../trunk-jpl/src/m/miscellaneous/PythonFuncs.py
2231===================================================================
2232--- ../trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 26300)
2233+++ ../trunk-jpl/src/m/miscellaneous/PythonFuncs.py (revision 26301)
2234@@ -1,25 +1,22 @@
2235 import numpy as np
2236
2237
2238-def logical_and_n(*arg):
2239-
2240+def logical_and_n(*arg): #{{{
2241 if len(arg):
2242 result = arg[0]
2243 for item in arg[1:]:
2244 result = np.logical_and(result, item)
2245 return result
2246-
2247 else:
2248 return None
2249+#}}}
2250
2251-
2252-def logical_or_n(*arg):
2253-
2254+def logical_or_n(*arg): #{{{
2255 if len(arg):
2256 result = arg[0]
2257 for item in arg[1:]:
2258 result = np.logical_or(result, item)
2259 return result
2260-
2261 else:
2262 return None
2263+#}}}
2264Index: ../trunk-jpl/src/m/solve/marshall.m
2265===================================================================
2266--- ../trunk-jpl/src/m/solve/marshall.m (revision 26300)
2267+++ ../trunk-jpl/src/m/solve/marshall.m (revision 26301)
2268@@ -44,7 +44,7 @@
2269 st=fclose(fid);
2270
2271 % Uncomment the following to make a copy of the binary input file for debugging
2272-% purposes (can be fed into scripts/ReadBin.py).
2273+% purposes (can be fed into scripts/BinRead.py).
2274 % copyfile([md.miscellaneous.name '.bin'], [md.miscellaneous.name '.m.bin'])
2275
2276 if st==-1,
2277Index: ../trunk-jpl/src/m/solve/marshall.py
2278===================================================================
2279--- ../trunk-jpl/src/m/solve/marshall.py (revision 26300)
2280+++ ../trunk-jpl/src/m/solve/marshall.py (revision 26301)
2281@@ -45,7 +45,7 @@
2282 fid.close()
2283
2284 # Uncomment the following to make a copy of the binary input file for
2285- # debugging purposes (can be fed into scripts/ReadBin.py).
2286+ # debugging purposes (can be fed into scripts/BinRead.py).
2287 # copy_cmd = "cp {}.bin {}.py.bin".format(md.miscellaneous.name, md.miscellaneous.name)
2288 # subprocess.call(copy_cmd, shell=True)
2289
2290Index: ../trunk-jpl/src/m/classes/dsl.py
2291===================================================================
2292--- ../trunk-jpl/src/m/classes/dsl.py (revision 26300)
2293+++ ../trunk-jpl/src/m/classes/dsl.py (revision 26301)
2294@@ -7,16 +7,21 @@
2295
2296
2297 class dsl(object):
2298- """DSL - slass definition
2299+ """DSL - class definition
2300
2301 Usage:
2302- dsl = dsl()
2303+ dsl = dsl() # dynamic sea level class, based on CMIP5 outputs
2304 """
2305
2306- def __init__(self): #{{{
2307+ def __init__(self, *args): #{{{
2308 self.global_average_thermosteric_sea_level = np.nan # Corresponds to zostoga field in CMIP5 archives. Specified as a temporally variable quantity (in m).
2309 self.sea_surface_height_above_geoid = np.nan # Corresponds to zos field in CMIP5 archives. Spatial average is 0. Specified as a spatio-temporally variable quantity (in m).
2310 self.sea_water_pressure_at_sea_floor = np.nan # Corresponds to bpo field in CMIP5 archives. Specified as a spatio-temporally variable quantity (in m equivalent, not in Pa!).
2311+
2312+ if len(args) == 0:
2313+ self.setdefaultparameters()
2314+ else:
2315+ raise Exception('constructor not supported')
2316 #}}}
2317
2318 def __repr__(self): #{{{
2319@@ -27,27 +32,15 @@
2320 return s
2321 #}}}
2322
2323- def defaultoutputs(self, md): #{{{
2324- return []
2325+ def setdefaultparameters(self): #{{{
2326+ self.global_average_thermosteric_sea_level = np.nan
2327+ self.sea_surface_height_above_geoid = np.nan
2328+ self.sea_water_pressure_at_sea_floor = np.nan
2329 #}}}
2330
2331- def initialize(self, md): #{{{
2332- if np.isnan(self.global_average_thermosteric_sea_level):
2333- self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
2334- print(' no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
2335-
2336- if np.isnan(self.sea_surface_height_above_geoid):
2337- self.sea_surface_height_above_geoid = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
2338- print(' no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
2339-
2340- if np.isnan(self.sea_water_pressure_at_sea_floor):
2341- self.sea_water_pressure_at_sea_floor = np.array(np.zeros(md.mesh.numberofvertices)).reshape(-1, 1)
2342- print(' no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
2343- #}}}
2344-
2345 def checkconsistency(self, md, solution, analyses): #{{{
2346 # Early return
2347- if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
2348+ if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
2349 return md
2350
2351 md = checkfield(md, 'fieldname', 'dsl.global_average_thermosteric_sea_level', 'NaN', 1, 'Inf', 1)
2352@@ -60,16 +53,31 @@
2353 return md
2354 # }}}
2355
2356+ def marshall(self, prefix, md, fid): #{{{
2357+ yts = md.constants.yts
2358+ WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
2359+ WriteData(fid, prefix, 'object', self, 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeseries', 1, 'yts', yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
2360+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify DSL at vertex locations.
2361+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts) # mattype 1 because we specify bottom pressure at vertex locations.
2362+ # }}}
2363+
2364 def extrude(self, md): #{{{
2365- self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node')
2366- self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node')
2367+ self.sea_surface_height_above_geoid = project3d(md, 'vector', self.sea_surface_height_above_geoid, 'type', 'node', 'layer', 1)
2368+ self.sea_water_pressure_at_sea_floor = project3d(md, 'vector', self.sea_water_pressure_at_sea_floor, 'type', 'node', 'layer', 1)
2369 return self
2370 #}}}
2371
2372- def marshall(self, prefix, md, fid): #{{{
2373- yts = md.constants.yts
2374- WriteData(fid, prefix, 'name', 'md.dsl.model', 'data', 1, 'format', 'Integer')
2375- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'global_average_thermosteric_sea_level', 'format', 'DoubleMat', 'mattype', 2, 'timeserieslength', 1+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 2, because we are sending a GMSL value identical everywhere on each element.
2376- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_water_pressure_at_sea_floor', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts, 'scale', 1e-3/md.constants.yts) # mattype 1 because we specify DSL at vertex locations.
2377- WriteData(fid, prefix, 'object', self, 'class', 'dsl', 'fieldname', 'sea_surface_height_above_geoid', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices+1, 'yts', md.constants.yts) # mattype 1 because we specify bottom pressure at vertex locations.
2378- # }}}
2379+
2380+ def initialize(self, md): #{{{
2381+ if np.isnan(self.global_average_thermosteric_sea_level):
2382+ self.global_average_thermosteric_sea_level = np.array([0, 0]).reshape(-1, 1)
2383+ print(' no dsl.global_average_thermosteric_sea_level specified: transient values set at zero')
2384+
2385+ if np.isnan(self.sea_surface_height_above_geoid):
2386+ self.sea_surface_height_above_geoid = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
2387+ print(' no dsl.sea_surface_height_above_geoid specified: transient values set at zero')
2388+
2389+ if np.isnan(self.sea_water_pressure_at_sea_floor):
2390+ self.sea_water_pressure_at_sea_floor = np.append(np.zeros((md.mesh.numberofvertices, 1)), 0).reshape(-1, 1)
2391+ print(' no dsl.sea_water_pressure_at_sea_floor specified: transient values set at zero')
2392+ #}}}
2393Index: ../trunk-jpl/src/m/classes/hydrologyshreve.m
2394===================================================================
2395--- ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26300)
2396+++ ../trunk-jpl/src/m/classes/hydrologyshreve.m (revision 26301)
2397@@ -10,8 +10,6 @@
2398 requested_outputs = {};
2399 end
2400 methods
2401- function self = extrude(self,md) % {{{
2402- end % }}}
2403 function self = hydrologyshreve(varargin) % {{{
2404 switch nargin
2405 case 0
2406@@ -22,10 +20,13 @@
2407 error('constructor not supported');
2408 end
2409 end % }}}
2410- function list = defaultoutputs(self,md) % {{{
2411- list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
2412- end % }}}
2413+ function disp(self) % {{{
2414+ disp(sprintf(' hydrologyshreve solution parameters:'));
2415+ fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
2416+ fielddisplay(self,'stabilization','artificial diffusivity (default: 1). can be more than 1 to increase diffusivity.');
2417+ fielddisplay(self,'requested_outputs','additional outputs requested');
2418
2419+ end % }}}
2420 function self = setdefaultparameters(self) % {{{
2421
2422 %Type of stabilization to use 0:nothing 1:artificial_diffusivity
2423@@ -33,7 +34,7 @@
2424 self.requested_outputs = {'default'};
2425 end % }}}
2426 function md = checkconsistency(self,md,solution,analyses) % {{{
2427-
2428+
2429 %Early return
2430 if ~ismember('HydrologyShreveAnalysis',analyses) | (strcmp(solution,'TransientSolution') & md.transient.ishydrology==0), return; end
2431
2432@@ -40,25 +41,23 @@
2433 md = checkfield(md,'fieldname','hydrology.spcwatercolumn','Inf',1,'timeseries',1);
2434 md = checkfield(md,'fieldname','hydrology.stabilization','>=',0);
2435 end % }}}
2436- function disp(self) % {{{
2437- disp(sprintf(' hydrologyshreve solution parameters:'));
2438- fielddisplay(self,'spcwatercolumn','water thickness constraints (NaN means no constraint) [m]');
2439- fielddisplay(self,'stabilization','artificial diffusivity (default is 1). can be more than 1 to increase diffusivity.');
2440- fielddisplay(self,'requested_outputs','additional outputs requested');
2441-
2442+ function list = defaultoutputs(self,md) % {{{
2443+ list = {'Watercolumn','HydrologyWaterVx','HydrologyWaterVy'};
2444 end % }}}
2445 function marshall(self,prefix,md,fid) % {{{
2446 WriteData(fid,prefix,'name','md.hydrology.model','data',2,'format','Integer');
2447 WriteData(fid,prefix,'object',self,'fieldname','spcwatercolumn','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
2448 WriteData(fid,prefix,'object',self,'fieldname','stabilization','format','Double');
2449- outputs = self.requested_outputs;
2450- pos = find(ismember(outputs,'default'));
2451- if ~isempty(pos),
2452- outputs(pos) = []; %remove 'default' from outputs
2453- outputs = [outputs defaultoutputs(self,md)]; %add defaults
2454- end
2455- WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
2456+ outputs = self.requested_outputs;
2457+ pos = find(ismember(outputs,'default'));
2458+ if ~isempty(pos),
2459+ outputs(pos) = []; %remove 'default' from outputs
2460+ outputs = [outputs defaultoutputs(self,md)]; %add defaults
2461+ end
2462+ WriteData(fid,prefix,'data',outputs,'name','md.hydrology.requested_outputs','format','StringArray');
2463 end % }}}
2464+ function self = extrude(self,md) % {{{
2465+ end % }}}
2466 function savemodeljs(self,fid,modelname) % {{{
2467
2468 writejs1Darray(fid,[modelname '.hydrology.spcwatercolumn'],self.spcwatercolumn);
2469Index: ../trunk-jpl/src/m/classes/initialization.py
2470===================================================================
2471--- ../trunk-jpl/src/m/classes/initialization.py (revision 26300)
2472+++ ../trunk-jpl/src/m/classes/initialization.py (revision 26301)
2473@@ -10,32 +10,33 @@
2474 """INITIALIZATION class definition
2475
2476 Usage:
2477- initialization = initialization()
2478+ initialization = initialization()
2479 """
2480
2481- def __init__(self): # {{{
2482+ def __init__(self): #{{{
2483 self.vx = np.nan
2484 self.vy = np.nan
2485 self.vz = np.nan
2486 self.vel = np.nan
2487- self.enthalpy = np.nan
2488 self.pressure = np.nan
2489 self.temperature = np.nan
2490+ self.enthalpy = np.nan
2491 self.waterfraction = np.nan
2492- self.watercolumn = np.nan
2493 self.sediment_head = np.nan
2494 self.epl_head = np.nan
2495 self.epl_thickness = np.nan
2496+ self.watercolumn = np.nan
2497 self.hydraulic_potential = np.nan
2498 self.channelarea = np.nan
2499 self.sealevel = np.nan
2500 self.bottompressure = np.nan
2501+ self.dsl = np.nan
2502+ self.str = np.nan
2503 self.sample = np.nan
2504
2505- #set defaults
2506 self.setdefaultparameters()
2507 #}}}
2508- def __repr__(self): # {{{
2509+ def __repr__(self): #{{{
2510 s = ' initial field values:\n'
2511 s += '{}\n'.format(fielddisplay(self, 'vx', 'x component of velocity [m / yr]'))
2512 s += '{}\n'.format(fielddisplay(self, 'vy', 'y component of velocity [m / yr]'))
2513@@ -51,37 +52,13 @@
2514 s += '{}\n'.format(fielddisplay(self, 'epl_thickness', 'thickness of the epl [m]'))
2515 s += '{}\n'.format(fielddisplay(self, 'hydraulic_potential', 'Hydraulic potential (for GlaDS) [Pa]'))
2516 s += '{}\n'.format(fielddisplay(self, 'channelarea', 'subglaciale water channel area (for GlaDS) [m2]'))
2517- s += '{}\n'.format(fielddisplay(self,'sample','Realization of a Gaussian random field'))
2518+ s += '{}\n'.format(fielddisplay(self,'sample', 'Realization of a Gaussian random field'))
2519 return s
2520 #}}}
2521- def extrude(self, md): # {{{
2522- self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
2523- self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
2524- self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
2525- self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
2526- self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
2527- self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
2528- self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
2529- self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
2530- self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
2531- self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
2532- self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
2533- self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
2534- self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
2535-
2536- #Lithostatic pressure by default
2537- if np.ndim(md.geometry.surface) == 2:
2538- print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
2539- self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, ) - md.mesh.z)
2540- else:
2541- self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
2542-
2543- return self
2544+ def setdefaultparameters(self): #{{{
2545+ return
2546 #}}}
2547- def setdefaultparameters(self): # {{{
2548- return self
2549- #}}}
2550- def checkconsistency(self, md, solution, analyses): # {{{
2551+ def checkconsistency(self, md, solution, analyses): #{{{
2552 if 'StressbalanceAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.isstressbalance:
2553 if not np.any(np.logical_or(np.isnan(md.initialization.vx), np.isnan(md.initialization.vy))):
2554 md = checkfield(md, 'fieldname', 'initialization.vx', 'NaN', 1, 'Inf', 1, 'size', [md.mesh.numberofvertices])
2555@@ -129,17 +106,20 @@
2556 md = checkfield(md, 'fieldname', 'initialization.channelarea', 'NaN', 1, 'Inf', 1, '>=', 0, 'size', [md.mesh.numberofelements])
2557 if 'SamplingAnalysis' in analyses and not solution == 'TransientSolution' and not md.transient.issampling:
2558 if np.any(np.isnan(md.initialization.sample)):
2559- md = checkfield(md,'fieldname','initialization.sample','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
2560+ md = checkfield(md, 'fieldname', 'initialization.sample', 'NaN', 1,'Inf', 1, 'size', [md.mesh.numberofvertices])
2561 return md
2562- # }}}
2563- def marshall(self, prefix, md, fid): # {{{
2564+ #}}}
2565+ def marshall(self, prefix, md, fid): #{{{
2566 yts = md.constants.yts
2567- WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
2568- WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
2569- WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1. / yts)
2570+
2571+ WriteData(fid, prefix, 'object', self, 'fieldname', 'vx', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
2572+ WriteData(fid, prefix, 'object', self, 'fieldname', 'vy', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
2573+ WriteData(fid, prefix, 'object', self, 'fieldname', 'vz', 'format', 'DoubleMat', 'mattype', 1, 'scale', 1 / yts)
2574 WriteData(fid, prefix, 'object', self, 'fieldname', 'pressure', 'format', 'DoubleMat', 'mattype', 1)
2575- WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1, 'timeserieslength', md.mesh.numberofvertices + 1, 'yts', md.constants.yts)
2576+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevel', 'format', 'DoubleMat', 'mattype', 1,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts)
2577 WriteData(fid, prefix, 'object', self, 'fieldname', 'bottompressure', 'format', 'DoubleMat', 'mattype', 1)
2578+ WriteData(fid, prefix, 'object', self, 'fieldname', 'str', 'format', 'DoubleMat', 'mattype', 1)
2579+ WriteData(fid, prefix, 'object', self, 'fieldname', 'dsl', 'format', 'DoubleMat', 'mattype', 1)
2580 WriteData(fid, prefix, 'object', self, 'fieldname', 'temperature', 'format', 'DoubleMat', 'mattype', 1)
2581 WriteData(fid, prefix, 'object', self, 'fieldname', 'waterfraction', 'format', 'DoubleMat', 'mattype', 1)
2582 WriteData(fid, prefix, 'object', self, 'fieldname', 'sediment_head', 'format', 'DoubleMat', 'mattype', 1)
2583@@ -148,13 +128,39 @@
2584 WriteData(fid, prefix, 'object', self, 'fieldname', 'watercolumn', 'format', 'DoubleMat', 'mattype', 1)
2585 WriteData(fid, prefix, 'object', self, 'fieldname', 'channelarea', 'format', 'DoubleMat', 'mattype', 1)
2586 WriteData(fid, prefix, 'object', self, 'fieldname', 'hydraulic_potential', 'format', 'DoubleMat', 'mattype', 1)
2587- WriteData(fid,prefix,'object',self,'fieldname','sample','format','DoubleMat','mattype',1)
2588+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sample', 'format', 'DoubleMat', 'mattype', 1)
2589+
2590 if md.thermal.isenthalpy:
2591 if (np.size(self.enthalpy) <= 1):
2592+ # Reconstruct enthalpy
2593 tpmp = md.materials.meltingpoint - md.materials.beta * md.initialization.pressure
2594- pos = np.nonzero(md.initialization.waterfraction > 0.)[0]
2595+ pos = np.where(md.initialization.waterfraction > 0)[0]
2596 self.enthalpy = md.materials.heatcapacity * (md.initialization.temperature - md.constants.referencetemperature)
2597- self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, ) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, )
2598+ self.enthalpy[pos] = md.materials.heatcapacity * (tpmp[pos].reshape(-1, 1) - md.constants.referencetemperature) + md.materials.latentheat * md.initialization.waterfraction[pos].reshape(-1, 1)
2599
2600 WriteData(fid, prefix, 'data', self.enthalpy, 'format', 'DoubleMat', 'mattype', 1, 'name', 'md.initialization.enthalpy')
2601- # }}}
2602+ #}}}
2603+ def extrude(self, md): #{{{
2604+ self.vx = project3d(md, 'vector', self.vx, 'type', 'node')
2605+ self.vy = project3d(md, 'vector', self.vy, 'type', 'node')
2606+ self.vz = project3d(md, 'vector', self.vz, 'type', 'node')
2607+ self.vel = project3d(md, 'vector', self.vel, 'type', 'node')
2608+ self.temperature = project3d(md, 'vector', self.temperature, 'type', 'node')
2609+ self.enthalpy = project3d(md, 'vector', self.enthalpy, 'type', 'node')
2610+ self.waterfraction = project3d(md, 'vector', self.waterfraction, 'type', 'node')
2611+ self.watercolumn = project3d(md, 'vector', self.watercolumn, 'type', 'node')
2612+ self.sediment_head = project3d(md, 'vector', self.sediment_head, 'type', 'node', 'layer', 1)
2613+ self.epl_head = project3d(md, 'vector', self.epl_head, 'type', 'node', 'layer', 1)
2614+ self.epl_thickness = project3d(md, 'vector', self.epl_thickness, 'type', 'node', 'layer', 1)
2615+ self.sealevel = project3d(md, 'vector', self.sealevel, 'type', 'node', 'layer', 1)
2616+ self.bottompressure = project3d(md, 'vector', self.bottompressure, 'type', 'node', 'layer', 1)
2617+
2618+ # Lithostatic pressure by default
2619+ if np.ndim(md.geometry.surface) == 2:
2620+ print('Reshaping md.geometry.surface for your convenience but you should fix it in your model set up')
2621+ self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface.reshape(-1, 1) - md.mesh.z)
2622+ else:
2623+ self.pressure = md.constants.g * md.materials.rho_ice * (md.geometry.surface - md.mesh.z)
2624+
2625+ return self
2626+ #}}}
2627Index: ../trunk-jpl/src/m/classes/lovenumbers.py
2628===================================================================
2629--- ../trunk-jpl/src/m/classes/lovenumbers.py (revision 26300)
2630+++ ../trunk-jpl/src/m/classes/lovenumbers.py (revision 26301)
2631@@ -9,8 +9,11 @@
2632 """LOVENUMBERS class definition
2633
2634 Usage:
2635- lovenumbers = lovenumbers() #will setup love numbers deg 1001 by default
2636- lovenumbers = lovenumbers('maxdeg', 10001); #supply numbers of degrees required (here, 10001)
2637+ lovenumbers = lovenumbers()
2638+ lovenumbers = lovenumbers('maxdeg', 10000, 'referenceframe', 'CF');
2639+
2640+ Choose numbers of degrees required (1000 by default) and reference frame
2641+ (between CF and CM; CM by default)
2642 """
2643
2644 def __init__(self, *args): #{{{
2645@@ -25,6 +28,10 @@
2646 self.tl = []
2647 self.tk2secular = 0 # deg 2 secular number
2648
2649+ # Time/frequency for visco-elastic love numbers
2650+ self.timefreq = []
2651+ self.istime = 1
2652+
2653 options = pairoptions(*args)
2654 maxdeg = options.getfieldvalue('maxdeg', 1000)
2655 referenceframe = options.getfieldvalue('referenceframe', 'CM')
2656@@ -31,6 +38,20 @@
2657 self.setdefaultparameters(maxdeg, referenceframe)
2658 #}}}
2659
2660+ def __repr__(self): #{{{
2661+ s = ' lovenumbers parameters:\n'
2662+ s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
2663+ s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
2664+ s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
2665+ s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
2666+ s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
2667+ s += '{}\n'.format(fielddisplay(self, 'tl', 'tidal load Love number (deg 2)'))
2668+ s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
2669+ s += '{}\n'.format(fielddisplay(self, 'istime', 'time (default: 1) or frequency love numbers (0)'))
2670+ s += '{}\n'.format(fielddisplay(self, 'timefreq', 'time/frequency vector (yr or 1/yr)'))
2671+ return s
2672+ #}}}
2673+
2674 def setdefaultparameters(self, maxdeg, referenceframe): #{{{
2675 # Initialize love numbers
2676 self.h = getlovenumbers('type', 'loadingverticaldisplacement', 'referenceframe', referenceframe, 'maxdeg', maxdeg)
2677@@ -42,11 +63,15 @@
2678
2679 # Secular fluid love number
2680 self.tk2secular = 0.942
2681+
2682+ # Time
2683+ self.istime = 1 # Temporal love numbers by default
2684+ self.timefreq = 0 # Elastic case by default
2685 return self
2686 #}}}
2687
2688 def checkconsistency(self, md, solution, analyses): #{{{
2689- if ('SealevelriseAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
2690+ if ('SealevelchangeAnalysis' not in analyses) or (solution == 'TransientSolution' and not md.transient.isslc):
2691 return
2692 md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.h', 'NaN', 1, 'Inf', 1)
2693 md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.k', 'NaN', 1, 'Inf', 1)
2694@@ -55,9 +80,17 @@
2695 md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.th', 'NaN', 1, 'Inf', 1)
2696 md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk', 'NaN', 1, 'Inf', 1)
2697 md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.tk2secular', 'NaN', 1, 'Inf', 1)
2698+ md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.timefreq', 'NaN', 1, 'Inf', 1)
2699+ md = checkfield(md, 'fieldname', 'solidearth.lovenumbers.istime', 'NaN', 1, 'Inf', 1, 'values', [0, 1])
2700+
2701 # Check that love numbers are provided at the same level of accuracy
2702 if (self.h.shape[0] != self.k.shape[0]) or (self.h.shape[0] != self.l.shape[0]):
2703 raise ValueError('lovenumbers error message: love numbers should be provided at the same level of accuracy')
2704+
2705+ ntf = len(self.timefreq)
2706+ if (self.h.shape[1] != ntf or self.k.shape[1] != ntf or self.l.shape[1] != ntf or self.th.shape[1] != ntf or self.tk.shape[1] != ntf or self.tl.shape[1] != ntf):
2707+ raise ValueError('lovenumbers error message: love numbers should have as many time/frequency steps as the time/frequency vector')
2708+
2709 return md
2710 #}}}
2711
2712@@ -65,17 +98,6 @@
2713 return[]
2714 #}}}
2715
2716- def __repr__(self): #{{{
2717- s = ' lovenumbers parameters:\n'
2718- s += '{}\n'.format(fielddisplay(self, 'h', 'load Love number for radial displacement'))
2719- s += '{}\n'.format(fielddisplay(self, 'k', 'load Love number for gravitational potential perturbation'))
2720- s += '{}\n'.format(fielddisplay(self, 'l', 'load Love number for horizontal displacements'))
2721- s += '{}\n'.format(fielddisplay(self, 'th', 'tidal load Love number (deg 2)'))
2722- s += '{}\n'.format(fielddisplay(self, 'tk', 'tidal load Love number (deg 2)'))
2723- s += '{}\n'.format(fielddisplay(self, 'tk2secular', 'secular fluid Love number'))
2724- return s
2725- #}}}
2726-
2727 def marshall(self, prefix, md, fid): #{{{
2728 WriteData(fid, prefix, 'object', self, 'fieldname', 'h', 'name', 'md.solidearth.lovenumbers.h', 'format', 'DoubleMat', 'mattype', 1)
2729 WriteData(fid, prefix, 'object', self, 'fieldname', 'k', 'name', 'md.solidearth.lovenumbers.k', 'format', 'DoubleMat', 'mattype', 1)
2730@@ -85,6 +107,13 @@
2731 WriteData(fid, prefix, 'object', self, 'fieldname', 'tk', 'name', 'md.solidearth.lovenumbers.tk', 'format', 'DoubleMat', 'mattype', 1)
2732 WriteData(fid, prefix, 'object', self, 'fieldname', 'tl', 'name', 'md.solidearth.lovenumbers.tl', 'format', 'DoubleMat', 'mattype', 1)
2733 WriteData(fid, prefix, 'object', self, 'data', self.tk2secular, 'fieldname', 'lovenumbers.tk2secular', 'format', 'Double')
2734+
2735+ if (self.istime):
2736+ scale = md.constants.yts
2737+ else:
2738+ scale = 1.0 / md.constants.yts
2739+ WriteData(fid, prefix, 'object', self, 'fieldname', 'istime', 'name', 'md.solidearth.lovenumbers.istime', 'format', 'Boolean')
2740+ WriteData(fid, prefix, 'object', self, 'fieldname', 'timefreq', 'name', 'md.solidearth.lovenumbers.timefreq', 'format', 'DoubleMat', 'mattype', 1, 'scale', scale);
2741 #}}}
2742
2743 def extrude(self, md): #{{{
2744Index: ../trunk-jpl/src/m/classes/materials.py
2745===================================================================
2746--- ../trunk-jpl/src/m/classes/materials.py (revision 26300)
2747+++ ../trunk-jpl/src/m/classes/materials.py (revision 26301)
2748@@ -24,6 +24,7 @@
2749 if not(self.nature[i] == 'litho' or self.nature[i] == 'ice' or self.nature[i] == 'hydro'):
2750 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
2751
2752+ # Start filling in the dynamic fields (not truly dynamic under Python)
2753 for i in range(len(self.nature)):
2754 nat = self.nature[i]
2755 if nat == 'ice':
2756@@ -66,7 +67,6 @@
2757 raise RuntimeError("materials constructor error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
2758 self.earth_density = 0
2759
2760- # Set default parameters
2761 self.setdefaultparameters()
2762 #}}}
2763
2764@@ -104,11 +104,11 @@
2765 s += '{}\n'.format(fielddisplay(self, 'ebm_alpha', 'array describing each layer\'s exponent parameter controlling the shape of shear modulus curve between taul and tauh, only for EBM rheology (numlayers)'))
2766 s += '{}\n'.format(fielddisplay(self, 'ebm_delta', 'array describing each layer\'s amplitude of the transient relaxation (ratio between elastic rigity to pre-maxwell relaxation rigity), only for EBM rheology (numlayers)'))
2767 s += '{}\n'.format(fielddisplay(self, 'ebm_taul', 'array describing each layer\'s starting period for transient relaxation, only for EBM rheology (numlayers) [s]'))
2768- s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'))
2769+ s += '{}\n'.format(fielddisplay(self, 'ebm_tauh', 'array describing each layer''s array describing each layer\'s end period for transient relaxation, only for Burgers rheology (numlayers) [s]'))
2770
2771 s += '{}\n'.format(fielddisplay(self, 'rheologymodel', 'array describing whether we adopt a Maxwell (0), Burgers (1) or EBM (2) rheology (default: 0)'))
2772 s += '{}\n'.format(fielddisplay(self, 'density', 'array describing each layer\'s density (numlayers) [kg/m^3]'))
2773- s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default 1) (numlayers)'))
2774+ s += '{}\n'.format(fielddisplay(self, 'issolid', 'array describing whether the layer is solid or liquid (default: 1) (numlayers)'))
2775 elif nat == 'hydro':
2776 s += 'Hydro:\n'
2777 s += '{}\n'.format(fielddisplay(self, 'rho_ice', 'ice density [kg/m^3]'))
2778@@ -126,19 +126,19 @@
2779 nat = self.nature[i]
2780 if nat == 'ice':
2781 # Ice density (kg/m^3)
2782- self.rho_ice = 917.
2783+ self.rho_ice = 917
2784
2785 # Ocean water density (kg/m^3)
2786- self.rho_water = 1023.
2787+ self.rho_water = 1023
2788
2789 # Fresh water density (kg/m^3)
2790- self.rho_freshwater = 1000.
2791+ self.rho_freshwater = 1000
2792
2793 # Water viscosity (N.s/m^2)
2794 self.mu_water = 0.001787
2795
2796 # Ice heat capacity cp (J/kg/K)
2797- self.heatcapacity = 2093.
2798+ self.heatcapacity = 2093
2799
2800 # Ice latent heat of fusion L (J/kg)
2801 self.latentheat = 3.34 * 1e5
2802@@ -159,7 +159,7 @@
2803 self.beta = 9.8 * 1e-8
2804
2805 # Mixed layer (ice-water interface) heat capacity (J/kg/K)
2806- self.mixed_layer_capacity = 3974.
2807+ self.mixed_layer_capacity = 3974
2808
2809 # Thermal exchange velocity (ice-water interface) (m/s)
2810 self.thermal_exchange_velocity = 1.00 * 1e-4
2811@@ -192,17 +192,17 @@
2812 self.ebm_taul = [np.nan, np.nan]
2813 self.ebm_tauh = [np.nan, np.nan]
2814 self.rheologymodel = [0, 0]
2815- self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3]
2816- self.issolid = [1, 1] # Is layer solid or liquid?
2817+ self.density = [5.51e3, 5.50e3] # (Pa) # Mantle and lithosphere density [kg/m^3]
2818+ self.issolid = [1, 1] # Is layer solid or liquid?
2819 elif nat == 'hydro':
2820 # Ice density (kg/m^3)
2821- self.rho_ice = 917.
2822+ self.rho_ice = 917
2823
2824 # Ocean water density (kg/m^3)
2825- self.rho_water = 1023.
2826+ self.rho_water = 1023
2827
2828 # Fresh water density (kg/m^3)
2829- self.rho_freshwater = 1000.
2830+ self.rho_freshwater = 1000
2831 else:
2832 raise RuntimeError("materials setdefaultparameters error message: nature of the material not supported yet! ('ice' or 'litho' or 'hydro')")
2833
2834@@ -250,7 +250,7 @@
2835 raise RuntimeError('First layer must be solid (issolid[0] > 0 AND lame_mu[0] > 0). Add a weak inner core if necessary.')
2836 ind = np.where(md.materials.issolid == 0)[0]
2837 if np.sum(np.in1d(np.diff(ind),1) >= 1): # If there are at least two consecutive indices that contain issolid = 0
2838- raise RuntimeError('Fluid layers detected at layers #{} but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.'.format(i))
2839+ raise RuntimeError('Fluid layers detected at layers #' + indices + ', but having 2 or more adjacent fluid layers is not supported yet. Consider merging them.')
2840
2841 elif nat == 'hydro':
2842 md = checkfield(md, 'fieldname', 'materials.rho_ice', '>', 0)
2843@@ -266,7 +266,7 @@
2844 def marshall(self, prefix, md, fid): #{{{
2845 #1: MatdamageiceEnum 2: MatestarEnum 3: MaticeEnum 4: MatenhancediceEnum 5: MaterialsEnum
2846 WriteData(fid, prefix, 'name', 'md.materials.nature', 'data', naturetointeger(self.nature), 'format', 'IntMat', 'mattype', 3)
2847- WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER, this can evolve if you have classes
2848+ WriteData(fid, prefix, 'name', 'md.materials.type', 'data', 5, 'format', 'Integer') #DANGER: this can evolve if you have classes
2849 for i in range(len(self.nature)):
2850 nat = self.nature[i]
2851 if nat == 'ice':
2852@@ -306,6 +306,7 @@
2853 for i in range(self.numlayers):
2854 earth_density = earth_density + (pow(self.radius[i + 1], 3) - pow(self.radius[i], 3)) * self.density[i]
2855 earth_density = earth_density / pow(self.radius[self.numlayers], 3)
2856+ print(earth_density)
2857 self.earth_density = earth_density
2858 elif nat == 'hydro':
2859 WriteData(fid, prefix, 'object', self, 'class', 'materials', 'fieldname', 'rho_ice', 'format', 'Double')
2860Index: ../trunk-jpl/src/m/classes/solidearth.m
2861===================================================================
2862--- ../trunk-jpl/src/m/classes/solidearth.m (revision 26300)
2863+++ ../trunk-jpl/src/m/classes/solidearth.m (revision 26301)
2864@@ -44,6 +44,24 @@
2865 error('solidearth constructor error message: zero or one argument only!');
2866 end
2867 end % }}}
2868+ function disp(self) % {{{
2869+ disp(sprintf(' solidearth inputs, forcings and settings:'));
2870+
2871+ fielddisplay(self,'planetradius','planet radius [m]');
2872+ fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
2873+ fielddisplay(self,'requested_outputs','additional outputs requested');
2874+ fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
2875+ fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
2876+ fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
2877+ if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
2878+ self.settings.disp();
2879+ self.lovenumbers.disp();
2880+ self.rotational.disp();
2881+ if ~isempty(self.external),
2882+ self.external.disp();
2883+ end
2884+
2885+ end % }}}
2886 function self = setdefaultparameters(self,planet) % {{{
2887
2888 %output default:
2889@@ -86,24 +104,6 @@
2890 function list=defaultoutputs(self,md) % {{{
2891 list = {'Sealevel'};
2892 end % }}}
2893- function disp(self) % {{{
2894- disp(sprintf(' solidearth inputs, forcings and settings:'));
2895-
2896- fielddisplay(self,'planetradius','planet radius [m]');
2897- fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps');
2898- fielddisplay(self,'requested_outputs','additional outputs requested');
2899- fielddisplay(self,'partitionice','ice partition vector for barystatic contribution');
2900- fielddisplay(self,'partitionhydro','hydro partition vector for barystatic contribution');
2901- fielddisplay(self,'partitionocean','ocean partition vector for barystatic contribution');
2902- if isempty(self.external), fielddisplay(self,'external','external solution, of the type solidearthsolution'); end
2903- self.settings.disp();
2904- self.lovenumbers.disp();
2905- self.rotational.disp();
2906- if ~isempty(self.external),
2907- self.external.disp();
2908- end
2909-
2910- end % }}}
2911 function marshall(self,prefix,md,fid) % {{{
2912
2913 WriteData(fid,prefix,'object',self,'fieldname','planetradius','format','Double');
2914Index: ../trunk-jpl/src/m/classes/solidearthsettings.py
2915===================================================================
2916--- ../trunk-jpl/src/m/classes/solidearthsettings.py (revision 26300)
2917+++ ../trunk-jpl/src/m/classes/solidearthsettings.py (revision 26301)
2918@@ -120,9 +120,9 @@
2919 if md.mesh.__class__.__name__ is 'mesh3dsurface':
2920 if self.grdmodel == 2:
2921 raise RuntimeException('model requires a 2D mesh to run gia Ivins computations (change mesh from mesh3dsurface to mesh2d)')
2922- else:
2923- if self.grdmodel == 1:
2924- raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
2925+ else:
2926+ if self.grdmodel == 1:
2927+ raise RuntimeException('model requires a 3D surface mesh to run GRD computations (change mesh from mesh2d to mesh3dsurface)')
2928 if self.sealevelloading and not self.grdocean:
2929 raise RuntimeException('solidearthsettings checkconsistency error message: need grdocean on if sealevelloading flag is set')
2930
2931@@ -133,24 +133,24 @@
2932 #}}}
2933
2934 def marshall(self, prefix, md, fid): #{{{
2935- WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double')
2936- WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double')
2937- WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer')
2938- WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean')
2939- WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean')
2940- WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean')
2941- WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean')
2942- WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean')
2943- WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Integer')
2944- WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer')
2945- WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double')
2946- WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale', md.constants.yts)
2947- WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer')
2948- WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer')
2949- WriteData(fid, prefix, 'object', self, 'fieldname','isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer')
2950- WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer')
2951- WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer')
2952- WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer')
2953+ WriteData(fid, prefix, 'object', self, 'fieldname', 'reltol', 'name', 'md.solidearth.settings.reltol', 'format', 'Double');
2954+ WriteData(fid, prefix, 'object', self, 'fieldname', 'abstol', 'name', 'md.solidearth.settings.abstol', 'format', 'Double');
2955+ WriteData(fid, prefix, 'object', self, 'fieldname', 'maxiter', 'name', 'md.solidearth.settings.maxiter', 'format', 'Integer');
2956+ WriteData(fid, prefix, 'object', self, 'fieldname', 'selfattraction', 'name', 'md.solidearth.settings.selfattraction', 'format', 'Boolean');
2957+ WriteData(fid, prefix, 'object', self, 'fieldname', 'elastic', 'name', 'md.solidearth.settings.elastic', 'format', 'Boolean');
2958+ WriteData(fid, prefix, 'object', self, 'fieldname', 'viscous', 'name', 'md.solidearth.settings.viscous', 'format', 'Boolean');
2959+ WriteData(fid, prefix, 'object', self, 'fieldname', 'rotation', 'name', 'md.solidearth.settings.rotation', 'format', 'Boolean');
2960+ WriteData(fid, prefix, 'object', self, 'fieldname', 'grdocean', 'name', 'md.solidearth.settings.grdocean', 'format', 'Boolean');
2961+ WriteData(fid, prefix, 'object', self, 'fieldname', 'ocean_area_scaling', 'name', 'md.solidearth.settings.ocean_area_scaling', 'format', 'Boolean');
2962+ WriteData(fid, prefix, 'object', self, 'fieldname', 'runfrequency', 'name', 'md.solidearth.settings.runfrequency', 'format', 'Integer');
2963+ WriteData(fid, prefix, 'object', self, 'fieldname', 'degacc', 'name', 'md.solidearth.settings.degacc', 'format', 'Double');
2964+ WriteData(fid, prefix, 'object', self, 'fieldname', 'timeacc', 'name', 'md.solidearth.settings.timeacc', 'format', 'Double', 'scale',md.constants.yts);
2965+ WriteData(fid, prefix, 'object', self, 'fieldname', 'horiz', 'name', 'md.solidearth.settings.horiz', 'format', 'Integer');
2966+ WriteData(fid, prefix, 'object', self, 'fieldname', 'sealevelloading', 'name', 'md.solidearth.settings.sealevelloading', 'format', 'Integer');
2967+ WriteData(fid, prefix, 'object', self, 'fieldname', 'isgrd', 'name', 'md.solidearth.settings.isgrd', 'format', 'Integer');
2968+ WriteData(fid, prefix, 'object', self, 'fieldname', 'compute_bp_grd', 'name', 'md.solidearth.settings.compute_bp_grd', 'format', 'Integer');
2969+ WriteData(fid, prefix, 'object', self, 'fieldname', 'grdmodel', 'name', 'md.solidearth.settings.grdmodel', 'format', 'Integer');
2970+ WriteData(fid, prefix, 'object', self, 'fieldname', 'cross_section_shape', 'name', 'md.solidearth.settings.cross_section_shape', 'format', 'Integer');
2971 #}}}
2972
2973 def extrude(self, md): #{{{
2974Index: ../trunk-jpl/src/m/classes/solidearthsolution.py
2975===================================================================
2976--- ../trunk-jpl/src/m/classes/solidearthsolution.py (nonexistent)
2977+++ ../trunk-jpl/src/m/classes/solidearthsolution.py (revision 26301)
2978@@ -0,0 +1,80 @@
2979+import numpy as np
2980+
2981+from checkfield import checkfield
2982+from fielddisplay import fielddisplay
2983+from WriteData import WriteData
2984+
2985+
2986+class solidearthsolution(object):
2987+ """SOLIDEARTHSOLUTION class definition
2988+
2989+ Usage:
2990+ solidearthsolution = solidearthsolution()
2991+ """
2992+
2993+ def __init__(self, *args): #{{{
2994+ self.displacementeast = None
2995+ self.displacementnorth = None
2996+ self.displacementup = None
2997+ self.geoid = None
2998+
2999+ if len(args) == 0:
3000+ self.setdefaultparameters()
3001+ else:
3002+ raise RuntimeError('constructor not supported')
3003+ #}}}
3004+
3005+ def __repr__(self): #{{{
3006+ s = ' units for time series is (yr)\n'
3007+ s += '{}\n'.format(fielddisplay(self, 'displacementeast', 'solid-Earth Eastwards bedrock displacement series (m)'))
3008+ s += '{}\n'.format(fielddisplay(self, 'displacementnorth', 'solid-Earth Northwards bedrock displacement time series (m)'))
3009+ s += '{}\n'.format(fielddisplay(self, 'displacementup', 'solid-Earth bedrock uplift time series (m)'))
3010+ s += '{}\n'.format(fielddisplay(self, 'geoid', 'solid-Earth geoid time series (m)'))
3011+
3012+ return s
3013+ #}}}
3014+
3015+ def setdefaultparameters(self): #{{{
3016+ self.displacementeast = []
3017+ self.displacementnorth = []
3018+ self.displacementup = []
3019+ self.geoid = []
3020+ #}}}
3021+
3022+ def checkconsistency(self, md, solution, analyses): #{{{
3023+ md = checkfield(md, 'fieldname', 'solidearth.external.displacementeast', 'Inf', 1, 'timeseries', 1)
3024+ md = checkfield(md, 'fieldname', 'solidearth.external.displacementnorth', 'Inf', 1, 'timeseries', 1)
3025+ md = checkfield(md, 'fieldname', 'solidearth.external.displacementup', 'Inf', 1, 'timeseries', 1)
3026+ md = checkfield(md, 'fieldname', 'solidearth.external.geoid', 'Inf', 1, 'timeseries', 1)
3027+ #}}}
3028+
3029+ def marshall(self, prefix, md, fid): #{{{
3030+ yts = md.constants.yts
3031+
3032+ # Transform our time series into time series rates
3033+ if np.shape(self.displacementeast, 1) == 1:
3034+ print('External solidearthsolution warning: only one time step provided, assuming the values are rates per year')
3035+ displacementeast_rate = np.append(np.array(self.displacementeast).reshape(-1, 1), 0)
3036+ displacementnorth_rate = np.append(np.array(self.displacementnorth).reshape(-1, 1), 0)
3037+ displacementup_rate = np.append(np.array(self.displacementup).reshape(-1, 1), 0)
3038+ geoid_rate = np.append(np.array(self.geoid).reshape(-1, 1), 0)
3039+ else:
3040+ time = self.displacementeast[-1, :]
3041+ dt = np.diff(time, 1, 1)
3042+ displacementeast_rate = np.diff(self.displacementeast[0:-2, :], 1, 1) / dt
3043+ displacementeast_rate.append(time[0:-2])
3044+ displacementnorth_rate = np.diff(self.displacementnorth[0:-2, :], 1, 1) / dt
3045+ displacementnorth_rate.append(time[0:-2])
3046+ displacementup_rate = np.diff(self.displacementup[0:-2, :], 1, 1) / dt
3047+ displacementup_rate.append(time[0:-2])
3048+ geoid_rate = np.diff(self.geoid[0:-2, :], 1, 1) / dt
3049+ geoid_rate.append(time[0:-2])
3050+ WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementeast', 'data', displacementeast_rate, 'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementeast', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
3051+ WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementup', 'data', displacementup_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementup', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
3052+ WriteData(fid, prefix, 'object', self, 'fieldname', 'displacementnorth', 'data', displacementnorth_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.displacementnorth', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
3053+ WriteData(fid, prefix, 'object', self, 'fieldname', 'geoid', 'data', geoid_rate,'format', 'DoubleMat', 'name', 'md.solidearth.external.geoid', 'mattype', 1, 'scale', 1 / yts,'timeserieslength', md.mesh.numberofvertices + 1, 'yts', yts);
3054+ #}}}
3055+
3056+ def extrude(self, md): #{{{
3057+ return self
3058+ #}}}
3059Index: ../trunk-jpl/test/NightlyRun/test2002.m
3060===================================================================
3061--- ../trunk-jpl/test/NightlyRun/test2002.m (revision 26300)
3062+++ ../trunk-jpl/test/NightlyRun/test2002.m (revision 26301)
3063@@ -14,6 +14,7 @@
3064 %solidearth loading: {{{
3065 md.masstransport.spcthickness=[md.geometry.thickness;0];
3066 md.smb.mass_balance=zeros(md.mesh.numberofvertices,1);
3067+
3068 %antarctica
3069 xe=md.mesh.x(md.mesh.elements)*[1;1;1]/3;
3070 ye=md.mesh.y(md.mesh.elements)*[1;1;1]/3;
3071@@ -25,6 +26,7 @@
3072 pos=find(late < -80);
3073 md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
3074 posant=pos;
3075+
3076 %greenland
3077 pos=find(late>60 & late<90 & longe>-75 & longe<-15);
3078 md.masstransport.spcthickness(md.mesh.elements(pos,:))= md.masstransport.spcthickness(md.mesh.elements(pos,:))-100;
3079@@ -95,6 +97,7 @@
3080 md=solve(md,'Transient');
3081 Seustatic=md.results.TransientSolution.Sealevel;
3082 Beustatic=md.results.TransientSolution.Bed;
3083+pause
3084
3085 %eustatic + selfattraction run:
3086 md.solidearth.settings.selfattraction=1;
3087Index: ../trunk-jpl/test/NightlyRun/test2002.py
3088===================================================================
3089--- ../trunk-jpl/test/NightlyRun/test2002.py (revision 26300)
3090+++ ../trunk-jpl/test/NightlyRun/test2002.py (revision 26301)
3091@@ -1,6 +1,8 @@
3092 #Test Name: EarthSlc
3093 import numpy as np
3094
3095+from MatlabFuncs import *
3096+
3097 from gmshplanet import *
3098 from gmtmask import *
3099 from lovenumbers import *
3100@@ -11,98 +13,137 @@
3101 from solve import *
3102
3103
3104-#mesh earth:
3105+# Mesh earth
3106+#
3107+# NOTE: In MATLAB, we currently use cached mesh to account for differences in
3108+# mesh generated under Linux versus under macOS
3109+#
3110 md = model()
3111-md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) #700 km resolution mesh
3112+md.mesh = gmshplanet('radius', 6.371012 * 1e3, 'resolution', 700.) # 700 km resolution mesh
3113
3114-#parameterize solidearth solution:
3115-#solidearth loading:
3116-md.solidearth.surfaceload.icethicknesschange = np.zeros((md.mesh.numberofelements, 1))
3117-md.solidearth.initialsealevel = np.zeros((md.mesh.numberofvertices, 1))
3118-md.dsl.global_average_thermosteric_sea_level_change = np.zeros((2, 1))
3119-md.dsl.sea_surface_height_change_above_geoid = np.zeros((md.mesh.numberofvertices + 1, 1))
3120-md.dsl.sea_water_pressure_change_at_sea_floor = np.zeros((md.mesh.numberofvertices + 1, 1))
3121+# Geometry for the bed, arbitrary thickness of 100
3122+md.geometry.bed = np.zeros((md.mesh.numberofvertices, 1))
3123+md.geometry.base = md.geometry.bed
3124+md.geometry.thickness = 100 * np.ones((md.mesh.numberofvertices, 1))
3125+md.geometry.surface = md.geometry.bed + md.geometry.thickness
3126
3127-#antarctica
3128-late = md.mesh.lat[md.mesh.elements - 1].sum(axis=1) / 3
3129-longe = md.mesh.long[md.mesh.elements - 1].sum(axis=1) / 3
3130+# Solidearth loading #{{{
3131+md.masstransport.spcthickness = np.append(md.geometry.thickness, 0)
3132+md.smb.mass_balance = np.zeros((md.mesh.numberofvertices, 1))
3133+
3134+# Antarctica
3135+xe = md.mesh.x[md.mesh.elements - 1].sum(axis=1) / 3
3136+ye = md.mesh.y[md.mesh.elements - 1].sum(axis=1) / 3
3137+ze = md.mesh.z[md.mesh.elements - 1].sum(axis=1) / 3
3138+re = pow((pow(xe, 2) + pow(ye, 2) + pow(ze, 2)), 0.5)
3139+
3140+late = asind(ze / re)
3141+longe = atan2d(ye, xe)
3142 pos = np.where(late < -80)[0]
3143-md.solidearth.surfaceload.icethicknesschange[pos] = -100
3144-#greenland
3145-pos = np.where(np.logical_and.reduce((late > 70, late < 80, longe > -60, longe < -30)))[0]
3146-md.solidearth.surfaceload.icethicknesschange[pos] = -100
3147+md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
3148+posant = pos
3149
3150-#elastic loading from love numbers:
3151-md.solidearth.lovenumbers = lovenumbers('maxdeg', 100)
3152+# Greenland
3153+pos = np.where(np.logical_and.reduce((late > 60, late < 90, longe > -75, longe < -15)))[0]
3154+md.masstransport.spcthickness[md.mesh.elements[pos]] = md.masstransport.spcthickness[md.mesh.elements[pos]] - 100
3155+posgre = pos
3156+
3157+# Elastic loading from love numbers:
3158+md.solidearth.lovenumbers = lovenumbers('maxdeg', 1000)
3159 #}}}
3160
3161-#mask: {{{
3162+# Mask: {{{
3163 mask = gmtmask(md.mesh.lat, md.mesh.long)
3164-icemask = np.ones((md.mesh.numberofvertices, 1))
3165+oceanmask = -1 * np.ones((md.mesh.numberofvertices, 1))
3166 pos = np.where(mask == 0)[0]
3167-icemask[pos] = -1
3168-pos = np.where(mask[md.mesh.elements - 1].sum(axis=1) < 3)[0]
3169-icemask[md.mesh.elements[pos, :] - 1] = -1
3170-md.mask.ice_levelset = icemask
3171-md.mask.ocean_levelset = -icemask
3172+oceanmask[pos] = 1
3173
3174-#make sure that the elements that have loads are fully grounded
3175-pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0]
3176-md.mask.ocean_levelset[md.mesh.elements[pos, :] - 1] = 1
3177+icemask = np.ones((md.mesh.numberofvertices, 1))
3178+icemask[md.mesh.elements[posant]] = -1
3179+icemask[md.mesh.elements[posgre]] = -1
3180
3181-#make sure wherever there is an ice load, that the mask is set to ice:
3182-#pos = np.nonzero(md.solidearth.surfaceload.icethicknesschange)[0] # TODO: Do we need to do this twice?
3183-md.mask.ice_levelset[md.mesh.elements[pos, :] - 1] = -1
3184-# }}}
3185+md.mask.ice_levelset = icemask
3186+md.mask.ocean_levelset = oceanmask
3187+#}}}
3188
3189-md.solidearth.settings.ocean_area_scaling = 0
3190+# Time stepping {{{
3191+md.timestepping.start_time = 0
3192+md.timestepping.time_step = 1
3193+md.timestepping.final_time = 1
3194+#}}}
3195
3196-#geometry for the bed, arbitrary
3197-md.geometry.bed = -np.ones((md.mesh.numberofvertices, 1))
3198+# Masstransport
3199+md.basalforcings.groundedice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
3200+md.basalforcings.floatingice_melting_rate = np.zeros((md.mesh.numberofvertices, 1))
3201+md.initialization.vx = np.zeros((md.mesh.numberofvertices, 1))
3202+md.initialization.vy = np.zeros((md.mesh.numberofvertices, 1))
3203+md.initialization.sealevel = np.zeros((md.mesh.numberofvertices, 1))
3204+md.initialization.str = 0
3205
3206-#materials
3207-md.materials=materials('hydro')
3208+# Materials
3209+md.materials = materials('hydro')
3210
3211-#Miscellaneous
3212+# Miscellaneous
3213 md.miscellaneous.name = 'test2002'
3214
3215-#Solution parameters
3216+# Solution parameters
3217+md.cluster.np = 3
3218 md.solidearth.settings.reltol = np.nan
3219 md.solidearth.settings.abstol = 1e-3
3220-md.solidearth.settings.computesealevelchange = 1
3221+md.solidearth.settings.sealevelloading = 1
3222+md.solidearth.settings.isgrd = 1
3223+md.solidearth.settings.ocean_area_scaling = 0
3224+md.solidearth.settings.grdmodel = 1
3225
3226-#max number of iterations reverted back to 10 (i.e., the original default value)
3227+# Physics
3228+md.transient.issmb = 0
3229+md.transient.isstressbalance = 0
3230+md.transient.isthermal = 0
3231+md.transient.ismasstransport = 1
3232+md.transient.isslc = 1
3233+md.solidearth.requested_outputs = ['Sealevel', 'Bed']
3234+
3235+# Max number of iterations reverted back to 10 (i.e., the original default value)
3236 md.solidearth.settings.maxiter = 10
3237
3238-#eustatic run:
3239-md.solidearth.settings.rigid = 0
3240+# Eustatic run
3241+md.solidearth.settings.selfattraction = 0
3242 md.solidearth.settings.elastic = 0
3243 md.solidearth.settings.rotation = 0
3244-md = solve(md, 'Sealevelrise')
3245-Seustatic = md.results.SealevelriseSolution.Sealevel
3246+md.solidearth.settings.viscous = 0
3247
3248-#eustatic + rigid run:
3249-md.solidearth.settings.rigid = 1
3250+md = solve(md, 'Transient')
3251+Seustatic = md.results.TransientSolution.Sealevel
3252+Beustatic = md.results.TransientSolution.Bed
3253+
3254+# Eustatic + selfattraction run
3255+md.solidearth.settings.selfattraction = 1
3256 md.solidearth.settings.elastic = 0
3257 md.solidearth.settings.rotation = 0
3258-md = solve(md, 'Sealevelrise')
3259-Srigid = md.results.SealevelriseSolution.Sealevel
3260+md.solidearth.settings.viscous = 0
3261+md = solve(md, 'tr')
3262+Sselfattraction = md.results.TransientSolution.Sealevel
3263+Bselfattraction = md.results.TransientSolution.Bed
3264
3265-#eustatic + rigid + elastic run:
3266-md.solidearth.settings.rigid = 1
3267+# Eustatic + selfattraction + elastic run
3268+md.solidearth.settings.selfattraction = 1
3269 md.solidearth.settings.elastic = 1
3270 md.solidearth.settings.rotation = 0
3271-md = solve(md, 'Sealevelrise')
3272-Selastic = md.results.SealevelriseSolution.Sealevel
3273+md.solidearth.settings.viscous = 0
3274+md = solve(md, 'tr')
3275+Selastic = md.results.TransientSolution.Sealevel
3276+Belastic = md.results.TransientSolution.Bed
3277
3278-#eustatic + rigid + elastic + rotation run:
3279-md.solidearth.settings.rigid = 1
3280+# Eustatic + selfattraction + elastic + rotation run
3281+md.solidearth.settings.selfattraction = 1
3282 md.solidearth.settings.elastic = 1
3283 md.solidearth.settings.rotation = 1
3284-md = solve(md, 'Sealevelrise')
3285-Srotation = md.results.SealevelriseSolution.Sealevel
3286+md.solidearth.settings.viscous = 0
3287+md = solve(md, 'tr')
3288+Srotation = md.results.TransientSolution.Sealevel
3289+Brotation = md.results.TransientSolution.Bed
3290
3291-#Fields and tolerances to track changes
3292-field_names = ['Eustatic', 'Rigid', 'Elastic', 'Rotation']
3293-field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13]
3294-field_values = [Seustatic, Srigid, Selastic, Srotation]
3295+# Fields and tolerances to track changes
3296+field_names = ['Seustatic', 'Sselfattraction', 'Selastic', 'Srotation', 'Beustatic', 'Bselfattraction', 'Belastic', 'Brotation']
3297+field_tolerances = [1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13, 1e-13]
3298+field_values = [Seustatic, Sselfattraction, Selastic, Srotation, Beustatic, Bselfattraction, Belastic, Brotation]
Note: See TracBrowser for help on using the repository browser.