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

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

do not check evaluation_concurrency if qmu_analysis
File size: 13.7 KB
Line 
1function bool=ismodelselfconsistent(md,package),
2%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
3%
4% Usage:
5% bool=ismodelselfconsistent(md,package),
6
7%tolerance we use in litmus tests for the consistency of the model
8tolerance=10^-12;
9if (nargin~=2 )
10 help ismodelselfconsistent
11 error('ismodelselfconsistent error message: wrong usage');
12end
13
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;
20end
21
22%NAME
23if isempty(md.name),
24 disp(['model is not correctly configured: missing name!']);
25 bool=0;return;
26end
27
28%MESH
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
38%ELEMENTSTYPE
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
57
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
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(md.(fields{i})),
70 if any(isnan(md.(fields{i}))),
71 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
72 bool=0; return;
73 end
74 end
75end
76
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(md.(fields{i})),
83 if any(md.(fields{i})<0),
84 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
85 bool=0; return;
86 end
87 end
88end
89if any(md.p<=0),
90 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
91 bool=0; return;
92end
93
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(md.(fields{i})),
100 if any(md.(fields{i})==0),
101 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
102 bool=0; return;
103 end
104 end
105end
106
107%SIZE NUMBEROFELEMENTS
108fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
109for i=1:size(fields,2),
110 if (size(md.(fields{i}),1)~=md.numberofelements),
111 disp(['model ' md.name ' field ' 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 length(md.(fields{i}))~=md.numberofgrids,
120 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
121 bool=0; return;
122 end
123end
124
125%THICKNESS = SURFACE - BED
126if any((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!']);
135 bool=0;return;
136 end
137end
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
145
146%ARTIFICIAL DIFFUSIVITY
147if ~isscalar(md.artificial_diffusivity),
148 disp('artificial_diffusivity should be a scalar (1 or 0)');
149 bool=0;return;
150end
151
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
158%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
159%QMU
160if md.qmu_analysis,
161 if md.qmu_params.evaluation_concurrency~=1,
162 disp(['concurrency should be set to 1 when running dakota in library mode']);
163 bool=0;return;
164 end
165end
166
167%DIAGNOSTIC
168if strcmpi(md.analysis_type,'diagnostic')
169
170 %HUTTER ON ICESHELF WARNING
171 if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
172 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
173 end
174
175 %SINGULAR
176 if ~any(md.gridondirichlet_diag),
177 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
178 bool=0;return;
179 end
180
181 %DIRICHLET VALUES
182 if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
183 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
184 bool=0;return;
185 end
186
187 %DIRICHLET IF THICKNESS <= 0
188 if any(md.thickness<=0),
189 pos=find(md.thickness<=0);
190 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
191 disp(['model ' md.name ' has some grids with 0 thickness']);
192 bool=0; return;
193 end
194 end
195end
196
197%PROGNOSTIC
198if strcmp(md.analysis_type,'prognostic'),
199 %old testing
200 % if isempty(find(md.gridondirichlet_diag)),
201 % disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
202 % bool=0;return;
203 % end
204 % if isempty(find(md.gridondirichlet_prog)),
205 % disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
206 % bool=0;return;
207 % end
208
209 %VELOCITIES
210 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
211 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
212 bool=0; return;
213 end
214end
215
216%THERMAL TRANSIENT
217if strcmp(md.analysis_type,'thermal')
218 if strcmp(md.sub_analysis_type,'transient')
219
220 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
221 if isempty(md.temperature),
222 disp(['An initial temperature is needed for a transient thermal computation'])
223 bool=0;return;
224 end
225 if isstruct(md.temperature) | isstruct(md.melting) | isstruct(md.accumulation),
226 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
227 bool=0;return;
228 end
229 end
230end
231
232%THERMAL STEADY AND THERMAL TRANSIENT
233if strcmpi(md.analysis_type,'thermal'),
234
235 %EXTRUSION
236 if strcmp(md.type,'2d'),
237 disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
238 bool=0;return;
239 end
240
241 %VELOCITIES AND PRESSURE
242 if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy) | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
243 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
244 bool=0;return;
245 end
246 if (isempty(md.pressure) | isnan(md.pressure)),
247 disp(['pressure is required. Run ''diagnostic'' solution first!'])
248 bool=0;return;
249 end
250end
251
252%THERMAL TRANSIENT AND TRANSIENT
253if strcmp(md.analysis_type,'thermal'),
254
255 %DT and NDT
256 fields={'dt','ndt'};
257 for i=1:length(fields),
258 if any(md.(fields{i})<0),
259 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
260 bool=0; return;
261 end
262 end
263
264 %INITIAL TEMPERATURE
265 if isstruct(md.temperature),
266 disp(['The initial temperature should be empty or a list but not a structure'])
267 bool=0;return;
268 end
269end
270
271%PARAMETERS
272if strcmp(md.analysis_type,'parameters')
273
274 %PACKAGE
275 if ~strcmpi(package,'ice'),
276 disp('parameter solution only supported by package ice yet');
277 bool=0;return;
278 end
279
280 %OUTPUT
281 if ~iscell(md.parameteroutput)
282 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
283 bool=0; return;
284 end
285 for i=1:length(md.parameteroutput)
286 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
287 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
288 disp(['one of the parameteroutput is not supported yet']);
289 bool=0; return;
290 end
291 end
292 %VELOCITY
293 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
294 disp(['velocities are required!']);
295 bool=0;return;
296 end
297
298 %HUTTER
299 if any(md.elements_type(:,1)==hutterenum);
300 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
301 bool=0;return;
302 end
303end
304
305%CONTROL
306if strcmpi(md.analysis_type,'control'),
307
308 %CONTROL TYPE
309 if ~ischar(md.control_type),
310 disp('control_type should be a string');
311 bool=0;return;
312 end
313
314 %LENGTH CONTROL FIELDS
315 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
316 disp('maxiter, optscal and fit must have the length specified by nsteps')
317 bool=0;return;
318 end
319
320 %FIT
321 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
322 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
323 bool=0;return;
324 end
325
326 %OBSERVED VELOCITIES
327 fields={'vx_obs','vy_obs'};
328 for i=1:length(fields),
329 if any(length(md.(fields{i}))~=md.numberofgrids),
330 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
331 bool=0; return;
332 end
333 end
334end
335
336%QMU
337if strcmpi(md.analysis_type,'qmu'),
338 if ~strcmpi(md.cluster,'none'),
339 if md.waitonlock==0,
340 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
341 bool=0;return;
342 end
343 end
344end
345
346%MESH
347if strcmpi(md.analysis_type,'mesh'),
348 %this solution is a little special. It should come right after the md=model; operation. So a lot less checks!
349
350 bool=1;
351 return;
352end
353
354%MESH2GRID
355if strcmpi(md.analysis_type,'mesh2grid'),
356 if ~strcmpi(md.cluster,'none'),
357 disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
358 bool=0;return;
359 end
360end
361
362
363%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
364
365%CIELO
366if strcmpi(package,'cielo'),
367
368 %NAN VALUES
369 fields={'sparsity'};
370 for i=1:length(fields),
371 if ~isempty(md.(fields{i})),
372 if any(isnan(md.(fields{i}))),
373 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
374 bool=0; return;
375 end
376 end
377 end
378
379 %FIELD > 0
380 fields={'sparsity'};
381 for i=1:length(fields),
382 if ~isempty(md.(fields{i})),
383 if any(md.(fields{i})<0),
384 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
385 bool=0; return;
386 end
387 end
388 end
389
390 %FIELD ~= 0
391 fields={'sparsity'};
392 for i=1:length(fields),
393 if ~isempty(md.(fields{i})),
394 if any(md.(fields{i})==0),
395 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
396 bool=0; return;
397 end
398 end
399 end
400
401 %SPARSITY BETWEEN 0 AND 1
402 if ( (md.sparsity<=0) | (md.sparsity>1)),
403 disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
404 bool=0; return;
405 end
406
407 %RIFTS
408 if md.numrifts,
409 if ~strcmpi(md.cluster,'none')
410 if isempty(findstr('aijmumps',md.solverstring)),
411 disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
412 bool=0;return;
413 end
414 end
415 end
416
417 %CONNECTIVITY
418 if strcmpi(md.type,'2d'),
419 if md.connectivity<9,
420 disp('connectivity should be at least 9 for 2d models');
421 bool=0;return;
422 end
423 end
424 if strcmpi(md.type,'3d'),
425 if md.connectivity<24,
426 disp('connectivity should be at least 24 for 3d models');
427 bool=0;return;
428 end
429 end
430
431 %LOWMEM = 0 or 1
432 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
433 disp(['model ' md.name ' lowmem field should be 0 or 1']);
434 bool=0; return;
435 end
436
437%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
438
439 if ~strcmpi(md.cluster,'none'),
440
441 %NAN VALUES
442 fields={'time','np'};
443 for i=1:length(fields),
444 if ~isempty(md.(fields{i})),
445 if any(isnan(md.(fields{i}))),
446 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
447 bool=0; return;
448 end
449 end
450 end
451
452 %FIELD > 0
453 fields={'time','np'};
454 for i=1:length(fields),
455 if ~isempty(md.(fields{i})),
456 if any(md.(fields{i})<0),
457 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
458 bool=0; return;
459 end
460 end
461 end
462
463 %FIELD ~= 0
464 fields={'time','np'};
465 for i=1:length(fields),
466 if ~isempty(md.(fields{i})),
467 if any(md.(fields{i})==0),
468 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
469 bool=0; return;
470 end
471 end
472 end
473
474 end
475
476end
477
478%No problems, just return 1;
479bool=1;
480return;
Note: See TracBrowser for help on using the repository browser.