source: issm/trunk/src/py3/classes/taoinversion.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
  • Property svn:executable set to *
File size: 11.0 KB
Line 
1import numpy
2from project3d import project3d
3from WriteData import WriteData
4from checkfield import checkfield
5from fielddisplay import fielddisplay
6from IssmConfig import IssmConfig
7from EnumDefinitions import *
8from marshallcostfunctions import marshallcostfunctions
9
10
11class taoinversion:
12 def __init__(self):
13 iscontrol = 0
14 incomplete_adjoint = 0
15 control_parameters = float('NaN')
16 maxsteps = 0
17 maxiter = 0
18 fatol = 0
19 frtol = 0
20 gatol = 0
21 grtol = 0
22 gttol = 0
23 algorithm = ''
24 cost_functions = float('NaN')
25 cost_functions_coefficients = float('NaN')
26 min_parameters = float('NaN')
27 max_parameters = float('NaN')
28 vx_obs = float('NaN')
29 vy_obs = float('NaN')
30 vz_obs = float('NaN')
31 vel_obs = float('NaN')
32 thickness_obs = float('NaN')
33 surface_obs = float('NaN')
34
35 def __repr__(self):
36 string = ' taoinversion parameters:'
37 string = "%s\n\%s"%(string, fieldstring(self,'iscontrol','is inversion activated?'))
38 string="%s\n%s"%(string,fieldstring(self,'mantle_viscosity','mantle viscosity constraints (NaN means no constraint) (Pa s)'))
39 string="%s\n%s"%(string,fieldstring(self,'lithosphere_thickness','lithosphere thickness constraints (NaN means no constraint) (m)'))
40 string="%s\n%s"%(string,fieldstring(self,'cross_section_shape',"1: square-edged, 2: elliptical-edged surface"))
41 string="%s\n%s"%(string,fieldstring(self,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity'))
42 string="%s\n%s"%(string,fieldstring(self,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}'))
43 string="%s\n%s"%(string,fieldstring(self,'maxsteps','maximum number of iterations (gradient computation)'))
44 string="%s\n%s"%(string,fieldstring(self,'maxiter','maximum number of Function evaluation (forward run)'))
45 string="%s\n%s"%(string,fieldstring(self,'fatol','convergence criterion: f(X)-f(X*) (X: current iteration, X*: "true" solution, f: cost function)'))
46 string="%s\n%s"%(string,fieldstring(self,'frtol','convergence criterion: |f(X)-f(X*)|/|f(X*)|'))
47 string="%s\n%s"%(string,fieldstring(self,'gatol','convergence criterion: ||g(X)|| (g: gradient of the cost function)'))
48 string="%s\n%s"%(string,fieldstring(self,'grtol','convergence criterion: ||g(X)||/|f(X)|'))
49 string="%s\n%s"%(string,fieldstring(self,'gttol','convergence criterion: ||g(X)||/||g(X0)|| (g(X0): gradient at initial guess X0)'))
50 string="%s\n%s"%(string,fieldstring(self,'algorithm','minimization algorithm: ''tao_blmvm'', ''tao_cg'', ''tao_lmvm'''))
51 string="%s\n%s"%(string,fieldstring(self,'cost_functions','indicate the type of response for each optimization step'))
52 string="%s\n%s"%(string,fieldstring(self,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter'))
53 string="%s\n%s"%(string,fieldstring(self,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex'))
54 string="%s\n%s"%(string,fieldstring(self,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex'))
55 string="%s\n%s"%(string,fieldstring(self,'vx_obs','observed velocity x component [m/yr]'))
56 string="%s\n%s"%(string,fieldstring(self,'vy_obs','observed velocity y component [m/yr]'))
57 string="%s\n%s"%(string,fieldstring(self,'vel_obs','observed velocity magnitude [m/yr]'))
58 string="%s\n%s"%(string,fieldstring(self,'thickness_obs','observed thickness [m]'))
59 string="%s\n%s"%(string,fieldstring(self,'surface_obs','observed surface elevation [m]'))
60 string="%s\n%s"%(string,'Available cost functions:')
61 string="%s\n%s"%(string, ' 101: SurfaceAbsVelMisfit')
62 string="%s\n%s"%(string, ' 102: SurfaceRelVelMisfit')
63 string="%s\n%s"%(string, ' 103: SurfaceLogVelMisfit')
64 string="%s\n%s"%(string, ' 104: SurfaceLogVxVyMisfit')
65 string="%s\n%s"%(string, ' 105: SurfaceAverageVelMisfit')
66 string="%s\n%s"%(string, ' 201: ThicknessAbsMisfit')
67 string="%s\n%s"%(string, ' 501: DragCoefficientAbsGradient')
68 string="%s\n%s"%(string, ' 502: RheologyBbarAbsGradient')
69 string="%s\n%s"%(string, ' 503: ThicknessAbsGradient')
70 return string
71 def setdefaultparameters(self):
72
73 #default is incomplete adjoint for now
74 self.incomplete_adjoint=1
75
76 #parameter to be inferred by control methods (only
77 #drag and B are supported yet)
78 self.control_parameters=['FrictionCoefficient']
79
80 #number of iterations and steps
81 self.maxsteps=20;
82 self.maxiter =30;
83
84 #default tolerances
85 self.fatol = 0;
86 self.frtol = 0;
87 self.gatol = 0;
88 self.grtol = 0;
89 self.gttol = 1e-4;
90
91 #minimization algorithm
92 PETSCMAJOR = IssmConfig('_PETSC_MAJOR_')
93 PETSCMINOR = IssmConfig('_PETSC_MINOR_')
94 if(PETSCMAJOR>3 or (PETSCMAJOR==3 and PETSCMINOR>=5)):
95 self.algorithm = 'blmvm';
96 else:
97 self.algorithm = 'tao_blmvm';
98
99 #several responses can be used:
100 self.cost_functions=101;
101
102 return self
103
104 def extrude(self,md):
105 self.vx_obs=project3d(md,'vector',self.vx_obs,'type','node')
106 self.vy_obs=project3d(md,'vector',self.vy_obs,'type','node')
107 self.vel_obs=project3d(md,'vector',self.vel_obs,'type','node')
108 self.thickness_obs=project3d(md,'vector',self.thickness_obs,'type','node')
109
110 if numel(self.cost_functions_coefficients) > 1:
111 self.cost_functions_coefficients=project3d(md,'vector',self.cost_functions_coefficients,'type','node')
112
113 if numel(self.min_parameters) > 1:
114 self.min_parameters=project3d(md,'vector',self.min_parameters,'type','node')
115
116 if numel(self.max_parameters)>1:
117 self.max_parameters=project3d(md,'vector',self.max_parameters,'type','node')
118
119 return self
120
121 def checkconsistency(self,md,solution,analyses):
122 if not self.control:
123 return md
124 if not IssmConfig('_HAVE_TAO_'):
125 md = checkmessage(md,['TAO has not been installed, ISSM needs to be reconfigured and recompiled with TAO'])
126
127
128 num_controls= numpy.numel(md.inversion.control_parameters)
129 num_costfunc= numpy.size(md.inversion.cost_functions,2)
130
131 md = checkfield(md,'fieldname','inversion.iscontrol','values',[0, 1])
132 md = checkfield(md,'fieldname','inversion.incomplete_adjoint','values',[0, 1])
133 md = checkfield(md,'fieldname','inversion.control_parameters','cell',1,'values',supportedcontrols())
134 md = checkfield(md,'fieldname','inversion.maxsteps','numel',1,'>=',0)
135 md = checkfield(md,'fieldname','inversion.maxiter','numel',1,'>=',0)
136 md = checkfield(md,'fieldname','inversion.fatol','numel',1,'>=',0)
137 md = checkfield(md,'fieldname','inversion.frtol','numel',1,'>=',0)
138 md = checkfield(md,'fieldname','inversion.gatol','numel',1,'>=',0)
139 md = checkfield(md,'fieldname','inversion.grtol','numel',1,'>=',0)
140 md = checkfield(md,'fieldname','inversion.gttol','numel',1,'>=',0)
141
142
143 PETSCMAJOR = IssmConfig('_PETSC_MAJOR_')
144 PETSCMINOR = IssmConfig('_PETSC_MINOR_')
145 if(PETSCMAJOR>3 or (PETSCMAJOR==3 and PETSCMINOR>=5)):
146 md = checkfield(md,'fieldname','inversion.algorithm','values',{'blmvm','cg','lmvm'})
147 else:
148 md = checkfield(md,'fieldname','inversion.algorithm','values',{'tao_blmvm','tao_cg','tao_lmvm'})
149
150
151 md = checkfield(md,'fieldname','inversion.cost_functions','size',[1, num_costfunc],'values',supportedcostfunctions())
152 md = checkfield(md,'fieldname','inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices, num_costfunc],'>=',0)
153 md = checkfield(md,'fieldname','inversion.min_parameters','size',[md.mesh.numberofvertices, num_controls])
154 md = checkfield(md,'fieldname','inversion.max_parameters','size',[md.mesh.numberofvertices, num_controls])
155
156
157 if solution==BalancethicknessSolutionEnum():
158 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
159 elif solution==BalancethicknessSoftSolutionEnum():
160 md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
161 else:
162 md = checkfield(md,'fieldname','inversion.vx_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
163 md = checkfield(md,'fieldname','inversion.vy_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
164
165 def marshall(self, md, fid):
166
167 yts=365.0*24.0*3600.0;
168 WriteData(fid,'object',self,'class','inversion','fieldname','iscontrol','format','Boolean')
169 WriteData(fid,'enum',InversionTypeEnum(),'data',1,'format','Integer')
170 if not self.iscontrol:
171 return
172 WriteData(fid,'object',self,'class','inversion','fieldname','incomplete_adjoint','format','Boolean')
173 WriteData(fid,'object',self,'class','inversion','fieldname','maxsteps','format','Integer')
174 WriteData(fid,'object',self,'class','inversion','fieldname','maxiter','format','Integer')
175 WriteData(fid,'object',self,'class','inversion','fieldname','fatol','format','Double')
176 WriteData(fid,'object',self,'class','inversion','fieldname','frtol','format','Double')
177 WriteData(fid,'object',self,'class','inversion','fieldname','gatol','format','Double')
178 WriteData(fid,'object',self,'class','inversion','fieldname','grtol','format','Double')
179 WriteData(fid,'object',self,'class','inversion','fieldname','gttol','format','Double')
180 WriteData(fid,'object',self,'class','inversion','fieldname','algorithm','format','String')
181 WriteData(fid,'object',self,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1)
182 WriteData(fid,'object',self,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3)
183 WriteData(fid,'object',self,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3)
184 WriteData(fid,'object',self,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts)
185 WriteData(fid,'object',self,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts)
186 WriteData(fid,'object',self,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts)
187 WriteData(fid,'object',self,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1)
188 WriteData(fid,'object',self,'class','inversion','fieldname','surface_obs','format','DoubleMat','mattype',1)
189
190 #process control parameters
191 num_control_parameters = numpy.numel(self.control_parameters)
192 data = numpy.array([StringToEnum(self.control_parameter[0]) for control_parameter in self.control_parameters]).reshape(1,-1)
193 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3)
194 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer')
195
196 #process cost functions
197 num_cost_functions = numpy.size(self.cost_functions,2)
198 data= marshallcostfunctions(self.cost_functions)
199 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3)
200 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer')
Note: See TracBrowser for help on using the repository browser.