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

Last change on this file since 23189 was 23189, checked in by Mathieu Morlighem, 7 years ago

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