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

Last change on this file since 2991 was 2991, checked in by seroussi, 15 years ago

check size and not length in ismodelselfconsistent
File size: 15.1 KB
Line 
1function ismodelselfconsistent(md),
2%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
3%
4% Usage:
5% ismodelselfconsistent(md),
6
7%tolerance we use in litmus checks for the consistency of the model
8tolerance=10^-12;
9
10%check usage
11if nargin~=1,
12 help ismodelselfconsistent
13 error('ismodelselfconsistent error message: wrong usage');
14end
15
16%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSIENT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
17
18%check analysis
19if md.analysis_type==TransientAnalysisEnum,
20 if md.dt<=0,
21 error('model not consistent: field dt must be positive for a transient run')
22 end
23
24 %recursive call to ismodelselfconsistent
25 if strcmpi(md.type,'2d'),
26 analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum];
27 else
28 analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum ThermalAnalysisEnum];
29 end
30
31 for i=1:length(analysis),
32 md.analysis_type=analysis(i);
33 ismodelselfconsistent(md);
34 end
35end
36
37
38%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMON CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39
40%COUNTER
41if md.counter<3,
42 error(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize,setelementstype)!']);
43end
44
45%NAME
46if isempty(md.name),
47 error(['model is not correctly configured: missing name!']);
48end
49
50%ELEMENTSTYPE
51if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
52 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
53end
54if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum) + (md.elements_type(:,1)==PattynFormulationEnum)))
55 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
56end
57if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
58 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
59end
60if strcmpi(md.type,'2d'),
61 if (ismember(PattynFormulationEnum,md.elements_type(:,1)) | ismember(StokesFormulationEnum,md.elements_type(:,2))),
62 error(['model not consistent: for a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
63 end
64end
65if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
66 error(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
67end
68if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
69 error(['model not consistent: coupling MacAyeal/Pattyn not implemented yet']);
70end
71if (md.isstokes & md.analysis_type==TransientAnalysisEnum());
72 error(['model not consistent: Stokes transient not implemented yet']);
73end
74
75%ICEFRONT
76if strcmpi(md.type,'2d'),
77 fields={'pressureload'};
78 checksize(md,fields,[NaN 3]);
79elseif strcmpi(md.type,'3d'),
80 fields={'pressureload'};
81 checksize(md,fields,[NaN 5]);
82end
83
84%NO NAN
85fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
86 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
87 'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
88checknan(md,fields);
89
90%FIELDS >= 0
91fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
92 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
93 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
94checkgreater(md,fields,0);
95
96%FIELDS > 0
97fields={'numberofelements','numberofgrids','elements','drag_type','p',...
98 'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
99 'sparsity','deltaH','DeltaH','timeacc','timedec'};
100checkgreaterstrict(md,fields,0);
101
102%SIZE NUMBEROFELEMENTS
103fields={'p','q','elementoniceshelf','n','elementonbed'};
104checksize(md,fields,[md.numberofelements 1]);
105
106%SIZE NUMBEROFELEMENTS
107fields={'elements'};
108checklength(md,fields,md.numberofelements);
109
110%SIZE NUMBEROFGRIDS
111fields={'x','y','z','B','drag','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
112checksize(md,fields,[md.numberofgrids 1]);
113
114%SIZE 6
115fields={'spcvelocity'};
116checksize(md,fields,[md.numberofgrids 6]);
117
118%THICKNESS = SURFACE - BED
119if any((md.thickness-md.surface+md.bed)>tolerance),
120 error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
121end
122
123%RIFTS
124if md.numrifts,
125 if ~strcmpi(md.type,'2d'),
126 error(['model not consistent: models with rifts are only supported in 2d for now!']);
127 end
128 if ~isstruct(md.rifts),
129 error(['model not consistent: md.rifts should be a structure!']);
130 end
131 if ~isempty(find(md.segmentmarkers>=2)),
132 %We have segments with rift markers, but no rift structure!
133 error(['model not consistent: model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
134 end
135 %Check that rifts are filled with proper material
136 checkvalues(md,{'rifts.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
137else
138 if ~isnans(md.rifts),
139 error(['model not consistent: md.rifts shoud be NaN since md.numrifts is 0!']);
140 end
141end
142
143%ARTIFICIAL DIFFUSIVITY
144if ~ismember(md.artificial_diffusivity,[0 1]),
145 error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
146end
147
148%PARAMETEROUTPUT
149if md.numoutput~=length(md.parameteroutput),
150 error('model not consistent: numoutput should be the same size as parameteroutput');
151end
152
153%CONNECTIVITY
154if strcmpi(md.type,'2d'),
155 if md.connectivity<9,
156 error('model not consistent: connectivity should be at least 9 for 2d models');
157 end
158end
159if strcmpi(md.type,'3d'),
160 if md.connectivity<24,
161 error('model not consistent: connectivity should be at least 24 for 3d models');
162 end
163end
164
165%LOWMEM = 0 or 1
166if ((md.lowmem ~= 1) & (md.lowmem~=0)),
167 error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
168end
169
170%PARALLEL
171if ~strcmpi(md.cluster,'none'),
172
173 %NAN VALUES
174 fields={'time','np'};
175 checknan(md,fields);
176
177 %FIELD > 0
178 fields={'time','np'};
179 checkgreaterstrict(md,fields,0);
180
181end
182
183%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184
185%QMU
186if md.qmu_analysis,
187 if md.qmu_params.evaluation_concurrency~=1,
188 error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
189 end
190 if ~isempty(md.part),
191 if numel(md.part)~=md.numberofgrids,
192 error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
193 end
194 if find(md.part)>=md.numberofgrids,
195 error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
196 end
197 if md.npart~=md.numberofgrids,
198 error(['model not consistent: user supplied partition should have same size as md.npart']);
199 end
200 end
201 if md.eps_rel>10^-3,
202 error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
203 end
204end
205
206%DIAGNOSTIC
207if md.analysis_type==DiagnosticAnalysisEnum,
208
209 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
210 if strcmpi(md.type,'2d'),
211 if isempty(find(~md.spcvelocity(:,1:2))),
212 error(['model not consistent: model ' md.name ' is totally constrained, no need to solve!']);
213 end
214 end
215
216 %HUTTER ON ICESHELF WARNING
217 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
218 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
219 end
220
221 %SINGULAR
222 if ~any(sum(md.spcvelocity(:,1:2),2)==2),
223 error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
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 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
231 end
232 end
233
234 %INITIAL VELOCITY
235 if length(md.vx)==md.numberofgrids & length(md.vy)==md.numberofgrids,
236 fields={'vx','vy'};
237 checknan(md,fields);
238 end
239end
240
241%PROGNOSTIC
242if md.analysis_type==PrognosticAnalysisEnum,
243
244 %INITIAL VELOCITIES
245 fields={'vx','vy'};
246 checksize(md,fields,[md.numberofgrids 1]);
247 checknan(md,fields);
248
249 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
250 if strcmpi(md.type,'2d'),
251 if isempty(find(~md.spcthickness(:,1))),
252 error(['model not consistent: model ' md.name ' is totally constrained for prognostic, no need to solve!']);
253 end
254 end
255
256end
257
258%STEADYSTATE
259if md.analysis_type==SteadystateAnalysisEnum,
260
261 %NDT
262 if md.dt~=0,
263 error(['model not consistent: for a steadystate computation, dt must be zero.']);
264 end
265
266 %PRESSURE
267 if isnans(md.pressure),
268 error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
269 end
270
271 %eps:
272 if isnan(md.eps_rel),
273 error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
274 end
275
276 %dim:
277 if strcmpi(md.type,'2d'),
278 error(['model not consistent: for a steadystate computation, model needs to be 3d']);
279 end
280end
281
282%THERMAL STEADY AND THERMAL TRANSIENT
283if md.analysis_type==ThermalAnalysisEnum,
284
285 %EXTRUSION
286 if strcmp(md.type,'2d'),
287 error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
288 end
289
290 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
291 if strcmpi(md.type,'2d'),
292 if isempty(find(~md.spctemperature(:,1))),
293 error(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
294 end
295 end
296
297 %VELOCITIES AND PRESSURE
298 fields={'vx','vy','vz','pressure'};
299 checksize(md,fields,[md.numberofgrids 1]);
300 checknan(md,fields);
301
302end
303
304%THERMAL TRANSIENT
305if md.analysis_type==ThermalAnalysisEnum & md.dt~=0,
306
307 %DT and NDT
308 fields={'dt','ndt'};
309 checkgreaterstrict(md,fields,0);
310
311 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
312 fields={'temperature','accumulation','melting'};
313 checksize(md,fields,[md.numberofgrids 1]);
314 checknan(md,fields);
315
316 %INITIAL TEMPERATURE
317 fields={'temperature','spctemperature(:,2)','observed_temperature'};
318 checkgreater(md,fields,0)
319
320end
321
322%BALANCEDTHICKNESS
323if md.analysis_type==BalancedthicknessAnalysisEnum
324
325 %VELOCITIES MELTING AND ACCUMULATION
326 fields={'vx','vy','accumulation','melting'};
327 checksize(md,fields,[md.numberofgrids 1]);
328 checknan(md,fields);
329
330 %SPC
331 if any(md.spcthickness(find(md.gridonboundary))~=1),
332 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);
333 end
334end
335
336%BALANCEDVELOCITIES
337if md.analysis_type==BalancedvelocitiesAnalysisEnum
338
339 %VELOCITIES MELTING AND ACCUMULATION
340 fields={'vx','vy','accumulation','melting'};
341 checksize(md,fields,[md.numberofgrids 1]);
342 checknan(md,fields);
343
344 %SPC
345 if any(md.spcvelocity(find(md.gridonboundary),[1:2])~=1),
346 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcvelocity']);
347 end
348end
349
350%CONTROL
351if md.control_analysis,
352
353 %CONTROL TYPE
354 if ~ischar(md.control_type),
355 error('model not consistent: control_type should be a string');
356 end
357
358 %LENGTH CONTROL FIELDS
359 fields={'maxiter','optscal','fit','cm_jump'};
360 checksize(md,fields,[md.nsteps 1]);
361
362 %FIT
363 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
364 error('model not consistent: wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
365 end
366
367 %OBSERVED VELOCITIES
368 fields={'vx_obs','vy_obs'};
369 checksize(md,fields,[md.numberofgrids 1]);
370 checknan(md,fields);
371
372 %DIRICHLET IF THICKNESS <= 0
373 if any(md.thickness<=0),
374 pos=find(md.thickness<=0);
375 if any(find(md.spcthickness(pos,1)==0)),
376 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
377 end
378 end
379
380 %parameters
381 fields={'cm_noisedmp'};
382 checknan(md,fields);
383end
384
385%QMU
386if strcmpi(md.analysis_type,'qmu'),
387 if ~strcmpi(md.cluster,'none'),
388 if md.waitonlock==0,
389 error(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
390 end
391 end
392end
393
394end %end function
395
396%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
397
398function checklength(md,fields,fieldlength)
399 %CHECKSIZE - check length of a field
400 for i=1:length(fields),
401 if length(eval(['md.' fields{i}]))~=fieldlength,
402 error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
403 end
404 end
405end
406
407function checksize(md,fields,fieldsize)
408 %CHECKSIZE - check size of a field
409 for i=1:length(fields),
410 if isnan(fieldsize(1)),
411 if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
412 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
413 end
414 elseif isnan(fieldsize(2)),
415 if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
416 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
417 end
418 else
419 if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) | (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
420 error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
421 end
422 end
423 end
424end
425
426function checknan(md,fields)
427 %CHECKNAN - check nan values of a field
428 for i=1:length(fields),
429 if any(isnan(eval(['md.' fields{i}]))),
430 error(['model not consistent: NaN values in field ' fields{i}]);
431 end
432 end
433end
434
435function checkreal(md,fields)
436 %CHECKREAL - check real values of a field
437 for i=1:length(fields),
438 if any(eval(['~isreal(md.' fields{i} ')'])),
439 error(['model not consistent: complex values in field ' fields{i}]);
440 end
441 end
442end
443
444function checkgreaterstrict(md,fields,lowerbound)
445 %CHECKGREATERSTRICT - check values of a field
446 for i=1:length(fields),
447 if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
448 error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
449 end
450 end
451end
452
453function checkgreater(md,fields,lowerbound)
454 %CHECKGREATER - check values of a field
455 for i=1:length(fields),
456 if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
457 error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
458 end
459 end
460end
461
462function checklessstrict(md,fields,upperbound)
463 %CHECKLESSSTRICT - check values of a field
464 for i=1:length(fields),
465 if any(eval(['md.' fields{i} '>=' num2str(upperbound) ])),
466 error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
467 end
468 end
469end
470
471function checkless(md,fields,upperbound)
472 %CHECKLESS - check values of a field
473 for i=1:length(fields),
474 if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
475 error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
476 end
477 end
478end
479
480function checkvalues(md,fields,values)
481 %CHECKVALUE - check that a field has a certain value
482 for i=1:length(fields),
483 if eval(['~ismember( md.' fields{i} ',values)']),
484 error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
485 end
486 end
487end
Note: See TracBrowser for help on using the repository browser.