ismodelselfconsistent

PURPOSE ^

ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.

SYNOPSIS ^

function bool=ismodelselfconsistent(md,solutiontype,package),

DESCRIPTION ^

ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.

   Usage:
      bool=ismodelselfconsistent(md,solutiontype,package),

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function bool=ismodelselfconsistent(md,solutiontype,package),
0002 %ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
0003 %
0004 %   Usage:
0005 %      bool=ismodelselfconsistent(md,solutiontype,package),
0006 
0007 if strcmpi(solutiontype,'mesh'),
0008     %this solution is a little special. It should come right after the md=model;  operation. So a lot less checks!
0009 
0010     bool=1;
0011     return;
0012 end
0013 
0014 %tolerance we use in litmus tests for the consistency of the model
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 %check elementstype
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 %check that if there are rifts, model is 2d
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 %If there are rifts, check that the solver is set to mumps in parallel
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 %check elementstype in diagnostic
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 %Check that mesh is 3d
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 %control_type should be a cell
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 %check artificial_diffusivity is an integer
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 %check on connectivity
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 %Check previous computation
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 %Check that problem is not singular
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 %Check we have values for dirichletvalues_diag
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     %old testing
0209 %    if isempty(find(md.gridondirichlet_diag)),
0210 %        disp(['model ' md.name ' diagnostic is not well posed (singular). You need at least one grid with fixed velocity!'])
0211 %        bool=0;return;
0212 %    end
0213 %    if isempty(find(md.gridondirichlet_prog)),
0214 %        disp(['model ' md.name ' prognostic is not well posed (singular). You need at least one grid with fixed thickness!'])
0215 %        bool=0;return;
0216 %    end
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 %some checks for control methods
0223 if strcmp(solutiontype,'control')
0224     %control_type should be a cell
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     %check fields used by control
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 %Check that no NaN values popped up:
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 %Check that these fields are > 0
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 %Check that these fields are ~=0
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 %Check that these fields are of size md.numberofelements
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 %Check that these fields are of size md.numberofgrids
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 %Make sure lowmem is either 0 or 1
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 %Make sure sparsity is between 0 and 1:
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 %Make sure thickness=surface-bed;
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 %Check that rifts were processed correctly:
0318 if ~isstruct(md.rifts),
0319     if ~isempty(find(md.segmentmarkers>=2)),
0320         %We have segments with rift markers, but no rift structure!
0321         disp(['model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
0322         bool=0; return;
0323     end
0324 end
0325 
0326 
0327 %Check that thickness is not null. If is it, then the corresponding grids should be spc'd.
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 %Check that p>0 in the sliding law u = k * sigma ^ p * Neff ^ q :
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 %Check parameteroutput
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 %Put as many checks as you want here:
0359 
0360 
0361 %No problems, just return 1;
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

Generated on Sun 29-Mar-2009 20:22:55 by m2html © 2003