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

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

fields name
File size: 16.9 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 4]);
79elseif strcmpi(md.type,'3d'),
80 fields={'pressureload'};
81 checksize(md,fields,[NaN 6]);
82else
83 error('type should be either ''2d'' or ''3d''');
84end
85
86%ELEMENTS
87fields={'elements'};
88if strcmpi(md.type,'2d'),
89 checksize(md,fields,[md.numberofelements 3]);
90else
91 checksize(md,fields,[md.numberofelements 6]);
92end
93
94%NO NAN
95fields={'numberofelements','numberofgrids','x','y','z','drag_coefficient','drag_type','drag_p','drag_q',...
96 'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
97 'tolx','np','eps_res','max_nonlinear_iterations','exclusive','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
98checknan(md,fields);
99
100%FIELDS >= 0
101fields={'numberofelements','numberofgrids','elements','drag_coefficient','drag_type','drag_p','drag_q',...
102 'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
103 'sparsity','lowmem','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
104checkgreater(md,fields,0);
105
106%FIELDS > 0
107fields={'numberofelements','numberofgrids','elements','drag_type','drag_p',...
108 'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter','tolx',...
109 'sparsity','deltaH','DeltaH','timeacc','timedec'};
110checkgreaterstrict(md,fields,0);
111
112%SIZE NUMBEROFELEMENTS
113fields={'drag_p','drag_q','elementoniceshelf','n','elementonbed'};
114checksize(md,fields,[md.numberofelements 1]);
115
116%SIZE NUMBEROFGRIDS
117fields={'x','y','z','rheology_B','drag_coefficient','melting_rate','accumulation_rate','surface','thickness','bed','gridonbed','gridonsurface'};
118checksize(md,fields,[md.numberofgrids 1]);
119
120%SIZE 6
121fields={'spcvelocity'};
122checksize(md,fields,[md.numberofgrids 6]);
123
124%THICKNESS = SURFACE - BED
125if any((md.thickness-md.surface+md.bed)>tolerance),
126 error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
127end
128
129%RIFTS
130if md.numrifts,
131 if ~strcmpi(md.type,'2d'),
132 error(['model not consistent: models with rifts are only supported in 2d for now!']);
133 end
134 if ~isstruct(md.rifts),
135 error(['model not consistent: md.rifts should be a structure!']);
136 end
137 if ~isempty(find(md.segmentmarkers>=2)),
138 %We have segments with rift markers, but no rift structure!
139 error(['model not consistent: model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
140 end
141 %Check that rifts are filled with proper material
142 checkvalues(md,{'rifts.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
143else
144 if ~isnans(md.rifts),
145 error(['model not consistent: md.rifts shoud be NaN since md.numrifts is 0!']);
146 end
147end
148
149%ARTIFICIAL DIFFUSIVITY
150if ~ismember(md.artificial_diffusivity,[0 1]),
151 error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
152end
153
154%PARAMETEROUTPUT
155if md.numoutput~=length(md.parameteroutput),
156 error('model not consistent: numoutput should be the same size as parameteroutput');
157end
158
159%CONNECTIVITY
160if strcmpi(md.type,'2d'),
161 if md.connectivity<9,
162 error('model not consistent: connectivity should be at least 9 for 2d models');
163 end
164end
165if strcmpi(md.type,'3d'),
166 if md.connectivity<24,
167 error('model not consistent: connectivity should be at least 24 for 3d models');
168 end
169end
170
171%LOWMEM = 0 or 1
172if ((md.lowmem ~= 1) & (md.lowmem~=0)),
173 error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
174end
175
176%PARALLEL
177if ~strcmpi(md.cluster,'none'),
178
179 %NAN VALUES
180 fields={'time','np'};
181 checknan(md,fields);
182
183 %FIELD > 0
184 fields={'time','np'};
185 checkgreaterstrict(md,fields,0);
186
187end
188
189%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
190
191%QMU
192if md.qmu_analysis,
193 if md.qmu_params.evaluation_concurrency~=1,
194 error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
195 end
196 if ~isempty(md.part),
197 if numel(md.part)~=md.numberofgrids,
198 error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
199 end
200 if find(md.part)>=md.numberofgrids,
201 error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
202 end
203 if min(md.part)~=0,
204 error(['model not consistent: partition vector not indexed from 0 on']);
205 end
206 if max(md.part)>=md.numberofgrids,
207 error(['model not consistent: partition vector cannot have maximum index larger than number of grids']);
208 end
209 if ~isempty(find(md.part<0)),
210 error(['model not consistent: partition vector cannot have values less than 0']);
211 end
212 if ~isempty(find(md.part>=md.npart)),
213 error(['model not consistent: partition vector cannot have values more than md.npart-1']);
214 end
215
216 end
217 if md.eps_rel>1.1*10^-5,
218 error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
219 end
220end
221
222%DIAGNOSTIC
223if md.analysis_type==DiagnosticAnalysisEnum,
224
225 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
226 if strcmpi(md.type,'2d'),
227 if isempty(find(~md.spcvelocity(:,1:2))),
228 error(['model not consistent: model ' md.name ' is totally constrained, no need to solve!']);
229 end
230 end
231
232 %HUTTER ON ICESHELF WARNING
233 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
234 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
235 end
236
237 %SINGULAR
238 if ~any(sum(md.spcvelocity(:,1:2),2)==2),
239 error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
240 end
241
242 %DIRICHLET IF THICKNESS <= 0
243 if any(md.thickness<=0),
244 pos=find(md.thickness<=0);
245 if any(find(md.spcthickness(pos,1)==0)),
246 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
247 end
248 end
249
250 %INITIAL VELOCITY
251 if length(md.vx)==md.numberofgrids & length(md.vy)==md.numberofgrids,
252 fields={'vx','vy'};
253 checknan(md,fields);
254 end
255end
256
257%PROGNOSTIC
258if md.analysis_type==PrognosticAnalysisEnum,
259
260 %INITIAL VELOCITIES
261 fields={'vx','vy'};
262 checksize(md,fields,[md.numberofgrids 1]);
263 checknan(md,fields);
264
265 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
266 if strcmpi(md.type,'2d'),
267 if isempty(find(~md.spcthickness(:,1))),
268 error(['model not consistent: model ' md.name ' is totally constrained for prognostic, no need to solve!']);
269 end
270 end
271
272end
273
274%STEADYSTATE
275if md.analysis_type==SteadystateAnalysisEnum,
276
277 %NDT
278 if md.dt~=0,
279 error(['model not consistent: for a steadystate computation, dt must be zero.']);
280 end
281
282 %PRESSURE
283 if isnans(md.pressure),
284 error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
285 end
286
287 %eps:
288 if isnan(md.eps_rel),
289 error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
290 end
291
292 %dim:
293 if strcmpi(md.type,'2d'),
294 error(['model not consistent: for a steadystate computation, model needs to be 3d']);
295 end
296end
297
298%THERMAL STEADY AND THERMAL TRANSIENT
299if md.analysis_type==ThermalAnalysisEnum,
300
301 %EXTRUSION
302 if strcmp(md.type,'2d'),
303 error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
304 end
305
306 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
307 if strcmpi(md.type,'2d'),
308 if isempty(find(~md.spctemperature(:,1))),
309 error(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
310 end
311 end
312
313 %VELOCITIES AND PRESSURE
314 fields={'vx','vy','vz','pressure'};
315 checksize(md,fields,[md.numberofgrids 1]);
316 checknan(md,fields);
317
318end
319
320%THERMAL TRANSIENT
321if md.analysis_type==ThermalAnalysisEnum & md.dt~=0,
322
323 %DT and NDT
324 fields={'dt','ndt'};
325 checkgreaterstrict(md,fields,0);
326
327 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
328 fields={'temperature','accumulation_rate','melting_rate'};
329 checksize(md,fields,[md.numberofgrids 1]);
330 checknan(md,fields);
331
332 %INITIAL TEMPERATURE
333 fields={'temperature','spctemperature(:,2)','observed_temperature'};
334 checkgreater(md,fields,0)
335
336end
337
338%BALANCEDTHICKNESS
339if md.analysis_type==BalancedthicknessAnalysisEnum
340
341 %VELOCITIES MELTING AND ACCUMULATION
342 fields={'vx','vy','accumulation_rate','melting_rate'};
343 checksize(md,fields,[md.numberofgrids 1]);
344 checknan(md,fields);
345
346 %SPC
347 if any(md.spcthickness(find(md.gridonboundary))~=1),
348 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);
349 end
350end
351
352%BALANCEDTHICKNESS2
353if md.analysis_type==Balancedthickness2AnalysisEnum
354
355 %VELOCITIES MELTING AND ACCUMULATION
356 fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
357 checksize(md,fields,[md.numberofgrids 1]);
358 checknan(md,fields);
359
360 %CHECK THE COMPATIBILITY OF THE EDGES
361 fields={'edges'};
362 checksize(md,fields,[NaN 4]);
363 checknan(md,fields);
364end
365
366%BALANCEDVELOCITIES
367if md.analysis_type==BalancedvelocitiesAnalysisEnum
368
369 %VELOCITIES MELTING AND ACCUMULATION
370 fields={'vx','vy','accumulation_rate','melting_rate'};
371 checksize(md,fields,[md.numberofgrids 1]);
372 checknan(md,fields);
373
374 %SPC
375 if any(md.spcvelocity(find(md.gridonboundary),[1:2])~=1),
376 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcvelocity']);
377 end
378end
379
380%CONTROL
381if md.control_analysis,
382
383 %CONTROL TYPE
384 if ~ischar(md.control_type),
385 error('model not consistent: control_type should be a string');
386 end
387
388 %LENGTH CONTROL FIELDS
389 fields={'maxiter','optscal','fit','cm_jump'};
390 checksize(md,fields,[md.nsteps 1]);
391
392 %FIT
393 checkvalues(md,{'fit'},[0 1 2 3 4]);
394
395 %WEIGHTS
396 fields={'weights'};
397 checksize(md,fields,[md.numberofgrids 1]);
398 checkgreater(md,fields,0);
399
400 %OBSERVED VELOCITIES
401 fields={'vx_obs','vy_obs'};
402 checksize(md,fields,[md.numberofgrids 1]);
403 checknan(md,fields);
404
405 %DIRICHLET IF THICKNESS <= 0
406 if any(md.thickness<=0),
407 pos=find(md.thickness<=0);
408 if any(find(md.spcthickness(pos,1)==0)),
409 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
410 end
411 end
412
413 %parameters
414 fields={'cm_noisedmp'};
415 checknan(md,fields);
416end
417
418%QMU
419if strcmpi(md.analysis_type,'qmu'),
420 if ~strcmpi(md.cluster,'none'),
421 if md.waitonlock==0,
422 error(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
423 end
424 end
425end
426
427end %end function
428
429%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
430
431function checklength(md,fields,fieldlength)
432 %CHECKSIZE - check length of a field
433 for i=1:length(fields),
434 if length(eval(['md.' fields{i}]))~=fieldlength,
435 error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
436 end
437 end
438end
439
440function checksize(md,fields,fieldsize)
441 %CHECKSIZE - check size of a field
442 for i=1:length(fields),
443 if isnan(fieldsize(1)),
444 if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
445 %LOG warnings:
446 if strcmpi(fields{i},'pressureload'),
447 disp(' ');
448 disp(' If pressureload does not have the right size (one column missing), this is due to a change in md.pressure.');
449 disp(' To update your model, use the following line');
450 disp(' md.pressureload=[md.pressureload WaterEnum*md.elementoniceshelf(md.pressureload(:,end))+AirEnum*md.elementonicesheet(md.pressureload(:,end))];');
451 disp(' ');
452 end
453 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
454 end
455 elseif isnan(fieldsize(2)),
456 if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
457 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
458 end
459 else
460 if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) | (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
461 %LOG warnings:
462 if strcmpi(fields{i},'weights'),
463 disp(' ');
464 disp(' ''weights'' is a new field that should be of length numberofgrids.');
465 disp(' To update your model, use the following line');
466 disp(' md.weights=ones(md.numberofgrids,1)');
467 disp(' ');
468 end
469 error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
470 end
471 end
472 end
473end
474
475function checknan(md,fields)
476 %CHECKNAN - check nan values of a field
477 for i=1:length(fields),
478 if any(isnan(eval(['md.' fields{i}]))),
479 error(['model not consistent: NaN values in field ' fields{i}]);
480 end
481 end
482end
483
484function checkreal(md,fields)
485 %CHECKREAL - check real values of a field
486 for i=1:length(fields),
487 if any(eval(['~isreal(md.' fields{i} ')'])),
488 error(['model not consistent: complex values in field ' fields{i}]);
489 end
490 end
491end
492
493function checkgreaterstrict(md,fields,lowerbound)
494 %CHECKGREATERSTRICT - check values of a field
495 for i=1:length(fields),
496 if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
497 error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
498 end
499 end
500end
501
502function checkgreater(md,fields,lowerbound)
503 %CHECKGREATER - check values of a field
504 for i=1:length(fields),
505 if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
506 error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
507 end
508 end
509end
510
511function checklessstrict(md,fields,upperbound)
512 %CHECKLESSSTRICT - check values of a field
513 for i=1:length(fields),
514 if any(eval(['md.' fields{i} '>=' num2str(upperbound) ])),
515 error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
516 end
517 end
518end
519
520function checkless(md,fields,upperbound)
521 %CHECKLESS - check values of a field
522 for i=1:length(fields),
523 if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
524 error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
525 end
526 end
527end
528
529function checkvalues(md,fields,values)
530 %CHECKVALUE - check that a field has a certain value
531 for i=1:length(fields),
532 if eval(['~ismember( md.' fields{i} ',values)']),
533 error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
534 end
535 end
536end
Note: See TracBrowser for help on using the repository browser.