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

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

updated enums and some fixing to launch parallel solution

File size: 16.8 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;
[2326]9
10%check usage
11if nargin~=1,
[27]12 help ismodelselfconsistent
13 error('ismodelselfconsistent error message: wrong usage');
[1]14end
15
[2326]16%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% TRANSIENT %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[1670]17
[2326]18%check analysis
[4089]19if (md.analysis_type==Transient2DAnalysisEnum | md.analysis_type==Transient3DAnalysisEnum),
[2326]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
[3994]25 if (md.dim==2),
[2528]26 analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum];
27 else
28 analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum ThermalAnalysisEnum];
29 end
30
[2326]31 for i=1:length(analysis),
[2520]32 md.analysis_type=analysis(i);
[2326]33 ismodelselfconsistent(md);
34 end
35end
36
37
[27]38%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMMON CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
39
40%COUNTER
41if md.counter<3,
[2326]42 error(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize,setelementstype)!']);
[1]43end
44
[27]45%NAME
[1]46if isempty(md.name),
[2326]47 error(['model is not correctly configured: missing name!']);
[1]48end
49
[27]50%ELEMENTSTYPE
[1]51if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
[2326]52 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
[1]53end
[1651]54if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum) + (md.elements_type(:,1)==PattynFormulationEnum)))
[2326]55 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
[1]56end
[1651]57if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
[2326]58 error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
[1]59end
[3994]60if (md.dim==2),
[1651]61 if (ismember(PattynFormulationEnum,md.elements_type(:,1)) | ismember(StokesFormulationEnum,md.elements_type(:,2))),
[2326]62 error(['model not consistent: for a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
[1]63 end
64end
[358]65if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
[2326]66 error(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
[358]67end
[2164]68if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
[2326]69 error(['model not consistent: coupling MacAyeal/Pattyn not implemented yet']);
[2164]70end
[4089]71if (md.isstokes & md.analysis_type==Transient3DAnalysisEnum());
[2326]72 error(['model not consistent: Stokes transient not implemented yet']);
[2164]73end
[358]74
[1796]75%ICEFRONT
[3994]76if (md.dim==2),
[2326]77 fields={'pressureload'};
[3099]78 checksize(md,fields,[NaN 4]);
[3994]79elseif md.dim==3,
[2326]80 fields={'pressureload'};
[3099]81 checksize(md,fields,[NaN 6]);
[3356]82else
[3357]83 error('type should be either ''2d'' or ''3d''');
[1796]84end
[358]85
[3356]86%ELEMENTS
87fields={'elements'};
[3994]88if (md.dim==2),
[3356]89 checksize(md,fields,[md.numberofelements 3]);
90else
91 checksize(md,fields,[md.numberofelements 6]);
92end
93
[27]94%NO NAN
[3758]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'};
[2366]98checknan(md,fields);
[1]99
[1666]100%FIELDS >= 0
[3758]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'};
[2366]104checkgreater(md,fields,0);
[1]105
[2326]106%FIELDS > 0
[3758]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',...
[27]109 'sparsity','deltaH','DeltaH','timeacc','timedec'};
[2366]110checkgreaterstrict(md,fields,0);
[1]111
[27]112%SIZE NUMBEROFELEMENTS
[3759]113fields={'drag_p','drag_q','elementoniceshelf','rheology_n','elementonbed'};
[2991]114checksize(md,fields,[md.numberofelements 1]);
115
[27]116%SIZE NUMBEROFGRIDS
[3758]117fields={'x','y','z','rheology_B','drag_coefficient','melting_rate','accumulation_rate','surface','thickness','bed','gridonbed','gridonsurface'};
[2991]118checksize(md,fields,[md.numberofgrids 1]);
[27]119
[2405]120%SIZE 6
121fields={'spcvelocity'};
122checksize(md,fields,[md.numberofgrids 6]);
123
[27]124%THICKNESS = SURFACE - BED
[681]125if any((md.thickness-md.surface+md.bed)>tolerance),
[2326]126 error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
[27]127end
128
129%RIFTS
130if md.numrifts,
[3994]131 if ~(md.dim==2),
[2326]132 error(['model not consistent: models with rifts are only supported in 2d for now!']);
[1]133 end
[2660]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
[2748]144 if ~isnans(md.rifts),
[2660]145 error(['model not consistent: md.rifts shoud be NaN since md.numrifts is 0!']);
[27]146 end
147end
[1]148
[27]149%ARTIFICIAL DIFFUSIVITY
[2326]150if ~ismember(md.artificial_diffusivity,[0 1]),
151 error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
[1]152end
153
[2307]154%PARAMETEROUTPUT
155if md.numoutput~=length(md.parameteroutput),
[2326]156 error('model not consistent: numoutput should be the same size as parameteroutput');
[2307]157end
158
[2326]159%CONNECTIVITY
[3994]160if (md.dim==2),
[2326]161 if md.connectivity<9,
162 error('model not consistent: connectivity should be at least 9 for 2d models');
163 end
164end
[3994]165if md.dim==3,
[2326]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'};
[2366]181 checknan(md,fields);
[2326]182
183 %FIELD > 0
184 fields={'time','np'};
[2366]185 checkgreaterstrict(md,fields,0);
[2326]186
187end
188
[27]189%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SOLUTION CHECKS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2326]190
[804]191%QMU
[966]192if md.qmu_analysis,
193 if md.qmu_params.evaluation_concurrency~=1,
[2326]194 error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
[966]195 end
[1111]196 if ~isempty(md.part),
197 if numel(md.part)~=md.numberofgrids,
[2326]198 error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
[1111]199 end
200 if find(md.part)>=md.numberofgrids,
[2326]201 error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
[1111]202 end
[3037]203 if min(md.part)~=0,
204 error(['model not consistent: partition vector not indexed from 0 on']);
[1111]205 end
[3037]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
[3055]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
[1111]216 end
[3037]217 if md.eps_rel>1.1*10^-5,
[2326]218 error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
[2052]219 end
[804]220end
[27]221
222%DIAGNOSTIC
[1781]223if md.analysis_type==DiagnosticAnalysisEnum,
[27]224
[2208]225 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
[3994]226 if (md.dim==2),
[2208]227 if isempty(find(~md.spcvelocity(:,1:2))),
[2326]228 error(['model not consistent: model ' md.name ' is totally constrained, no need to solve!']);
[2208]229 end
230 end
231
[27]232 %HUTTER ON ICESHELF WARNING
[1651]233 if any(md.elements_type(:,1)==HutterFormulationEnum & md.elementoniceshelf),
[27]234 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
235 end
236
237 %SINGULAR
[2282]238 if ~any(sum(md.spcvelocity(:,1:2),2)==2),
[2326]239 error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
[1]240 end
241
[27]242 %DIRICHLET IF THICKNESS <= 0
[681]243 if any(md.thickness<=0),
[27]244 pos=find(md.thickness<=0);
[1755]245 if any(find(md.spcthickness(pos,1)==0)),
[2326]246 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
[27]247 end
248 end
[2357]249
[2367]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
[1]255end
256
[27]257%PROGNOSTIC
[1781]258if md.analysis_type==PrognosticAnalysisEnum,
[1]259
[2357]260 %INITIAL VELOCITIES
261 fields={'vx','vy'};
[2366]262 checksize(md,fields,[md.numberofgrids 1]);
263 checknan(md,fields);
[2357]264
[2519]265 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
[3994]266 if (md.dim==2),
[2519]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
[1]272end
273
[1867]274%STEADYSTATE
275if md.analysis_type==SteadystateAnalysisEnum,
[2326]276
[1901]277 %NDT
278 if md.dt~=0,
[2326]279 error(['model not consistent: for a steadystate computation, dt must be zero.']);
[1901]280 end
[2326]281
[1787]282 %PRESSURE
283 if isnans(md.pressure),
[2326]284 error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
[1787]285 end
286
287 %eps:
288 if isnan(md.eps_rel),
[2326]289 error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
[1787]290 end
291
292 %dim:
[3994]293 if (md.dim==2),
[2326]294 error(['model not consistent: for a steadystate computation, model needs to be 3d']);
[1787]295 end
296end
297
[27]298%THERMAL STEADY AND THERMAL TRANSIENT
[1781]299if md.analysis_type==ThermalAnalysisEnum,
[27]300
301 %EXTRUSION
[3994]302 if (md.dim==2),
[2326]303 error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
[1]304 end
[27]305
[2519]306 %CHECK THAT WE ARE NOT FULLY CONSTRAINED
[3994]307 if isempty(find(~md.spctemperature(:,1))),
308 error(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
[2519]309 end
310
[27]311 %VELOCITIES AND PRESSURE
[2357]312 fields={'vx','vy','vz','pressure'};
[2366]313 checksize(md,fields,[md.numberofgrids 1]);
314 checknan(md,fields);
[2357]315
[1]316end
317
[1311]318%THERMAL TRANSIENT
[2520]319if md.analysis_type==ThermalAnalysisEnum & md.dt~=0,
[27]320
321 %DT and NDT
322 fields={'dt','ndt'};
[2366]323 checkgreaterstrict(md,fields,0);
[27]324
[1311]325 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
[3758]326 fields={'temperature','accumulation_rate','melting_rate'};
[2366]327 checksize(md,fields,[md.numberofgrids 1]);
[2440]328 checknan(md,fields);
[1]329
[2714]330 %INITIAL TEMPERATURE
[2669]331 fields={'temperature','spctemperature(:,2)','observed_temperature'};
332 checkgreater(md,fields,0)
333
[1]334end
335
[2714]336%BALANCEDTHICKNESS
337if md.analysis_type==BalancedthicknessAnalysisEnum
338
339 %VELOCITIES MELTING AND ACCUMULATION
[3758]340 fields={'vx','vy','accumulation_rate','melting_rate'};
[3589]341 checksize(md,fields,[md.numberofgrids 1]);
342 checknan(md,fields);
343
344 %SPC
345 if any(md.spcthickness(find(md.gridonboundary))~=1),
346 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);
347 end
348end
349
350%BALANCEDTHICKNESS2
351if md.analysis_type==Balancedthickness2AnalysisEnum
352
353 %VELOCITIES MELTING AND ACCUMULATION
[3758]354 fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
[2714]355 checksize(md,fields,[md.numberofgrids 1]);
356 checknan(md,fields);
357
[3653]358 %CHECK THE COMPATIBILITY OF THE EDGES
359 fields={'edges'};
360 checksize(md,fields,[NaN 4]);
361 checknan(md,fields);
[2714]362end
363
[2723]364%BALANCEDVELOCITIES
365if md.analysis_type==BalancedvelocitiesAnalysisEnum
366
367 %VELOCITIES MELTING AND ACCUMULATION
[3758]368 fields={'vx','vy','accumulation_rate','melting_rate'};
[2723]369 checksize(md,fields,[md.numberofgrids 1]);
370 checknan(md,fields);
371
372 %SPC
373 if any(md.spcvelocity(find(md.gridonboundary),[1:2])~=1),
374 error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcvelocity']);
375 end
376end
377
[27]378%CONTROL
[1907]379if md.control_analysis,
[1]380
[27]381 %CONTROL TYPE
[45]382 if ~ischar(md.control_type),
[2326]383 error('model not consistent: control_type should be a string');
[1]384 end
385
[27]386 %LENGTH CONTROL FIELDS
[2326]387 fields={'maxiter','optscal','fit','cm_jump'};
[2991]388 checksize(md,fields,[md.nsteps 1]);
[27]389
390 %FIT
[3189]391 checkvalues(md,{'fit'},[0 1 2 3 4]);
[1]392
[3203]393 %WEIGHTS
394 fields={'weights'};
395 checksize(md,fields,[md.numberofgrids 1]);
396 checkgreater(md,fields,0);
397
[27]398 %OBSERVED VELOCITIES
399 fields={'vx_obs','vy_obs'};
[2366]400 checksize(md,fields,[md.numberofgrids 1]);
401 checknan(md,fields);
[1361]402
403 %DIRICHLET IF THICKNESS <= 0
404 if any(md.thickness<=0),
405 pos=find(md.thickness<=0);
[1755]406 if any(find(md.spcthickness(pos,1)==0)),
[2326]407 error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
[1361]408 end
409 end
[2686]410
411 %parameters
412 fields={'cm_noisedmp'};
413 checknan(md,fields);
[1]414end
415
[27]416%QMU
[465]417if strcmpi(md.analysis_type,'qmu'),
[27]418 if ~strcmpi(md.cluster,'none'),
419 if md.waitonlock==0,
[2326]420 error(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
[1]421 end
422 end
423end
424
[2326]425end %end function
[1]426
[2366]427%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2326]428
[2366]429function checklength(md,fields,fieldlength)
430 %CHECKSIZE - check length of a field
[2326]431 for i=1:length(fields),
432 if length(eval(['md.' fields{i}]))~=fieldlength,
433 error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
[2164]434 end
435 end
436end
[1]437
[2366]438function checksize(md,fields,fieldsize)
439 %CHECKSIZE - check size of a field
[2326]440 for i=1:length(fields),
441 if isnan(fieldsize(1)),
442 if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
[3203]443 %LOG warnings:
[3104]444 if strcmpi(fields{i},'pressureload'),
445 disp(' ');
446 disp(' If pressureload does not have the right size (one column missing), this is due to a change in md.pressure.');
447 disp(' To update your model, use the following line');
448 disp(' md.pressureload=[md.pressureload WaterEnum*md.elementoniceshelf(md.pressureload(:,end))+AirEnum*md.elementonicesheet(md.pressureload(:,end))];');
449 disp(' ');
450 end
[2326]451 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
452 end
453 elseif isnan(fieldsize(2)),
454 if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
455 error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
456 end
457 else
458 if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) | (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
[3219]459 %LOG warnings:
460 if strcmpi(fields{i},'weights'),
461 disp(' ');
462 disp(' ''weights'' is a new field that should be of length numberofgrids.');
463 disp(' To update your model, use the following line');
464 disp(' md.weights=ones(md.numberofgrids,1)');
465 disp(' ');
466 end
[2326]467 error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
468 end
[2164]469 end
470 end
471end
472
[2366]473function checknan(md,fields)
474 %CHECKNAN - check nan values of a field
[2326]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}]);
[2164]478 end
479 end
480end
481
[2366]482function checkreal(md,fields)
483 %CHECKREAL - check real values of a field
[2326]484 for i=1:length(fields),
485 if any(eval(['~isreal(md.' fields{i} ')'])),
486 error(['model not consistent: complex values in field ' fields{i}]);
487 end
488 end
[2164]489end
490
[2366]491function checkgreaterstrict(md,fields,lowerbound)
492 %CHECKGREATERSTRICT - check values of a field
[2326]493 for i=1:length(fields),
494 if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
495 error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
496 end
[2164]497 end
498end
499
[2366]500function checkgreater(md,fields,lowerbound)
501 %CHECKGREATER - check values of a field
[27]502 for i=1:length(fields),
[2326]503 if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
504 error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
[27]505 end
506 end
[2326]507end
[1]508
[2366]509function checklessstrict(md,fields,upperbound)
510 %CHECKLESSSTRICT - 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 stricly below ' num2str(upperbound)]);
[27]514 end
515 end
[2326]516end
[1]517
[2366]518function checkless(md,fields,upperbound)
519 %CHECKLESS - check values of a field
[27]520 for i=1:length(fields),
[2326]521 if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
522 error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
[27]523 end
[1]524 end
525end
[2660]526
527function checkvalues(md,fields,values)
528 %CHECKVALUE - check that a field has a certain value
529 for i=1:length(fields),
530 if eval(['~ismember( md.' fields{i} ',values)']),
531 error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
532 end
533 end
534end
Note: See TracBrowser for help on using the repository browser.