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

Last change on this file since 4756 was 4756, checked in by Mathieu Morlighem, 15 years ago

Added folds

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