source: issm/branches/trunk-jpl-damage/src/m/classes/diagnostic.py@ 13101

Last change on this file since 13101 was 13101, checked in by cborstad, 13 years ago

merged trunk-jpl through revision 13099 into branch
File size: 9.5 KB
Line 
1#module imports
2import numpy
3import sys
4from fielddisplay import fielddisplay
5from EnumDefinitions import *
6from checkfield import *
7from WriteData import *
8
9class diagnostic(object):
10 """
11 DIAGNOSTIC class definition
12
13 Usage:
14 diagnostic=diagnostic();
15 """
16
17 #properties
18 def __init__(self):
19 # {{{ Properties
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.stokesreconditioning = 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.requested_outputs = float('NaN')
38
39 #set defaults
40 self.setdefaultparameters()
41
42 #}}}
43 def __repr__(self):
44 # {{{ Display
45
46 string='\n Diagnostic solution parameters:'
47 string="%s\n\n%s"%(string,' Convergence criteria:')
48
49 string="%s\n%s"%(string,fielddisplay(self,'restol','mechanical equilibrium residual convergence criterion'))
50 string="%s\n%s"%(string,fielddisplay(self,'reltol','velocity relative convergence criterion, NaN -> not applied'))
51 string="%s\n%s"%(string,fielddisplay(self,'abstol','velocity absolute convergence criterion, NaN -> not applied'))
52 string="%s\n%s"%(string,fielddisplay(self,'isnewton',"Apply Newton's method instead of a Picard fixed point method"))
53 string="%s\n%s"%(string,fielddisplay(self,'maxiter','maximum number of nonlinear iterations'))
54 string="%s\n%s"%(string,fielddisplay(self,'viscosity_overshoot','over-shooting constant new=new+C*(new-old)'))
55
56 string="%s\n%s"%(string,' boundary conditions:')
57
58 string="%s\n%s"%(string,fielddisplay(self,'spcvx','x-axis velocity constraint (NaN means no constraint)'))
59 string="%s\n%s"%(string,fielddisplay(self,'spcvy','y-axis velocity constraint (NaN means no constraint)'))
60 string="%s\n%s"%(string,fielddisplay(self,'spcvz','z-axis velocity constraint (NaN means no constraint)'))
61 string="%s\n%s"%(string,fielddisplay(self,'icefront','segments on ice front list (last column 0-> Air, 1-> Water, 2->Ice'))
62
63 string="%s\n%s"%(string,' Rift options:')
64 string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_threshold','threshold for instability of mechanical constraints'))
65 string="%s\n%s"%(string,fielddisplay(self,'rift_penalty_lock','number of iterations before rift penalties are locked'))
66
67 string="%s\n%s"%(string,' Penalty options:')
68 string="%s\n%s"%(string,fielddisplay(self,'penalty_factor','offset used by penalties: penalty = Kmax*10^offset'))
69 string="%s\n%s"%(string,fielddisplay(self,'vertex_pairing','pairs of vertices that are penalized'))
70
71 string="%s\n%s"%(string,' Other:')
72 string="%s\n%s"%(string,fielddisplay(self,'shelf_dampening','use dampening for floating ice ? Only for Stokes model'))
73 string="%s\n%s"%(string,fielddisplay(self,'stokesreconditioning','multiplier for incompressibility equation. Only for Stokes model'))
74 string="%s\n%s"%(string,fielddisplay(self,'referential','local referential'))
75 string="%s\n%s"%(string,fielddisplay(self,'requested_outputs','additional outputs requested'))
76
77 return string
78 #}}}
79
80 def setdefaultparameters(self):
81 # {{{setdefaultparameters
82 #maximum of non-linear iterations.
83 self.maxiter=100
84
85 #Convergence criterion: absolute, relative and residual
86 self.restol=10**-4
87 self.reltol=0.01
88 self.abstol=10
89
90 self.stokesreconditioning=10**13
91 self.shelf_dampening=0
92
93 #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
94 self.penalty_factor=3
95
96 #coefficient to update the viscosity between each iteration of
97 #a diagnostic according to the following formula
98 #viscosity(n)=viscosity(n)+viscosity_overshoot(viscosity(n)-viscosity(n-1))
99 self.viscosity_overshoot=0
100
101 #Stop the iterations of rift if below a threshold
102 self.rift_penalty_threshold=0
103
104 #in some solutions, it might be needed to stop a run when only
105 #a few constraints remain unstable. For thermal computation, this
106 #parameter is often used.
107 self.rift_penalty_lock=10
108
109 return self
110 #}}}
111
112 def checkconsistency(self,md,solution,analyses): # {{{
113
114 #Early return
115 if DiagnosticHorizAnalysisEnum() not in analyses:
116 return md
117 #if (DiagnosticHorizAnalysisEnum() not in analyses) | (solution==TransientSolutionEnum() and not md.transient.isdiagnostic):
118 # return md
119
120 md = checkfield(md,'diagnostic.spcvx','forcing',1)
121 md = checkfield(md,'diagnostic.spcvy','forcing',1)
122 if md.mesh.dimension==3:
123 md = checkfield(md,'diagnostic.spcvz','forcing',1)
124 md = checkfield(md,'diagnostic.restol','size',[1],'>',0)
125 md = checkfield(md,'diagnostic.reltol','size',[1])
126 md = checkfield(md,'diagnostic.abstol','size',[1])
127 md = checkfield(md,'diagnostic.isnewton','numel',[1],'values',[0,1])
128 md = checkfield(md,'diagnostic.stokesreconditioning','size',[1],'NaN',1)
129 md = checkfield(md,'diagnostic.viscosity_overshoot','size',[1],'NaN',1)
130 if md.mesh.dimension==2:
131 md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),4],'NaN',1)
132 else:
133 md = checkfield(md,'diagnostic.icefront','size',[float('NaN'),6],'NaN',1)
134 md = checkfield(md,'diagnostic.icefront[:,-1]','values',[0,1,2])
135 md = checkfield(md,'diagnostic.maxiter','size',[1],'>=',1)
136 md = checkfield(md,'diagnostic.referential','size',[md.mesh.numberofvertices,6])
137 if not md.diagnostic.requested_outputs:
138 md = checkfield(md,'diagnostic.requested_outputs','size',[float('NaN'),1])
139
140 #singular solution
141# if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2),
142 if not numpy.any(numpy.logical_and(numpy.logical_not(numpy.isnan(md.diagnostic.spcvx)),numpy.logical_not(numpy.isnan(md.diagnostic.spcvy)))):
143 md.checkmessage("model is not well posed (singular). You need at least one node with fixed velocity!")
144 #CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
145# if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6),
146 if numpy.any(numpy.logical_and(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=0,numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)!=6)):
147 md.checkmessage("Each line of diagnostic.referential should contain either only NaN values or no NaN values")
148 #CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
149# if any(sum(isnan(md.diagnostic.referential),2)==0),
150 if numpy.any(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)==0):
151 pos=[i for i,item in enumerate(numpy.sum(numpy.isnan(md.diagnostic.referential),axis=1)) if item==0]
152# numpy.inner (and numpy.dot) calculate all the dot product permutations, resulting in a full matrix multiply
153# if numpy.any(numpy.abs(numpy.inner(md.diagnostic.referential[pos,0:2],md.diagnostic.referential[pos,3:5]).diagonal())>sys.float_info.epsilon):
154# md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
155 for item in md.diagnostic.referential[pos,:]:
156 if numpy.abs(numpy.inner(item[0:2],item[3:5]))>sys.float_info.epsilon:
157 md.checkmessage("Vectors in diagnostic.referential (columns 1 to 3 and 4 to 6) must be orthogonal")
158 #CHECK THAT NO rotation specified for FS Grounded ice at base
159# if md.mesh.dimension==3 & md.flowequation.isstokes,
160 if md.mesh.dimension==3 and md.flowequation.isstokes:
161 pos=numpy.nonzero(numpy.logical_and(md.mask.vertexongroundedice,md.mesh.vertexonbed))
162 if numpy.any(numpy.logical_not(numpy.isnan(md.diagnostic.referential[pos,:]))):
163 md.checkmessage("no referential should be specified for basal vertices of grounded ice")
164
165 return md
166 # }}}
167
168 def marshall(self,fid): # {{{
169 WriteData(fid,'object',self,'fieldname','spcvx','format','DoubleMat','mattype',1)
170 WriteData(fid,'object',self,'fieldname','spcvy','format','DoubleMat','mattype',1)
171 WriteData(fid,'object',self,'fieldname','spcvz','format','DoubleMat','mattype',1)
172 WriteData(fid,'object',self,'fieldname','restol','format','Double')
173 WriteData(fid,'object',self,'fieldname','reltol','format','Double')
174 WriteData(fid,'object',self,'fieldname','abstol','format','Double')
175 WriteData(fid,'object',self,'fieldname','isnewton','format','Boolean')
176 WriteData(fid,'object',self,'fieldname','stokesreconditioning','format','Double')
177 WriteData(fid,'object',self,'fieldname','viscosity_overshoot','format','Double')
178 WriteData(fid,'object',self,'fieldname','maxiter','format','Integer')
179 WriteData(fid,'object',self,'fieldname','shelf_dampening','format','Integer')
180 WriteData(fid,'object',self,'fieldname','vertex_pairing','format','DoubleMat','mattype',3)
181 WriteData(fid,'object',self,'fieldname','penalty_factor','format','Double')
182 WriteData(fid,'object',self,'fieldname','rift_penalty_lock','format','Integer')
183 WriteData(fid,'object',self,'fieldname','rift_penalty_threshold','format','Integer')
184 WriteData(fid,'object',self,'fieldname','referential','format','DoubleMat','mattype',1)
185 WriteData(fid,'object',self,'fieldname','requested_outputs','format','DoubleMat','mattype',3)
186
187 #marshall ice front
188 data=self.icefront
189 data[[i for i,item in enumerate(data[:,-1]) if item==0],-1]=AirEnum()
190 data[[i for i,item in enumerate(data[:,-1]) if item==1],-1]=WaterEnum()
191 data[[i for i,item in enumerate(data[:,-1]) if item==2],-1]=IceEnum()
192 WriteData(fid,'data',data,'enum',DiagnosticIcefrontEnum(),'format','DoubleMat','mattype',3)
193 # }}}
194
Note: See TracBrowser for help on using the repository browser.