| [19895] | 1 | import numpy
|
|---|
| 2 | from project3d import project3d
|
|---|
| 3 | from WriteData import WriteData
|
|---|
| 4 | from checkfield import checkfield
|
|---|
| 5 | from fielddisplay import fielddisplay
|
|---|
| 6 | from IssmConfig import IssmConfig
|
|---|
| 7 | from EnumDefinitions import *
|
|---|
| 8 | from marshallcostfunctions import marshallcostfunctions
|
|---|
| 9 |
|
|---|
| 10 |
|
|---|
| 11 | class 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():
|
|---|
| [19898] | 158 | md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
|
|---|
| [19895] | 159 | elif solution==BalancethicknessSoftSolutionEnum():
|
|---|
| [19898] | 160 | md = checkfield(md,'fieldname','inversion.thickness_obs','size',[md.mesh.numberofvertices, 1],'NaN',1,'Inf',1)
|
|---|
| [19895] | 161 | else:
|
|---|
| [19898] | 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)
|
|---|
| [19895] | 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')
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
