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

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

no dirichlet needed in balancedthickness DG

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