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

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

Changed antzoom to basinzoom.
marshall: added parameteroutput in marshalling.
extrude: added extrusion of balanced velocities.
ismodelselfconsistent: checked that numoutput and parameteroutput are same size.
File size: 14.4 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%PARAMETEROUTPUT
177if md.numoutput~=length(md.parameteroutput),
178 disp('numoutput should be the same size as parameteroutput');
179 bool=0;return;
180end
181
182%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
183%QMU
184if md.qmu_analysis,
185 if md.qmu_params.evaluation_concurrency~=1,
186 disp(['concurrency should be set to 1 when running dakota in library mode']);
187 bool=0;return;
188 end
189 if ~isempty(md.part),
190 if numel(md.part)~=md.numberofgrids,
191 disp(['user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
192 bool=0;return;
193 end
194 if find(md.part)>=md.numberofgrids,
195 disp(['user supplied partition should be indexed from 0 (c-convention)']);
196 bool=0;return;
197 end
198 if md.npart~=md.numberofgrids,
199 disp(['user supplied partition should have same size as md.npart']);
200 bool=0;return;
201 end
202
203 end
204 if md.eps_rel>10^-3,
205 disp(['for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
206 bool=0;return;
207 end
208end
209
210%DIAGNOSTIC
211if md.analysis_type==DiagnosticAnalysisEnum,
212
213 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
214 if strcmpi(md.type,'2d'),
215 if isempty(find(~md.spcvelocity(:,1:2))),
216 disp(['model ' md.name ' is totally constrained, no need to solve!']);
217 bool=0;return;
218 end
219 end
220
221 %HUTTER ON ICESHELF WARNING
222 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
223 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
224 end
225
226 %SINGULAR
227 if ~any(sum(md.spcvelocity(:,1:2),2)==2),
228 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
229 bool=0;return;
230 end
231
232 %DIRICHLET IF THICKNESS <= 0
233 if any(md.thickness<=0),
234 pos=find(md.thickness<=0);
235 if any(find(md.spcthickness(pos,1)==0)),
236 disp(['model ' md.name ' has some grids with 0 thickness']);
237 bool=0; return;
238 end
239 end
240end
241
242%PROGNOSTIC
243if md.analysis_type==PrognosticAnalysisEnum,
244
245 %VELOCITIES
246 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
247 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
248 bool=0; return;
249 end
250end
251
252%STEADYSTATE
253if md.analysis_type==SteadystateAnalysisEnum,
254 %NDT
255 if md.dt~=0,
256 disp(['For a steadystate computation, dt must be zero.']);
257 bool=0;return;
258 end
259 %PRESSURE
260 if isnans(md.pressure),
261 disp(['For a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
262 bool=0;return;
263 end
264
265 %eps:
266 if isnan(md.eps_rel),
267 disp(['For a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
268 bool=0;return;
269 end
270
271 %dim:
272 if strcmpi(md.type,'2d'),
273 disp(['For a steadystate computation, model needs to be 3d']);
274 bool=0;return;
275 end
276end
277
278%THERMAL STEADY AND THERMAL TRANSIENT
279if md.analysis_type==ThermalAnalysisEnum,
280
281 %EXTRUSION
282 if strcmp(md.type,'2d'),
283 disp(['For a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
284 bool=0;return;
285 end
286
287 %VELOCITIES AND PRESSURE
288 if (length(md.vx)~=md.numberofgrids | length(md.vy)~=md.numberofgrids | length(md.vz)~=md.numberofgrids),
289 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
290 bool=0;return;
291 end
292 if (length(md.pressure)~=md.numberofgrids),
293 disp(['pressure is required. Run ''diagnostic'' solution first!'])
294 bool=0;return;
295 end
296end
297
298%THERMAL TRANSIENT
299if md.analysis_type==ThermalAnalysisEnum & md.dt~=0
300
301 %DT and NDT
302 fields={'dt','ndt'};
303 for i=1:length(fields),
304 if any(md.(fields{i})<0),
305 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
306 bool=0; return;
307 end
308 end
309
310 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
311 if (length(md.temperature)~=md.numberofgrids),
312 disp(['An initial temperature is needed for a transient thermal computation'])
313 bool=0;return;
314 end
315 if (length(md.temperature)~=md.numberofgrids | length(md.accumulation)~=md.numberofgrids | length(md.melting)~=md.numberofgrids),
316 disp(['The initial temperature, melting or accumulation should be a list and not a structure'])
317 bool=0;return;
318 end
319end
320
321%PARAMETERS
322if md.analysis_type==ParametersAnalysisEnum,
323
324 %OUTPUT
325 if ~iscell(md.parameteroutput)
326 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
327 bool=0; return;
328 end
329 for i=1:length(md.parameteroutput)
330 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
331 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
332 disp(['one of the parameteroutput is not supported yet']);
333 bool=0; return;
334 end
335 end
336 %VELOCITY
337 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
338 disp(['velocities are required!']);
339 bool=0;return;
340 end
341
342 %HUTTER
343 if any(md.elements_type(:,1)==HutterFormulationEnum);
344 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
345 bool=0;return;
346 end
347end
348
349%CONTROL
350if md.control_analysis,
351
352 %CONTROL TYPE
353 if ~ischar(md.control_type),
354 disp('control_type should be a string');
355 bool=0;return;
356 end
357
358 %LENGTH CONTROL FIELDS
359 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps | length(md.cm_jump)~=md.nsteps)
360 disp('maxiter, optscal, fit and cm_jump must have the length specified by nsteps')
361 bool=0;return;
362 end
363
364 %FIT
365 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
366 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
367 bool=0;return;
368 end
369
370 %OBSERVED VELOCITIES
371 fields={'vx_obs','vy_obs'};
372 for i=1:length(fields),
373 if any(length(md.(fields{i}))~=md.numberofgrids),
374 disp(['model ' md.name ' field ' fields{i} ' should be of size ' num2str(md.numberofgrids) '!']);
375 bool=0; return;
376 end
377 end
378
379 %SINGULAR
380 if ~any(sum(md.spcvelocity(:,1:3),2)==3),
381 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
382 bool=0;return;
383 end
384
385 %DIRICHLET IF THICKNESS <= 0
386 if any(md.thickness<=0),
387 pos=find(md.thickness<=0);
388 if any(find(md.spcthickness(pos,1)==0)),
389 disp(['model ' md.name ' has some grids with 0 thickness']);
390 bool=0; return;
391 end
392 end
393end
394
395%QMU
396if strcmpi(md.analysis_type,'qmu'),
397 if ~strcmpi(md.cluster,'none'),
398 if md.waitonlock==0,
399 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
400 bool=0;return;
401 end
402 end
403end
404
405%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PACKAGE CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
406
407%NAN VALUES
408fields={'sparsity'};
409for i=1:length(fields),
410 if ~isempty(md.(fields{i})),
411 if any(isnan(md.(fields{i}))),
412 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
413 bool=0; return;
414 end
415 end
416end
417
418%FIELD > 0
419fields={'sparsity'};
420for i=1:length(fields),
421 if ~isempty(md.(fields{i})),
422 if any(md.(fields{i})<0),
423 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
424 bool=0; return;
425 end
426 end
427end
428
429%FIELD ~= 0
430fields={'sparsity'};
431for i=1:length(fields),
432 if ~isempty(md.(fields{i})),
433 if any(md.(fields{i})==0),
434 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
435 bool=0; return;
436 end
437 end
438end
439
440%SPARSITY BETWEEN 0 AND 1
441if ( (md.sparsity<=0) | (md.sparsity>1)),
442 disp(['model ' md.name ' sparsity should be inside the [0 1] range']);
443 bool=0; return;
444end
445
446%CONNECTIVITY
447if strcmpi(md.type,'2d'),
448 if md.connectivity<9,
449 disp('connectivity should be at least 9 for 2d models');
450 bool=0;return;
451 end
452end
453if strcmpi(md.type,'3d'),
454 if md.connectivity<24,
455 disp('connectivity should be at least 24 for 3d models');
456 bool=0;return;
457 end
458end
459
460%LOWMEM = 0 or 1
461if ((md.lowmem ~= 1) & (md.lowmem~=0)),
462 disp(['model ' md.name ' lowmem field should be 0 or 1']);
463 bool=0; return;
464end
465
466%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARALLEL CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
467
468if ~strcmpi(md.cluster,'none'),
469
470 %NAN VALUES
471 fields={'time','np'};
472 for i=1:length(fields),
473 if ~isempty(md.(fields{i})),
474 if any(isnan(md.(fields{i}))),
475 disp(['model ' md.name ' has an NaN value in field ' fields{i} '!']);
476 bool=0; return;
477 end
478 end
479 end
480
481 %FIELD > 0
482 fields={'time','np'};
483 for i=1:length(fields),
484 if ~isempty(md.(fields{i})),
485 if any(md.(fields{i})<0),
486 disp(['model ' md.name ' has a <0 value in field ' fields{i} '!']);
487 bool=0; return;
488 end
489 end
490 end
491
492 %FIELD ~= 0
493 fields={'time','np'};
494 for i=1:length(fields),
495 if ~isempty(md.(fields{i})),
496 if any(md.(fields{i})==0),
497 disp(['model ' md.name ' has a =0 value in field ' fields{i} '!']);
498 bool=0; return;
499 end
500 end
501 end
502
503end
Note: See TracBrowser for help on using the repository browser.