source: issm/trunk/src/m/classes/taoinversion.m@ 16560

Last change on this file since 16560 was 16560, checked in by Mathieu Morlighem, 11 years ago

merged trunk-jpl and trunk for revision 16554
File size: 8.0 KB
RevLine 
[15396]1%TAOINVERSION class definition
[14165]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 nsteps = 0
12 cost_functions = NaN
13 cost_functions_coefficients = NaN
14 min_parameters = NaN
15 max_parameters = NaN
16 vx_obs = NaN
17 vy_obs = NaN
18 vz_obs = NaN
19 vel_obs = NaN
20 thickness_obs = NaN
21 end
22 methods
23 function obj = taoinversion(varargin) % {{{
24 switch nargin
25 case 0
26 obj=setdefaultparameters(obj);
[14170]27 case 1
28 if isa(varargin{1},'inversion'),
29 disp('converting inversion to taoinversion');
30 in=varargin{1};
31 obj.iscontrol = in.iscontrol;
32 obj.incomplete_adjoint = in.incomplete_adjoint;
33 obj.control_parameters = in.control_parameters;
34 obj.nsteps = in.nsteps;
35 obj.cost_functions = in.cost_functions(1,:); %Keep first line only
36 obj.cost_functions_coefficients = in.cost_functions_coefficients;
37 obj.min_parameters = in.min_parameters;
38 obj.max_parameters = in.max_parameters;
39 obj.vx_obs = in.vx_obs;
40 obj.vy_obs = in.vy_obs;
41 obj.vz_obs = in.vz_obs;
42 obj.vel_obs = in.vel_obs;
43 obj.thickness_obs = in.thickness_obs;
44 end
[14165]45 otherwise
46 error('constructor not supported');
47 end
48 end % }}}
49 function obj = setdefaultparameters(obj) % {{{
50
51 %default is incomplete adjoint for now
52 obj.incomplete_adjoint=1;
53
54 %parameter to be inferred by control methods (only
55 %drag and B are supported yet)
56 obj.control_parameters={'FrictionCoefficient'};
57
58 %number of steps in the control methods
59 obj.nsteps=20;
60
61 %several responses can be used:
62 obj.cost_functions=101*ones(obj.nsteps,1);
63
64 end % }}}
65 function md = checkconsistency(obj,md,solution,analyses) % {{{
66
67 %Early return
68 if ~obj.iscontrol, return; end
69
[16560]70 if ~IssmConfig('_HAVE_TAO_'),
71 md = checkmessage(md,['TAO has not been installed, ISSM needs to be reconfigured and recompiled with TAO']);
72 end
73
[14165]74 num_controls=numel(md.inversion.control_parameters);
75 num_costfunc=size(md.inversion.cost_functions,2);
76
77 md = checkfield(md,'inversion.iscontrol','values',[0 1]);
78 md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
79 md = checkfield(md,'inversion.control_parameters','cell',1,'values',...
80 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
[16560]81 md = checkfield(md,'inversion.nsteps','numel',1,'>=',0);
[14165]82 md = checkfield(md,'inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
83 md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
84 md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
85 md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
86
87 if solution==BalancethicknessSolutionEnum()
88 md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
[15396]89 elseif solution==BalancethicknessSoftSolutionEnum()
[14165]90 md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
91 else
92 md = checkfield(md,'inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
93 md = checkfield(md,'inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
94 end
95 end % }}}
96 function disp(obj) % {{{
97 disp(sprintf(' taoinversion parameters:'));
98 fielddisplay(obj,'iscontrol','is inversion activated?');
[15396]99 fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity');
100 fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
[14165]101 fielddisplay(obj,'nsteps','number of optimization searches');
102 fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step');
103 fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
104 fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
105 fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
[15396]106 fielddisplay(obj,'vx_obs','observed velocity x component [m/yr]');
107 fielddisplay(obj,'vy_obs','observed velocity y component [m/yr]');
108 fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]');
[14165]109 fielddisplay(obj,'thickness_obs','observed thickness [m]');
110 disp('Available cost functions:');
111 disp(' 101: SurfaceAbsVelMisfit');
112 disp(' 102: SurfaceRelVelMisfit');
113 disp(' 103: SurfaceLogVelMisfit');
114 disp(' 104: SurfaceLogVxVyMisfit');
115 disp(' 105: SurfaceAverageVelMisfit');
116 disp(' 201: ThicknessAbsMisfit');
117 disp(' 501: DragCoefficientAbsGradient');
118 disp(' 502: RheologyBbarAbsGradient');
119 disp(' 503: ThicknessAbsGradient');
120 end % }}}
[15396]121 function marshall(obj,md,fid) % {{{
[14165]122
[16137]123 yts=365.0*24.0*3600.0;
124
[14165]125 WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean');
126 WriteData(fid,'enum',InversionTaoEnum(),'data',true,'format','Boolean');
127 if ~obj.iscontrol, return; end
128 WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
129 WriteData(fid,'object',obj,'class','inversion','fieldname','nsteps','format','Integer');
130 WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
131 WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
132 WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
[16137]133 WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
134 WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
135 WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
[14165]136 WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
137
138 %process control parameters
139 num_control_parameters=numel(obj.control_parameters);
140 data=zeros(1,num_control_parameters);
141 for i=1:num_control_parameters,
142 data(i)=StringToEnum(obj.control_parameters{i});
143 end
144 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3);
145 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer');
146
147 %process cost functions
148 num_cost_functions=size(obj.cost_functions,2);
149 data=obj.cost_functions;
150 pos=find(data==101); data(pos)=SurfaceAbsVelMisfitEnum();
151 pos=find(data==102); data(pos)=SurfaceRelVelMisfitEnum();
152 pos=find(data==103); data(pos)=SurfaceLogVelMisfitEnum();
153 pos=find(data==104); data(pos)=SurfaceLogVxVyMisfitEnum();
154 pos=find(data==105); data(pos)=SurfaceAverageVelMisfitEnum();
155 pos=find(data==201); data(pos)=ThicknessAbsMisfitEnum();
156 pos=find(data==501); data(pos)=DragCoefficientAbsGradientEnum();
157 pos=find(data==502); data(pos)=RheologyBbarAbsGradientEnum();
158 pos=find(data==503); data(pos)=ThicknessAbsGradientEnum();
159 pos=find(data==504); data(pos)=ThicknessAlongGradientEnum();
160 pos=find(data==505); data(pos)=ThicknessAcrossGradientEnum();
161 pos=find(data==506); data(pos)=BalancethicknessMisfitEnum();
162 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
163 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
164 end % }}}
165 end
166end
Note: See TracBrowser for help on using the repository browser.