source: issm/trunk-jpl/src/py3/classes/slr.py@ 23677

Last change on this file since 23677 was 23677, checked in by bdef, 6 years ago

CHG: adding missing directories and cleaning code
File size: 11.9 KB
Line 
1from fielddisplay import fielddisplay
2from MatlabFuncs import *
3from model import *
4import numpy as np
5from checkfield import checkfield
6from WriteData import WriteData
7
8class slr(object):
9 """
10 SLR class definition
11
12 Usage:
13 slr=slr()
14 """
15 def __init__(self): # {{{
16 self.deltathickness = float('NaN')
17 self.sealevel = float('NaN')
18 self.spcthickness = float('NaN')
19 self.maxiter = 0
20 self.reltol = 0
21 self.abstol = 0
22 self.love_h = 0 #provided by PREM model()
23 self.love_k = 0 #ideam
24 self.love_l = 0 #ideam
25 self.tide_love_k = 0 #ideam
26 self.tide_love_h = 0 #ideam
27 self.fluid_love = 0
28 self.equatorial_moi = 0
29 self.polar_moi = 0
30 self.angular_velocity = 0
31 self.rigid = 0
32 self.elastic = 0
33 self.rotation = 0
34 self.ocean_area_scaling = 0
35 self.steric_rate = 0 #rate of ocean expansion from steric effects.
36 self.geodetic_run_frequency = 1 #how many time steps we skip before we run the geodetic part of the solver during transient
37 self.geodetic = 0 #compute geodetic SLR? (in addition to steric?)
38 self.degacc = 0
39 self.loop_increment = 0
40 self.horiz = 0
41 self.Ngia = float('NaN')
42 self.Ugia = float('NaN')
43 self.requested_outputs = []
44 self.transitions = []
45
46 #set defaults
47 self.setdefaultparameters()
48 #}}}
49
50 def __repr__(self): # {{{
51 string=' slr parameters:'
52 string="%s\n%s"%(string,fielddisplay(self,'deltathickness','thickness change: ice height equivalent [m]'))
53 string="%s\n%s"%(string,fielddisplay(self,'sealevel','current sea level (prior to computation) [m]'))
54 string="%s\n%s"%(string,fielddisplay(self,'spcthickness','thickness constraints (NaN means no constraint) [m]'))
55 string="%s\n%s"%(string,fielddisplay(self,'reltol','sea level rise relative convergence criterion, (NaN: not applied)'))
56 string="%s\n%s"%(string,fielddisplay(self,'abstol','sea level rise absolute convergence criterion, (default, NaN: not applied)'))
57 string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations'))
58 string="%s\n%s"%(string,fielddisplay(self,'love_h','load Love number for radial displacement'))
59 string="%s\n%s"%(string,fielddisplay(self,'love_k','load Love number for gravitational potential perturbation'))
60 string="%s\n%s"%(string,fielddisplay(self,'love_l','load Love number for horizontal displaements'))
61 string="%s\n%s"%(string,fielddisplay(self,'tide_love_k','tidal load Love number (degree 2)'))
62 string="%s\n%s"%(string,fielddisplay(self,'tide_love_h','tidal load Love number (degree 2)'))
63 string="%s\n%s"%(string,fielddisplay(self,'fluid_love','secular fluid Love number'))
64 string="%s\n%s"%(string,fielddisplay(self,'equatorial_moi','mean equatorial moment of inertia [kg m^2]'))
65 string="%s\n%s"%(string,fielddisplay(self,'polar_moi','polar moment of inertia [kg m^2]'))
66 string="%s\n%s"%(string,fielddisplay(self,'angular_velocity','mean rotational velocity of earth [per second]'))
67 string="%s\n%s"%(string,fielddisplay(self,'ocean_area_scaling','correction for model representation of ocean area [default: No correction]'))
68 string="%s\n%s"%(string,fielddisplay(self,'steric_rate','rate of steric ocean expansion [mm/yr]'))
69 string="%s\n%s"%(string,fielddisplay(self,'Ngia','rate of viscous (GIA) geoid expansion (in mm/yr)'))
70 string="%s\n%s"%(string,fielddisplay(self,'Ugia','rate of viscous (GIA) bedrock uplift (in mm/yr)'))
71 string="%s\n%s"%(string,fielddisplay(self,'loop_increment','vector assembly (in the convolution) framentation'))
72 string="%s\n%s"%(string,fielddisplay(self,'geodetic','compute geodetic SLR? ( in addition to steric?) default 0'))
73 string="%s\n%s"%(string,fielddisplay(self,'geodetic_run_frequency','how many time steps we skip before we run SLR solver during transient (default: 1)'))
74 string="%s\n%s"%(string,fielddisplay(self,'rigid','rigid earth graviational potential perturbation'))
75 string="%s\n%s"%(string,fielddisplay(self,'elastic','elastic earth graviational potential perturbation'))
76 string="%s\n%s"%(string,fielddisplay(self,'rotation','earth rotational potential perturbation'))
77 string="%s\n%s"%(string,fielddisplay(self,'degacc','accuracy (default .01 deg) for numerical discretization of the Green''s functions'))
78 string="%s\n%s"%(string,fielddisplay(self,'transitions','indices into parts of the mesh that will be icecaps'))
79 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
80
81 return string
82 # }}}
83
84 def setdefaultparameters(self): # {{{
85 #Convergence criterion: absolute, relative and residual
86 self.reltol = 0.01 #default
87 self.abstol = float('NaN') #1 mm of sea level rise
88
89 #maximum of non-linear iterations.
90 self.maxiter = 5
91 self.loop_increment = 200
92
93 #computational flags:
94 self.geodetic = 0
95 self.rigid = 1
96 self.elastic = 1
97 self.ocean_area_scaling = 0
98 self.rotation = 1
99
100 #tidal love numbers:
101 self.tide_love_h = 0.6149 #degree 2
102 self.tide_love_k = 0.3055 #degree 2
103
104 #secular fluid love number:
105 self.fluid_love = 0.942
106
107 #moment of inertia:
108 self.equatorial_moi = 8.0077*10**37 # [kg m^2]
109 self.polar_moi = 8.0345*10**37 # [kg m^2]
110
111 #mean rotational velocity of earth
112 self.angular_velocity = 7.2921*10**-5 # [s^-1]
113
114 #numerical discretization accuracy
115 self.degacc = .01
116
117 #steric:
118 self.steric_rate = 0
119
120 #how many time steps we skip before we run SLR solver during transient
121 self.geodetic_run_frequency = 1
122
123 #output default:
124 self.requested_outputs = ['default']
125
126 #transitions should be a cell array of vectors:
127 self.transitions = []
128
129 #horizontal displacement? (not by default)
130 self.horiz = 0
131
132 return self
133 #}}}
134
135 def checkconsistency(self,md,solution,analyses): # {{{
136 #Early return
137 if (solution!='SealevelriseAnalysis'):
138 return md
139
140 md = checkfield(md,'fieldname','slr.deltathickness','NaN',1,'Inf',1,'size',[md.mesh.numberofelements])
141 md = checkfield(md,'fieldname','slr.sealevel','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
142 md = checkfield(md,'fieldname','slr.spcthickness','Inf',1,'timeseries',1)
143 md = checkfield(md,'fieldname','slr.love_h','NaN',1,'Inf',1)
144 md = checkfield(md,'fieldname','slr.love_k','NaN',1,'Inf',1)
145 md = checkfield(md,'fieldname','slr.love_l','NaN',1,'Inf',1)
146 md = checkfield(md,'fieldname','slr.tide_love_h','NaN',1,'Inf',1)
147 md = checkfield(md,'fieldname','slr.tide_love_k','NaN',1,'Inf',1)
148 md = checkfield(md,'fieldname','slr.fluid_love','NaN',1,'Inf',1)
149 md = checkfield(md,'fieldname','slr.equatorial_moi','NaN',1,'Inf',1)
150 md = checkfield(md,'fieldname','slr.polar_moi','NaN',1,'Inf',1)
151 md = checkfield(md,'fieldname','slr.angular_velocity','NaN',1,'Inf',1)
152 md = checkfield(md,'fieldname','slr.reltol','size',[1,1])
153 md = checkfield(md,'fieldname','slr.abstol','size',[1,1])
154 md = checkfield(md,'fieldname','slr.maxiter','size',[1,1],'>=',1)
155 md = checkfield(md,'fieldname','slr.geodetic_run_frequency','size',[1,1],'>=',1)
156 md = checkfield(md,'fieldname','slr.steric_rate','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
157 md = checkfield(md,'fieldname','slr.degacc','size',[1,1],'>=',1e-10)
158 md = checkfield(md,'fieldname','slr.requested_outputs','stringrow',1)
159 md = checkfield(md,'fieldname','slr.loop_increment','NaN',1,'Inf',1,'>=',1)
160 md = checkfield(md,'fieldname','slr.horiz','NaN',1,'Inf',1,'values',[0,1])
161 md = checkfield(md,'fieldname','slr.Ngia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
162 md = checkfield(md,'fieldname','slr.Ugia','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices])
163
164 #check that love numbers are provided at the same level of accuracy:
165 if (size(self.love_h,0) != size(self.love_k,0) | size(self.love_h,0) != size(self.love_l,0)):
166 error('slr error message: love numbers should be provided at the same level of accuracy')
167
168 #cross check that whereever we have an ice load, the mask is <0 on each vertex:
169 pos=np.where(self.deltathickness)
170 maskpos=md.mask.ice_levelset[md.mesh.elements[pos,:]]
171 els=np.where(maskpos>0)
172 if len(els[0])>0:
173 warnings.warn('slr checkconsistency fail: there are elements with ice loads where some vertices are not on the ice!')
174
175 #check that if geodetic is requested, we are a mesh3dsurface model (planet), or if we are not,
176 #a coupler to a planet model is provided.
177 if self.geodetic and not md.transient.iscoupler and domaintype(md.mesh)!='mesh3dsurface':
178 error('model is requesting geodetic computations without being a mesh3dsurface, or being coupled to one!')
179 return md
180 # }}}
181
182 def defaultoutputs(self,md): # {{{
183 return ['Sealevel']
184 # }}}
185
186 def marshall(self,prefix,md,fid): # {{{
187 WriteData(fid,prefix,'object',self,'fieldname','deltathickness','format','DoubleMat','mattype',2)
188 WriteData(fid,prefix,'object',self,'fieldname','sealevel','mattype',1,'format','DoubleMat','timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
189 WriteData(fid,prefix,'object',self,'fieldname','spcthickness','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
190 WriteData(fid,prefix,'object',self,'fieldname','reltol','format','Double')
191 WriteData(fid,prefix,'object',self,'fieldname','abstol','format','Double')
192 WriteData(fid,prefix,'object',self,'fieldname','maxiter','format','Integer')
193 WriteData(fid,prefix,'object',self,'fieldname','love_h','format','DoubleMat','mattype',1)
194 WriteData(fid,prefix,'object',self,'fieldname','love_k','format','DoubleMat','mattype',1)
195 WriteData(fid,prefix,'object',self,'fieldname','love_l','format','DoubleMat','mattype',1)
196 WriteData(fid,prefix,'object',self,'fieldname','tide_love_h','format','Double')
197 WriteData(fid,prefix,'object',self,'fieldname','tide_love_k','format','Double')
198 WriteData(fid,prefix,'object',self,'fieldname','fluid_love','format','Double')
199 WriteData(fid,prefix,'object',self,'fieldname','equatorial_moi','format','Double')
200 WriteData(fid,prefix,'object',self,'fieldname','polar_moi','format','Double')
201 WriteData(fid,prefix,'object',self,'fieldname','angular_velocity','format','Double')
202 WriteData(fid,prefix,'object',self,'fieldname','rigid','format','Boolean')
203 WriteData(fid,prefix,'object',self,'fieldname','elastic','format','Boolean')
204 WriteData(fid,prefix,'object',self,'fieldname','rotation','format','Boolean')
205 WriteData(fid,prefix,'object',self,'fieldname','ocean_area_scaling','format','Boolean')
206 WriteData(fid,prefix,'object',self,'fieldname','geodetic_run_frequency','format','Integer')
207 WriteData(fid,prefix,'object',self,'fieldname','steric_rate','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts)
208 WriteData(fid,prefix,'object',self,'fieldname','Ngia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts)
209 WriteData(fid,prefix,'object',self,'fieldname','Ugia','format','DoubleMat','mattype',1,'scale',1e-3/md.constants.yts)
210 WriteData(fid,prefix,'object',self,'fieldname','degacc','format','Double')
211 WriteData(fid,prefix,'object',self,'fieldname','transitions','format','MatArray')
212 WriteData(fid,prefix,'object',self,'fieldname','loop_increment','format','Integer')
213 WriteData(fid,prefix,'object',self,'fieldname','horiz','format','Integer')
214 WriteData(fid,prefix,'object',self,'fieldname','geodetic','format','Integer')
215
216 #process requested outputs
217 outputs = self.requested_outputs
218 indices = [i for i, x in enumerate(outputs) if x == 'default']
219 if len(indices) > 0:
220 outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
221 outputs =outputscopy
222 WriteData(fid,prefix,'data',outputs,'name','md.slr.requested_outputs','format','StringArray')
223 # }}}
Note: See TracBrowser for help on using the repository browser.