| [12038] | 1 | #module imports
|
|---|
| [13101] | 2 | import numpy
|
|---|
| 3 | import sys
|
|---|
| [12038] | 4 | from fielddisplay import fielddisplay
|
|---|
| [13101] | 5 | from EnumDefinitions import *
|
|---|
| 6 | from checkfield import *
|
|---|
| 7 | from WriteData import *
|
|---|
| [12038] | 8 |
|
|---|
| [13101] | 9 | class diagnostic(object):
|
|---|
| 10 | """
|
|---|
| 11 | DIAGNOSTIC class definition
|
|---|
| 12 |
|
|---|
| 13 | Usage:
|
|---|
| 14 | diagnostic=diagnostic();
|
|---|
| 15 | """
|
|---|
| 16 |
|
|---|
| [12038] | 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')
|
|---|
| [12123] | 38 |
|
|---|
| 39 | #set defaults
|
|---|
| 40 | self.setdefaultparameters()
|
|---|
| 41 |
|
|---|
| [12038] | 42 | #}}}
|
|---|
| [13101] | 43 | def __repr__(self):
|
|---|
| [12038] | 44 | # {{{ Display
|
|---|
| 45 |
|
|---|
| 46 | string='\n Diagnostic solution parameters:'
|
|---|
| 47 | string="%s\n\n%s"%(string,' Convergence criteria:')
|
|---|
| 48 |
|
|---|
| [13101] | 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)'))
|
|---|
| [12038] | 55 |
|
|---|
| 56 | string="%s\n%s"%(string,' boundary conditions:')
|
|---|
| 57 |
|
|---|
| [13101] | 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'))
|
|---|
| [12038] | 62 |
|
|---|
| 63 | string="%s\n%s"%(string,' Rift options:')
|
|---|
| [13101] | 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'))
|
|---|
| [12038] | 66 |
|
|---|
| 67 | string="%s\n%s"%(string,' Penalty options:')
|
|---|
| [13101] | 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'))
|
|---|
| [12038] | 70 |
|
|---|
| 71 | string="%s\n%s"%(string,' Other:')
|
|---|
| [13101] | 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'))
|
|---|
| [12038] | 76 |
|
|---|
| 77 | return string
|
|---|
| 78 | #}}}
|
|---|
| [12123] | 79 |
|
|---|
| [13101] | 80 | def setdefaultparameters(self):
|
|---|
| [12123] | 81 | # {{{setdefaultparameters
|
|---|
| 82 | #maximum of non-linear iterations.
|
|---|
| [13101] | 83 | self.maxiter=100
|
|---|
| [12123] | 84 |
|
|---|
| 85 | #Convergence criterion: absolute, relative and residual
|
|---|
| [13101] | 86 | self.restol=10**-4
|
|---|
| 87 | self.reltol=0.01
|
|---|
| 88 | self.abstol=10
|
|---|
| [12123] | 89 |
|
|---|
| [13101] | 90 | self.stokesreconditioning=10**13
|
|---|
| 91 | self.shelf_dampening=0
|
|---|
| [12123] | 92 |
|
|---|
| 93 | #Penalty factor applied kappa=max(stiffness matrix)*10^penalty_factor
|
|---|
| [13101] | 94 | self.penalty_factor=3
|
|---|
| [12123] | 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))
|
|---|
| [13101] | 99 | self.viscosity_overshoot=0
|
|---|
| [12123] | 100 |
|
|---|
| 101 | #Stop the iterations of rift if below a threshold
|
|---|
| [13101] | 102 | self.rift_penalty_threshold=0
|
|---|
| [12123] | 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.
|
|---|
| [13101] | 107 | self.rift_penalty_lock=10
|
|---|
| [12123] | 108 |
|
|---|
| [13101] | 109 | return self
|
|---|
| [12123] | 110 | #}}}
|
|---|
| [13101] | 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 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
