source: issm/trunk/src/m/solutions/macayeal/diagnostic.m@ 119

Last change on this file since 119 was 119, checked in by Mathieu Morlighem, 16 years ago

change the way viscosity is computed -> same as ICE and CIELO
File size: 12.8 KB
Line 
1function md=diagnostic(md)
2%DIAGNOSTIC - compute the velocity field of a 2d model using MacAyeal solution
3%
4% this routine solves the problem using MacAyeal's model. It calculates the velocity
5% field corresponding to the parameters and the geometry given by the model md
6%
7% Usage:
8% md=diagnostic(md)
9
10%First check we do have the correct argument number
11if ((nargin~=1) || (nargout~=1)),
12 macayealdiagnosticusage();
13 error('macayealdiagnostic error message: incorrect number of input and output arguments');
14end
15
16if ~isa(md,'model'),
17 macayealdiagnosticusage();
18 error('macayealdiagnostic error message: input argument is not a @model object');
19end
20
21%start timing
22t1=clock;
23
24%Transfer model fields into matlab variables
25x=md.x;
26y=md.y;
27
28index=md.elements;index=sort(index,2); %necessary
29nods=md.numberofgrids;
30nel=md.numberofelements;
31z_thick=md.thickness;
32z_surf=md.surface;
33z_bed=md.bed;
34z_thick_bar=(z_thick(index(:,1))+z_thick(index(:,2))+z_thick(index(:,3)))/3;
35rho_ice=md.rho_ice;
36rho_water=md.rho_water;
37g=md.g;
38viscosity_overshoot=md.viscosity_overshoot;
39index_icefront=md.segmentonneumann_diag; index_icefront=index_icefront(:,1:2); %we strip the last column, which holds the element number for the boundary segment
40nodes_on_boundary=md.gridonboundary;
41nodes_on_icefront=zeros(nods,1); nodes_on_icefront(index_icefront)=1;
42node_on_dirichlet=md.gridondirichlet_diag;
43nodes_on_icesheet=md.gridonicesheet;
44element_on_icesheet=md.elementonicesheet;
45qcoeff=md.q;
46pcoeff=md.p;
47drag_type=md.drag_type;
48drag=md.drag;
49
50criterion_rel=md.eps_rel;
51criterion_abs=md.eps_abs;
52yts=md.yts;
53B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
54glen_coeff=md.n;
55
56%average of p and q over the grids (size nel->nods)
57pcoeff_grid=zeros(nods,1);
58qcoeff_grid=zeros(nods,1);
59for i=1:nods
60 %1: find the elements that contain the grid i
61 neighbors_gridi=[];
62 for j=1:3
63 neighbors_gridi=[neighbors_gridi find(index(:,j)==i)'];
64 end
65 numberofneighbors_gridi=length(neighbors_gridi);
66 %2 retrieve the value of p and q over each of these elements. The average is
67 %plugged into dbx_grid
68 qcoeff_grid(i)=sum(qcoeff(neighbors_gridi))/numberofneighbors_gridi;
69 pcoeff_grid(i)=sum(pcoeff(neighbors_gridi))/numberofneighbors_gridi;
70end
71
72%Build length_icefront and normal_icefront:
73[length_icefront,normal_icefront]=buildicefrontnormal(x,y,index_icefront);
74
75%Start building areas
76aire=zeros(md.numberofelements,1);
77
78for n=1:nel
79 aire(n)=1/2 * det([1 1 1;x(index(n,:))';y(index(n,:))']);
80end
81
82aire=abs(aire); % if index is sorted from its original value, then aire could be negative
83
84alpha=zeros(nel,3);
85beta=zeros(nel,3);
86gamma=zeros(nel,3);
87
88for n=1:nel
89 X=inv([x(index(n,:)) y(index(n,:)) ones(3,1)]);
90 alpha(n,:)=X(1,:);
91 beta(n,:)=X(2,:);
92 gamma(n,:)=X(3,:);
93end
94
95clear X
96
97
98%Do once and for all the initial computation of matrix-locations:
99row_location=zeros(nel*3*3,1);
100col_location=zeros(nel*3*3,1);
101row_location_AD=zeros(nel*6,1);
102col_location_AD=zeros(nel*6,1);
103
104count=-nel+1;
105
106for i=1:3
107 for j=1:3
108 count=count+nel;
109 row_location(count:count+nel-1)=index(:,i);
110 col_location(count:count+nel-1)=index(:,j);
111 end
112end
113
114count=-nel+1;
115
116for i=1:3
117 for j=i:3
118 count=count+nel;
119 row_location_AD(count:count+nel-1)=index(:,i);
120 col_location_AD(count:count+nel-1)=index(:,j);
121 end
122end
123
124permanent_pieces_of_A=zeros(nel*6,1);
125permanent_pieces_of_B=zeros(nel*3*3,1);
126permanent_pieces_of_C=zeros(nel*3*3,1);
127permanent_pieces_of_D=zeros(nel*6,1);
128
129count=-nel+1;
130
131for i=1:3
132 for j=i:3
133 count=count+nel;
134 permanent_pieces_of_A(count:count+nel-1)= z_thick_bar .* aire ...
135 .*(2*alpha(:,i).*alpha(:,j) + 1/2*beta(:,i).*beta(:,j));
136 % This loop structure works only when index is sorted
137 permanent_pieces_of_D(count:count+nel-1)= z_thick_bar .* aire ...
138 .*(2*beta(:,i).*beta(:,j) + 1/2*alpha(:,i).*alpha(:,j));
139 end
140end
141
142count=-nel+1;
143
144for i=1:3
145 for j=1:3
146 count=count+nel;
147
148 permanent_pieces_of_B(count:count+nel-1)= z_thick_bar .* aire ...
149 .*(alpha(:,j).*beta(:,i) + 1/2*beta(:,j).*alpha(:,i));
150
151 permanent_pieces_of_C(count:count+nel-1)= z_thick_bar .* aire ...
152 .*(beta(:,j).*alpha(:,i) + 1/2*alpha(:,j).*beta(:,i));
153 end
154end
155
156% Step 3 -- Set up right-hand side of the problem.
157
158% (Note to myself: to avoid vector dependency problem, yet still vectorize,
159% I treat the Rhs vector as a sparse matrix
160
161Rhs_x=zeros(nel*27,1);
162Rhs_y=zeros(nel*27,1);
163Rhs_y=ones(nel*27,1);
164row_rhs=zeros(nel*27,1);
165
166count=-nel+1;
167for i=1:3
168 for n=1:3
169 for m=1:3
170 count=count+nel;
171
172 Rhs_x(count:count+nel-1)= -rho_ice * g * ...
173 z_thick(index(:,n)).*z_surf(index(:,m)) ...
174 .* aire(:) .* alpha(:,m) * ( (n==i)/6 + (n~=i)/12 );
175
176 Rhs_y(count:count+nel-1)= -rho_ice * g * ...
177 z_thick(index(:,n)).*z_surf(index(:,m)) ...
178 .* aire(:) .* beta(:,m) .* ( (n==i)/6 + (n~=i)/12 );
179
180 row_rhs(count:count+nel-1)=index(:,i);
181
182 end
183 end
184end
185
186Rhs=full([sparse(row_rhs,ones(nel*27,1),Rhs_x,nods,1)
187 sparse(row_rhs,ones(nel*27,1),Rhs_y,nods,1)]);
188
189 for k=1:2
190 for l=1:2
191 for j=1:2
192 Rhs(index_icefront(:,k))=Rhs(index_icefront(:,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) );
200
201 Rhs(index_icefront(:,k)+nods)=Rhs(index_icefront(:,k)+nods) + ...
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) );
209 end
210 end
211end
212
213clear Rhs_x Rhs_y row_rhs
214
215% Step 4 -- Create the parsing matrix to "wring out"
216% the known boundary conditions
217
218% kinematic condition (u,v=0) specified at non-ice-front boundaries;
219num_specified= 2*sum(node_on_dirichlet);
220row=zeros(2*nods-num_specified,1);
221col=zeros(2*nods-num_specified,1);
222value=zeros(2*nods-num_specified,1);
223count=0;
224for n=1:nods
225 if(~node_on_dirichlet(n))
226 count=count+1;
227 row(count)=count;
228 col(count)=n;
229 value(count)=1;
230 end
231end
232for n=1:nods
233 if (~node_on_dirichlet(n))
234 count=count+1;
235 row(count)=count;
236 col(count)=nods+n;
237 value(count)=1;
238 end
239end
240
241P=sparse(row,col,value,2*nods-num_specified,2*nods);
242
243specified_velocity=zeros(2*nods,1);
244pos=find(node_on_dirichlet);
245specified_velocity(pos)=md.dirichletvalues_diag(pos,1)/md.yts;
246specified_velocity(nods+pos)=md.dirichletvalues_diag(pos,2)/md.yts;
247
248
249% ######## an attention grabbing break ########
250% Iterate on flow law till converged
251
252converged_yet=0;
253loop=0;
254
255u_old=zeros(nods,1);
256v_old=zeros(nods,1);
257convergence_count=1;
258
259while (~converged_yet)
260
261 if (loop>(100))
262 warning('Maximum Viscosity Iterations Reached.')
263 break;
264 end;
265 loop=loop+1;
266
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
284
285 % Step 4 -- Set up stress-balance matrix (whos solution is the velocity
286 % field):
287
288 matrix_value_A=zeros(nel*6,1);
289 matrix_value_B=zeros(nel*3*3,1);
290 matrix_value_C=zeros(nel*3*3,1);
291 matrix_value_D=zeros(nel*6,1);
292
293 count=-nel+1;
294
295 for i=1:3
296 for j=i:3
297 count=count+nel;
298
299 matrix_value_A(count:count+nel-1)=nu_bar .* ...
300 permanent_pieces_of_A(count:count+nel-1);
301
302 matrix_value_D(count:count+nel-1)=nu_bar .* ...
303 permanent_pieces_of_D(count:count+nel-1);
304 end
305 end
306
307 count=-nel+1;
308
309 for i=1:3
310 for j=1:3
311 count=count+nel;
312
313 matrix_value_B(count:count+nel-1)=nu_bar .* ...
314 permanent_pieces_of_B(count:count+nel-1);
315
316 matrix_value_C(count:count+nel-1)=nu_bar .* ...
317 permanent_pieces_of_C(count:count+nel-1);
318 end
319 end
320
321 A=sparse(row_location_AD,col_location_AD,matrix_value_A,nods,nods);
322 A=A+triu(A,1)';
323 B=sparse(row_location,col_location,matrix_value_B,nods,nods);
324 C=sparse(row_location,col_location,matrix_value_C,nods,nods);
325 D=sparse(row_location_AD,col_location_AD,matrix_value_D,nods,nods);
326 D=D+triu(D,1)';
327
328 F=[A C
329 B D];
330
331 %Now, take care of the basal friction if there is any: to make things easier, we translate u = k *Neff^(-q)* sigma^p into
332 % sigma= drag^2 * Neff ^(r) * u^s with r=q/p and s=1/p : */
333
334 if (drag_type==2),
335 %compute coeffs:
336 rcoeff=qcoeff_grid./pcoeff_grid;
337 scoeff=1./pcoeff_grid;
338
339 %initialization of basal drag stiffness
340 Dragoperator=spalloc(2*nods,2*nods,0);
341
342 if loop~=1,
343
344 %retrieve the velocity magnitude
345 velocity_mag=sqrt(solution(1:nods).^2+solution(nods+1:2*nods).^2);
346
347 %Computation of the effective pressure
348 Neff=g*(rho_ice*z_thick+rho_water*z_bed);
349
350 %If effective pressure becomes negative, sliding becomes unstable (Paterson 4th edition p 148). This is because
351 %the water pressure is so high, the ice sheet elevates over its ice bumps and slides. But the limit behaviour
352 %for this should be an ice shelf sliding (no basal drag). Therefore, for any effective pressure Neff < 0, we should
353 %replace it by Neff=0 (ie, equival it to an ice shelf)*/
354 pos=find(Neff<0);
355 Neff(pos)=0;
356
357 %Basal drag coefficient: Tau_x=-alpha^2 u, Tau_y=-alpha^2 v (See
358 %MacAyeal)
359 alpha2=(drag.^2).*(Neff.^rcoeff).*(velocity_mag.^(scoeff-1));
360
361 %stiffness due to basal drag
362 %initialization
363 count=0;
364 value=zeros(nel,27);
365 row=zeros(nel,27);
366 col=zeros(nel,27);
367
368 for m=1:3
369 for k=1:3
370 for l=1:3
371 if ( (m==k) + (m==l) + (l==k) )==3
372 fac=1/10;
373 elseif ( (m==k) + (m==l) + (l==k) )==1
374 fac=1/30;
375 else
376 fac=1/60;
377 end
378
379 count=count+1;
380 row(:,count)=index(:,m);
381 col(:,count)=index(:,k);
382 value(:,count)=fac*aire(:).* ...
383 (alpha2(index(:,l)));%.*element_on_icesheet(index(:,l));
384
385 end
386 end
387 end
388
389 Dragoperator=[sparse(row,col,value,nods,nods) spalloc(nods,nods,0)
390 spalloc(nods,nods,0) sparse(row,col,value,nods,nods)];
391 end
392
393 F=F+Dragoperator; %plug into the global stiffness matrix
394
395 end
396
397 Rhs_parsed=P*(Rhs - F*specified_velocity);
398 F=P*F*P';
399
400 % Step 5 -- Solve the problem:
401
402 % Digression: if need be clear up some memory:
403
404 clear matrix_value_A ...
405 matrix_value_B matrix_value_C matrix_value_D A B C D
406
407 % We can use either the LU or the Cholesky decomposition, but the
408 % Cholesky decomposition is twice as efficient as LU for symmetric
409 % definite positive matrix
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);
414
415 %Add spcs to the calculated solution
416 solution=P'*solution + specified_velocity;
417
418 %Recover solution vector
419 u=solution(1:nods);
420 v=solution(nods+1:2*nods);
421
422 %Test for direct shooting convergence
423 if convergence_count>1,
424
425 ug=[u_old;v_old];
426 nug=norm(ug,2);
427 dug=[u; v]-[u_old; v_old];
428 ndug=norm(dug,2);
429 relative_change=ndug/nug;
430
431 %Figure out if viscosity converged
432 if relative_change<criterion_rel,
433 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' < ',criterion_rel,' m/yr'));end
434 converged_yet=1;
435 else
436 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: norm(du)/norm(u)=',relative_change,' > ',criterion_rel,' m/yr'));end
437 converged_yet=0;
438 end
439
440 if ~isnan(criterion_abs),
441 change=max(abs(dug))*yts;
442 if change<criterion_abs
443 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' < ',criterion_abs,' m/yr'));end
444 else
445 if md.debug, disp(sprintf('%s %g %s %g %s',' Convergence criterion: max(du)=',change,' > ',criterion_abs,' m/yr'));end
446 converged_yet=0;
447 end
448 end
449 else
450 converged_yet=0;
451 end
452
453 convergence_count=convergence_count+1;
454
455end % This end statement terminates the "while" command way above
456
457%Load results onto md:
458md.vx=u*yts;
459md.vy=v*yts;
460md.vel=sqrt(md.vx.^2+md.vy.^2);
461
462%stop timing
463t2=clock;
464
465disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']));
466
467function macayealdiagnosticusage();
468disp('md=macayealdiagnostic(md)');
469disp(' where md is a structure of class @model');
Note: See TracBrowser for help on using the repository browser.