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

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

merged trunk-jpl and trunk for revision 16135
File size: 7.9 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
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);
[15396]85 elseif solution==BalancethicknessSoftSolutionEnum()
[14165]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?');
[15396]95 fielddisplay(obj,'incomplete_adjoint','1: linear viscosity, 0: non-linear viscosity');
96 fielddisplay(obj,'control_parameters','ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''}');
[14165]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');
[15396]102 fielddisplay(obj,'vx_obs','observed velocity x component [m/yr]');
103 fielddisplay(obj,'vy_obs','observed velocity y component [m/yr]');
104 fielddisplay(obj,'vel_obs','observed velocity magnitude [m/yr]');
[14165]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 % }}}
[15396]117 function marshall(obj,md,fid) % {{{
[14165]118
[16137]119 yts=365.0*24.0*3600.0;
120
[14165]121 WriteData(fid,'object',obj,'class','inversion','fieldname','iscontrol','format','Boolean');
122 WriteData(fid,'enum',InversionTaoEnum(),'data',true,'format','Boolean');
123 if ~obj.iscontrol, return; end
124 WriteData(fid,'object',obj,'class','inversion','fieldname','incomplete_adjoint','format','Boolean');
125 WriteData(fid,'object',obj,'class','inversion','fieldname','nsteps','format','Integer');
126 WriteData(fid,'object',obj,'class','inversion','fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
127 WriteData(fid,'object',obj,'class','inversion','fieldname','min_parameters','format','DoubleMat','mattype',3);
128 WriteData(fid,'object',obj,'class','inversion','fieldname','max_parameters','format','DoubleMat','mattype',3);
[16137]129 WriteData(fid,'object',obj,'class','inversion','fieldname','vx_obs','format','DoubleMat','mattype',1,'scale',1./yts);
130 WriteData(fid,'object',obj,'class','inversion','fieldname','vy_obs','format','DoubleMat','mattype',1,'scale',1./yts);
131 WriteData(fid,'object',obj,'class','inversion','fieldname','vz_obs','format','DoubleMat','mattype',1,'scale',1./yts);
[14165]132 WriteData(fid,'object',obj,'class','inversion','fieldname','thickness_obs','format','DoubleMat','mattype',1);
133
134 %process control parameters
135 num_control_parameters=numel(obj.control_parameters);
136 data=zeros(1,num_control_parameters);
137 for i=1:num_control_parameters,
138 data(i)=StringToEnum(obj.control_parameters{i});
139 end
140 WriteData(fid,'data',data,'enum',InversionControlParametersEnum(),'format','DoubleMat','mattype',3);
141 WriteData(fid,'data',num_control_parameters,'enum',InversionNumControlParametersEnum(),'format','Integer');
142
143 %process cost functions
144 num_cost_functions=size(obj.cost_functions,2);
145 data=obj.cost_functions;
146 pos=find(data==101); data(pos)=SurfaceAbsVelMisfitEnum();
147 pos=find(data==102); data(pos)=SurfaceRelVelMisfitEnum();
148 pos=find(data==103); data(pos)=SurfaceLogVelMisfitEnum();
149 pos=find(data==104); data(pos)=SurfaceLogVxVyMisfitEnum();
150 pos=find(data==105); data(pos)=SurfaceAverageVelMisfitEnum();
151 pos=find(data==201); data(pos)=ThicknessAbsMisfitEnum();
152 pos=find(data==501); data(pos)=DragCoefficientAbsGradientEnum();
153 pos=find(data==502); data(pos)=RheologyBbarAbsGradientEnum();
154 pos=find(data==503); data(pos)=ThicknessAbsGradientEnum();
155 pos=find(data==504); data(pos)=ThicknessAlongGradientEnum();
156 pos=find(data==505); data(pos)=ThicknessAcrossGradientEnum();
157 pos=find(data==506); data(pos)=BalancethicknessMisfitEnum();
158 WriteData(fid,'data',data,'enum',InversionCostFunctionsEnum(),'format','DoubleMat','mattype',3);
159 WriteData(fid,'data',num_cost_functions,'enum',InversionNumCostFunctionsEnum(),'format','Integer');
160 end % }}}
161 end
162end
Note: See TracBrowser for help on using the repository browser.