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

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

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