Changeset 408

Timestamp:

05/13/09 17:24:05 (16 years ago)

Author:

Mathieu Morlighem

Message:

minor

Location:

issm/trunk/src/m/solutions/macayeal

Files:

• control.m (modified) (7 diffs)
• grad_J_flow.m (modified) (1 diff)

issm/trunk/src/m/solutions/macayeal/control.m

@@ -20 +20 @@
 
 %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 +48 @@
 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 +71 @@
 create_el2nod_matrices
 old_direction=zeros(nods,1);
-
 
 %setup some fake distributions.
@@ -81 +83 @@
 z_bed_ybar=matrix_ybar*z_bed;
 
-
 %data_elements and data_nodes are the elements and nodes where 
 %we have observations
@@ -104 +105 @@
 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 +246 @@
     md.cont_parameter=B;
 
-elseif strcmp(md.control_type,'basal'),
+elseif strcmp(md.control_type,'drag'),
     for niteration=1:nsteps,
 
@@ -264 +271 @@
         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

issm/trunk/src/m/solutions/macayeal/grad_J_flow.m

@@ -56 +56 @@
            (  2*beta(:,k).*v(index(:,k))  +  alpha(:,k).*u(index(:,k))  ).*2.*  beta(:,l).*adjointv(index(:,l))...
            );
-        
 
       end