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

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

neumann boundary conditions gone, corrected displaybc.
added check for diagnostic horiz, when all grids are spc'd.
File size: 14.5 KB
Line 
1function bool=ismodelselfconsistent(md),
2%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
3%
4% Usage:
5% bool=ismodelselfconsistent(md),
6
7%tolerance we use in litmus tests for the consistency of the model
8tolerance=10^-12;
9if (nargin~=1 )
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
60if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
61 disp(['no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
62 bool=0;return;
63end
64if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
65 disp(['coupling MacAyeal/Pattyn not implemented yet']);
66 bool=0;return;
67end
68if (md.isstokes & md.analysis_type==TransientAnalysisEnum());
69 disp(['Stokes transient not implemented yet']);
70 bool=0;return;
71end
72
73%ICEFRONT
74if strcmpi(md.type,'2d'),
75 if size(md.pressureload,2)~=3 ,
76 disp(['pressureload in ' md.name ' shoud have three columns']);
77 bool=0;return;
78 end
79elseif strcmpi(md.type,'3d'),
80 if size(md.pressureload,2)~=5 ,
81 disp(['pressureload in ' md.name ' shoud have five columns']);
82 bool=0;return;
83 end
84end
85
86%NO NAN
87fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
88 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
89 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
90for i=1:length(fields),
91 if ~isempty(md.(fields{i})),
92 if any(isnan(md.(fields{i}))),
93 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
94 bool=0; return;
95 end
96 end
97end
98
99%FIELDS >= 0
100fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
101 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
102 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
103for i=1:length(fields),
104 if ~isempty(md.(fields{i})),
105 if any(md.(fields{i})<0),
106 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
107 bool=0; return;
108 end
109 end
110end
111if any(md.p<=0),
112 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
113 bool=0; return;
114end
115
116%FIELDS ~=0
117fields={'numberofelements','numberofgrids','elements','drag_type',...
118 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
119 'sparsity','deltaH','DeltaH','timeacc','timedec'};
120for i=1:length(fields),
121 if ~isempty(md.(fields{i})),
122 if any(md.(fields{i})==0),
123 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
124 bool=0; return;
125 end
126 end
127end
128
129%SIZE NUMBEROFELEMENTS
130fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
131for i=1:size(fields,2),
132 if (size(md.(fields{i}),1)~=md.numberofelements),
133 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofelements) '!']);
134 bool=0; return;
135 end
136end
137
138%SIZE NUMBEROFGRIDS
139fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
140for i=1:length(fields),
141 if length(md.(fields{i}))~=md.numberofgrids,
142 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
143 bool=0; return;
144 end
145end
146
147%THICKNESS = SURFACE - BED
148if any((md.thickness-md.surface+md.bed)>tolerance),
149 disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
150 bool=0; return;
151end
152
153%RIFTS
154if md.numrifts,
155 if ~strcmpi(md.type,'2d'),
156 disp(['Models with rifts are only supported in 2d for now!']);
157 bool=0;return;
158 end
159end
160if ~isstruct(md.rifts),
161 if ~isnan(md.rifts),
162 if ~isempty(find(md.segmentmarkers>=2)),
163 %We have segments with rift markers, but no rift structure!
164 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
165 bool=0; return;
166 end
167 end
168end
169
170%ARTIFICIAL DIFFUSIVITY
171if ~isscalar(md.artificial_diffusivity),
172 disp('artificial_diffusivity should be a scalar (1 or 0)');
173 bool=0;return;
174end
175
176%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
177%QMU
178if md.qmu_analysis,
179 if md.qmu_params.evaluation_concurrency~=1,
180 disp(['concurrency should be set to 1 when running dakota in library mode']);
181 bool=0;return;
182 end
183 if ~isempty(md.part),
184 if numel(md.part)~=md.numberofgrids,
185 disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
186 bool=0;return;
187 end
188 if find(md.part)>=md.numberofgrids,
189 disp(['user supplied partition should be indexed from 0 (c-convention)']);
190 bool=0;return;
191 end
192 if md.npart~=md.numberofgrids,
193 disp(['user supplied partition should have same size as md.npart']);
194 bool=0;return;
195 end
196
197 end
198 if md.eps_rel>10^-3,
199 disp(['for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
200 bool=0;return;
201 end
202end
203
204%DIAGNOSTIC
205if md.analysis_type==DiagnosticAnalysisEnum,
206
207 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
208 if strcmpi(md.type,'2d'),
209 if isempty(find(~md.spcvelocity(:,1:2))),
210 disp(['model ' md.name ' is totally constrained, no need to solve!']);
211 bool=0;return;
212 end
213 end
214
215 %HUTTER ON ICESHELF WARNING
216 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
217 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
218 end
219
220 %SINGULAR
221 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
222 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
223 bool=0;return;
224 end
225
226 %DIRICHLET IF THICKNESS <= 0
227 if any(md.thickness<=0),
228 pos=find(md.thickness<=0);
229 if any(find(md.spcthickness(pos,1)==0)),
230 disp(['model ' md.name ' has some grids with 0 thickness']);
231 bool=0; return;
232 end
233 end
234end
235
236%PROGNOSTIC
237if md.analysis_type==PrognosticAnalysisEnum,
238
239 %VELOCITIES
240 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
241 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
242 bool=0; return;
243 end
244end
245
246%STEADYSTATE
247if md.analysis_type==SteadystateAnalysisEnum,
248 %NDT
249 if md.dt~=0,
250 disp(['For a steadystate computation, dt must be zero.']);
251 bool=0;return;
252 end
253 %PRESSURE
254 if isnans(md.pressure),
255 disp(['For a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
256 bool=0;return;
257 end
258
259 %eps:
260 if isnan(md.eps_rel),
261 disp(['For a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
262 bool=0;return;
263 end
264
265 %dim:
266 if strcmpi(md.type,'2d'),
267 disp(['For a steadystate computation, model needs to be 3d']);
268 bool=0;return;
269 end
270end
271
272%THERMAL STEADY AND THERMAL TRANSIENT
273if md.analysis_type==ThermalAnalysisEnum,
274
275 %EXTRUSION
276 if strcmp(md.type,'2d'),
277 disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
278 bool=0;return;
279 end
280
281 %VELOCITIES AND PRESSURE
282 if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
283 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
284 bool=0;return;
285 end
286 if (length(md.pressure)~=md.numberofgrids),
287 disp(['pressure is required. Run ''diagnostic'' solution first!'])
288 bool=0;return;
289 end
290end
291
292%THERMAL TRANSIENT
293if md.analysis_type==ThermalAnalysisEnum & md.dt~=0
294
295 %DT and NDT
296 fields={'dt','ndt'};
297 for i=1:length(fields),
298 if any(md.(fields{i})<0),
299 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
300 bool=0; return;
301 end
302 end
303
304 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
305 if (length(md.temperature)~=md.numberofgrids),
306 disp(['An initial temperature is needed for a transient thermal computation'])
307 bool=0;return;
308 end
309 if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
310 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
311 bool=0;return;
312 end
313end
314
315%PARAMETERS
316if md.analysis_type==ParametersAnalysisEnum,
317
318 %OUTPUT
319 if ~iscell(md.parameteroutput)
320 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
321 bool=0; return;
322 end
323 for i=1:length(md.parameteroutput)
324 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
325 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
326 disp(['one of the parameteroutput is not supported yet']);
327 bool=0; return;
328 end
329 end
330 %VELOCITY
331 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
332 disp(['velocities are required!']);
333 bool=0;return;
334 end
335
336 %HUTTER
337 if any(md.elements_type(:,1)==HutterFormulationEnum);
338 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
339 bool=0;return;
340 end
341end
342
343%CONTROL
344if md.control_analysis,
345
346 %CONTROL TYPE
347 if ~ischar(md.control_type),
348 disp('control_type should be a string');
349 bool=0;return;
350 end
351
352 %LENGTH CONTROL FIELDS
353 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
354 disp('maxiter, optscal and fit must have the length specified by nsteps')
355 bool=0;return;
356 end
357
358 %FIT
359 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
360 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
361 bool=0;return;
362 end
363
364 %OBSERVED VELOCITIES
365 fields={'vx_obs','vy_obs'};
366 for i=1:length(fields),
367 if any(length(md.(fields{i}))~=md.numberofgrids),
368 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
369 bool=0; return;
370 end
371 end
372
373 %SINGULAR
374 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
375 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
376 bool=0;return;
377 end
378
379 %DIRICHLET IF THICKNESS <= 0
380 if any(md.thickness<=0),
381 pos=find(md.thickness<=0);
382 if any(find(md.spcthickness(pos,1)==0)),
383 disp(['model ' md.name ' has some grids with 0 thickness']);
384 bool=0; return;
385 end
386 end
387end
388
389%QMU
390if strcmpi(md.analysis_type,'qmu'),
391 if ~strcmpi(md.cluster,'none'),
392 if md.waitonlock==0,
393 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
394 bool=0;return;
395 end
396 end
397end
398
399%MESH2GRID
400if md.analysis_type==Mesh2gridAnalysisEnum,
401 if ~strcmpi(md.cluster,'none'),
402 disp(['model is not correctly configured: mesh2grid not supported in parallel yet!']);
403 bool=0;return;
404 end
405end
406
407
408%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
409
410%NAN VALUES
411fields={'sparsity'};
412for 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
419end
420
421%FIELD > 0
422fields={'sparsity'};
423for 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
430end
431
432%FIELD ~= 0
433fields={'sparsity'};
434for 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
441end
442
443%SPARSITY BETWEEN 0 AND 1
444if ( (md.sparsity<=0) | (md.sparsity>1)),
445 disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
446 bool=0; return;
447end
448
449%CONNECTIVITY
450if 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
455end
456if 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
461end
462
463%LOWMEM = 0 or 1
464if ((md.lowmem ~= 1) & (md.lowmem~=0)),
465 disp(['model ' md.name ' lowmem field should be 0 or 1']);
466 bool=0; return;
467end
468
469%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
470
471if ~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
506end
Note: See TracBrowser for help on using the repository browser.