Index: /issm/trunk/src/m/solutions/macayeal/control.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/control.m (revision 407) +++ /issm/trunk/src/m/solutions/macayeal/control.m (revision 408) @@ -20,7 +20,7 @@ %Check that control is done on flow law, and not drag (not supported yet): -if ~strcmp(md.control_type,'B') - error('macayealcontrol error message: only ''flow'' law inversion supported yet'); -end +%if ~strcmp(md.control_type,'B') +% error('macayealcontrol error message: only ''B'' inversion supported yet'); +%end %Transfer model fields into matlab variables @@ -48,11 +48,14 @@ tolx=md.tolx; maxiter=md.maxiter(1); -B_ini=md.B; -drag_coeff_ini=md.drag; +if strcmp(md.control_type,'B') + B_ini=md.B; + B_to_p=10^-3*yts^(-1/3); +else + drag_coeff_ini=md.drag; +end glen_coeff=md.n; vx_obs=md.vx_obs/md.yts; %From m/a to m/s vy_obs=md.vy_obs/md.yts; nsteps=md.nsteps; -B_to_p=10^-3*yts^(-1/3); %Build length_icefront and normal_icefront: @@ -68,5 +71,4 @@ create_el2nod_matrices old_direction=zeros(nods,1); - %setup some fake distributions. @@ -81,5 +83,4 @@ z_bed_ybar=matrix_ybar*z_bed; - %data_elements and data_nodes are the elements and nodes where %we have observations @@ -104,20 +105,26 @@ setup_control -%setup initial flow law paramter distribution to startup the optimization. -B_bar_ini=(B_ini(index(:,1))+B_ini(index(:,2))+B_ini(index(:,3)))/3; -B=B_ini; %variables used in optimziation are B and B_bar. -B_bar=B_bar_ini; - %setup initial drag paramter distribution to startup the optimization. drag_type=md.drag_type; qcoeff=md.q; pcoeff=md.p; -if strcmp(md.control_type,'basal') & (drag_type~=2), error('md.drag_type must be 2 for control methods'); end -drag_coeff_bar_ini=(drag_coeff_ini(index(:,1))+drag_coeff_ini(index(:,2))+drag_coeff_ini(index(:,3)))/3; -drag_coeff=drag_coeff_ini; %variables used in optimziation are drag_coeff and drag_coeff_bar. -drag_coeff_bar=drag_coeff_bar_ini; + +if strcmp(md.control_type,'B') + %setup initial flow law paramter distribution to startup the optimization. + B_bar_ini=(B_ini(index(:,1))+B_ini(index(:,2))+B_ini(index(:,3)))/3; + B=B_ini; %variables used in optimziation are B and B_bar. + B_bar=B_bar_ini; +else + B_bar=(md.B(index(:,1))+md.B(index(:,2))+md.B(index(:,3)))/3; + if (drag_type~=2), error('md.drag_type must be 2 for control methods'); end + drag_coeff_bar_ini=(drag_coeff_ini(index(:,1))+drag_coeff_ini(index(:,2))+drag_coeff_ini(index(:,3)))/3; + drag_coeff=drag_coeff_ini; %variables used in optimziation are drag_coeff and drag_coeff_bar. + drag_coeff_bar=drag_coeff_bar_ini; +end + + %check that the model is not a pure ice shelf (no friction on ice shelves) -if strcmp(md.control_type,'basal') & (length(find(drag_coeff))==0), +if strcmp(md.control_type,'drag') & (length(find(drag_coeff))==0), disp(sprintf('\n No drag pure for ice shelves => STOP')) return @@ -239,5 +246,5 @@ md.cont_parameter=B; -elseif strcmp(md.control_type,'basal'), +elseif strcmp(md.control_type,'drag'), for niteration=1:nsteps, @@ -264,75 +271,70 @@ v_objective=v; - %Orthogonalize direction - direction=real(direction); - direction=direction/sqrt(direction'*direction); - pos=find(nodes_on_iceshelf); - direction(pos)=0; - - - % rough orthagonalization - direction=direction-(direction'*old_direction)*old_direction; - old_direction=direction; - - %visualize direction - if md.plot==1, - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow; - end - - - % normalize direction to 10^4, so that when variations on drag are computed - %they will be significant. - if abs(max(direction))>abs(min(direction)) - direction=10^4*direction/abs(max(direction)); - else - direction=10^4*direction/abs(min(direction)); - end - - %during optimization, bounds on drag variations can vary. Here, they are less - %strict in the first iteration. - if niteration<=2, - upperbound=2; - lowerbound=-2; - else - upperbound=1; - lowerbound=-1; - end - - %search the multiplicative constant to the direction, that will - %be used to modify drag. - - search_constant=fminbnd('objectivefunction_C_drag',lowerbound,upperbound, ... - options,B,glen_coeff,drag_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ... - vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion); - - - %update value of drag - drag_coeff_old=drag_coeff; - drag_coeff_new=drag_coeff+search_constant*direction; - drag_coeff=drag_coeff_new; - - if length(find(drag_coeff<0))~=0, - disp(sprintf('\n Some basal drag coefficient negative => STOP')) - break - end - - %Average drag over elements: - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) - end - drag_coeff_bar=(drag_coeff(index(:,1))+drag_coeff(index(:,2))+drag_coeff(index(:,3)))/3; - - %visualize new distribution. - if md.plot==1, - plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);drawnow; - end - - %evaluate new misfit. - %@@@AK load F_file %this file was created in objectivefunction_C, and is reloaded + %Orthogonalize direction + direction=real(direction); + direction=direction/sqrt(direction'*direction); + pos=find(nodes_on_iceshelf); + direction(pos)=0; + + + % rough orthagonalization + direction=direction-(direction'*old_direction)*old_direction; + old_direction=direction; + + %visualize direction + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[10 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',direction,'title','Orthogonal direction','figure',2);drawnow; + end + + + % normalize direction to 50, so that when variations on drag are computed + %they will be significant. + direction=50*direction/max(abs(direction)); + + %during optimization, bounds on drag variations can vary. Here, they are less + %strict in the first iteration. + if niteration<=2, + upperbound=2; + lowerbound=0; + else + upperbound=1; + lowerbound=0; + end + + %search the multiplicative constant to the direction, that will + %be used to modify drag. + + search_constant=fminbnd('objectivefunction_C_drag',lowerbound,upperbound, ... + options,B,glen_coeff,drag_coeff,direction,Rhs,S,F,P,specified_velocity,nods,nel,vx_obs, ... + vy_obs,u_objective, v_objective, index,alpha, beta, area,weighting,rowD,colD,rowDshort,colDshort,valueuu,valuevv,valueuv,matrix_bar,criterion); + + %update value of drag + drag_coeff_old=drag_coeff; + drag_coeff_new=drag_coeff+search_constant*direction; + drag_coeff=drag_coeff_new; + + if length(find(drag_coeff<0))~=0, + disp(sprintf('\n Some basal drag coefficient negative => STOP')) + break + end + + %Average drag over elements: + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)*1.97/3 scrsz(4)*1/3 scrsz(3)/3 scrsz(4)/3]) + end + drag_coeff_bar=(drag_coeff(index(:,1))+drag_coeff(index(:,2))+drag_coeff(index(:,3)))/3; + + %visualize new distribution. + if md.plot==1, + plotmodel(md,'data',drag_coeff,'title',['Drag at iteration ' num2str(niteration)],'caxis',[0 1000],'colorbar','on','figure',3);drawnow; + end + + %evaluate new misfit. + %@@@AK load F_file %this file was created in objectivefunction_C, and is reloaded %here to win some computation time. load F_file Index: /issm/trunk/src/m/solutions/macayeal/grad_J_flow.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/grad_J_flow.m (revision 407) +++ /issm/trunk/src/m/solutions/macayeal/grad_J_flow.m (revision 408) @@ -56,5 +56,4 @@ ( 2*beta(:,k).*v(index(:,k)) + alpha(:,k).*u(index(:,k)) ).*2.* beta(:,l).*adjointv(index(:,l))... ); - end