![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 119
- Timestamp:
- 04/29/09 08:36:46 (16 years ago)
- File:
-
- 1 edited
-
issm/trunk/src/m/solutions/macayeal/diagnostic.m (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/solutions/macayeal/diagnostic.m
r67 r119 95 95 clear X 96 96 97 %Initialize viscosity 98 nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff); 99 100 101 % Do once and for all the initial computation of matrix-locations: 102 97 98 %Do once and for all the initial computation of matrix-locations: 103 99 row_location=zeros(nel*3*3,1); 104 100 col_location=zeros(nel*3*3,1); … … 116 112 end 117 113 118 119 114 count=-nel+1; 120 115 … … 126 121 end 127 122 end 128 129 130 123 131 124 permanent_pieces_of_A=zeros(nel*6,1); … … 144 137 permanent_pieces_of_D(count:count+nel-1)= z_thick_bar .* aire ... 145 138 .*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j)); 146 147 139 end 148 140 end … … 159 151 permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ... 160 152 .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)); 161 162 163 153 end 164 154 end … … 173 163 Rhs_y=ones(nel*27,1); 174 164 row_rhs=zeros(nel*27,1); 175 176 165 177 166 count=-nel+1; … … 198 187 sparse(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]); 199 188 200 for k=1:2201 for l=1:2202 for j=1:2189 for k=1:2 190 for l=1:2 191 for j=1:2 203 192 Rhs(index_icefront(:,k))=Rhs(index_icefront(:,k)) + ... 204 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ...205 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ...206 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ...207 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))...208 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))...209 .*normal_icefront(:,1).*length_icefront(:) ...210 /(4*(l==k & j==k) + 12*(l~=k | j~=k) );193 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... 194 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... 195 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... 196 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... 197 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... 198 .*normal_icefront(:,1).*length_icefront(:) ... 199 /(4*(l==k & j==k) + 12*(l~=k | j~=k) ); 211 200 212 201 Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ... 213 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... 214 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... 215 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... 216 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... 217 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... 218 .*normal_icefront(:,2).*length_icefront(:) ... 219 /(4*(l==k & j==k) + 12*(l~=k | j~=k) ); 220 202 (rho_ice*g/2*z_thick(index_icefront(:,l)).*z_thick(index_icefront(:,j)) ... 203 + rho_water*g/2*(md.gridoniceshelf(index_icefront(:,k))) ... 204 .*(-min(0,z_bed(index_icefront(:,l))).*min(0,z_bed(index_icefront(:,j))) ... 205 +min(0,z_thick(index_icefront(:,l))+z_bed(index_icefront(:,l)))... 206 .*min(0,z_thick(index_icefront(:,j))+z_bed(index_icefront(:,j)))))... 207 .*normal_icefront(:,2).*length_icefront(:) ... 208 /(4*(l==k & j==k) + 12*(l~=k | j~=k) ); 221 209 end 222 210 end … … 253 241 P=sparse(row,col,value,2*nods-num_specified,2*nods); 254 242 255 256 243 specified_velocity=zeros(2*nods,1); 257 244 pos=find(node_on_dirichlet); … … 278 265 loop=loop+1; 279 266 280 nu_bar_oldvalue=nu_bar; 267 %Compute viscosity (as in ICE and CIELO) 268 if convergence_count==1; 269 %Initialize viscosity 270 nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff); 271 elseif convergence_count==2, 272 nu_bar_oldvalue=nu_bar; 273 nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); 274 nu_bar=nu_bar + viscosity_overshoot*(nu_bar-nu_bar_oldvalue); 275 nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0)); 276 u_old=u; v_old=v; 277 else 278 nu_bar_oldvalue=viscosity(index,nel,alpha,beta,u_old,v_old,B_bar,glen_coeff); 279 nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); 280 nu_bar=nu_bar + viscosity_overshoot*(nu_bar-nu_bar_oldvalue); 281 nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0)); 282 u_old=u; v_old=v; 283 end 281 284 282 285 % Step 4 -- Set up stress-balance matrix (whos solution is the velocity … … 299 302 matrix_value_D(count:count+nel-1)=nu_bar .* ... 300 303 permanent_pieces_of_D(count:count+nel-1); 301 302 303 304 end 304 305 end … … 315 316 matrix_value_C(count:count+nel-1)=nu_bar .* ... 316 317 permanent_pieces_of_C(count:count+nel-1); 317 318 319 end 320 end 321 322 318 end 319 end 323 320 324 321 A=sparse(row_location_AD,col_location_AD,matrix_value_A,nods,nods); … … 328 325 D=sparse(row_location_AD,col_location_AD,matrix_value_D,nods,nods); 329 326 D=D+triu(D,1)'; 330 331 327 332 328 F=[A C … … 399 395 end 400 396 401 402 397 Rhs_parsed=P*(Rhs - F*specified_velocity); 403 398 F=P*F*P'; 404 405 399 406 400 % Step 5 -- Solve the problem: … … 414 408 % Cholesky decomposition is twice as efficient as LU for symmetric 415 409 % definite positive matrix 416 417 if strcmpi(md.solver_type,'lu'), 418 % Solve by LU decomposition. 419 [L,U] = lu(F); 420 solution= U\(L\Rhs_parsed); 421 elseif strcmpi(md.solver_type,'cholesky'), 422 % Solve by Choleski decomposition. 423 L = chol(F); 424 solution= L\(L'\Rhs_parsed); 425 else 426 % use matlab's generic solver 427 solution = F\Rhs_parsed; 428 end 429 430 410 if md.debug, 411 disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(F))); 412 end 413 solution=Solver(F,Rhs_parsed,md.solver_type); 431 414 432 415 %Add spcs to the calculated solution … … 436 419 u=solution(1:nods); 437 420 v=solution(nods+1:2*nods); 438 439 %Compute viscosity440 nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);441 change=1 - nu_bar_oldvalue./nu_bar;442 443 %Update viscosity444 nu_bar=nu_bar.*(1+viscosity_overshoot*change);445 location=find(nu_bar<=0);446 nu_bar(location)=-nu_bar(location);447 421 448 422 %Test for direct shooting convergence … … 455 429 relative_change=ndug/nug; 456 430 457 458 431 %Figure out if viscosity converged 459 432 if relative_change<criterion_rel, 460 disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' < ',criterion_rel,' m/yr'));433 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' < ',criterion_rel,' m/yr'));end 461 434 converged_yet=1; 462 435 else 463 disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' > ',criterion_rel,' m/yr'));436 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' > ',criterion_rel,' m/yr'));end 464 437 converged_yet=0; 465 438 end … … 468 441 change=max(abs(dug))*yts; 469 442 if change<criterion_abs 470 disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' < ',criterion_abs,' m/yr'));443 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' < ',criterion_abs,' m/yr'));end 471 444 else 472 disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' > ',criterion_abs,' m/yr'));445 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' > ',criterion_abs,' m/yr'));end 473 446 converged_yet=0; 474 447 end … … 477 450 converged_yet=0; 478 451 end 479 480 u_old=u;481 v_old=v;482 452 483 453 convergence_count=convergence_count+1;
Note:
See TracChangeset
for help on using the changeset viewer.
