source: issm/trunk/src/m/classes/public/ismodelselfconsistent.m@ 358

Last change on this file since 358 was 358, checked in by Eric.Larour, 16 years ago

New slope extrude module. New slope computation in Penta and Tria
File size: 13.8 KB
RevLine 
[1]1function bool=ismodelselfconsistent(md,solutiontype,package),
2%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
3%
4% Usage:
5% bool=ismodelselfconsistent(md,solutiontype,package),
6
7%tolerance we use in litmus tests for the consistency of the model
8tolerance=10^-12;
9if (nargin~=3 )
[27]10 help ismodelselfconsistent
11 error('ismodelselfconsistent error message: wrong usage');
[1]12end
13
[27]14%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMON CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
15
16%COUNTER
17if md.counter<3,
18 disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
19 bool=0;return;
[1]20end
21
[27]22%NAME
[1]23if isempty(md.name),
24 disp(['model is not correctly configured: missing name!']);
25 bool=0;return;
26end
27
[27]28%MESH
[1]29if md.numberofelements<=0,
30 disp(['model ' md.name ' does not have any elements!']);
31 bool=0; return;
32end
33if md.numberofgrids<=0,
34 disp(['model ' md.name ' does not have any grids!']);
35 bool=0; return;
36end
37
[27]38%ELEMENTSTYPE
[1]39if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
40 disp(['Types of elements have not been set properly, run setelementstype first'])
41 bool=0;return;
42end
43if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum) + (md.elements_type(:,1)==pattynenum)))
44 disp(['Types of elements have not been set properly, run setelementstype first'])
45 bool=0;return;
46end
47if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum)))
48 disp(['Types of elements have not been set properly, run setelementstype first'])
49 bool=0;return;
50end
51if strcmpi(md.type,'2d'),
52 if (ismember(pattynenum,md.elements_type(:,1)) | ismember(stokesenum,md.elements_type(:,2))),
53 disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
54 bool=0;return;
55 end
56end
[27]57
[358]58if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
59 disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
60 bool=0;return;
61end
62
63
[27]64%NO NAN
65fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
66 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
67 'tolx','np','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
68for i=1:length(fields),
69 if ~isempty(eval(['md.' char(fields(i))])),
70 if find(isnan(eval(['md.' char(fields(i))]))),
71 disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
72 bool=0; return;
73 end
[1]74 end
75end
76
[27]77%FIELDS > 0
78fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
79 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
80 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
81for i=1:length(fields),
82 if ~isempty(eval(['md.' char(fields(i))])),
83 if find((eval(['md.' char(fields(i))]))<0),
84 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
85 bool=0; return;
[1]86 end
87 end
88end
[27]89if any(md.p<=0),
90 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
91 bool=0; return;
92end
[1]93
[27]94%FIELDS ~=0
95fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
96 'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx',...
97 'sparsity','deltaH','DeltaH','timeacc','timedec'};
98for i=1:length(fields),
99 if ~isempty(eval(['md.' char(fields(i))])),
100 if find((eval(['md.' char(fields(i))]))==0),
101 disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
102 bool=0; return;
103 end
[1]104 end
105end
106
[27]107%SIZE NUMBEROFELEMENTS
108fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
109for i=1:size(fields,2),
110 if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
111 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
112 bool=0; return;
113 end
114end
115
116%SIZE NUMBEROFGRIDS
117fields={'x','y','z','B','drag','gridondirichlet_diag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
118for i=1:length(fields),
119 if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
120 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
121 bool=0; return;
122 end
123end
124
125%THICKNESS = SURFACE - BED
126if find((md.thickness-md.surface+md.bed)>tolerance),
127 disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
128 bool=0; return;
129end
130
131%RIFTS
132if md.numrifts,
133 if ~strcmpi(md.type,'2d'),
134 disp(['Models with rifts are only supported in 2d for now!']);
[1]135 bool=0;return;
136 end
137end
[27]138if ~isstruct(md.rifts),
139 if ~isempty(find(md.segmentmarkers>=2)),
140 %We have segments with rift markers, but no rift structure!
141 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
142 bool=0; return;
143 end
144end
[1]145
[27]146%ARTIFICIAL DIFFUSIVITY
[1]147if ~isscalar(md.artificial_diffusivity),
148 disp('artificial_diffusivity should be a scalar (1 or 0)');
149 bool=0;return;
150end
151
[58]152%CLUSTER
153if ~strcmpi(package,'cielo') & ~strcmpi(md.cluster,'none')
154 disp(['parallel solution not supported by package ' package '. Use cluster=''none'' ']);
155 bool=0;return;
156end
157
[27]158%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159
160%DIAGNOSTIC
161if strcmpi(solutiontype,'diagnostic')
162
163 %HUTTER ON ICESHELF WARNING
164 if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
165 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
166 end
167
168 %SINGULAR
169 if ~any(md.gridondirichlet_diag),
170 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
[1]171 bool=0;return;
172 end
173
[27]174 %DIRICHLET VALUES
175 if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
176 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
[1]177 bool=0;return;
178 end
[27]179
180 %DIRICHLET IF THICKNESS <= 0
181 if ~isempty(find(md.thickness<=0)),
182 pos=find(md.thickness<=0);
183 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
184 disp(['model ' md.name ' has some grids with 0 thickness']);
185 bool=0; return;
186 end
187 end
[1]188end
189
[27]190%PROGNOSTIC
191if strcmp(solutiontype,'prognostic'),
192 %old testing
193 % if isempty(find(md.gridondirichlet_diag)),
194 % disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
195 % bool=0;return;
196 % end
197 % if isempty(find(md.gridondirichlet_prog)),
198 % disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
199 % bool=0;return;
200 % end
[1]201
[27]202 %VELOCITIES
203 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
[1]204 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
[27]205 bool=0; return;
[1]206 end
207end
208
[27]209%THERMAL TRANSIENT
210if strcmp(solutiontype,'thermaltransient')
211
212 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
213 if isempty(md.temperature),
214 disp(['An initial temperature is needed for a transient thermal computation'])
[1]215 bool=0;return;
216 end
[27]217 if isstruct(md.temperature) | isstruct(md.melting) | isstruct(md.accumulation),
218 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
[1]219 bool=0;return;
220 end
221end
222
[27]223%THERMAL STEADY AND THERMAL TRANSIENT
224if strcmpi(solutiontype,'thermalsteady') | strcmp(solutiontype,'thermaltransient'),
225
226 %EXTRUSION
227 if strcmp(md.type,'2d'),
228 disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
[1]229 bool=0;return;
230 end
[27]231
232 %VELOCITIES AND PRESSURE
233 if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy) | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
234 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
[1]235 bool=0;return;
236 end
[27]237 if (isempty(md.pressure) | isnan(md.pressure)),
238 disp(['pressure is required. Run ''diagnostic'' solution first!'])
239 bool=0;return;
240 end
[1]241end
242
[27]243%THERMAL TRANSIENT AND TRANSIENT
244if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'transient'),
245
246 %DT and NDT
247 fields={'dt','ndt'};
248 for i=1:length(fields),
249 if find((eval(['md.' char(fields(i))]))<0),
250 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
251 bool=0; return;
252 end
253 end
254
255 %INITIAL TEMPERATURE
[1]256 if isstruct(md.temperature),
257 disp(['The initial temperature should be empty or a list but not a structure'])
258 bool=0;return;
259 end
260end
261
[27]262%PARAMETERS
[1]263if strcmp(solutiontype,'parameters')
[27]264
265 %OUTPUT
266 if ~iscell(md.parameteroutput)
267 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
268 bool=0; return;
269 end
270 for i=1:length(md.parameteroutput)
271 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
272 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
273 disp(['one of the parameteroutput is not supported yet']);
274 bool=0; return;
275 end
276 end
277 %VELOCITY
[1]278 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
279 disp(['velocities are required!']);
280 bool=0;return;
281 end
[27]282
283 %HUTTER
[1]284 if any(md.elements_type(:,1)==hutterenum);
285 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
286 bool=0;return;
287 end
288end
289
[27]290%CONTROL
291if strcmpi(solutiontype,'control'),
[1]292
[27]293 %CONTROL TYPE
[45]294 if ~ischar(md.control_type),
295 disp('control_type should be a string');
[1]296 bool=0;return;
297 end
298
[27]299 %LENGTH CONTROL FIELDS
[1]300 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
301 disp('maxiter, optscal and fit must have the length specified by nsteps')
302 bool=0;return;
303 end
[27]304
305 %FIT
[1]306 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
307 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
308 bool=0;return;
309 end
310
[27]311 %OBSERVED VELOCITIES
312 fields={'vx_obs','vy_obs'};
313 for i=1:length(fields),
314 if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
315 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
[1]316 bool=0; return;
317 end
318 end
319end
320
[27]321%QMU
322if strcmpi(solutiontype,'qmu'),
323 if ~strcmpi(md.cluster,'none'),
324 if md.waitonlock==0,
325 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
326 bool=0;return;
[1]327 end
328 end
329end
330
[27]331%MESH
332if strcmpi(solutiontype,'mesh'),
333 %this solution is a little special. It should come right after the md=model; operation. So a lot less checks!
[1]334
[27]335 bool=1;
336 return;
[1]337end
338
[132]339%MESH2GRID
340if strcmpi(solutiontype,'mesh2grid'),
341 if ~strcmpi(md.cluster,'none'),
342 disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
343 bool=0;return;
344 end
345end
346
347
[27]348%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[1]349
[27]350%CIELO
351if strcmpi(package,'cielo'),
[1]352
[27]353 %NAN VALUES
[64]354 fields={'sparsity'};
[27]355 for i=1:length(fields),
356 if ~isempty(eval(['md.' char(fields(i))])),
357 if find(isnan(eval(['md.' char(fields(i))]))),
358 disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
359 bool=0; return;
360 end
361 end
362 end
[1]363
[27]364 %FIELD > 0
[64]365 fields={'sparsity'};
[27]366 for i=1:length(fields),
367 if ~isempty(eval(['md.' char(fields(i))])),
368 if find((eval(['md.' char(fields(i))]))<0),
369 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
370 bool=0; return;
371 end
372 end
373 end
[1]374
[27]375 %FIELD ~= 0
[64]376 fields={'sparsity'};
[27]377 for i=1:length(fields),
378 if ~isempty(eval(['md.' char(fields(i))])),
379 if find((eval(['md.' char(fields(i))]))==0),
380 disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
381 bool=0; return;
382 end
383 end
[1]384 end
385
[27]386 %SPARSITY BETWEEN 0 AND 1
387 if ( (md.sparsity<=0) | (md.sparsity>1)),
[39]388 disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
[1]389 bool=0; return;
390 end
391
[27]392 %RIFTS
393 if md.numrifts,
394 if ~strcmpi(md.cluster,'none')
395 if isempty(findstr('aijmumps',md.solverstring)),
396 disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
397 bool=0;return;
398 end
399 end
400 end
[1]401
[27]402 %CONNECTIVITY
403 if strcmpi(md.type,'2d'),
404 if md.connectivity<9,
405 disp('connectivity should be at least 9 for 2d models');
406 bool=0;return;
407 end
408 end
409 if strcmpi(md.type,'3d'),
410 if md.connectivity<24,
411 disp('connectivity should be at least 24 for 3d models');
412 bool=0;return;
413 end
414 end
[1]415
[27]416 %LOWMEM = 0 or 1
417 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
418 disp(['model ' md.name ' lowmem field should be 0 or 1']);
[1]419 bool=0; return;
420 end
[64]421
422%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
423
424 if ~strcmpi(md.cluster,'none'),
425
426 %NAN VALUES
427 fields={'time','np'};
428 for i=1:length(fields),
429 if ~isempty(eval(['md.' char(fields(i))])),
430 if find(isnan(eval(['md.' char(fields(i))]))),
431 disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
432 bool=0; return;
433 end
434 end
435 end
436
437 %FIELD > 0
438 fields={'time','np'};
439 for i=1:length(fields),
440 if ~isempty(eval(['md.' char(fields(i))])),
441 if find((eval(['md.' char(fields(i))]))<0),
442 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
443 bool=0; return;
444 end
445 end
446 end
447
448 %FIELD ~= 0
449 fields={'time','np'};
450 for i=1:length(fields),
451 if ~isempty(eval(['md.' char(fields(i))])),
452 if find((eval(['md.' char(fields(i))]))==0),
453 disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
454 bool=0; return;
455 end
456 end
457 end
458
459 end
460
[1]461end
462
463%No problems, just return 1;
464bool=1;
465return;
Note: See TracBrowser for help on using the repository browser.