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

Last change on this file since 14310 was 14310, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 14308
File size: 7.8 KB
Line 
1%INVERSION 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 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);
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
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
70 num_controls=numel(md.inversion.control_parameters);
71 num_costfunc=size(md.inversion.cost_functions,2);
72
73 md = checkfield(md,'inversion.iscontrol','values',[0 1]);
74 md = checkfield(md,'inversion.incomplete_adjoint','values',[0 1]);
75 md = checkfield(md,'inversion.control_parameters','cell',1,'values',...
76 {'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'MaterialsRheologyZbar' 'Vx' 'Vy' 'Thickness'});
77 md = checkfield(md,'inversion.nsteps','numel',1,'>=',1);
78 md = checkfield(md,'inversion.cost_functions','size',[1 num_costfunc],'values',[101:105 201 501:506]);
79 md = checkfield(md,'inversion.cost_functions_coefficients','size',[md.mesh.numberofvertices num_costfunc],'>=',0);
80 md = checkfield(md,'inversion.min_parameters','size',[md.mesh.numberofvertices num_controls]);
81 md = checkfield(md,'inversion.max_parameters','size',[md.mesh.numberofvertices num_controls]);
82
83 if solution==BalancethicknessSolutionEnum()
84 md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
85 elseif solution==WeakBalancethicknessSolutionEnum()
86 md = checkfield(md,'inversion.thickness_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
87 else
88 md = checkfield(md,'inversion.vx_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
89 md = checkfield(md,'inversion.vy_obs','size',[md.mesh.numberofvertices 1],'NaN',1);
90 end
91 end % }}}
92 function disp(obj) % {{{
93 disp(sprintf(' taoinversion parameters:'));
94 fielddisplay(obj,'iscontrol','is inversion activated?');
95 fielddisplay(obj,'incomplete_adjoint','do we assume linear viscosity?');
96 fielddisplay(obj,'control_parameters','parameter where inverse control is carried out; ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
97 fielddisplay(obj,'nsteps','number of optimization searches');
98 fielddisplay(obj,'cost_functions','indicate the type of response for each optimization step');
99 fielddisplay(obj,'cost_functions_coefficients','cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter');
100 fielddisplay(obj,'min_parameters','absolute minimum acceptable value of the inversed parameter on each vertex');
101 fielddisplay(obj,'max_parameters','absolute maximum acceptable value of the inversed parameter on each vertex');
102 fielddisplay(obj,'vx_obs','observed velocity x component [m/a]');
103 fielddisplay(obj,'vy_obs','observed velocity y component [m/a]');
104 fielddisplay(obj,'vel_obs','observed velocity magnitude [m/a]');
105 fielddisplay(obj,'thickness_obs','observed thickness [m]');
106 disp('Available cost functions:');
107 disp(' 101: SurfaceAbsVelMisfit');
108 disp(' 102: SurfaceRelVelMisfit');
109 disp(' 103: SurfaceLogVelMisfit');
110 disp(' 104: SurfaceLogVxVyMisfit');
111 disp(' 105: SurfaceAverageVelMisfit');
112 disp(' 201: ThicknessAbsMisfit');
113 disp(' 501: DragCoefficientAbsGradient');
114 disp(' 502: RheologyBbarAbsGradient');
115 disp(' 503: ThicknessAbsGradient');
116 end % }}}
117 function marshall(obj,fid) % {{{
118
119 WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean');
120 WriteData(fid,'enum',InversionTaoEnum(),'data',true,'format','Boolean');
121 if ~obj.iscontrol, return; end
122 WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
123 WriteData(fid,'object',obj,'class','inversion','fieldname','nsteps','format','Integer');
124 WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
125 WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
126 WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
127 WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1);
128 WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1);
129 WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1);
130 WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
131
132 %process control parameters
133 num_control_parameters=numel(obj.control_parameters);
134 data=zeros(1,num_control_parameters);
135 for i=1:num_control_parameters,
136 data(i)=StringToEnum(obj.control_parameters{i});
137 end
138 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3);
139 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer');
140
141 %process cost functions
142 num_cost_functions=size(obj.cost_functions,2);
143 data=obj.cost_functions;
144 pos=find(data==101); data(pos)=SurfaceAbsVelMisfitEnum();
145 pos=find(data==102); data(pos)=SurfaceRelVelMisfitEnum();
146 pos=find(data==103); data(pos)=SurfaceLogVelMisfitEnum();
147 pos=find(data==104); data(pos)=SurfaceLogVxVyMisfitEnum();
148 pos=find(data==105); data(pos)=SurfaceAverageVelMisfitEnum();
149 pos=find(data==201); data(pos)=ThicknessAbsMisfitEnum();
150 pos=find(data==501); data(pos)=DragCoefficientAbsGradientEnum();
151 pos=find(data==502); data(pos)=RheologyBbarAbsGradientEnum();
152 pos=find(data==503); data(pos)=ThicknessAbsGradientEnum();
153 pos=find(data==504); data(pos)=ThicknessAlongGradientEnum();
154 pos=find(data==505); data(pos)=ThicknessAcrossGradientEnum();
155 pos=find(data==506); data(pos)=BalancethicknessMisfitEnum();
156 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
157 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
158 end % }}}
159 end
160end
Note: See TracBrowser for help on using the repository browser.