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

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

ismodelselfconsistent: for qmu analysis to be efficient, eps_rel should be at least <10-3, better 10-15. yeeh ahh!
updated dakota nightly run with that in mind.
File size: 14.6 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%initialize output
15bool=1;
16
17%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMON CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18
19%COUNTER
20if md.counter<3,
21 disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
22 bool=0;return;
23end
24
25%NAME
26if isempty(md.name),
27 disp(['model is not correctly configured: missing name!']);
28 bool=0;return;
29end
30
31%MESH
32if md.numberofelements<=0,
33 disp(['model ' md.name ' does not have any elements!']);
34 bool=0; return;
35end
36if md.numberofgrids<=0,
37 disp(['model ' md.name ' does not have any grids!']);
38 bool=0; return;
39end
40
41%ELEMENTSTYPE
42if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
43 disp(['Types of elements have not been set properly, run setelementstype first'])
44 bool=0;return;
45end
46if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum) + (md.elements_type(:,1)==PattynFormulationEnum)))
47 disp(['Types of elements have not been set properly, run setelementstype first'])
48 bool=0;return;
49end
50if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
51 disp(['Types of elements have not been set properly, run setelementstype first'])
52 bool=0;return;
53end
54if strcmpi(md.type,'2d'),
55 if (ismember(PattynFormulationEnum,md.elements_type(:,1)) | ismember(StokesFormulationEnum,md.elements_type(:,2))),
56 disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
57 bool=0;return;
58 end
59end
60
61if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
62 disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
63 bool=0;return;
64end
65
66%ICEFRONT
67if strcmpi(md.type,'2d'),
68 if size(md.pressureload,2)~=3 ,
69 disp(['pressureload in ' md.name ' shoud have three columns']);
70 bool=0;return;
71 end
72elseif strcmpi(md.type,'3d'),
73 if size(md.pressureload,2)~=5 ,
74 disp(['pressureload in ' md.name ' shoud have five columns']);
75 bool=0;return;
76 end
77end
78
79%NO NAN
80fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
81 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
82 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
83for i=1:length(fields),
84 if ~isempty(md.(fields{i})),
85 if any(isnan(md.(fields{i}))),
86 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
87 bool=0; return;
88 end
89 end
90end
91
92%FIELDS >= 0
93fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
94 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
95 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
96for i=1:length(fields),
97 if ~isempty(md.(fields{i})),
98 if any(md.(fields{i})<0),
99 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
100 bool=0; return;
101 end
102 end
103end
104if any(md.p<=0),
105 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
106 bool=0; return;
107end
108
109%FIELDS ~=0
110fields={'numberofelements','numberofgrids','elements','drag_type',...
111 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
112 'sparsity','deltaH','DeltaH','timeacc','timedec'};
113for i=1:length(fields),
114 if ~isempty(md.(fields{i})),
115 if any(md.(fields{i})==0),
116 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
117 bool=0; return;
118 end
119 end
120end
121
122%SIZE NUMBEROFELEMENTS
123fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
124for i=1:size(fields,2),
125 if (size(md.(fields{i}),1)~=md.numberofelements),
126 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
127 bool=0; return;
128 end
129end
130
131%SIZE NUMBEROFGRIDS
132fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
133for i=1:length(fields),
134 if length(md.(fields{i}))~=md.numberofgrids,
135 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
136 bool=0; return;
137 end
138end
139
140%THICKNESS = SURFACE - BED
141if any((md.thickness-md.surface+md.bed)>tolerance),
142 disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
143 bool=0; return;
144end
145
146%RIFTS
147if md.numrifts,
148 if ~strcmpi(md.type,'2d'),
149 disp(['Models with rifts are only supported in 2d for now!']);
150 bool=0;return;
151 end
152end
153if ~isstruct(md.rifts),
154 if ~isnan(md.rifts),
155 if ~isempty(find(md.segmentmarkers>=2)),
156 %We have segments with rift markers, but no rift structure!
157 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
158 bool=0; return;
159 end
160 end
161end
162
163%ARTIFICIAL DIFFUSIVITY
164if ~isscalar(md.artificial_diffusivity),
165 disp('artificial_diffusivity should be a scalar (1 or 0)');
166 bool=0;return;
167end
168
169%CLUSTER
170if ~strcmpi(package,'cielo') & ~strcmpi(md.cluster,'none')
171 disp(['parallel solution not supported by package ' package '. Use cluster=''none'' ']);
172 bool=0;return;
173end
174
175%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
176%QMU
177if md.qmu_analysis,
178 if md.qmu_params.evaluation_concurrency~=1,
179 disp(['concurrency should be set to 1 when running dakota in library mode']);
180 bool=0;return;
181 end
182 if ~isempty(md.part),
183 if numel(md.part)~=md.numberofgrids,
184 disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
185 bool=0;return;
186 end
187 if find(md.part)>=md.numberofgrids,
188 disp(['user supplied partition should be indexed from 0 (c-convention)']);
189 bool=0;return;
190 end
191 if md.npart~=md.numberofgrids,
192 disp(['user supplied partition should have same size as md.npart']);
193 bool=0;return;
194 end
195
196 end
197 if md.eps_rel>10^-3,
198 disp(['for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
199 bool=0;return;
200 end
201end
202
203%DIAGNOSTIC
204if md.analysis_type==DiagnosticAnalysisEnum,
205
206 %HUTTER ON ICESHELF WARNING
207 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
208 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
209 end
210
211 %SINGULAR
212 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
213 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
214 bool=0;return;
215 end
216
217 %DIRICHLET IF THICKNESS <= 0
218 if any(md.thickness<=0),
219 pos=find(md.thickness<=0);
220 if any(find(md.spcthickness(pos,1)==0)),
221 disp(['model ' md.name ' has some grids with 0 thickness']);
222 bool=0; return;
223 end
224 end
225end
226
227%PROGNOSTIC
228if md.analysis_type==PrognosticAnalysisEnum,
229
230 %VELOCITIES
231 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
232 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
233 bool=0; return;
234 end
235end
236
237%STEADYSTATE
238if md.analysis_type==SteadystateAnalysisEnum,
239 %NDT
240 if md.dt~=0,
241 disp(['For a steadystate computation, dt must be zero.']);
242 bool=0;return;
243 end
244 %PRESSURE
245 if isnans(md.pressure),
246 disp(['For a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
247 bool=0;return;
248 end
249
250 %eps:
251 if isnan(md.eps_rel),
252 disp(['For a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
253 bool=0;return;
254 end
255
256 %dim:
257 if strcmpi(md.type,'2d'),
258 disp(['For a steadystate computation, model needs to be 3d']);
259 bool=0;return;
260 end
261end
262
263%THERMAL STEADY AND THERMAL TRANSIENT
264if md.analysis_type==ThermalAnalysisEnum,
265
266 %EXTRUSION
267 if strcmp(md.type,'2d'),
268 disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
269 bool=0;return;
270 end
271
272 %VELOCITIES AND PRESSURE
273 if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
274 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
275 bool=0;return;
276 end
277 if (length(md.pressure)~=md.numberofgrids),
278 disp(['pressure is required. Run ''diagnostic'' solution first!'])
279 bool=0;return;
280 end
281end
282
283%THERMAL TRANSIENT
284if md.analysis_type==ThermalAnalysisEnum & md.dt~=0
285
286 %DT and NDT
287 fields={'dt','ndt'};
288 for i=1:length(fields),
289 if any(md.(fields{i})<0),
290 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
291 bool=0; return;
292 end
293 end
294
295 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
296 if (length(md.temperature)~=md.numberofgrids),
297 disp(['An initial temperature is needed for a transient thermal computation'])
298 bool=0;return;
299 end
300 if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
301 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
302 bool=0;return;
303 end
304end
305
306%PARAMETERS
307if md.analysis_type==ParametersAnalysisEnum,
308
309 %PACKAGE
310 if ~strcmpi(package,'ice'),
311 disp('parameter solution only supported by package ice yet');
312 bool=0;return;
313 end
314
315 %OUTPUT
316 if ~iscell(md.parameteroutput)
317 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
318 bool=0; return;
319 end
320 for i=1:length(md.parameteroutput)
321 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
322 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
323 disp(['one of the parameteroutput is not supported yet']);
324 bool=0; return;
325 end
326 end
327 %VELOCITY
328 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
329 disp(['velocities are required!']);
330 bool=0;return;
331 end
332
333 %HUTTER
334 if any(md.elements_type(:,1)==HutterFormulationEnum);
335 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
336 bool=0;return;
337 end
338end
339
340%CONTROL
341if md.control_analysis,
342
343 %CONTROL TYPE
344 if ~ischar(md.control_type),
345 disp('control_type should be a string');
346 bool=0;return;
347 end
348
349 %LENGTH CONTROL FIELDS
350 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
351 disp('maxiter, optscal and fit must have the length specified by nsteps')
352 bool=0;return;
353 end
354
355 %FIT
356 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
357 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
358 bool=0;return;
359 end
360
361 %OBSERVED VELOCITIES
362 fields={'vx_obs','vy_obs'};
363 for i=1:length(fields),
364 if any(length(md.(fields{i}))~=md.numberofgrids),
365 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
366 bool=0; return;
367 end
368 end
369
370 %SINGULAR
371 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
372 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
373 bool=0;return;
374 end
375
376 %DIRICHLET IF THICKNESS <= 0
377 if any(md.thickness<=0),
378 pos=find(md.thickness<=0);
379 if any(find(md.spcthickness(pos,1)==0)),
380 disp(['model ' md.name ' has some grids with 0 thickness']);
381 bool=0; return;
382 end
383 end
384end
385
386%QMU
387if strcmpi(md.analysis_type,'qmu'),
388 if ~strcmpi(md.cluster,'none'),
389 if md.waitonlock==0,
390 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
391 bool=0;return;
392 end
393 end
394end
395
396%MESH2GRID
397if md.analysis_type==Mesh2gridAnalysisEnum,
398 if ~strcmpi(md.cluster,'none'),
399 disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
400 bool=0;return;
401 end
402end
403
404
405%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406
407%CIELO
408if strcmpi(package,'cielo'),
409
410 %NAN VALUES
411 fields={'sparsity'};
412 for i=1:length(fields),
413 if ~isempty(md.(fields{i})),
414 if any(isnan(md.(fields{i}))),
415 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
416 bool=0; return;
417 end
418 end
419 end
420
421 %FIELD > 0
422 fields={'sparsity'};
423 for i=1:length(fields),
424 if ~isempty(md.(fields{i})),
425 if any(md.(fields{i})<0),
426 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
427 bool=0; return;
428 end
429 end
430 end
431
432 %FIELD ~= 0
433 fields={'sparsity'};
434 for i=1:length(fields),
435 if ~isempty(md.(fields{i})),
436 if any(md.(fields{i})==0),
437 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
438 bool=0; return;
439 end
440 end
441 end
442
443 %SPARSITY BETWEEN 0 AND 1
444 if ( (md.sparsity<=0) | (md.sparsity>1)),
445 disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
446 bool=0; return;
447 end
448
449 %CONNECTIVITY
450 if strcmpi(md.type,'2d'),
451 if md.connectivity<9,
452 disp('connectivity should be at least 9 for 2d models');
453 bool=0;return;
454 end
455 end
456 if strcmpi(md.type,'3d'),
457 if md.connectivity<24,
458 disp('connectivity should be at least 24 for 3d models');
459 bool=0;return;
460 end
461 end
462
463 %LOWMEM = 0 or 1
464 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
465 disp(['model ' md.name ' lowmem field should be 0 or 1']);
466 bool=0; return;
467 end
468
469%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470
471 if ~strcmpi(md.cluster,'none'),
472
473 %NAN VALUES
474 fields={'time','np'};
475 for i=1:length(fields),
476 if ~isempty(md.(fields{i})),
477 if any(isnan(md.(fields{i}))),
478 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
479 bool=0; return;
480 end
481 end
482 end
483
484 %FIELD > 0
485 fields={'time','np'};
486 for i=1:length(fields),
487 if ~isempty(md.(fields{i})),
488 if any(md.(fields{i})<0),
489 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
490 bool=0; return;
491 end
492 end
493 end
494
495 %FIELD ~= 0
496 fields={'time','np'};
497 for i=1:length(fields),
498 if ~isempty(md.(fields{i})),
499 if any(md.(fields{i})==0),
500 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
501 bool=0; return;
502 end
503 end
504 end
505
506 end
507
508end
509
510%ICE and MACAYEAL
511if (strcmpi(package,'ice') | strcmpi(package,'macayeal')),
512
513 if (isnan(md.eps_rel) & isnan(md.eps_abs)),
514 disp(['At least eps_rel must be positive since eps_res not supported by package ' package '!']);
515 bool=0; return;
516 end
517
518end
Note: See TracBrowser for help on using the repository browser.