source: issm/trunk/src/m/classes/stressbalance.py@ 21341

Last change on this file since 21341 was 21341, checked in by Mathieu Morlighem, 8 years ago

merged trunk-jpl and trunk for revision 21337
File size: 11.3 KB
Line 
1import numpy as np
2import sys
3import copy
4from project3d import project3d
5from fielddisplay import fielddisplay
6from checkfield import checkfield
7from WriteData import WriteData
8import MatlabFuncs as m
9
10class stressbalance(object):
11 """
12 STRESSBALANCE class definition
13
14 Usage:
15 stressbalance=stressbalance();
16 """
17
18 def __init__(self): # {{{
19 self.spcvx = float('NaN')
20 self.spcvy = float('NaN')
21 self.spcvz = float('NaN')
22 self.restol = 0
23 self.reltol = 0
24 self.abstol = 0
25 self.isnewton = 0
26 self.FSreconditioning = 0
27 self.viscosity_overshoot = 0
28 self.icefront = float('NaN')
29 self.maxiter = 0
30 self.shelf_dampening = 0
31 self.vertex_pairing = float('NaN')
32 self.penalty_factor = float('NaN')
33 self.rift_penalty_lock = float('NaN')
34 self.rift_penalty_threshold = 0
35 self.referential = float('NaN')
36 self.loadingforce = float('NaN')
37 self.requested_outputs = []
38
39 #set defaults
40 self.setdefaultparameters()
41
42 #}}}
43 def __repr__(self): # {{{
44
45 string=' StressBalance solution parameters:'
46 string="%s\n%s"%(string,' Convergence criteria:')
47 string="%s\n%s"%(string,fielddisplay(self,'restol','mechanical equilibrium residual convergence criterion'))
48 string="%s\n%s"%(string,fielddisplay(self,'reltol','velocity relative convergence criterion, NaN: not applied'))
49 string="%s\n%s"%(string,fielddisplay(self,'abstol','velocity absolute convergence criterion, NaN: not applied'))
50 string="%s\n%s"%(string,fielddisplay(self,'isnewton',"0: Picard's fixed point, 1: Newton's method, 2: hybrid"))
51 string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations'))
52 string="%s\n%s"%(string,fielddisplay(self,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'))
53
54 string="%s\n%s"%(string,'\n boundary conditions:')
55 string="%s\n%s"%(string,fielddisplay(self,'spcvx','x-axis velocity constraint (NaN means no constraint) [m/yr]'))
56 string="%s\n%s"%(string,fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint) [m/yr]'))
57 string="%s\n%s"%(string,fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint) [m/yr]'))
58 string="%s\n%s"%(string,fielddisplay(self,'icefront','segments on ice front list (last column 0: Air, 1: Water, 2: Ice'))
59
60 string="%s\n%s"%(string,'\n Rift options:')
61 string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints'))
62 string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked'))
63
64 string="%s\n%s"%(string,'\n Penalty options:')
65 string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'))
66 string="%s\n%s"%(string,fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized'))
67
68 string="%s\n%s"%(string,'\n Other:')
69 string="%s\n%s"%(string,fielddisplay(self,'shelf_dampening','use dampening for floating ice ? Only for FS model'))
70 string="%s\n%s"%(string,fielddisplay(self,'FSreconditioning','multiplier for incompressibility equation. Only for FS model'))
71 string="%s\n%s"%(string,fielddisplay(self,'referential','local referential'))
72 string="%s\n%s"%(string,fielddisplay(self,'loadingforce','loading force applied on each point [N/m^3]'))
73 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
74
75 return string
76 #}}}
77 def extrude(self,md): # {{{
78 self.spcvx=project3d(md,'vector',self.spcvx,'type','node')
79 self.spcvy=project3d(md,'vector',self.spcvy,'type','node')
80 self.spcvz=project3d(md,'vector',self.spcvz,'type','node')
81 self.referential=project3d(md,'vector',self.referential,'type','node')
82 self.loadingforce=project3d(md,'vector',self.loadingforce,'type','node')
83
84 return self
85 #}}}
86 def setdefaultparameters(self): # {{{
87 #maximum of non-linear iterations.
88 self.maxiter=100
89
90 #Convergence criterion: absolute, relative and residual
91 self.restol=10**-4
92 self.reltol=0.01
93 self.abstol=10
94
95 self.FSreconditioning=10**13
96 self.shelf_dampening=0
97
98 #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
99 self.penalty_factor=3
100
101 #coefficient to update the viscosity between each iteration of
102 #a stressbalance according to the following formula
103 #viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
104 self.viscosity_overshoot=0
105
106 #Stop the iterations of rift if below a threshold
107 self.rift_penalty_threshold=0
108
109 #in some solutions, it might be needed to stop a run when only
110 #a few constraints remain unstable. For thermal computation, this
111 #parameter is often used.
112 self.rift_penalty_lock=10
113
114 #output default:
115 self.requested_outputs=['default']
116
117 return self
118 #}}}
119 def defaultoutputs(self,md): # {{{
120
121 if md.mesh.dimension()==3:
122 list = ['Vx','Vy','Vz','Vel','Pressure']
123 else:
124 list = ['Vx','Vy','Vel','Pressure']
125 return list
126
127 #}}}
128 def checkconsistency(self,md,solution,analyses): # {{{
129
130 #Early return
131 if 'StressbalanceAnalysis' not in analyses:
132 return md
133
134 md = checkfield(md,'fieldname','stressbalance.spcvx','Inf',1,'timeseries',1)
135 md = checkfield(md,'fieldname','stressbalance.spcvy','Inf',1,'timeseries',1)
136 if m.strcmp(md.mesh.domaintype(),'3D'):
137 md = checkfield(md,'fieldname','stressbalance.spcvz','Inf',1,'timeseries',1)
138 md = checkfield(md,'fieldname','stressbalance.restol','size',[1],'>',0)
139 md = checkfield(md,'fieldname','stressbalance.reltol','size',[1])
140 md = checkfield(md,'fieldname','stressbalance.abstol','size',[1])
141 md = checkfield(md,'fieldname','stressbalance.isnewton','numel',[1],'values',[0,1,2])
142 md = checkfield(md,'fieldname','stressbalance.FSreconditioning','size',[1],'NaN',1,'Inf',1)
143 md = checkfield(md,'fieldname','stressbalance.viscosity_overshoot','size',[1],'NaN',1,'Inf',1)
144 md = checkfield(md,'fieldname','stressbalance.maxiter','size',[1],'>=',1)
145 md = checkfield(md,'fieldname','stressbalance.referential','size',[md.mesh.numberofvertices,6])
146 md = checkfield(md,'fieldname','stressbalance.loadingforce','size',[md.mesh.numberofvertices,3])
147 md = checkfield(md,'fieldname','stressbalance.requested_outputs','stringrow',1);
148
149 #singular solution
150# if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2),
151 if not np.any(np.logical_and(np.logical_not(np.isnan(md.stressbalance.spcvx)),np.logical_not(np.isnan(md.stressbalance.spcvy)))):
152 print "\n !!! Warning: no spc applied, model might not be well posed if no basal friction is applied, check for solution crash\n"
153 #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
154# if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6),
155 if np.any(np.logical_and(np.sum(np.isnan(md.stressbalance.referential),axis=1)!=0,np.sum(np.isnan(md.stressbalance.referential),axis=1)!=6)):
156 md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values")
157 #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
158# if any(sum(isnan(md.stressbalance.referential),2)==0),
159 if np.any(np.sum(np.isnan(md.stressbalance.referential),axis=1)==0):
160 pos=[i for i,item in enumerate(np.sum(np.isnan(md.stressbalance.referential),axis=1)) if item==0]
161# np.inner (and np.dot) calculate all the dot product permutations, resulting in a full matrix multiply
162# if np.any(np.abs(np.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):
163# md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
164 for item in md.stressbalance.referential[pos,:]:
165 if np.abs(np.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:
166 md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
167 #CHECK THAT NO rotation specified for FS Grounded ice at base
168 if m.strcmp(md.mesh.domaintype(),'3D') and md.flowequation.isFS:
169 pos=np.nonzero(np.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbase))
170 if np.any(np.logical_not(np.isnan(md.stressbalance.referential[pos,:]))):
171 md.checkmessage("no referential should be specified for basal vertices of grounded ice")
172
173 return md
174 # }}}
175 def marshall(self,prefix,md,fid): # {{{
176
177 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3)
178
179 yts=md.constants.yts
180
181 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
182 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
183 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts)
184 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','restol','format','Double')
185 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','reltol','format','Double')
186 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','abstol','format','Double','scale',1./yts)
187 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','isnewton','format','Integer')
188 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','FSreconditioning','format','Double')
189 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','viscosity_overshoot','format','Double')
190 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','maxiter','format','Integer')
191 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','shelf_dampening','format','Integer')
192 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','penalty_factor','format','Double')
193 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer')
194 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer')
195 WriteData(fid,prefix,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1)
196
197 if isinstance(self.loadingforce, (list, tuple, np.ndarray)):
198 lx=self.loadingforce[:,0];
199 ly=self.loadingforce[:,1];
200 lz=self.loadingforce[:,2];
201 else:
202 lx=float('NaN'); ly=float('NaN'); lz=float('NaN');
203
204 WriteData(fid,prefix,'data',lx,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcex')
205 WriteData(fid,prefix,'data',ly,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcey')
206 WriteData(fid,prefix,'data',lz,'format','DoubleMat','mattype',1,'name','md.stressbalance.loadingforcez')
207
208 #process requested outputs
209 outputs = self.requested_outputs
210 indices = [i for i, x in enumerate(outputs) if x == 'default']
211 if len(indices) > 0:
212 outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
213 outputs =outputscopy
214 WriteData(fid,prefix,'data',outputs,'name','md.stressbalance.requested_outputs','format','StringArray')
215 # }}}
Note: See TracBrowser for help on using the repository browser.