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

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

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