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

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 9 years ago

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