0001 function md=macayealcontrol(md)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 if ((nargin~=1) || (nargout~=1)),
0012 velfinderusage();
0013 error('macayealcontrol error message: incorrect number of input and output arguments');
0014 end
0015
0016 if ~isa(md,'model'),
0017 macayealcontrolusage();
0018 error('macayealcontrol error message: input argument is not a @model object');
0019 end
0020
0021
0022 if ~strcmp(md.control_type,'B')
0023 error('macayealcontrol error message: only ''flow'' law inversion supported yet');
0024 end
0025
0026
0027 x=md.x;
0028 y=md.y;
0029 index=md.elements;
0030 index=sort(index,2);
0031 nods=md.numberofgrids;
0032 nel=md.numberofelements;
0033 z_thick=md.thickness;
0034 z_surf=md.surface;
0035 z_bed=md.bed;
0036 z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3;
0037 rho_ice=md.rho_ice;
0038 rho_water=md.rho_water;
0039 g=md.g;
0040 index_icefront=md.segmentonneumann_diag; index_icefront=index_icefront(:,1:2);
0041 nodes_on_boundary=md.gridonboundary;
0042 nodes_on_dirichlet=md.gridondirichlet_diag;
0043 nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1;
0044 nodes_on_iceshelf=md.gridoniceshelf;
0045
0046 criterion=md.eps_rel;
0047 yts=md.yts;
0048 tolx=md.tolx;
0049 maxiter=md.maxiter(1);
0050 B_ini=md.B;
0051 drag_coeff_ini=md.drag;
0052 glen_coeff=md.n;
0053 vx_obs=md.vx_obs/md.yts;
0054 vy_obs=md.vy_obs/md.yts;
0055 nsteps=md.nsteps;
0056 B_to_p=10^-3*yts^(-1/3);
0057
0058
0059 [length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront);
0060
0061
0062 [alpha, beta, gamma, area]=shape(index,x,y,nel,nods);
0063
0064 [matrix_bar, matrix_xbar, matrix_ybar]=...
0065 bar_maker(nel,nods,index,alpha,beta);
0066
0067
0068 create_el2nod_matrices
0069 old_direction=zeros(nods,1);
0070
0071
0072
0073 z_thick_bar=matrix_bar*z_thick;
0074 z_surf_bar=matrix_bar*z_surf;
0075 z_surf_xbar=matrix_xbar*z_surf;
0076 z_surf_ybar=matrix_ybar*z_surf;
0077 z_surf_x=el2nod\(el2nodRhs*z_surf_xbar);
0078 z_surf_y=el2nod\(el2nodRhs*z_surf_ybar);
0079 z_bed_bar=matrix_bar*z_bed;
0080 z_bed_xbar=matrix_xbar*z_bed;
0081 z_bed_ybar=matrix_ybar*z_bed;
0082
0083
0084
0085
0086 data_elements=[1:nel]';
0087 data_nodes=ones(nods,1);
0088
0089
0090 J=zeros(nsteps,1);
0091
0092
0093
0094 weighting=ones(nods,1);
0095
0096
0097 options=optimset('fminbnd');
0098 options=optimset(options,'Display','iter');
0099 options=optimset(options,'MaxFunEvals',maxiter);
0100 options=optimset(options,'MaxIter',100);
0101 options=optimset(options,'TolX',tolx);
0102
0103
0104 setup_control
0105
0106
0107 B_bar_ini=(B_ini(index(:,1))+B_ini(index(:,2))+B_ini(index(:,3)))/3;
0108 B=B_ini;
0109 B_bar=B_bar_ini;
0110
0111
0112 drag_type=md.drag_type;
0113 qcoeff=md.q;
0114 pcoeff=md.p;
0115 if strcmp(md.control_type,'basal') & (drag_type~=2), error('md.drag_type must be 2 for control methods'); end
0116 drag_coeff_bar_ini=(drag_coeff_ini(index(:,1))+drag_coeff_ini(index(:,2))+drag_coeff_ini(index(:,3)))/3;
0117 drag_coeff=drag_coeff_ini;
0118 drag_coeff_bar=drag_coeff_bar_ini;
0119
0120
0121 if strcmp(md.control_type,'basal') & (length(find(drag_coeff))==0),
0122 disp(sprintf('\n No drag pure for ice shelves => STOP'))
0123 return
0124 end
0125
0126
0127 nu_bar=viscosity(index,nel,alpha,beta,vx_obs,vy_obs,B_bar,glen_coeff);
0128
0129
0130
0131 disp('calculating the velocity with the initial parameters')
0132 c_velfinder;
0133
0134 if strcmp(md.control_type{1},'B'),
0135 for niteration=1:nsteps,
0136
0137 disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]);
0138
0139
0140 [u,v,adjointu,adjointv,direction]= ...
0141 grad_J_flow(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,...
0142 index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar);
0143
0144 if niteration==1
0145 scrsz = get(0,'ScreenSize');
0146 figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
0147 end
0148 plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',...
0149 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',...
0150 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',...
0151 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#4',[0 100],'figure',1);pause(1);
0152
0153
0154 u_objective=u;
0155 v_objective=v;
0156
0157
0158 direction=real(direction);
0159 direction=direction/sqrt(direction'*direction);
0160
0161
0162
0163 direction=direction-(direction'*old_direction)*old_direction;
0164 old_direction=direction;
0165
0166
0167 if niteration==1
0168 scrsz = get(0,'ScreenSize');
0169 figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
0170 end
0171 plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1);
0172
0173
0174
0175
0176 if abs(max(direction))>abs(min(direction))
0177 direction=10^7*direction/abs(max(direction));
0178 else
0179 direction=10^7*direction/abs(min(direction));
0180 end
0181
0182
0183
0184 if niteration<=2,
0185 upperbound=20;
0186 lowerbound=0;
0187 else
0188 upperbound=10;
0189 lowerbound=0;
0190 end
0191
0192
0193
0194
0195 search_constant=fminbnd('objectivefunction_C_flow',lowerbound,upperbound, ...
0196 options,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ...
0197 vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion);
0198
0199
0200
0201 B_old=B;
0202 B_new=B+search_constant*direction;
0203 B=B_new;
0204 pos=find(B<0);B(pos)=-B(pos);
0205
0206
0207 B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
0208
0209
0210 if niteration==1
0211 scrsz = get(0,'ScreenSize');
0212 figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
0213 end
0214 plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'colorbar','on','figure',3);pause(1);
0215
0216
0217
0218
0219 load F_file
0220
0221 J(niteration)=objectivefunction_C_flow(0,B,glen_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ...
0222 vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion);
0223 disp(J(niteration));
0224
0225
0226 if(mod(niteration,5)==0),
0227 save temporary_control_results_flow B B_bar u v J direction
0228 end
0229 end
0230
0231
0232 md.cont_J=J;
0233 md.cont_parameters=B;
0234
0235 elseif strcmp(md.control_type,'basal'),
0236 for niteration=1:nsteps,
0237
0238 disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]);
0239
0240
0241 [u,v,adjointu,adjointv,direction]= ...
0242 grad_J_drag(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vx_obs,...
0243 index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,drag_coeff,drag_coeff_bar);
0244
0245 if niteration==1
0246 scrsz = get(0,'ScreenSize');
0247 figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3])
0248 end
0249 plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',...
0250 'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',...
0251 'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',...
0252 'data',100*abs(sqrt(vx_obs.^2+vy_obs.^2)-sqrt(u.^2+v.^2))./sqrt(vx_obs.^2+vy_obs.^2),'title','Relative misfit','caxis#3',[0 100],'figure',1);pause(1);
0253
0254
0255 u_objective=u;
0256 v_objective=v;
0257
0258
0259 direction=real(direction);
0260 direction=direction/sqrt(direction'*direction);
0261 pos=find(nodes_on_iceshelf);
0262 direction(pos)=0;
0263
0264
0265
0266 direction=direction-(direction'*old_direction)*old_direction;
0267 old_direction=direction;
0268
0269
0270 if niteration==1
0271 scrsz = get(0,'ScreenSize');
0272 figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
0273 end
0274 plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);pause(1);
0275
0276
0277
0278
0279 if abs(max(direction))>abs(min(direction))
0280 direction=10^4*direction/abs(max(direction));
0281 else
0282 direction=10^4*direction/abs(min(direction));
0283 end
0284
0285
0286
0287 if niteration<=2,
0288 upperbound=2;
0289 lowerbound=-2;
0290 else
0291 upperbound=1;
0292 lowerbound=-1;
0293 end
0294
0295
0296
0297
0298 search_constant=fminbnd('objectivefunction_C_drag',lowerbound,upperbound, ...
0299 options,B,glen_coeff,drag_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ...
0300 vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion);
0301
0302
0303
0304 drag_coeff_old=drag_coeff;
0305 drag_coeff_new=drag_coeff+search_constant*direction;
0306 drag_coeff=drag_coeff_new;
0307
0308 if length(find(drag_coeff<0))~=0,
0309 disp(sprintf('\n Some basal drag coefficient negative => STOP'))
0310 break
0311 end
0312
0313
0314 if niteration==1
0315 scrsz = get(0,'ScreenSize');
0316 figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3])
0317 end
0318 drag_coeff_bar=(drag_coeff(index(:,1))+drag_coeff(index(:,2))+drag_coeff(index(:,3)))/3;
0319
0320
0321 plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);pause(1);
0322
0323
0324
0325
0326 load F_file
0327
0328 J(niteration)=objectivefunction_C_drag(0,B,glen_coeff,drag_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ...
0329 vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion);
0330 disp(J(niteration));
0331
0332
0333 if(mod(niteration,5)==0),
0334 save temporary_control_results_drag drag_coeff drag_coeff_bar u v J direction
0335 end
0336 end
0337 end
0338
0339 function macayealcontrolusage();
0340 disp('md=macayealcontrol(md)');
0341 disp(' where md is a structure of class @model');