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

Last change on this file since 1781 was 1781, checked in by seroussi, 16 years ago

analysis_type are Enums now
File size: 13.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
67%NO NAN
68fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
69 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
70 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
71for i=1:length(fields),
72 if ~isempty(md.(fields{i})),
73 if any(isnan(md.(fields{i}))),
74 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
75 bool=0; return;
76 end
77 end
78end
79
80%FIELDS >= 0
81fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
82 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
83 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
84for i=1:length(fields),
85 if ~isempty(md.(fields{i})),
86 if any(md.(fields{i})<0),
87 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
88 bool=0; return;
89 end
90 end
91end
92if any(md.p<=0),
93 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
94 bool=0; return;
95end
96
97%FIELDS ~=0
98fields={'numberofelements','numberofgrids','elements','drag_type',...
99 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
100 'sparsity','deltaH','DeltaH','timeacc','timedec'};
101for i=1:length(fields),
102 if ~isempty(md.(fields{i})),
103 if any(md.(fields{i})==0),
104 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
105 bool=0; return;
106 end
107 end
108end
109
110%SIZE NUMBEROFELEMENTS
111fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
112for i=1:size(fields,2),
113 if (size(md.(fields{i}),1)~=md.numberofelements),
114 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
115 bool=0; return;
116 end
117end
118
119%SIZE NUMBEROFGRIDS
120fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
121for i=1:length(fields),
122 if length(md.(fields{i}))~=md.numberofgrids,
123 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
124 bool=0; return;
125 end
126end
127
128%THICKNESS = SURFACE - BED
129if any((md.thickness-md.surface+md.bed)>tolerance),
130 disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
131 bool=0; return;
132end
133
134%RIFTS
135if md.numrifts,
136 if ~strcmpi(md.type,'2d'),
137 disp(['Models with rifts are only supported in 2d for now!']);
138 bool=0;return;
139 end
140end
141if ~isstruct(md.rifts),
142 if ~isnan(md.rifts),
143 if ~isempty(find(md.segmentmarkers>=2)),
144 %We have segments with rift markers, but no rift structure!
145 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
146 bool=0; return;
147 end
148 end
149end
150
151%ARTIFICIAL DIFFUSIVITY
152if ~isscalar(md.artificial_diffusivity),
153 disp('artificial_diffusivity should be a scalar (1 or 0)');
154 bool=0;return;
155end
156
157%CLUSTER
158if ~strcmpi(package,'cielo') & ~strcmpi(md.cluster,'none')
159 disp(['parallel solution not supported by package ' package '. Use cluster=''none'' ']);
160 bool=0;return;
161end
162
163%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
164%QMU
165if md.qmu_analysis,
166 if md.qmu_params.evaluation_concurrency~=1,
167 disp(['concurrency should be set to 1 when running dakota in library mode']);
168 bool=0;return;
169 end
170 if ~isempty(md.part),
171 if numel(md.part)~=md.numberofgrids,
172 disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
173 bool=0;return;
174 end
175 if find(md.part)>=md.numberofgrids,
176 disp(['user supplied partition should be indexed from 0 (c-convention)']);
177 bool=0;return;
178 end
179 if md.npart~=md.numberofgrids,
180 disp(['user supplied partition should have same size as md.npart']);
181 bool=0;return;
182 end
183
184 end
185end
186
187%DIAGNOSTIC
188if md.analysis_type==DiagnosticAnalysisEnum,
189
190 %HUTTER ON ICESHELF WARNING
191 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
192 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
193 end
194
195 %SINGULAR
196 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
197 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
198 bool=0;return;
199 end
200
201 %DIRICHLET IF THICKNESS <= 0
202 if any(md.thickness<=0),
203 pos=find(md.thickness<=0);
204 if any(find(md.spcthickness(pos,1)==0)),
205 disp(['model ' md.name ' has some grids with 0 thickness']);
206 bool=0; return;
207 end
208 end
209end
210
211%PROGNOSTIC
212if md.analysis_type==PrognosticAnalysisEnum,
213
214 %VELOCITIES
215 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
216 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
217 bool=0; return;
218 end
219end
220
221%THERMAL STEADY AND THERMAL TRANSIENT
222if md.analysis_type==ThermalAnalysisEnum,
223
224 %EXTRUSION
225 if strcmp(md.type,'2d'),
226 disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
227 bool=0;return;
228 end
229
230 %VELOCITIES AND PRESSURE
231 if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
232 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
233 bool=0;return;
234 end
235 if (length(md.pressure)~=md.numberofgrids),
236 disp(['pressure is required. Run ''diagnostic'' solution first!'])
237 bool=0;return;
238 end
239end
240
241%THERMAL TRANSIENT
242if md.analysis_type==ThermalAnalysisEnum & md.sub_analysis_type==TransientAnalysisEnum,
243
244 %DT and NDT
245 fields={'dt','ndt'};
246 for i=1:length(fields),
247 if any(md.(fields{i})<0),
248 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
249 bool=0; return;
250 end
251 end
252
253 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
254 if (length(md.temperature)~=md.numberofgrids),
255 disp(['An initial temperature is needed for a transient thermal computation'])
256 bool=0;return;
257 end
258 if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
259 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
260 bool=0;return;
261 end
262end
263
264%PARAMETERS
265if md.analysis_type==ParametersAnalysisEnum,
266
267 %PACKAGE
268 if ~strcmpi(package,'ice'),
269 disp('parameter solution only supported by package ice yet');
270 bool=0;return;
271 end
272
273 %OUTPUT
274 if ~iscell(md.parameteroutput)
275 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
276 bool=0; return;
277 end
278 for i=1:length(md.parameteroutput)
279 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
280 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
281 disp(['one of the parameteroutput is not supported yet']);
282 bool=0; return;
283 end
284 end
285 %VELOCITY
286 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
287 disp(['velocities are required!']);
288 bool=0;return;
289 end
290
291 %HUTTER
292 if any(md.elements_type(:,1)==HutterFormulationEnum);
293 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
294 bool=0;return;
295 end
296end
297
298%CONTROL
299if md.analysis_type==ControlAnalysisEnum,
300
301 %CONTROL TYPE
302 if ~ischar(md.control_type),
303 disp('control_type should be a string');
304 bool=0;return;
305 end
306
307 %LENGTH CONTROL FIELDS
308 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
309 disp('maxiter, optscal and fit must have the length specified by nsteps')
310 bool=0;return;
311 end
312
313 %FIT
314 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
315 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
316 bool=0;return;
317 end
318
319 %OBSERVED VELOCITIES
320 fields={'vx_obs','vy_obs'};
321 for i=1:length(fields),
322 if any(length(md.(fields{i}))~=md.numberofgrids),
323 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
324 bool=0; return;
325 end
326 end
327
328 %SINGULAR
329 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
330 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
331 bool=0;return;
332 end
333
334 %DIRICHLET IF THICKNESS <= 0
335 if any(md.thickness<=0),
336 pos=find(md.thickness<=0);
337 if any(find(md.spcthickness(pos,1)==0)),
338 disp(['model ' md.name ' has some grids with 0 thickness']);
339 bool=0; return;
340 end
341 end
342end
343
344%QMU
345if strcmpi(md.analysis_type,'qmu'),
346 if ~strcmpi(md.cluster,'none'),
347 if md.waitonlock==0,
348 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
349 bool=0;return;
350 end
351 end
352end
353
354%MESH2GRID
355if md.analysis_type==Mesh2gridAnalysisEnum,
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 %CONNECTIVITY
408 if strcmpi(md.type,'2d'),
409 if md.connectivity<9,
410 disp('connectivity should be at least 9 for 2d models');
411 bool=0;return;
412 end
413 end
414 if strcmpi(md.type,'3d'),
415 if md.connectivity<24,
416 disp('connectivity should be at least 24 for 3d models');
417 bool=0;return;
418 end
419 end
420
421 %LOWMEM = 0 or 1
422 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
423 disp(['model ' md.name ' lowmem field should be 0 or 1']);
424 bool=0; return;
425 end
426
427%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
428
429 if ~strcmpi(md.cluster,'none'),
430
431 %NAN VALUES
432 fields={'time','np'};
433 for i=1:length(fields),
434 if ~isempty(md.(fields{i})),
435 if any(isnan(md.(fields{i}))),
436 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
437 bool=0; return;
438 end
439 end
440 end
441
442 %FIELD > 0
443 fields={'time','np'};
444 for i=1:length(fields),
445 if ~isempty(md.(fields{i})),
446 if any(md.(fields{i})<0),
447 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
448 bool=0; return;
449 end
450 end
451 end
452
453 %FIELD ~= 0
454 fields={'time','np'};
455 for i=1:length(fields),
456 if ~isempty(md.(fields{i})),
457 if any(md.(fields{i})==0),
458 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
459 bool=0; return;
460 end
461 end
462 end
463
464 end
465
466end
467
468%ICE and MACAYEAL
469if (strcmpi(package,'ice') | strcmpi(package,'macayeal')),
470
471 if (isnan(md.eps_rel) & isnan(md.eps_abs)),
472 disp(['At least eps_rel must be positive since eps_res not supported by package ' package '!']);
473 bool=0; return;
474 end
475
476end
Note: See TracBrowser for help on using the repository browser.