0001 function md=macayealdiagnostic(md)
0002
0003
0004
0005
0006
0007
0008
0009
0010
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
0022 t1=clock;
0023
0024
0025 x=md.x;
0026 y=md.y;
0027
0028 index=md.elements;index=sort(index,2);
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);
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
0057 pcoeff_grid=zeros(nods,1);
0058 qcoeff_grid=zeros(nods,1);
0059 for i=1:nods
0060
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
0067
0068 qcoeff_grid(i)=sum(qcoeff(neighbors_gridi))/numberofneighbors_gridi;
0069 pcoeff_grid(i)=sum(pcoeff(neighbors_gridi))/numberofneighbors_gridi;
0070 end
0071
0072
0073 [length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront);
0074
0075
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);
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
0098 nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
0099
0100
0101
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
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
0167
0168
0169
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
0228
0229
0230
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
0263
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
0283
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
0336
0337
0338 if (drag_type==2),
0339
0340 rcoeff=qcoeff_grid./pcoeff_grid;
0341 scoeff=1./pcoeff_grid;
0342
0343
0344 Dragoperator=spalloc(2*nods,2*nods,0);
0345
0346 if loop~=1,
0347
0348
0349 velocity_mag=sqrt(solution(1:nods).^2+solution(nods+1:2*nods).^2);
0350
0351
0352 Neff=g*(rho_ice*z_thick+rho_water*z_bed);
0353
0354
0355
0356
0357
0358 pos=find(Neff<0);
0359 Neff(pos)=0;
0360
0361
0362
0363 alpha2=(drag.^2).*(Neff.^rcoeff).*(velocity_mag.^(scoeff-1));
0364
0365
0366
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)));
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;
0398
0399 end
0400
0401
0402 Rhs_parsed=P*(Rhs - F*specified_velocity);
0403 F=P*F*P';
0404
0405
0406
0407
0408
0409
0410 clear matrix_value_A ...
0411 matrix_value_B matrix_value_C matrix_value_D A B C D
0412
0413
0414
0415
0416
0417 if strcmpi(md.solver_type,'lu'),
0418
0419 [L,U] = lu(F);
0420 solution= U\(L\Rhs_parsed);
0421 elseif strcmpi(md.solver_type,'cholesky'),
0422
0423 L = chol(F);
0424 solution= L\(L'\Rhs_parsed);
0425 else
0426
0427 solution = F\Rhs_parsed;
0428 end
0429
0430
0431
0432
0433 solution=P'*solution + specified_velocity;
0434
0435
0436 u=solution(1:nods);
0437 v=solution(nods+1:2*nods);
0438
0439
0440 nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
0441 change=1 - nu_bar_oldvalue./nu_bar;
0442
0443
0444 nu_bar=nu_bar.*(1+viscosity_overshoot*change);
0445 location=find(nu_bar<=0);
0446 nu_bar(location)=-nu_bar(location);
0447
0448
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
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
0487
0488
0489 md.vx=u*yts;
0490 md.vy=v*yts;
0491 md.vel=sqrt(md.vx.^2+md.vy.^2);
0492
0493
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');