macayealdiagnostic

PURPOSE ^

MACAYEALDIAGNOSTIC - compute the velocity field of a 2d model using MacAyeal solution

SYNOPSIS ^

function md=macayealdiagnostic(md)

DESCRIPTION ^

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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function md=macayealdiagnostic(md)
0002 %MACAYEALDIAGNOSTIC - compute the velocity field of a 2d model using MacAyeal solution
0003 %
0004 %   this routine solves the problem using MacAyeal's model. It calculates the velocity
0005 %   field corresponding to the parameters and the geometry given by the model md
0006 %
0007 %   Usage:
0008 %      md=macayealdiagnostic(md)
0009 
0010 %First check we do have the correct argument number
0011 if ((nargin~=1) || (nargout~=1)),
0012     macayealdiagnosticusage();
0013     error('macayealdiagnostic error message: incorrect number of input and output arguments');
0014 end
0015 
0016 if ~isa(md,'model'),
0017     macayealdiagnosticusage();
0018     error('macayealdiagnostic error message: input argument is not a @model object');
0019 end
0020 
0021 %start timing
0022 t1=clock;
0023 
0024 %Transfer model fields into matlab variables
0025 x=md.x;
0026 y=md.y;
0027 
0028 index=md.elements;index=sort(index,2); %necessary
0029 nods=md.numberofgrids;
0030 nel=md.numberofelements;
0031 z_thick=md.thickness;
0032 z_surf=md.surface;
0033 z_bed=md.bed;
0034 z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3;
0035 rho_ice=md.rho_ice;
0036 rho_water=md.rho_water;
0037 g=md.g;
0038 viscosity_overshoot=md.viscosity_overshoot;
0039 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
0040 nodes_on_boundary=md.gridonboundary;
0041 nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1;
0042 node_on_dirichlet=md.gridondirichlet_diag;
0043 nodes_on_icesheet=md.gridonicesheet;
0044 element_on_icesheet=md.elementonicesheet;
0045 qcoeff=md.q;
0046 pcoeff=md.p;
0047 drag_type=md.drag_type;
0048 drag=md.drag;
0049 
0050 criterion_rel=md.eps_rel;
0051 criterion_abs=md.eps_abs;
0052 yts=md.yts;
0053 B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
0054 glen_coeff=md.n;
0055 
0056 %average of p and q over the grids (size nel->nods)
0057 pcoeff_grid=zeros(nods,1);
0058 qcoeff_grid=zeros(nods,1);
0059 for i=1:nods
0060     %1: find the elements that contain the grid i
0061     neighbors_gridi=[];
0062     for j=1:3
0063         neighbors_gridi=[neighbors_gridi find(index(:,j)==i)'];
0064     end
0065     numberofneighbors_gridi=length(neighbors_gridi);
0066     %2 retrieve the value of p and q over each of these elements. The average is
0067     %plugged into dbx_grid
0068     qcoeff_grid(i)=sum(qcoeff(neighbors_gridi))/numberofneighbors_gridi;
0069     pcoeff_grid(i)=sum(pcoeff(neighbors_gridi))/numberofneighbors_gridi;
0070 end
0071 
0072 %Build length_icefront and normal_icefront:
0073 [length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront);
0074 
0075 %Start building areas
0076 aire=zeros(md.numberofelements,1);
0077 
0078 for n=1:nel
0079     aire(n)=1/2 * det([1 1 1;x(index(n,:))';y(index(n,:))']);
0080 end
0081 
0082 aire=abs(aire); % if index is sorted from its original value, then aire could be negative
0083 
0084 alpha=zeros(nel,3);
0085 beta=zeros(nel,3);
0086 gamma=zeros(nel,3);
0087 
0088 for n=1:nel
0089     X=inv([x(index(n,:)) y(index(n,:)) ones(3,1)]);
0090     alpha(n,:)=X(1,:);
0091     beta(n,:)=X(2,:);
0092     gamma(n,:)=X(3,:);
0093 end
0094 
0095 clear X
0096 
0097 %Initialize viscosity
0098 nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
0099 
0100 
0101  %  Do once and for all the initial computation of matrix-locations:
0102 
0103 row_location=zeros(nel*3*3,1);
0104 col_location=zeros(nel*3*3,1);
0105 row_location_AD=zeros(nel*6,1);
0106 col_location_AD=zeros(nel*6,1);
0107 
0108 count=-nel+1;
0109 
0110 for i=1:3
0111     for j=1:3
0112             count=count+nel;
0113             row_location(count:count+nel-1)=index(:,i);
0114             col_location(count:count+nel-1)=index(:,j);
0115     end
0116 end
0117 
0118 
0119 count=-nel+1;
0120 
0121 for i=1:3
0122     for j=i:3
0123             count=count+nel;
0124             row_location_AD(count:count+nel-1)=index(:,i);
0125             col_location_AD(count:count+nel-1)=index(:,j);
0126     end
0127 end
0128 
0129 
0130 
0131 permanent_pieces_of_A=zeros(nel*6,1);
0132 permanent_pieces_of_B=zeros(nel*3*3,1);
0133 permanent_pieces_of_C=zeros(nel*3*3,1);
0134 permanent_pieces_of_D=zeros(nel*6,1);
0135 
0136 count=-nel+1;
0137 
0138 for i=1:3
0139     for j=i:3
0140             count=count+nel;
0141             permanent_pieces_of_A(count:count+nel-1)= z_thick_bar .* aire ...
0142                 .*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j));                
0143                 % This loop structure works only when index is sorted
0144             permanent_pieces_of_D(count:count+nel-1)= z_thick_bar .* aire ...
0145                 .*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j));
0146 
0147     end
0148 end
0149 
0150 count=-nel+1;
0151 
0152 for i=1:3
0153     for j=1:3
0154             count=count+nel;
0155 
0156             permanent_pieces_of_B(count:count+nel-1)= z_thick_bar .* aire ...
0157             .*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i));
0158 
0159             permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ...
0160             .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i));
0161 
0162 
0163     end
0164 end
0165 
0166 % Step 3 -- Set up right-hand side of the problem.
0167 
0168 % (Note to myself: to avoid vector dependency problem, yet still vectorize,
0169 % I treat the Rhs vector as a sparse2 matrix
0170 
0171 Rhs_x=zeros(nel*27,1);
0172 Rhs_y=zeros(nel*27,1);
0173 Rhs_y=ones(nel*27,1);
0174 row_rhs=zeros(nel*27,1);
0175 
0176 
0177 count=-nel+1;
0178 for i=1:3
0179     for n=1:3
0180         for m=1:3
0181             count=count+nel;
0182 
0183             Rhs_x(count:count+nel-1)= -rho_ice * g  * ...
0184              z_thick(index(:,n)).*z_surf(index(:,m)) ...
0185             .* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 );
0186 
0187             Rhs_y(count:count+nel-1)= -rho_ice * g  * ...
0188              z_thick(index(:,n)).*z_surf(index(:,m)) ...
0189             .* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 );
0190 
0191             row_rhs(count:count+nel-1)=index(:,i);
0192 
0193         end
0194     end
0195 end
0196 
0197 Rhs=full([sparse2(row_rhs,ones(nel*27,1),Rhs_x,nods,1)
0198    sparse2(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]);
0199 
0200 for k=1:2
0201     for l=1:2
0202         for j=1:2
0203                 Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ...
0204 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
0205 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ...
0206 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ...
0207 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))...
0208 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))...
0209 .*normal_icefront(:,1).*length_icefront(:) ...
0210 /(4*(l==k & j==k) + 12*(l~=k | j~=k) );
0211 
0212                 Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ...
0213 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...
0214 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ...
0215 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ...
0216 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))...
0217 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))...
0218 .*normal_icefront(:,2).*length_icefront(:) ...
0219 /(4*(l==k & j==k) + 12*(l~=k | j~=k) );
0220 
0221         end
0222     end
0223 end
0224 
0225 clear Rhs_x Rhs_y row_rhs
0226 
0227 % Step 4 -- Create the parsing matrix to "wring out"
0228 %           the known boundary conditions
0229 
0230 % kinematic condition (u,v=0) specified at non-ice-front boundaries;
0231 num_specified= 2*sum(node_on_dirichlet);
0232 row=zeros(2*nods-num_specified,1);
0233 col=zeros(2*nods-num_specified,1);
0234 value=zeros(2*nods-num_specified,1);
0235 count=0;
0236 for n=1:nods
0237    if(~node_on_dirichlet(n))
0238            count=count+1;
0239         row(count)=count;
0240         col(count)=n;
0241         value(count)=1;
0242     end
0243 end
0244 for n=1:nods
0245    if (~node_on_dirichlet(n))
0246          count=count+1;
0247         row(count)=count;
0248         col(count)=nods+n;
0249         value(count)=1;
0250     end
0251 end
0252 
0253 P=sparse2(row,col,value,2*nods-num_specified,2*nods);
0254 
0255 
0256 specified_velocity=zeros(2*nods,1);
0257 pos=find(node_on_dirichlet);
0258 specified_velocity(pos)=md.dirichletvalues_diag(pos,1)/md.yts;
0259 specified_velocity(nods+pos)=md.dirichletvalues_diag(pos,2)/md.yts;
0260 
0261 
0262 % ######## an attention grabbing break ########
0263 %      Iterate on flow law till converged
0264 
0265 converged_yet=0;
0266 loop=0;
0267 
0268 u_old=zeros(nods,1);
0269 v_old=zeros(nods,1);
0270 convergence_count=1;
0271 
0272 while (~converged_yet)
0273 
0274     if (loop>(100))
0275         warning('Maximum Viscosity Iterations  Reached.')
0276         break;
0277     end;
0278     loop=loop+1;
0279 
0280     nu_bar_oldvalue=nu_bar;
0281 
0282     % Step 4 -- Set up stress-balance matrix  (whos solution is the velocity
0283     % field):
0284 
0285     matrix_value_A=zeros(nel*6,1);
0286     matrix_value_B=zeros(nel*3*3,1);
0287     matrix_value_C=zeros(nel*3*3,1);
0288     matrix_value_D=zeros(nel*6,1);
0289 
0290     count=-nel+1;
0291 
0292     for i=1:3
0293         for j=i:3
0294                 count=count+nel;
0295 
0296                 matrix_value_A(count:count+nel-1)=nu_bar .* ...
0297                   permanent_pieces_of_A(count:count+nel-1);
0298 
0299                 matrix_value_D(count:count+nel-1)=nu_bar .* ...
0300                   permanent_pieces_of_D(count:count+nel-1);
0301 
0302 
0303         end
0304     end
0305 
0306     count=-nel+1;
0307 
0308     for i=1:3
0309         for j=1:3
0310                 count=count+nel;
0311 
0312                 matrix_value_B(count:count+nel-1)=nu_bar .* ...
0313                   permanent_pieces_of_B(count:count+nel-1);
0314 
0315                 matrix_value_C(count:count+nel-1)=nu_bar .* ...
0316                   permanent_pieces_of_C(count:count+nel-1);
0317 
0318 
0319         end
0320     end
0321 
0322 
0323 
0324     A=sparse2(row_location_AD,col_location_AD,matrix_value_A,nods,nods);
0325     A=A+triu(A,1)';
0326     B=sparse2(row_location,col_location,matrix_value_B,nods,nods);
0327     C=sparse2(row_location,col_location,matrix_value_C,nods,nods);
0328     D=sparse2(row_location_AD,col_location_AD,matrix_value_D,nods,nods);
0329     D=D+triu(D,1)';
0330 
0331 
0332     F=[A C
0333     B D];
0334 
0335     %Now, take care of the basal friction if there is any: to make things easier, we translate u = k *Neff^(-q)* sigma^p into
0336     % sigma= drag^2 * Neff ^(r) * u^s with r=q/p and s=1/p : */
0337      
0338     if  (drag_type==2), 
0339         %compute coeffs:
0340         rcoeff=qcoeff_grid./pcoeff_grid;
0341         scoeff=1./pcoeff_grid;
0342         
0343         %initialization of basal drag stiffness
0344         Dragoperator=spalloc(2*nods,2*nods,0);
0345 
0346         if loop~=1,
0347             
0348             %retrieve the velocity magnitude
0349             velocity_mag=sqrt(solution(1:nods).^2+solution(nods+1:2*nods).^2);
0350         
0351             %Computation of the effective pressure
0352             Neff=g*(rho_ice*z_thick+rho_water*z_bed);
0353 
0354             %If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
0355             %the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
0356             %for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
0357             %replace it by Neff=0 (ie, equival it to an ice shelf)*/
0358             pos=find(Neff<0);
0359             Neff(pos)=0;
0360             
0361             %Basal drag coefficient: Tau_x=-alpha^2 u, Tau_y=-alpha^2 v (See
0362             %MacAyeal)
0363             alpha2=(drag.^2).*(Neff.^rcoeff).*(velocity_mag.^(scoeff-1));
0364 
0365             %stiffness due to basal drag
0366             %initialization
0367             count=0;
0368             value=zeros(nel,27);
0369             row=zeros(nel,27);
0370             col=zeros(nel,27);
0371 
0372             for m=1:3
0373                 for k=1:3
0374                 for l=1:3
0375                     if ( (m==k) + (m==l) + (l==k) )==3
0376                     fac=1/10;
0377                     elseif ( (m==k) + (m==l) + (l==k) )==1
0378                     fac=1/30;
0379                     else
0380                     fac=1/60;
0381                     end
0382 
0383                     count=count+1;
0384                     row(:,count)=index(:,m);
0385                     col(:,count)=index(:,k);
0386                     value(:,count)=fac*aire(:).* ...
0387                     (alpha2(index(:,l)));%.*element_on_icesheet(index(:,l));
0388 
0389                 end
0390                 end
0391             end
0392 
0393             Dragoperator=[sparse2(row,col,value,nods,nods) spalloc(nods,nods,0)
0394                 spalloc(nods,nods,0) sparse2(row,col,value,nods,nods)];
0395             end
0396 
0397         F=F+Dragoperator;  %plug into the global stiffness matrix
0398         
0399     end
0400 
0401 
0402     Rhs_parsed=P*(Rhs - F*specified_velocity);
0403     F=P*F*P';
0404 
0405 
0406     % Step 5 -- Solve the problem:
0407 
0408     % Digression: if need be clear up some memory:
0409 
0410     clear   matrix_value_A ...
0411     matrix_value_B matrix_value_C matrix_value_D A B C D 
0412 
0413     % We can use either the LU or the Cholesky decomposition, but the
0414     % Cholesky decomposition is twice as efficient as LU for symmetric
0415     % definite positive matrix
0416 
0417         if strcmpi(md.solver_type,'lu'),
0418             % Solve by LU decomposition.
0419             [L,U] = lu(F);
0420             solution= U\(L\Rhs_parsed);
0421         elseif strcmpi(md.solver_type,'cholesky'),
0422             % Solve by Choleski decomposition.
0423             L = chol(F);
0424             solution= L\(L'\Rhs_parsed);
0425         else
0426             % use matlab's generic solver
0427             solution = F\Rhs_parsed;
0428         end
0429        
0430 
0431 
0432     %Add spcs to the calculated solution
0433     solution=P'*solution + specified_velocity;
0434 
0435     %Recover solution vector
0436     u=solution(1:nods);
0437     v=solution(nods+1:2*nods);
0438 
0439     %Compute viscosity
0440     nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
0441     change=1 - nu_bar_oldvalue./nu_bar; 
0442 
0443     %Update viscosity
0444     nu_bar=nu_bar.*(1+viscosity_overshoot*change);
0445     location=find(nu_bar<=0);
0446     nu_bar(location)=-nu_bar(location);
0447 
0448     %Test for direct shooting convergence
0449     if convergence_count>1,
0450 
0451         ug=[u_old;v_old];
0452         nug=norm(ug,2);
0453         dug=[u; v]-[u_old; v_old];
0454         ndug=norm(dug,2);
0455         relative_change=ndug/nug;
0456 
0457 
0458         %Figure out if viscosity converged
0459         if relative_change<criterion_rel,
0460             disp(sprintf('%s %g %s %g %s','   Convergence criterion: norm(du)/norm(u)=',relative_change,' < ',criterion_rel,' m/yr'));
0461             converged_yet=1;
0462         else
0463             disp(sprintf('%s %g %s %g %s','   Convergence criterion: norm(du)/norm(u)=',relative_change,' > ',criterion_rel,' m/yr'));
0464             converged_yet=0;
0465         end
0466 
0467         if ~isnan(criterion_abs),
0468             change=max(abs(dug))*yts;
0469             if change<criterion_abs
0470                 disp(sprintf('%s %g %s %g %s','   Convergence criterion: max(du)=',change,' < ',criterion_abs,' m/yr'));
0471                 converged_yet=1;
0472             else
0473                 disp(sprintf('%s %g %s %g %s','   Convergence criterion: max(du)=',change,' > ',criterion_abs,' m/yr'));
0474                 converged_yet=0;
0475             end
0476         end
0477     else 
0478         converged_yet=0;
0479     end
0480 
0481     u_old=u;
0482     v_old=v;
0483 
0484     convergence_count=convergence_count+1;
0485         
0486 end % This end statement terminates the "while" command way above
0487 
0488 %Load results onto md:
0489 md.vx=u*yts;
0490 md.vy=v*yts;
0491 md.vel=sqrt(md.vx.^2+md.vy.^2);
0492 
0493 %stop timing
0494 t2=clock;
0495 
0496 disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']));
0497 
0498 function macayealdiagnosticusage();
0499 disp('md=macayealdiagnostic(md)');
0500 disp('   where md is a structure of class @model');

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