Index: /issm/trunk/src/m/solutions/macayeal/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/macayeal/diagnostic.m (revision 118) +++ /issm/trunk/src/m/solutions/macayeal/diagnostic.m (revision 119) @@ -95,10 +95,6 @@ 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: - + +%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); @@ -116,5 +112,4 @@ end - count=-nel+1; @@ -126,6 +121,4 @@ end end - - permanent_pieces_of_A=zeros(nel*6,1); @@ -144,5 +137,4 @@ 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 @@ -159,6 +151,4 @@ permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ... .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i)); - - end end @@ -173,5 +163,4 @@ Rhs_y=ones(nel*27,1); row_rhs=zeros(nel*27,1); - count=-nel+1; @@ -198,25 +187,24 @@ sparse(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]); -for k=1:2 - for l=1:2 - for j=1:2 + 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) ); + (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) ); - + (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 @@ -253,5 +241,4 @@ P=sparse(row,col,value,2*nods-num_specified,2*nods); - specified_velocity=zeros(2*nods,1); pos=find(node_on_dirichlet); @@ -278,5 +265,21 @@ loop=loop+1; - nu_bar_oldvalue=nu_bar; + %Compute viscosity (as in ICE and CIELO) + if convergence_count==1; + %Initialize viscosity + nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff); + elseif convergence_count==2, + nu_bar_oldvalue=nu_bar; + nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); + nu_bar=nu_bar + viscosity_overshoot*(nu_bar-nu_bar_oldvalue); + nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0)); + u_old=u; v_old=v; + else + nu_bar_oldvalue=viscosity(index,nel,alpha,beta,u_old,v_old,B_bar,glen_coeff); + nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff); + nu_bar=nu_bar + viscosity_overshoot*(nu_bar-nu_bar_oldvalue); + nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0)); + u_old=u; v_old=v; + end % Step 4 -- Set up stress-balance matrix (whos solution is the velocity @@ -299,6 +302,4 @@ matrix_value_D(count:count+nel-1)=nu_bar .* ... permanent_pieces_of_D(count:count+nel-1); - - end end @@ -315,10 +316,6 @@ matrix_value_C(count:count+nel-1)=nu_bar .* ... permanent_pieces_of_C(count:count+nel-1); - - - end - end - - + end + end A=sparse(row_location_AD,col_location_AD,matrix_value_A,nods,nods); @@ -328,5 +325,4 @@ D=sparse(row_location_AD,col_location_AD,matrix_value_D,nods,nods); D=D+triu(D,1)'; - F=[A C @@ -399,8 +395,6 @@ end - Rhs_parsed=P*(Rhs - F*specified_velocity); F=P*F*P'; - % Step 5 -- Solve the problem: @@ -414,19 +408,8 @@ % 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 - - + if md.debug, + disp(sprintf('%s%g',' condition number of stiffness matrix: ',condest(F))); + end + solution=Solver(F,Rhs_parsed,md.solver_type); %Add spcs to the calculated solution @@ -436,13 +419,4 @@ 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 @@ -455,11 +429,10 @@ relative_change=ndug/nug; - %Figure out if viscosity converged if relative_change ',criterion_rel,' m/yr')); + if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' > ',criterion_rel,' m/yr'));end converged_yet=0; end @@ -468,7 +441,7 @@ change=max(abs(dug))*yts; if change ',criterion_abs,' m/yr')); + if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' > ',criterion_abs,' m/yr'));end converged_yet=0; end @@ -477,7 +450,4 @@ converged_yet=0; end - - u_old=u; - v_old=v; convergence_count=convergence_count+1;