macayealcontrol

PURPOSE ^

MACAYEALCONTROL - launch a control method using MacAyeal solution

SYNOPSIS ^

function md=macayealcontrol(md)

DESCRIPTION ^

MACAYEALCONTROL - launch a control method using MacAyeal solution

   the routine is used for a control method. It determines the most adapted viscosity
   field so that the calculated velocity field is as close as possible to an observed velocity field

   Usage:
      md=macayealcontrol(md)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function md=macayealcontrol(md)
0002 %MACAYEALCONTROL - launch a control method using MacAyeal solution
0003 %
0004 %   the routine is used for a control method. It determines the most adapted viscosity
0005 %   field so that the calculated velocity field is as close as possible to an observed velocity field
0006 %
0007 %   Usage:
0008 %      md=macayealcontrol(md)
0009 
0010 %First check we do have the correct argument number
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 %Check that control is done on flow law, and not drag (not supported yet):
0022 if ~strcmp(md.control_type,'B')
0023     error('macayealcontrol error message: only ''flow'' law inversion supported yet');
0024 end
0025 
0026 %Transfer model fields into matlab variables
0027 x=md.x;
0028 y=md.y;
0029 index=md.elements;
0030 index=sort(index,2); %necessary
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); %we strip the last column, which holds the element number for the boundary segment
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; %From m/a to m/s
0054 vy_obs=md.vy_obs/md.yts;
0055 nsteps=md.nsteps;
0056 B_to_p=10^-3*yts^(-1/3);
0057 
0058 %Build length_icefront and normal_icefront:
0059 [length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront);
0060 
0061 %Building shape functions and derivative operators
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 %initialize some data
0068 create_el2nod_matrices
0069 old_direction=zeros(nods,1);
0070 
0071 
0072 %setup some fake distributions.
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 %data_elements and data_nodes are the elements and nodes where
0085 %we have observations
0086 data_elements=[1:nel]';
0087 data_nodes=ones(nods,1);
0088 
0089 %initialize misfit between model velocity and observed velocity
0090 J=zeros(nsteps,1);
0091 
0092 %Weighting: one areas where we don't want the control method to optimize
0093 %the misfit, we can set weighting to 0.
0094 weighting=ones(nods,1);
0095 
0096 %optimization parameters for the matlab functoin fminbnd
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 %build useful matrices.
0104 setup_control
0105 
0106 %setup initial flow law paramter distribution to startup the optimization.
0107 B_bar_ini=(B_ini(index(:,1))+B_ini(index(:,2))+B_ini(index(:,3)))/3; 
0108 B=B_ini;  %variables used in optimziation are B and B_bar.
0109 B_bar=B_bar_ini;
0110 
0111 %setup initial drag paramter distribution to startup the optimization.
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;  %variables used in optimziation are drag_coeff and drag_coeff_bar.
0118 drag_coeff_bar=drag_coeff_bar_ini;
0119 
0120 %check that the model is not a pure ice shelf (no friction on ice shelves)
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 %nu_bar=10^14*ones(nel,1);  %initial element viscosity distribution.
0127 nu_bar=viscosity(index,nel,alpha,beta,vx_obs,vy_obs,B_bar,glen_coeff);
0128 
0129 %%%AK velfinder;    %forward model that determines first velocity field used
0130 %to start optimization.
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        % Compute search direction -
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        %Keep track of u and v to use in objectivefunction_C:
0154        u_objective=u;
0155        v_objective=v;
0156 
0157        %Orthogonalize direction
0158        direction=real(direction);
0159        direction=direction/sqrt(direction'*direction);
0160 
0161 
0162        % rough orthagonalization
0163        direction=direction-(direction'*old_direction)*old_direction; 
0164        old_direction=direction;
0165 
0166        %visualize direction
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        % normalize direction to 10^7, so that when variations on B are computed
0175        %they will be significant.
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        %during optimization, bounds on B variations can vary. Here, they are less
0183        %strict in the first iteration.
0184        if niteration<=2,
0185           upperbound=20;
0186           lowerbound=0;
0187        else
0188           upperbound=10;
0189           lowerbound=0;
0190        end
0191 
0192        %search the multiplicative constant to the direction, that will
0193        %be used to modify B.
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         %update value of B
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        %Average B over elements:
0207        B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
0208 
0209        %visualize new distribution.
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        %evaluate new misfit.
0217        %@@@AK load F_file   %this file was created in objectivefunction_C, and is reloaded
0218        %here to win some computation time.
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       %do a backup every 5 iterations
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     %Load results onto md:
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         % Compute search direction -
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         %Keep track of u and v to use in objectivefunction_C:
0255         u_objective=u;
0256         v_objective=v;
0257 
0258         %Orthogonalize direction
0259        direction=real(direction);
0260        direction=direction/sqrt(direction'*direction);
0261        pos=find(nodes_on_iceshelf);
0262        direction(pos)=0;
0263 
0264 
0265        % rough orthagonalization
0266        direction=direction-(direction'*old_direction)*old_direction; 
0267        old_direction=direction;
0268 
0269        %visualize direction
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        % normalize direction to 10^4, so that when variations on drag are computed
0278        %they will be significant.
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        %during optimization, bounds on drag variations can vary. Here, they are less
0286        %strict in the first iteration.
0287        if niteration<=2,
0288           upperbound=2;
0289           lowerbound=-2;
0290        else
0291           upperbound=1;
0292           lowerbound=-1;
0293        end
0294 
0295        %search the multiplicative constant to the direction, that will
0296        %be used to modify drag.
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         %update value of drag
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        %Average drag over elements:
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        %visualize new distribution.
0321        plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);pause(1);
0322 
0323        %evaluate new misfit.
0324        %@@@AK load F_file   %this file was created in objectivefunction_C, and is reloaded
0325        %here to win some computation time.
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        %do a backup every 5 iterations
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');

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003