0001 function bool=ismodelselfconsistent(md,solutiontype,package),
0002
0003
0004
0005
0006
0007 if strcmpi(solutiontype,'mesh'),
0008
0009
0010 bool=1;
0011 return;
0012 end
0013
0014
0015 tolerance=10^-12;
0016 if (nargin~=3 )
0017 IsModelSelfConsistentUsage();
0018 error(' ');
0019 end
0020
0021 if strcmpi(solutiontype,'qmu'),
0022 if ~strcmpi(md.cluster,'none'),
0023 if md.waitonlock==0,
0024 disp(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
0025 bool=0;return;
0026 end
0027 end
0028 end
0029
0030 if isempty(md.name),
0031 disp(['model is not correctly configured: missing name!']);
0032 bool=0;return;
0033 end
0034
0035 if md.counter<3,
0036 disp(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize)!']);
0037 bool=0;return;
0038 end
0039
0040 if md.numberofelements<=0,
0041 disp(['model ' md.name ' does not have any elements!']);
0042 bool=0; return;
0043 end
0044
0045 if md.numberofgrids<=0,
0046 disp(['model ' md.name ' does not have any grids!']);
0047 bool=0; return;
0048 end
0049
0050
0051 if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
0052 disp(['Types of elements have not been set properly, run setelementstype first'])
0053 bool=0;return;
0054 end
0055 if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==hutterenum) + (md.elements_type(:,1)==macayealenum) + (md.elements_type(:,1)==pattynenum)))
0056 disp(['Types of elements have not been set properly, run setelementstype first'])
0057 bool=0;return;
0058 end
0059 if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==stokesenum) + (md.elements_type(:,2)==noneenum)))
0060 disp(['Types of elements have not been set properly, run setelementstype first'])
0061 bool=0;return;
0062 end
0063 if strcmpi(md.type,'2d'),
0064 if (ismember(pattynenum,md.elements_type(:,1)) | ismember(stokesenum,md.elements_type(:,2))),
0065 disp(['For a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
0066 bool=0;return;
0067 end
0068 end
0069
0070 if md.numrifts,
0071 if ~strcmpi(md.type,'2d'),
0072 disp(['Models with rifts are only supported in 2d for now!']);
0073 bool=0;return;
0074 end
0075 end
0076
0077
0078 if md.numrifts,
0079 if ~strcmpi(md.cluster,'none')
0080 if isempty(findstr('aijmumps',md.solverstring)),
0081 disp(['For a 2d model with rifts, running on a cluster, you should use a direct solver like MUMPS!']);
0082 bool=0;return;
0083 end
0084 end
0085 end
0086
0087
0088 if strcmpi(solutiontype,'diagnostic')
0089 if any(md.elements_type(:,1)==hutterenum & md.elementoniceshelf),
0090 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
0091 end
0092 end
0093
0094
0095 if strcmpi(solutiontype,'thermalsteady') | strcmp(solutiontype,'thermaltransient'),
0096 if strcmp(md.type,'2d'),
0097 disp(['For a ' solutiontype ' computation, the model must be 3d, extrude it first!'])
0098 bool=0;return;
0099 end
0100 end
0101
0102
0103 if ~iscell(md.control_type),
0104 disp('control_type should be a cell array of strings');
0105 bool=0;return;
0106 end
0107
0108
0109 if ~isscalar(md.artificial_diffusivity),
0110 disp('artificial_diffusivity should be a scalar (1 or 0)');
0111 bool=0;return;
0112 end
0113
0114
0115 if strcmpi(md.type,'2d'),
0116 if md.connectivity<9,
0117 disp('connectivity should be at least 9 for 2d models');
0118 bool=0;return;
0119 end
0120 end
0121
0122 if strcmpi(md.type,'3d'),
0123 if md.connectivity<24,
0124 disp('connectivity should be at least 24 for 3d models');
0125 bool=0;return;
0126 end
0127 end
0128
0129
0130
0131
0132
0133
0134 if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
0135 if (isempty(md.vx) | isempty(md.vy) | isempty(md.vz)),
0136 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
0137 bool=0;return;
0138 end
0139 end
0140
0141 if strcmp(solutiontype,'thermaltransient') | strcmp(solutiontype,'thermalsteady')
0142 if (isempty(md.pressure) | isnan(md.pressure)),
0143 disp(['pressure is required. Run ''diagnostic'' solution first!'])
0144 bool=0;return;
0145 end
0146
0147 if (isempty(md.vx) | isnan(md.vx) | isempty(md.vy) | isnan(md.vy) | isempty(md.vz) | isnan(md.vz)),
0148 disp(['velocities are required!']);
0149 bool=0;return;
0150 end
0151 end
0152
0153 if strcmp(solutiontype,'thermaltransient')
0154 if isempty(md.temperature),
0155 disp(['An initial temperature is needed for a transient thermal computation'])
0156 bool=0;return;
0157 end
0158 if isstruct(md.temperature) | isstruct(md.melting),
0159 disp(['The initial temperature or melting should be a list and not a structure'])
0160 bool=0;return;
0161 end
0162
0163 end
0164
0165 if strcmp(solutiontype,'transient')
0166 if isstruct(md.temperature),
0167 disp(['The initial temperature should be empty or a list but not a structure'])
0168 bool=0;return;
0169 end
0170
0171 end
0172
0173 if strcmp(solutiontype,'parameters')
0174 if ~(size(md.vx,1)==md.numberofgrids & size(md.vy,1)==md.numberofgrids & (size(md.vz,1)==md.numberofgrids | strcmpi(md.type,'2d')))
0175 disp(['velocities are required!']);
0176 bool=0;return;
0177 end
0178 if strcmpi(md.type,'2d') & md.acceleration
0179 disp(['solution requires acceleration=0']);
0180 bool=0;return;
0181 end
0182 if any(md.elements_type(:,1)==hutterenum);
0183 disp(['The model has Hutter''s elements. Impossible to compute parameters']);
0184 bool=0;return;
0185 end
0186 end
0187
0188
0189
0190 if strcmp(solutiontype,'diagnostic'),
0191 if isempty(find(md.gridondirichlet_diag)),
0192 disp(['model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
0193 bool=0;return;
0194 end
0195 end
0196
0197
0198 if strcmp(solutiontype,'diagnostic'),
0199 if isempty(md.gridondirichlet_diag) | isempty(md.dirichletvalues_diag),
0200 disp(['model ' md.name ' is not well posed. Missing dirichlet values for diagnostic run']);
0201 bool=0;return;
0202 end
0203 end
0204
0205
0206
0207 if strcmp(solutiontype,'prognostic'),
0208
0209
0210
0211
0212
0213
0214
0215
0216
0217 if (size(md.vx,1)~=md.numberofgrids | size(md.vy,1)~=md.numberofgrids),
0218 disp(['a 3d velocity is required. Run ''diagnostic'' solution first!'])
0219 bool=0; return;
0220 end
0221 end
0222
0223 if strcmp(solutiontype,'control')
0224
0225 if ~iscell(md.control_type),
0226 disp('control_type should be a cell array of strings');
0227 bool=0;return;
0228 end
0229
0230
0231 if (length(md.maxiter)~=md.nsteps | length(md.optscal)~=md.nsteps | length(md.fit)~=md.nsteps)
0232 disp('maxiter, optscal and fit must have the length specified by nsteps')
0233 bool=0;return;
0234 end
0235 if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
0236 disp('wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
0237 bool=0;return;
0238 end
0239 end
0240
0241
0242
0243 fields={'numberofelements','numberofgrids','x','y','z','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
0244 'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
0245 'tolx','np','time','eps_rel','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
0246 for i=1:length(fields),
0247 if ~isempty(eval(['md.' char(fields(i))])),
0248 if find(isnan(eval(['md.' char(fields(i))]))),
0249 disp(['model ' md.name ' has an NaN value in field ' char(fields(i)) '!']);
0250 bool=0; return;
0251 end
0252 end
0253 end
0254
0255
0256 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag','drag_type','p','q','gridondirichlet_diag',...
0257 'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_rel','eps_abs','nsteps','maxiter','tolx','np','time','exclusive',...
0258 'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
0259 for i=1:length(fields),
0260 if ~isempty(eval(['md.' char(fields(i))])),
0261 if find((eval(['md.' char(fields(i))]))<0),
0262 disp(['model ' md.name ' has a <0 value in field ' char(fields(i)) '!']);
0263 bool=0; return;
0264 end
0265 end
0266 end
0267
0268
0269 fields={'numberofelements','numberofgrids','elements','segmentonneumann_diag','drag_type',...
0270 'rho_ice','rho_water','B','thickness','g','eps_rel','eps_abs','maxiter','tolx','np','time',...
0271 'sparsity','deltaH','DeltaH','timeacc','timedec'};
0272 for i=1:length(fields),
0273 if ~isempty(eval(['md.' char(fields(i))])),
0274 if find((eval(['md.' char(fields(i))]))==0),
0275 disp(['model ' md.name ' has a =0 value in field ' char(fields(i)) '!']);
0276 bool=0; return;
0277 end
0278 end
0279 end
0280
0281
0282 fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
0283 for i=1:size(fields,2),
0284 if (size(eval(['md.' char(fields(i))]),1)~=md.numberofelements),
0285 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofelements) '!']);
0286 bool=0; return;
0287 end
0288 end
0289
0290
0291 fields={'melting','x','y','z','B','drag','gridondirichlet_diag','vx_obs','vy_obs','surface','thickness','bed','gridonbed','gridonsurface'};
0292 for i=1:length(fields),
0293 if find(length(eval(['md.' char(fields(i))]))~=md.numberofgrids),
0294 disp(['model ' md.name ' field ' char(fields(i)) ' should be of size ' num2str(md.numberofgrids) '!']);
0295 bool=0; return;
0296 end
0297 end
0298
0299
0300 if ((md.lowmem ~= 1) & (md.lowmem~=0)),
0301 disp(['model ' md.name ' lowmem field should be 0 or 1']);
0302 bool=0; return;
0303 end
0304
0305
0306 if ( (md.sparsity<=0) | (md.sparsity>1)),
0307 disp(['model ' md.name ' sparsity should be inside the ]0 1] range']);
0308 bool=0; return;
0309 end
0310
0311
0312 if find((md.thickness-md.surface+md.bed)>tolerance),
0313 disp(['model ' md.name ' violates the equality thickness=surface-bed!']);
0314 bool=0; return;
0315 end
0316
0317
0318 if ~isstruct(md.rifts),
0319 if ~isempty(find(md.segmentmarkers>=2)),
0320
0321 disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
0322 bool=0; return;
0323 end
0324 end
0325
0326
0327
0328 if ~isempty(find(md.thickness<=0)),
0329 pos=find(md.thickness<=0);
0330 if ~isempty(find(md.gridondirichlet_diag(pos)==0)),
0331 disp(['model ' md.name ' has some grids with 0 thickness']);
0332 bool=0; return;
0333 end
0334 end
0335
0336
0337 if ~isempty(find(md.p<=0)),
0338 disp(['model ' md.name ' has some p<0 friction coefficientsin sliding law']);
0339 bool=0; return;
0340 end
0341
0342
0343
0344 if ~iscell(md.parameteroutput)
0345 disp(['parameteroutput field must be a cell, example {''strainrate'',''stress'',''deviatoricstress'',''viscousheating''}']);
0346 bool=0; return;
0347 end
0348
0349
0350 for i=1:length(md.parameteroutput)
0351 if ~strcmpi(md.parameteroutput(i),'strainrate') & ~strcmpi(md.parameteroutput(i),'stress') & ~strcmpi(md.parameteroutput(i),'deviatoricstress') & ~strcmpi(md.parameteroutput(i),'viscousheating') ...
0352 & ~strcmpi(md.parameteroutput(i),'pressure_elem') & ~strcmpi(md.parameteroutput(i),'stress_bed') & ~strcmpi(md.parameteroutput(i),'stress_surface')
0353 disp(['one of the parameteroutput is not supported yet']);
0354 bool=0; return;
0355 end
0356 end
0357
0358
0359
0360
0361
0362 bool=1;
0363 return;
0364
0365 function IsModelSelfConsistentUsage()
0366 disp(' ');
0367 disp(' IsModelSelfConsistent usage: flag=IsModelSelfConsistent(model)');
0368 disp(' where model is an instance of the @model class, and flag a boolean');
0369
0370
0371
0372