Index: /issm/trunk/src/m/solutions/macayeal/control.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/control.m (revision 35) +++ /issm/trunk/src/m/solutions/macayeal/control.m (revision 35) @@ -0,0 +1,353 @@ +function md=control(md) +%CONTROL - 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=control(md) + +%First check we do have the correct argument number +if ((nargin~=1) || (nargout~=1)), + velfinderusage(); + error('macayealcontrol error message: incorrect number of input and output arguments'); +end + +if ~isa(md,'model'), + macayealcontrolusage(); + error('macayealcontrol error message: input argument is not a @model object'); +end + +%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 + +%Transfer model fields into matlab variables +x=md.x; +y=md.y; +index=md.elements; +index=sort(index,2); %necessary +nods=md.numberofgrids; +nel=md.numberofelements; +z_thick=md.thickness; +z_surf=md.surface; +z_bed=md.bed; +z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3; +rho_ice=md.rho_ice; +rho_water=md.rho_water; +g=md.g; +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 +nodes_on_boundary=md.gridonboundary; +nodes_on_dirichlet=md.gridondirichlet_diag; +nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1; +nodes_on_iceshelf=md.gridoniceshelf; + +criterion=md.eps_rel; +yts=md.yts; +tolx=md.tolx; +maxiter=md.maxiter(1); +B_ini=md.B; +drag_coeff_ini=md.drag; +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: +[length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront); + +%Building shape functions and derivative operators +[alpha, beta, gamma, area]=shape(index,x,y,nel,nods); + +[matrix_bar, matrix_xbar, matrix_ybar]=... + bar_maker(nel,nods,index,alpha,beta); + +%initialize some data +create_el2nod_matrices +old_direction=zeros(nods,1); + + +%setup some fake distributions. +z_thick_bar=matrix_bar*z_thick; +z_surf_bar=matrix_bar*z_surf; +z_surf_xbar=matrix_xbar*z_surf; +z_surf_ybar=matrix_ybar*z_surf; +z_surf_x=el2nod\(el2nodRhs*z_surf_xbar); +z_surf_y=el2nod\(el2nodRhs*z_surf_ybar); +z_bed_bar=matrix_bar*z_bed; +z_bed_xbar=matrix_xbar*z_bed; +z_bed_ybar=matrix_ybar*z_bed; + + +%data_elements and data_nodes are the elements and nodes where +%we have observations +data_elements=[1:nel]'; +data_nodes=ones(nods,1); + +%initialize misfit between model velocity and observed velocity +J=zeros(nsteps,1); + +%Weighting: one areas where we don't want the control method to optimize +%the misfit, we can set weighting to 0. +weighting=ones(nods,1); + +%optimization parameters for the matlab functoin fminbnd +options=optimset('fminbnd'); +options=optimset(options,'Display','iter'); +options=optimset(options,'MaxFunEvals',maxiter); +options=optimset(options,'MaxIter',100); +options=optimset(options,'TolX',tolx); + +%build useful matrices. +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; + +%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), + disp(sprintf('\n No drag pure for ice shelves => STOP')) + return +end + +%nu_bar=10^14*ones(nel,1); %initial element viscosity distribution. +nu_bar=viscosity(index,nel,alpha,beta,vx_obs,vy_obs,B_bar,glen_coeff); + +%%%AK velfinder; %forward model that determines first velocity field used +%to start optimization. +disp('calculating the velocity with the initial parameters') +c_velfinder; + +if strcmp(md.control_type{1},'B'), + for niteration=1:nsteps, + + disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]); + + % Compute search direction - + [u,v,adjointu,adjointv,direction]= ... + grad_J_flow(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,... + index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar); + + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',... + 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',... + 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',... + '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);drawnow; + end + + %Keep track of u and v to use in objectivefunction_C: + u_objective=u; + v_objective=v; + + %Orthogonalize direction + direction=real(direction); + direction=direction/sqrt(direction'*direction); + + + % 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^7, so that when variations on B are computed + %they will be significant. + if abs(max(direction))>abs(min(direction)) + direction=10^7*direction/abs(max(direction)); + else + direction=10^7*direction/abs(min(direction)); + end + + %during optimization, bounds on B variations can vary. Here, they are less + %strict in the first iteration. + if niteration<=2, + upperbound=20; + lowerbound=0; + else + upperbound=10; + lowerbound=0; + end + + %search the multiplicative constant to the direction, that will + %be used to modify B. + + search_constant=fminbnd('objectivefunction_C_flow',lowerbound,upperbound, ... + options,B,glen_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 B + B_old=B; + B_new=B+search_constant*direction; + B=B_new; + pos=find(B<0);B(pos)=-B(pos); + + %Average B over elements: + B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3; + + %visualize new distribution. + if md.plot==1, + 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 + plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'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 + + J(niteration)=objectivefunction_C_flow(0,B,glen_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); + disp(J(niteration)); + + %do a backup every 5 iterations + if(mod(niteration,5)==0), + save temporary_control_results_flow B B_bar u v J direction + end + end + + %Load results onto md: + md.cont_J=J; + md.cont_parameters=B; + +elseif strcmp(md.control_type,'basal'), + for niteration=1:nsteps, + + disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]); + + % Compute search direction - + [u,v,adjointu,adjointv,direction]= ... + grad_J_drag(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vx_obs,... + index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,drag_coeff,drag_coeff_bar); + + if md.plot==1, + if niteration==1 + scrsz = get(0,'ScreenSize'); + figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) + end + plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',... + 'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',... + 'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',... + '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);drawnow; + end + + %Keep track of u and v to use in objectivefunction_C: + u_objective=u; + 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 + %here to win some computation time. + load F_file + + J(niteration)=objectivefunction_C_drag(0,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); + disp(J(niteration)); + + %do a backup every 5 iterations + if(mod(niteration,5)==0), + save temporary_control_results_drag drag_coeff drag_coeff_bar u v J direction + end + end +end + +function macayealcontrolusage(); +disp('md=macayealcontrol(md)'); +disp(' where md is a structure of class @model'); Index: /issm/trunk/src/m/solutions/macayeal/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/diagnostic.m (revision 35) +++ /issm/trunk/src/m/solutions/macayeal/diagnostic.m (revision 35) @@ -0,0 +1,500 @@ +function md=diagnostic(md) +%DIAGNOSTIC - compute the velocity field of a 2d model using MacAyeal solution +% +% this routine solves the problem using MacAyeal's model. It calculates the velocity +% field corresponding to the parameters and the geometry given by the model md +% +% Usage: +% md=diagnostic(md) + +%First check we do have the correct argument number +if ((nargin~=1) || (nargout~=1)), + macayealdiagnosticusage(); + error('macayealdiagnostic error message: incorrect number of input and output arguments'); +end + +if ~isa(md,'model'), + macayealdiagnosticusage(); + error('macayealdiagnostic error message: input argument is not a @model object'); +end + +%start timing +t1=clock; + +%Transfer model fields into matlab variables +x=md.x; +y=md.y; + +index=md.elements;index=sort(index,2); %necessary +nods=md.numberofgrids; +nel=md.numberofelements; +z_thick=md.thickness; +z_surf=md.surface; +z_bed=md.bed; +z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3; +rho_ice=md.rho_ice; +rho_water=md.rho_water; +g=md.g; +viscosity_overshoot=md.viscosity_overshoot; +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 +nodes_on_boundary=md.gridonboundary; +nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1; +node_on_dirichlet=md.gridondirichlet_diag; +nodes_on_icesheet=md.gridonicesheet; +element_on_icesheet=md.elementonicesheet; +qcoeff=md.q; +pcoeff=md.p; +drag_type=md.drag_type; +drag=md.drag; + +criterion_rel=md.eps_rel; +criterion_abs=md.eps_abs; +yts=md.yts; +B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3; +glen_coeff=md.n; + +%average of p and q over the grids (size nel->nods) +pcoeff_grid=zeros(nods,1); +qcoeff_grid=zeros(nods,1); +for i=1:nods + %1: find the elements that contain the grid i + neighbors_gridi=[]; + for j=1:3 + neighbors_gridi=[neighbors_gridi find(index(:,j)==i)']; + end + numberofneighbors_gridi=length(neighbors_gridi); + %2 retrieve the value of p and q over each of these elements. The average is + %plugged into dbx_grid + qcoeff_grid(i)=sum(qcoeff(neighbors_gridi))/numberofneighbors_gridi; + pcoeff_grid(i)=sum(pcoeff(neighbors_gridi))/numberofneighbors_gridi; +end + +%Build length_icefront and normal_icefront: +[length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront); + +%Start building areas +aire=zeros(md.numberofelements,1); + +for n=1:nel + aire(n)=1/2 * det([1 1 1;x(index(n,:))';y(index(n,:))']); +end + +aire=abs(aire); % if index is sorted from its original value, then aire could be negative + +alpha=zeros(nel,3); +beta=zeros(nel,3); +gamma=zeros(nel,3); + +for n=1:nel + X=inv([x(index(n,:)) y(index(n,:)) ones(3,1)]); + alpha(n,:)=X(1,:); + beta(n,:)=X(2,:); + gamma(n,:)=X(3,:); +end + +clear X + +%Initialize viscosity +nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff); + + + % Do once and for all the initial computation of matrix-locations: + +row_location=zeros(nel*3*3,1); +col_location=zeros(nel*3*3,1); +row_location_AD=zeros(nel*6,1); +col_location_AD=zeros(nel*6,1); + +count=-nel+1; + +for i=1:3 + for j=1:3 + count=count+nel; + row_location(count:count+nel-1)=index(:,i); + col_location(count:count+nel-1)=index(:,j); + end +end + + +count=-nel+1; + +for i=1:3 + for j=i:3 + count=count+nel; + row_location_AD(count:count+nel-1)=index(:,i); + col_location_AD(count:count+nel-1)=index(:,j); + end +end + + + +permanent_pieces_of_A=zeros(nel*6,1); +permanent_pieces_of_B=zeros(nel*3*3,1); +permanent_pieces_of_C=zeros(nel*3*3,1); +permanent_pieces_of_D=zeros(nel*6,1); + +count=-nel+1; + +for i=1:3 + for j=i:3 + count=count+nel; + permanent_pieces_of_A(count:count+nel-1)= z_thick_bar .* aire ... + .*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j)); + % This loop structure works only when index is sorted + permanent_pieces_of_D(count:count+nel-1)= z_thick_bar .* aire ... + .*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j)); + + end +end + +count=-nel+1; + +for i=1:3 + for j=1:3 + count=count+nel; + + permanent_pieces_of_B(count:count+nel-1)= z_thick_bar .* aire ... + .*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i)); + + permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ... + .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)); + + + end +end + +% Step 3 -- Set up right-hand side of the problem. + +% (Note to myself: to avoid vector dependency problem, yet still vectorize, +% I treat the Rhs vector as a sparse matrix + +Rhs_x=zeros(nel*27,1); +Rhs_y=zeros(nel*27,1); +Rhs_y=ones(nel*27,1); +row_rhs=zeros(nel*27,1); + + +count=-nel+1; +for i=1:3 + for n=1:3 + for m=1:3 + count=count+nel; + + Rhs_x(count:count+nel-1)= -rho_ice * g * ... + z_thick(index(:,n)).*z_surf(index(:,m)) ... + .* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 ); + + Rhs_y(count:count+nel-1)= -rho_ice * g * ... + z_thick(index(:,n)).*z_surf(index(:,m)) ... + .* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 ); + + row_rhs(count:count+nel-1)=index(:,i); + + end + end +end + +Rhs=full([sparse(row_rhs,ones(nel*27,1),Rhs_x,nods,1) + sparse(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]); + +for k=1:2 + for l=1:2 + for j=1:2 + Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ... +(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... ++ rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... +.*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... ++min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... +.*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... +.*normal_icefront(:,1).*length_icefront(:) ... +/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); + + Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ... +(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... ++ rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... +.*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... ++min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... +.*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... +.*normal_icefront(:,2).*length_icefront(:) ... +/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); + + end + end +end + +clear Rhs_x Rhs_y row_rhs + +% Step 4 -- Create the parsing matrix to "wring out" +% the known boundary conditions + +% kinematic condition (u,v=0) specified at non-ice-front boundaries; +num_specified= 2*sum(node_on_dirichlet); +row=zeros(2*nods-num_specified,1); +col=zeros(2*nods-num_specified,1); +value=zeros(2*nods-num_specified,1); +count=0; +for n=1:nods + if(~node_on_dirichlet(n)) + count=count+1; + row(count)=count; + col(count)=n; + value(count)=1; + end +end +for n=1:nods + if (~node_on_dirichlet(n)) + count=count+1; + row(count)=count; + col(count)=nods+n; + value(count)=1; + end +end + +P=sparse(row,col,value,2*nods-num_specified,2*nods); + + +specified_velocity=zeros(2*nods,1); +pos=find(node_on_dirichlet); +specified_velocity(pos)=md.dirichletvalues_diag(pos,1)/md.yts; +specified_velocity(nods+pos)=md.dirichletvalues_diag(pos,2)/md.yts; + + +% ######## an attention grabbing break ######## +% Iterate on flow law till converged + +converged_yet=0; +loop=0; + +u_old=zeros(nods,1); +v_old=zeros(nods,1); +convergence_count=1; + +while (~converged_yet) + + if (loop>(100)) + warning('Maximum Viscosity Iterations Reached.') + break; + end; + loop=loop+1; + + nu_bar_oldvalue=nu_bar; + + % Step 4 -- Set up stress-balance matrix (whos solution is the velocity + % field): + + matrix_value_A=zeros(nel*6,1); + matrix_value_B=zeros(nel*3*3,1); + matrix_value_C=zeros(nel*3*3,1); + matrix_value_D=zeros(nel*6,1); + + count=-nel+1; + + for i=1:3 + for j=i:3 + count=count+nel; + + matrix_value_A(count:count+nel-1)=nu_bar .* ... + permanent_pieces_of_A(count:count+nel-1); + + matrix_value_D(count:count+nel-1)=nu_bar .* ... + permanent_pieces_of_D(count:count+nel-1); + + + end + end + + count=-nel+1; + + for i=1:3 + for j=1:3 + count=count+nel; + + matrix_value_B(count:count+nel-1)=nu_bar .* ... + permanent_pieces_of_B(count:count+nel-1); + + matrix_value_C(count:count+nel-1)=nu_bar .* ... + permanent_pieces_of_C(count:count+nel-1); + + + end + end + + + + A=sparse(row_location_AD,col_location_AD,matrix_value_A,nods,nods); + A=A+triu(A,1)'; + B=sparse(row_location,col_location,matrix_value_B,nods,nods); + C=sparse(row_location,col_location,matrix_value_C,nods,nods); + D=sparse(row_location_AD,col_location_AD,matrix_value_D,nods,nods); + D=D+triu(D,1)'; + + + F=[A C + B D]; + + %Now, take care of the basal friction if there is any: to make things easier, we translate u = k *Neff^(-q)* sigma^p into + % sigma= drag^2 * Neff ^(r) * u^s with r=q/p and s=1/p : */ + + if (drag_type==2), + %compute coeffs: + rcoeff=qcoeff_grid./pcoeff_grid; + scoeff=1./pcoeff_grid; + + %initialization of basal drag stiffness + Dragoperator=spalloc(2*nods,2*nods,0); + + if loop~=1, + + %retrieve the velocity magnitude + velocity_mag=sqrt(solution(1:nods).^2+solution(nods+1:2*nods).^2); + + %Computation of the effective pressure + Neff=g*(rho_ice*z_thick+rho_water*z_bed); + + %If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because + %the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour + %for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should + %replace it by Neff=0 (ie, equival it to an ice shelf)*/ + pos=find(Neff<0); + Neff(pos)=0; + + %Basal drag coefficient: Tau_x=-alpha^2 u, Tau_y=-alpha^2 v (See + %MacAyeal) + alpha2=(drag.^2).*(Neff.^rcoeff).*(velocity_mag.^(scoeff-1)); + + %stiffness due to basal drag + %initialization + count=0; + value=zeros(nel,27); + row=zeros(nel,27); + col=zeros(nel,27); + + for m=1:3 + for k=1:3 + for l=1:3 + if ( (m==k) + (m==l) + (l==k) )==3 + fac=1/10; + elseif ( (m==k) + (m==l) + (l==k) )==1 + fac=1/30; + else + fac=1/60; + end + + count=count+1; + row(:,count)=index(:,m); + col(:,count)=index(:,k); + value(:,count)=fac*aire(:).* ... + (alpha2(index(:,l)));%.*element_on_icesheet(index(:,l)); + + end + end + end + + Dragoperator=[sparse(row,col,value,nods,nods) spalloc(nods,nods,0) + spalloc(nods,nods,0) sparse(row,col,value,nods,nods)]; + end + + F=F+Dragoperator; %plug into the global stiffness matrix + + end + + + Rhs_parsed=P*(Rhs - F*specified_velocity); + F=P*F*P'; + + + % Step 5 -- Solve the problem: + + % Digression: if need be clear up some memory: + + clear matrix_value_A ... + matrix_value_B matrix_value_C matrix_value_D A B C D + + % We can use either the LU or the Cholesky decomposition, but the + % Cholesky decomposition is twice as efficient as LU for symmetric + % definite positive matrix + + if strcmpi(md.solver_type,'lu'), + % Solve by LU decomposition. + [L,U] = lu(F); + solution= U\(L\Rhs_parsed); + elseif strcmpi(md.solver_type,'cholesky'), + % Solve by Choleski decomposition. + L = chol(F); + solution= L\(L'\Rhs_parsed); + else + % use matlab's generic solver + solution = F\Rhs_parsed; + end + + + + %Add spcs to the calculated solution + solution=P'*solution + specified_velocity; + + %Recover solution vector + u=solution(1:nods); + v=solution(nods+1:2*nods); + + %Compute viscosity + nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); + change=1 - nu_bar_oldvalue./nu_bar; + + %Update viscosity + nu_bar=nu_bar.*(1+viscosity_overshoot*change); + location=find(nu_bar<=0); + nu_bar(location)=-nu_bar(location); + + %Test for direct shooting convergence + if convergence_count>1, + + ug=[u_old;v_old]; + nug=norm(ug,2); + dug=[u; v]-[u_old; v_old]; + ndug=norm(dug,2); + relative_change=ndug/nug; + + + %Figure out if viscosity converged + if relative_change ',criterion_rel,' m/yr')); + converged_yet=0; + end + + if ~isnan(criterion_abs), + change=max(abs(dug))*yts; + if change ',criterion_abs,' m/yr')); + converged_yet=0; + end + end + else + converged_yet=0; + end + + u_old=u; + v_old=v; + + convergence_count=convergence_count+1; + +end % This end statement terminates the "while" command way above + +%Load results onto md: +md.vx=u*yts; +md.vy=v*yts; +md.vel=sqrt(md.vx.^2+md.vy.^2); + +%stop timing +t2=clock; + +disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); + +function macayealdiagnosticusage(); +disp('md=macayealdiagnostic(md)'); +disp(' where md is a structure of class @model'); Index: sm/trunk/src/m/solutions/macayeal/macayealcontrol.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/macayealcontrol.m (revision 34) +++ (revision ) @@ -1,353 +1,0 @@ -function md=macayealcontrol(md) -%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) - -%First check we do have the correct argument number -if ((nargin~=1) || (nargout~=1)), - velfinderusage(); - error('macayealcontrol error message: incorrect number of input and output arguments'); -end - -if ~isa(md,'model'), - macayealcontrolusage(); - error('macayealcontrol error message: input argument is not a @model object'); -end - -%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 - -%Transfer model fields into matlab variables -x=md.x; -y=md.y; -index=md.elements; -index=sort(index,2); %necessary -nods=md.numberofgrids; -nel=md.numberofelements; -z_thick=md.thickness; -z_surf=md.surface; -z_bed=md.bed; -z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3; -rho_ice=md.rho_ice; -rho_water=md.rho_water; -g=md.g; -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 -nodes_on_boundary=md.gridonboundary; -nodes_on_dirichlet=md.gridondirichlet_diag; -nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1; -nodes_on_iceshelf=md.gridoniceshelf; - -criterion=md.eps_rel; -yts=md.yts; -tolx=md.tolx; -maxiter=md.maxiter(1); -B_ini=md.B; -drag_coeff_ini=md.drag; -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: -[length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront); - -%Building shape functions and derivative operators -[alpha, beta, gamma, area]=shape(index,x,y,nel,nods); - -[matrix_bar, matrix_xbar, matrix_ybar]=... - bar_maker(nel,nods,index,alpha,beta); - -%initialize some data -create_el2nod_matrices -old_direction=zeros(nods,1); - - -%setup some fake distributions. -z_thick_bar=matrix_bar*z_thick; -z_surf_bar=matrix_bar*z_surf; -z_surf_xbar=matrix_xbar*z_surf; -z_surf_ybar=matrix_ybar*z_surf; -z_surf_x=el2nod\(el2nodRhs*z_surf_xbar); -z_surf_y=el2nod\(el2nodRhs*z_surf_ybar); -z_bed_bar=matrix_bar*z_bed; -z_bed_xbar=matrix_xbar*z_bed; -z_bed_ybar=matrix_ybar*z_bed; - - -%data_elements and data_nodes are the elements and nodes where -%we have observations -data_elements=[1:nel]'; -data_nodes=ones(nods,1); - -%initialize misfit between model velocity and observed velocity -J=zeros(nsteps,1); - -%Weighting: one areas where we don't want the control method to optimize -%the misfit, we can set weighting to 0. -weighting=ones(nods,1); - -%optimization parameters for the matlab functoin fminbnd -options=optimset('fminbnd'); -options=optimset(options,'Display','iter'); -options=optimset(options,'MaxFunEvals',maxiter); -options=optimset(options,'MaxIter',100); -options=optimset(options,'TolX',tolx); - -%build useful matrices. -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; - -%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), - disp(sprintf('\n No drag pure for ice shelves => STOP')) - return -end - -%nu_bar=10^14*ones(nel,1); %initial element viscosity distribution. -nu_bar=viscosity(index,nel,alpha,beta,vx_obs,vy_obs,B_bar,glen_coeff); - -%%%AK velfinder; %forward model that determines first velocity field used -%to start optimization. -disp('calculating the velocity with the initial parameters') -c_velfinder; - -if strcmp(md.control_type{1},'B'), - for niteration=1:nsteps, - - disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]); - - % Compute search direction - - [u,v,adjointu,adjointv,direction]= ... - grad_J_flow(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vy_obs,... - index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,B_bar); - - if md.plot==1, - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',sqrt(u.^2+v.^2)*yts,'title','Modeled velocity',... - 'data',sqrt(adjointu.^2+adjointv.^2)*yts,'title','Adjoint vectors','colorbar#all', 'on',... - 'data',sqrt(vx_obs.^2+vy_obs.^2)*yts,'title','Observed velocity',... - '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);drawnow; - end - - %Keep track of u and v to use in objectivefunction_C: - u_objective=u; - v_objective=v; - - %Orthogonalize direction - direction=real(direction); - direction=direction/sqrt(direction'*direction); - - - % 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^7, so that when variations on B are computed - %they will be significant. - if abs(max(direction))>abs(min(direction)) - direction=10^7*direction/abs(max(direction)); - else - direction=10^7*direction/abs(min(direction)); - end - - %during optimization, bounds on B variations can vary. Here, they are less - %strict in the first iteration. - if niteration<=2, - upperbound=20; - lowerbound=0; - else - upperbound=10; - lowerbound=0; - end - - %search the multiplicative constant to the direction, that will - %be used to modify B. - - search_constant=fminbnd('objectivefunction_C_flow',lowerbound,upperbound, ... - options,B,glen_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 B - B_old=B; - B_new=B+search_constant*direction; - B=B_new; - pos=find(B<0);B(pos)=-B(pos); - - %Average B over elements: - B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3; - - %visualize new distribution. - if md.plot==1, - 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 - plotmodel(md,'data',B*B_to_p,'title',['B at iteration ' num2str(niteration)],'caxis',[200 900],'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 - - J(niteration)=objectivefunction_C_flow(0,B,glen_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); - disp(J(niteration)); - - %do a backup every 5 iterations - if(mod(niteration,5)==0), - save temporary_control_results_flow B B_bar u v J direction - end - end - - %Load results onto md: - md.cont_J=J; - md.cont_parameters=B; - -elseif strcmp(md.control_type,'basal'), - for niteration=1:nsteps, - - disp([' Step #' num2str(niteration) ' on ' num2str(nsteps)]); - - % Compute search direction - - [u,v,adjointu,adjointv,direction]= ... - grad_J_drag(Rhs,S,F,P,P0,area,specified_velocity,nods,vx_obs,vx_obs,... - index,x,y,nel,rho_ice,g,weighting,alpha,beta,z_thick_bar,drag_coeff,drag_coeff_bar); - - if md.plot==1, - if niteration==1 - scrsz = get(0,'ScreenSize'); - figure('Position',[scrsz(3)/3 scrsz(4)*1.8/3 scrsz(3)/3 scrsz(4)/3]) - end - plotmodel(md,'data',sqrt(u.^2+v.^2),'title','Modeled velocity',... - 'data',sqrt(adjointu.^2+adjointv.^2),'title','Adjoint vectors','colorbar', 'all',... - 'data',sqrt(vx_obs.^2+vy_obs.^2),'title','Observed velocity',... - '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);drawnow; - end - - %Keep track of u and v to use in objectivefunction_C: - u_objective=u; - 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 - %here to win some computation time. - load F_file - - J(niteration)=objectivefunction_C_drag(0,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); - disp(J(niteration)); - - %do a backup every 5 iterations - if(mod(niteration,5)==0), - save temporary_control_results_drag drag_coeff drag_coeff_bar u v J direction - end - end -end - -function macayealcontrolusage(); -disp('md=macayealcontrol(md)'); -disp(' where md is a structure of class @model'); Index: sm/trunk/src/m/solutions/macayeal/macayealdiagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/macayealdiagnostic.m (revision 34) +++ (revision ) @@ -1,500 +1,0 @@ -function md=macayealdiagnostic(md) -%MACAYEALDIAGNOSTIC - compute the velocity field of a 2d model using MacAyeal solution -% -% this routine solves the problem using MacAyeal's model. It calculates the velocity -% field corresponding to the parameters and the geometry given by the model md -% -% Usage: -% md=macayealdiagnostic(md) - -%First check we do have the correct argument number -if ((nargin~=1) || (nargout~=1)), - macayealdiagnosticusage(); - error('macayealdiagnostic error message: incorrect number of input and output arguments'); -end - -if ~isa(md,'model'), - macayealdiagnosticusage(); - error('macayealdiagnostic error message: input argument is not a @model object'); -end - -%start timing -t1=clock; - -%Transfer model fields into matlab variables -x=md.x; -y=md.y; - -index=md.elements;index=sort(index,2); %necessary -nods=md.numberofgrids; -nel=md.numberofelements; -z_thick=md.thickness; -z_surf=md.surface; -z_bed=md.bed; -z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3; -rho_ice=md.rho_ice; -rho_water=md.rho_water; -g=md.g; -viscosity_overshoot=md.viscosity_overshoot; -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 -nodes_on_boundary=md.gridonboundary; -nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1; -node_on_dirichlet=md.gridondirichlet_diag; -nodes_on_icesheet=md.gridonicesheet; -element_on_icesheet=md.elementonicesheet; -qcoeff=md.q; -pcoeff=md.p; -drag_type=md.drag_type; -drag=md.drag; - -criterion_rel=md.eps_rel; -criterion_abs=md.eps_abs; -yts=md.yts; -B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3; -glen_coeff=md.n; - -%average of p and q over the grids (size nel->nods) -pcoeff_grid=zeros(nods,1); -qcoeff_grid=zeros(nods,1); -for i=1:nods - %1: find the elements that contain the grid i - neighbors_gridi=[]; - for j=1:3 - neighbors_gridi=[neighbors_gridi find(index(:,j)==i)']; - end - numberofneighbors_gridi=length(neighbors_gridi); - %2 retrieve the value of p and q over each of these elements. The average is - %plugged into dbx_grid - qcoeff_grid(i)=sum(qcoeff(neighbors_gridi))/numberofneighbors_gridi; - pcoeff_grid(i)=sum(pcoeff(neighbors_gridi))/numberofneighbors_gridi; -end - -%Build length_icefront and normal_icefront: -[length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront); - -%Start building areas -aire=zeros(md.numberofelements,1); - -for n=1:nel - aire(n)=1/2 * det([1 1 1;x(index(n,:))';y(index(n,:))']); -end - -aire=abs(aire); % if index is sorted from its original value, then aire could be negative - -alpha=zeros(nel,3); -beta=zeros(nel,3); -gamma=zeros(nel,3); - -for n=1:nel - X=inv([x(index(n,:)) y(index(n,:)) ones(3,1)]); - alpha(n,:)=X(1,:); - beta(n,:)=X(2,:); - gamma(n,:)=X(3,:); -end - -clear X - -%Initialize viscosity -nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff); - - - % Do once and for all the initial computation of matrix-locations: - -row_location=zeros(nel*3*3,1); -col_location=zeros(nel*3*3,1); -row_location_AD=zeros(nel*6,1); -col_location_AD=zeros(nel*6,1); - -count=-nel+1; - -for i=1:3 - for j=1:3 - count=count+nel; - row_location(count:count+nel-1)=index(:,i); - col_location(count:count+nel-1)=index(:,j); - end -end - - -count=-nel+1; - -for i=1:3 - for j=i:3 - count=count+nel; - row_location_AD(count:count+nel-1)=index(:,i); - col_location_AD(count:count+nel-1)=index(:,j); - end -end - - - -permanent_pieces_of_A=zeros(nel*6,1); -permanent_pieces_of_B=zeros(nel*3*3,1); -permanent_pieces_of_C=zeros(nel*3*3,1); -permanent_pieces_of_D=zeros(nel*6,1); - -count=-nel+1; - -for i=1:3 - for j=i:3 - count=count+nel; - permanent_pieces_of_A(count:count+nel-1)= z_thick_bar .* aire ... - .*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j)); - % This loop structure works only when index is sorted - permanent_pieces_of_D(count:count+nel-1)= z_thick_bar .* aire ... - .*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j)); - - end -end - -count=-nel+1; - -for i=1:3 - for j=1:3 - count=count+nel; - - permanent_pieces_of_B(count:count+nel-1)= z_thick_bar .* aire ... - .*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i)); - - permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ... - .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)); - - - end -end - -% Step 3 -- Set up right-hand side of the problem. - -% (Note to myself: to avoid vector dependency problem, yet still vectorize, -% I treat the Rhs vector as a sparse2 matrix - -Rhs_x=zeros(nel*27,1); -Rhs_y=zeros(nel*27,1); -Rhs_y=ones(nel*27,1); -row_rhs=zeros(nel*27,1); - - -count=-nel+1; -for i=1:3 - for n=1:3 - for m=1:3 - count=count+nel; - - Rhs_x(count:count+nel-1)= -rho_ice * g * ... - z_thick(index(:,n)).*z_surf(index(:,m)) ... - .* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 ); - - Rhs_y(count:count+nel-1)= -rho_ice * g * ... - z_thick(index(:,n)).*z_surf(index(:,m)) ... - .* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 ); - - row_rhs(count:count+nel-1)=index(:,i); - - end - end -end - -Rhs=full([sparse2(row_rhs,ones(nel*27,1),Rhs_x,nods,1) - sparse2(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]); - -for k=1:2 - for l=1:2 - for j=1:2 - Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ... -(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... -+ rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... -.*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... -+min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... -.*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... -.*normal_icefront(:,1).*length_icefront(:) ... -/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); - - Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ... -(rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... -+ rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... -.*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... -+min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... -.*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... -.*normal_icefront(:,2).*length_icefront(:) ... -/(4*(l==k & j==k) + 12*(l~=k | j~=k) ); - - end - end -end - -clear Rhs_x Rhs_y row_rhs - -% Step 4 -- Create the parsing matrix to "wring out" -% the known boundary conditions - -% kinematic condition (u,v=0) specified at non-ice-front boundaries; -num_specified= 2*sum(node_on_dirichlet); -row=zeros(2*nods-num_specified,1); -col=zeros(2*nods-num_specified,1); -value=zeros(2*nods-num_specified,1); -count=0; -for n=1:nods - if(~node_on_dirichlet(n)) - count=count+1; - row(count)=count; - col(count)=n; - value(count)=1; - end -end -for n=1:nods - if (~node_on_dirichlet(n)) - count=count+1; - row(count)=count; - col(count)=nods+n; - value(count)=1; - end -end - -P=sparse2(row,col,value,2*nods-num_specified,2*nods); - - -specified_velocity=zeros(2*nods,1); -pos=find(node_on_dirichlet); -specified_velocity(pos)=md.dirichletvalues_diag(pos,1)/md.yts; -specified_velocity(nods+pos)=md.dirichletvalues_diag(pos,2)/md.yts; - - -% ######## an attention grabbing break ######## -% Iterate on flow law till converged - -converged_yet=0; -loop=0; - -u_old=zeros(nods,1); -v_old=zeros(nods,1); -convergence_count=1; - -while (~converged_yet) - - if (loop>(100)) - warning('Maximum Viscosity Iterations Reached.') - break; - end; - loop=loop+1; - - nu_bar_oldvalue=nu_bar; - - % Step 4 -- Set up stress-balance matrix (whos solution is the velocity - % field): - - matrix_value_A=zeros(nel*6,1); - matrix_value_B=zeros(nel*3*3,1); - matrix_value_C=zeros(nel*3*3,1); - matrix_value_D=zeros(nel*6,1); - - count=-nel+1; - - for i=1:3 - for j=i:3 - count=count+nel; - - matrix_value_A(count:count+nel-1)=nu_bar .* ... - permanent_pieces_of_A(count:count+nel-1); - - matrix_value_D(count:count+nel-1)=nu_bar .* ... - permanent_pieces_of_D(count:count+nel-1); - - - end - end - - count=-nel+1; - - for i=1:3 - for j=1:3 - count=count+nel; - - matrix_value_B(count:count+nel-1)=nu_bar .* ... - permanent_pieces_of_B(count:count+nel-1); - - matrix_value_C(count:count+nel-1)=nu_bar .* ... - permanent_pieces_of_C(count:count+nel-1); - - - end - end - - - - A=sparse2(row_location_AD,col_location_AD,matrix_value_A,nods,nods); - A=A+triu(A,1)'; - B=sparse2(row_location,col_location,matrix_value_B,nods,nods); - C=sparse2(row_location,col_location,matrix_value_C,nods,nods); - D=sparse2(row_location_AD,col_location_AD,matrix_value_D,nods,nods); - D=D+triu(D,1)'; - - - F=[A C - B D]; - - %Now, take care of the basal friction if there is any: to make things easier, we translate u = k *Neff^(-q)* sigma^p into - % sigma= drag^2 * Neff ^(r) * u^s with r=q/p and s=1/p : */ - - if (drag_type==2), - %compute coeffs: - rcoeff=qcoeff_grid./pcoeff_grid; - scoeff=1./pcoeff_grid; - - %initialization of basal drag stiffness - Dragoperator=spalloc(2*nods,2*nods,0); - - if loop~=1, - - %retrieve the velocity magnitude - velocity_mag=sqrt(solution(1:nods).^2+solution(nods+1:2*nods).^2); - - %Computation of the effective pressure - Neff=g*(rho_ice*z_thick+rho_water*z_bed); - - %If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because - %the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour - %for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should - %replace it by Neff=0 (ie, equival it to an ice shelf)*/ - pos=find(Neff<0); - Neff(pos)=0; - - %Basal drag coefficient: Tau_x=-alpha^2 u, Tau_y=-alpha^2 v (See - %MacAyeal) - alpha2=(drag.^2).*(Neff.^rcoeff).*(velocity_mag.^(scoeff-1)); - - %stiffness due to basal drag - %initialization - count=0; - value=zeros(nel,27); - row=zeros(nel,27); - col=zeros(nel,27); - - for m=1:3 - for k=1:3 - for l=1:3 - if ( (m==k) + (m==l) + (l==k) )==3 - fac=1/10; - elseif ( (m==k) + (m==l) + (l==k) )==1 - fac=1/30; - else - fac=1/60; - end - - count=count+1; - row(:,count)=index(:,m); - col(:,count)=index(:,k); - value(:,count)=fac*aire(:).* ... - (alpha2(index(:,l)));%.*element_on_icesheet(index(:,l)); - - end - end - end - - Dragoperator=[sparse2(row,col,value,nods,nods) spalloc(nods,nods,0) - spalloc(nods,nods,0) sparse2(row,col,value,nods,nods)]; - end - - F=F+Dragoperator; %plug into the global stiffness matrix - - end - - - Rhs_parsed=P*(Rhs - F*specified_velocity); - F=P*F*P'; - - - % Step 5 -- Solve the problem: - - % Digression: if need be clear up some memory: - - clear matrix_value_A ... - matrix_value_B matrix_value_C matrix_value_D A B C D - - % We can use either the LU or the Cholesky decomposition, but the - % Cholesky decomposition is twice as efficient as LU for symmetric - % definite positive matrix - - if strcmpi(md.solver_type,'lu'), - % Solve by LU decomposition. - [L,U] = lu(F); - solution= U\(L\Rhs_parsed); - elseif strcmpi(md.solver_type,'cholesky'), - % Solve by Choleski decomposition. - L = chol(F); - solution= L\(L'\Rhs_parsed); - else - % use matlab's generic solver - solution = F\Rhs_parsed; - end - - - - %Add spcs to the calculated solution - solution=P'*solution + specified_velocity; - - %Recover solution vector - u=solution(1:nods); - v=solution(nods+1:2*nods); - - %Compute viscosity - nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); - change=1 - nu_bar_oldvalue./nu_bar; - - %Update viscosity - nu_bar=nu_bar.*(1+viscosity_overshoot*change); - location=find(nu_bar<=0); - nu_bar(location)=-nu_bar(location); - - %Test for direct shooting convergence - if convergence_count>1, - - ug=[u_old;v_old]; - nug=norm(ug,2); - dug=[u; v]-[u_old; v_old]; - ndug=norm(dug,2); - relative_change=ndug/nug; - - - %Figure out if viscosity converged - if relative_change ',criterion_rel,' m/yr')); - converged_yet=0; - end - - if ~isnan(criterion_abs), - change=max(abs(dug))*yts; - if change ',criterion_abs,' m/yr')); - converged_yet=0; - end - end - else - converged_yet=0; - end - - u_old=u; - v_old=v; - - convergence_count=convergence_count+1; - -end % This end statement terminates the "while" command way above - -%Load results onto md: -md.vx=u*yts; -md.vy=v*yts; -md.vel=sqrt(md.vx.^2+md.vy.^2); - -%stop timing -t2=clock; - -disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); - -function macayealdiagnosticusage(); -disp('md=macayealdiagnostic(md)'); -disp(' where md is a structure of class @model');