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