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

Last change on this file since 16560 was 16560, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 16554
File size: 10.5 KB
Line 
1import numpy
2import sys
3import copy
4from fielddisplay import fielddisplay
5from EnumDefinitions import *
6from checkfield import *
7from WriteData import *
8from MatlabFuncs import *
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 setdefaultparameters(self): # {{{
78 #maximum of non-linear iterations.
79 self.maxiter=100
80
81 #Convergence criterion: absolute, relative and residual
82 self.restol=10**-4
83 self.reltol=0.01
84 self.abstol=10
85
86 self.FSreconditioning=10**13
87 self.shelf_dampening=0
88
89 #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
90 self.penalty_factor=3
91
92 #coefficient to update the viscosity between each iteration of
93 #a stressbalance according to the following formula
94 #viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
95 self.viscosity_overshoot=0
96
97 #Stop the iterations of rift if below a threshold
98 self.rift_penalty_threshold=0
99
100 #in some solutions, it might be needed to stop a run when only
101 #a few constraints remain unstable. For thermal computation, this
102 #parameter is often used.
103 self.rift_penalty_lock=10
104
105 #output default:
106 self.requested_outputs=['default']
107
108 return self
109 #}}}
110 def defaultoutputs(self,md): # {{{
111
112 if strcmp(md.mesh.meshtype(),'3D'):
113 list = ['Vx','Vy','Vz','Vel','Pressure']
114 elif strcmp(md.mesh.meshtype(),'2Dhorizontal'):
115 list = ['Vx','Vy','Vel','Pressure']
116 elif strcmp(md.mesh.meshtype(),'2Dvertical'):
117 list = ['Vx','Vy','Vel','Pressure']
118 else:
119 raise TypeError('mesh type not supported yet');
120 return list
121
122 #}}}
123 def checkconsistency(self,md,solution,analyses): # {{{
124
125 #Early return
126 if StressbalanceAnalysisEnum() not in analyses:
127 return md
128
129 md = checkfield(md,'stressbalance.spcvx','forcing',1)
130 md = checkfield(md,'stressbalance.spcvy','forcing',1)
131 if strcmp(md.mesh.meshtype(),'3D'):
132 md = checkfield(md,'stressbalance.spcvz','forcing',1)
133 md = checkfield(md,'stressbalance.restol','size',[1],'>',0)
134 md = checkfield(md,'stressbalance.reltol','size',[1])
135 md = checkfield(md,'stressbalance.abstol','size',[1])
136 md = checkfield(md,'stressbalance.isnewton','numel',[1],'values',[0,1,2])
137 md = checkfield(md,'stressbalance.FSreconditioning','size',[1],'NaN',1)
138 md = checkfield(md,'stressbalance.viscosity_overshoot','size',[1],'NaN',1)
139 md = checkfield(md,'stressbalance.maxiter','size',[1],'>=',1)
140 md = checkfield(md,'stressbalance.referential','size',[md.mesh.numberofvertices,6])
141 md = checkfield(md,'stressbalance.loadingforce','size',[md.mesh.numberofvertices,3])
142 md = checkfield(md,'stressbalance.requested_outputs','stringrow',1);
143
144 #singular solution
145# if ~any((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy))==2),
146 if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.stressbalance.spcvx)),numpy.logical_not(numpy.isnan(md.stressbalance.spcvy)))):
147 md.checkmessage("model is not well posed (singular). You need at least one node with fixed velocity!")
148 #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
149# if any(sum(isnan(md.stressbalance.referential),2)~=0 & sum(isnan(md.stressbalance.referential),2)~=6),
150 if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)!=6)):
151 md.checkmessage("Each line of stressbalance.referential should contain either only NaN values or no NaN values")
152 #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
153# if any(sum(isnan(md.stressbalance.referential),2)==0),
154 if numpy.any(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)==0):
155 pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.stressbalance.referential),axis=1)) if item==0]
156# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply
157# if numpy.any(numpy.abs(numpy.inner(md.stressbalance.referential[pos,0:2],md.stressbalance.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):
158# md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
159 for item in md.stressbalance.referential[pos,:]:
160 if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:
161 md.checkmessage("Vectors in stressbalance.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
162 #CHECK THAT NO rotation specified for FS Grounded ice at base
163 if strcmp(md.mesh.meshtype(),'3D') and md.flowequation.isFS:
164 pos=numpy.nonzero(numpy.logical_and(md.mask.groundedice_levelset,md.mesh.vertexonbed))
165 if numpy.any(numpy.logical_not(numpy.isnan(md.stressbalance.referential[pos,:]))):
166 md.checkmessage("no referential should be specified for basal vertices of grounded ice")
167
168 return md
169 # }}}
170 def marshall(self,md,fid): # {{{
171
172 yts=365.0*24.0*3600.0
173
174 WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvx','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
175 WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvy','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
176 WriteData(fid,'object',self,'class','stressbalance','fieldname','spcvz','format','DoubleMat','mattype',1,'scale',1./yts,'forcinglength',md.mesh.numberofvertices+1)
177 WriteData(fid,'object',self,'class','stressbalance','fieldname','restol','format','Double')
178 WriteData(fid,'object',self,'class','stressbalance','fieldname','reltol','format','Double')
179 WriteData(fid,'object',self,'class','stressbalance','fieldname','abstol','format','Double')
180 WriteData(fid,'object',self,'class','stressbalance','fieldname','isnewton','format','Integer')
181 WriteData(fid,'object',self,'class','stressbalance','fieldname','FSreconditioning','format','Double')
182 WriteData(fid,'object',self,'class','stressbalance','fieldname','viscosity_overshoot','format','Double')
183 WriteData(fid,'object',self,'class','stressbalance','fieldname','maxiter','format','Integer')
184 WriteData(fid,'object',self,'class','stressbalance','fieldname','shelf_dampening','format','Integer')
185 WriteData(fid,'object',self,'class','stressbalance','fieldname','vertex_pairing','format','DoubleMat','mattype',3)
186 WriteData(fid,'object',self,'class','stressbalance','fieldname','penalty_factor','format','Double')
187 WriteData(fid,'object',self,'class','stressbalance','fieldname','rift_penalty_lock','format','Integer')
188 WriteData(fid,'object',self,'class','stressbalance','fieldname','rift_penalty_threshold','format','Integer')
189 WriteData(fid,'object',self,'class','stressbalance','fieldname','referential','format','DoubleMat','mattype',1)
190
191 WriteData(fid,'data',self.loadingforce[:,0],'format','DoubleMat','mattype',1,'enum',LoadingforceXEnum())
192 WriteData(fid,'data',self.loadingforce[:,1],'format','DoubleMat','mattype',1,'enum',LoadingforceYEnum())
193 WriteData(fid,'data',self.loadingforce[:,2],'format','DoubleMat','mattype',1,'enum',LoadingforceZEnum())
194
195 #process requested outputs
196 outputs = self.requested_outputs
197 indices = [i for i, x in enumerate(outputs) if x == 'default']
198 if len(indices) > 0:
199 outputscopy=outputs[0:max(0,indices[0]-1)]+self.defaultoutputs(md)+outputs[indices[0]+1:]
200 outputs =outputscopy
201 WriteData(fid,'data',outputs,'enum',StressbalanceRequestedOutputsEnum(),'format','StringArray')
202 # }}}
Note: See TracBrowser for help on using the repository browser.