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

Last change on this file since 16137 was 16137, checked in by Mathieu Morlighem, 12 years ago

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