Changeset 9726


Ignore:
Timestamp:
09/09/11 08:44:42 (14 years ago)
Author:
Mathieu Morlighem
Message:

deleted ismodelselfconsistent because each object must check its own consistency

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/model/ismodelselfconsistent.m

    r9725 r9726  
    55%      ismodelselfconsistent(md),
    66
    7 %tolerance we use in litmus checks for the consistency of the model
    8 tolerance=10^-9;
    9 %check usage {{{1
    10 if nargin~=1,
    11         help ismodelselfconsistent
    12         message('ismodelselfconsistent message message: wrong usage');
    13 end
    14 %}}}
    15 
    16 %initialize consistency as true
    17 modelconsistency(true);
    18 
    19 % Flaim doesn't need the entire model, so take an early return
    20 %FLAIM {{{1
    21 if (md.private.solution == FlaimSolutionEnum),
    22         if ~exist(md.flaim.tracks,'file')
    23                 message(['model not consistent: flaim.tracks file ''' md.flaim.tracks ''' must exist.']);
    24         end
    25         %   probably going to need some checks on flaim.flightreqs here
    26         if (numel(md.flaim.criterion) ~= md.mesh.numberofvertices) && (numel(md.flaim.criterion) ~= md.mesh.numberofelements)
    27                 message(['model not consistent: flaim.criterion vector must have number of nodes (' int2str(md.mesh.numberofvertices) ') or elements (' int2str(md.mesh.numberofelements) ') values, not ' int2str(numel(md.flaim.criterion)) ' values.']);
    28         end
    29         return;
    30 end
    31 %}}}
    32 
    33 % Common checks
    34 %VERBOSE{{{1
    35 if ~strcmp('verbose',class(md.verbose))
    36         disp('WARNING: md.verbose must now be an object of class ''verbose''');
    37         disp('         To update your model use the following command:');
    38         disp(' ');
    39         disp('         md.verbose=verbose;');
    40         disp(' ');
    41         message(['field verbose should be of class ''verbose'' ']);
    42 end
    43 %}}}
    44 %NAME{{{1
    45 if isempty(md.miscellaneous.name),
    46         message(['model is not correctly configured: missing name!']);
    47 end
    48 %}}}
    49 %ELEMENTS{{{1
    50 fields={'elements'};
    51 if (md.mesh.dimension==2),
    52         checksize(md,fields,[md.mesh.numberofelements 3]);
    53 else
    54         checksize(md,fields,[md.mesh.numberofelements 6]);
    55 end
    56 if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.elements(:)))));
    57         message('orphan nodes have been found. Check the mesh');
    58 end
    59 %}}}
    60 
    61 %DG {{{1
    62 if md.prognostic.stabilization==3;
    63         %CHECK THE COMPATIBILITY OF THE EDGES
    64         fields={'edges'};
    65         checksize(md,fields,[NaN 4]);
    66         fields={'edges(:,1:3)'};
    67         checknan(md,fields);
    68         checkgreaterstrict(md,fields,0);
    69 end
    70 %}}}
    71 %PRESSURELOAD{{{1
    72 if (md.mesh.dimension==2),
    73         fields={'diagnostic.icefront'};
    74         checksize(md,fields,[NaN 4]);
    75 elseif md.mesh.dimension==3,
    76         fields={'diagnostic.icefront'};
    77         checksize(md,fields,[NaN 6]);
    78 else
    79         message('dim should be either 2 3');
    80 end
    81 checkvalues(md,{'diagnostic.icefront(:,end)'},[0 1 2]);
    82 %}}}
    83 %NO NAN {{{1
    84 fields={'mesh.numberofelements','mesh.numberofvertices','x','y','z','friction.coefficient','friction.p','friction.q',...
    85         'materials.rho_ice','materials.rho_water','materials.rheology_B','mask.elementonfloatingice','geometry.surface','geometry.thickness','geometry.bed','constants.g','settings.lowmem','inversion.nsteps','inversion.maxiter_per_step',...
    86         'diagnostic.restol','diagnostic.maxiter','materials.rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface','elementconnectivity'};
    87 checknan(md,fields);
    88 %}}}}
    89 %FIELDS >= 0 {{{1
    90 fields={'mesh.numberofelements','mesh.numberofvertices','elements','friction.coefficient','friction.p','friction.q',...
    91         'materials.rho_ice','materials.rho_water','materials.rheology_B','mask.elementonfloatingice','geometry.thickness','constants.g','diagnostic.restol','diagnostic.maxiter','diagnostic.reltol',...
    92         'diagnostic.abstol','inversion.nsteps','inversion.maxiter_per_step',...
    93         'settings.lowmem','materials.rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface'};
    94 checkgreater(md,fields,0);
    95 %}}}
    96 %FIELDS > 0 {{{1
    97 fields={'mesh.numberofelements','mesh.numberofvertices','elements','friction.p',...
    98         'materials.rho_ice','materials.rho_water','materials.rheology_B','geometry.thickness','constants.g','diagnostic.maxiter','diagnostic.restol','diagnostic.reltol','diagnostic.abstol','inversion.maxiter_per_step'};
    99 checkgreaterstrict(md,fields,0);
    100 %}}}
    101 %SIZE NUMBEROFELEMENTS {{{1
    102 fields={'friction.p','friction.q','mask.elementonfloatingice','materials.rheology_n','elementonbed'};
    103 checksize(md,fields,[md.mesh.numberofelements 1]);
    104 %}}}
    105 %SIZE NUMBEROFNODES {{{1
    106 fields={'x','y','z','materials.rheology_B','friction.coefficient','basalforcings.melting_rate','geometry.surface','geometry.thickness','geometry.bed','nodeonbed','nodeonsurface'};
    107 checksize(md,fields,[md.mesh.numberofvertices 1]);
    108 %}}}
    109 %OTHER SIZES {{{1
    110 fields={'diagnostic.referential'};
    111 checksize(md,fields,[md.mesh.numberofvertices 6]);
    112 if ~isempty(md.diagnostic.requested_outputs),
    113         if(size(md.diagnostic.requested_outputs,2)~=1),
    114                 message(['model ' md.miscellaneous.name ' requested outputs should be a column vector']);
    115         end
    116 end
    117 %}}}
    118 %THICKNESS = SURFACE - BED {{{1
    119 if any((md.geometry.thickness-md.geometry.surface+md.geometry.bed)>tolerance),
    120         message(['model not consistent: model ' md.miscellaneous.name ' violates the equality thickness=surface-bed!']);
    121 end
    122 %GROUNDING LINE MIGRATION {{{1
    123 if ~ismember({md.groundingline.migration},{'None' 'AgressiveMigration' 'SoftMigration'}),
    124         message(['model not consistent: model ' md.miscellaneous.name ' groundingline.migration field should be ''None'' ''AgressiveMigration'' or ''SoftMigration''']);
    125 end
    126 if ~strcmp(md.groundingline.migration,'None'),
    127         if (md.mesh.dimension==3 | strcmpi(md.cluster.name,'none')),
    128                 message(['model ' md.miscellaneous.name ' requesting grounding line migration, but grounding line module only implemented for 2d models and parallel runs!']);
    129         end
    130         if isnan(md.geometry.bathymetry),
    131                 message(['model not consistent: model ' md.miscellaneous.name ' requesting grounding line migration, but bathymetry is absent!']);
    132         end
    133         pos=find(md.mask.vertexongroundedice);
    134         if any(md.geometry.bed(pos)-md.geometry.bathymetry(pos)),
    135                 message(['model not consistent: model ' md.miscellaneous.name ' bathymetry not equal to bed on grounded ice !']);
    136         end
    137         pos=find(md.mask.vertexonfloatingice);
    138         if any(md.geometry.bathymetry(pos)-md.geometry.bed(pos)>tolerance),
    139                 message(['model not consistent: model ' md.miscellaneous.name ' bathymetry superior to bed on floating ice !']);
    140         end
    141 
    142 end
    143 
    144 %}}}
    145 %RIFTS{{{1
    146 if md.rifts.numrifts,
    147         if ~(md.mesh.dimension==2),
    148                 message(['model not consistent: models with rifts are only supported in 2d for now!']);
    149         end
    150         if ~isstruct(md.rifts.riftstruct),
    151                 message(['model not consistent: md.rifts should be a structure!']);
    152         end
    153         if ~isempty(find(md.mesh.segmentmarkers>=2)),
    154                 %We have segments with rift markers, but no rift structure!
    155                 message(['model not consistent: model ' md.miscellaneous.name ' should be processed for rifts (run meshprocessrifts)!']);
    156         end
    157         %Check that rifts are filled with proper material
    158         checkvalues(md,{'rifts.riftstruct.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
    159 else
    160         if ~isnans(md.rifts.riftstruct),
    161                 message(['model not consistent: md.rifts.riftstruct shoud be NaN since md.rifts.numrifts is 0!']);
    162         end
    163 end
    164 %}}}
    165 %FLAGS (0 or 1){{{1
    166 if ~ismember(md.prognostic.stabilization,[0 1 3]),
    167         message('model not consistent: prognostic.stabilization should be a scalar (0 or 1 or 3)');
    168 end
    169 if ~ismember(md.thermal.stabilization,[0 1 2]),
    170         message('model not consistent: thermal.stabilization should be a scalar (0 or 1 or 2)');
    171 end
    172 if ~ismember(md.balancethickness.stabilization,[0 1 3]),
    173         message('model not consistent: balancethickness.stabilization should be a scalar (0 or 1 or 3)');
    174 end
    175 if ~ismember(md.settings.lowmem,[0 1]),
    176         message(['model not consistent: model ' md.miscellaneous.name ' settings.lowmem field should be 0 or 1']);
    177 end
    178 if ~ismember(md.timestepping.time_adapt,[0 1]),
    179         message(['model not consistent: model ' md.miscellaneous.name ' time_adapt field should be 0 or 1']);
    180 end
    181 if ~ismember(md.prognostic.hydrostatic_adjustment,{'Absolute' 'Incremental'}),
    182         message(['model not consistent: model ' md.miscellaneous.name ' prognostic.hydrostatic_adjustment field should be AbsoluteEnum or IncrementalEnum']);
    183 end
    184 if ~ismember({md.materials.rheology_law},{'None' 'Paterson' 'Arrhenius'}),
    185         message(['model not consistent: model ' md.miscellaneous.name ' rheology_law field should be ''none'' ''paterson'' or ''arrhenius''']);
    186 end
    187 %}}}
    188 %CONNECTIVITY {{{1
    189 if (md.mesh.dimension==2),
    190         if md.mesh.average_vertex_connectivity<9,
    191                 message('model not consistent: connectivity should be at least 9 for 2d models');
    192         end
    193 end
    194 if md.mesh.dimension==3,
    195         if md.mesh.average_vertex_connectivity<24,
    196                 message('model not consistent: connectivity should be at least 24 for 3d models');
    197         end
    198 end
    199 %}}}
    200 %PARALLEL{{{1
    201 IsConsistent(md.cluster);
    202 %}}}
    203 %CONTROL{{{1
    204 if md.inversion.iscontrol,
    205 
    206         %CONTROL TYPE
    207         num_controls=numel(md.inversion.control_parameters);
    208         num_costfunc=size(md.inversion.cost_functions,2);
    209         if ~iscell(md.inversion.control_parameters)
    210                 message(['model not consistent: model ' md.miscellaneous.name ' inversion.control_parameters field should be a cell of strings']);
    211         end
    212         if ~ismember(md.inversion.control_parameters,{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'});
    213                 message(['model not consistent: model ' md.miscellaneous.name ' inversion.control_parameters field should be ''BalancethicknessThickeningRate'' ''FrictionCoefficient'' ''MaterialsRheologyBbar'' ''Vx'' ''Vy''']);
    214         end
    215 
    216         %LENGTH CONTROL FIELDS
    217         fields={'inversion.maxiter_per_step','inversion.step_threshold'};
    218         checksize(md,fields,[md.inversion.nsteps 1]);
    219         fields={'inversion.cost_functions'};
    220         checksize(md,fields,[md.inversion.nsteps num_costfunc]);
    221         fields={'inversion.gradient_scaling'};
    222         checksize(md,fields,[md.inversion.nsteps num_controls]);
    223         fields={'inversion.min_parameters','inversion.max_parameters'};
    224         checksize(md,fields,[md.mesh.numberofvertices num_controls]);
    225 
    226         %RESPONSES
    227         checkvalues(md,{'inversion.cost_functions'},[101:105 201 501:503]);
    228 
    229         %WEIGHTS
    230         fields={'inversion.cost_functions_coefficients'};
    231         checksize(md,fields,[md.mesh.numberofvertices num_costfunc]);
    232         checkgreater(md,fields,0);
    233 
    234         %OBSERVED VELOCITIES
    235         if md.private.solution==BalancethicknessSolutionEnum
    236                 fields={'inversion.thickness_obs'};
    237                 checksize(md,fields,[md.mesh.numberofvertices 1]);
    238                 checknan(md,fields);
    239         else
    240                 fields={'inversion.vx_obs','inversion.vy_obs'};
    241                 checksize(md,fields,[md.mesh.numberofvertices 1]);
    242                 checknan(md,fields);
    243         end
    244 
    245         %DIRICHLET IF THICKNESS <= 0
    246         if any(md.geometry.thickness<=0),
    247                 pos=find(md.geometry.thickness<=0);
    248                 if any(isnan(md.balancethickness.spcthickness(pos))),
    249                         message(['model not consistent: model ' md.miscellaneous.name ' has some nodes with 0 thickness']);
    250                 end
    251         end
    252 end
    253 %}}}
    254 %QMU {{{1
    255 if md.qmu.isdakota,
    256         if md.qmu.params.evaluation_concurrency~=1,
    257                 message(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
    258         end
    259         if ~isempty(md.qmu.partition),
    260                 if numel(md.qmu.partition)~=md.mesh.numberofvertices,
    261                         message(['model not consistent: user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1 ']);
    262                 end
    263                 if find(md.qmu.partition)>=md.mesh.numberofvertices,
    264                         message(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
    265                 end
    266                 if min(md.qmu.partition)~=0,
    267                         message(['model not consistent: partition vector not indexed from 0 on']);
    268                 end
    269                 if max(md.qmu.partition)>=md.mesh.numberofvertices,
    270                         message(['model not consistent: partition vector cannot have maximum index larger than number of nodes']);
    271                 end
    272                 if ~isempty(find(md.qmu.partition<0)),
    273                         message(['model not consistent: partition vector cannot have values less than 0']);
    274                 end
    275                 if ~isempty(find(md.qmu.partition>=md.qmu.numberofpartitions)),
    276                         message(['model not consistent: partition vector cannot have values more than md.qmu.numberofpartitions-1']);
    277                 end
    278                 if max(md.qmu.partition)>=md.qmu.numberofpartitions,
    279                         message(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']);
    280                 end
    281         end
    282 end
    283 
    284 if strcmpi(md.private.solution,'qmu'),
    285         if ~strcmpi(md.cluster.name,'none'),
    286                 if md.settings.waitonlock==0,
    287                         message(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
    288                 end
    289         end
    290 end
    291 %}}}
    292 
    293 %Solution specific check
    294 %TRANSIENT {{{1
    295 if (md.private.solution==TransientSolutionEnum),
    296 
    297         if md.timestepping.time_step<=0,
    298                 message('model not consistent: field timestepping.time_step must be positive for a transient run')
    299         end
    300         if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0),
    301                 message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between  0 and 1']);
    302         end
    303         if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0),
    304                 message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between  0 and 1']);
    305         end
    306         if ~ismember(md.transient.isdiagnostic,[0 1]),
    307                 message('model not consistent: transient.isdiagnostic should be a scalar (1 or 0)');
    308         end
    309         if ~ismember(md.transient.isprognostic,[0 1]),
    310                 message('model not consistent: transient.isprognostic should be a scalar (1 or 0)');
    311         end
    312         if ~ismember(md.transient.isthermal,[0 1]),
    313                 message('model not consistent: transient.isthermal should be a scalar (1 or 0)');
    314         end
    315         if ~ismember(md.transient.isgroundingline,[0 1]),
    316                 message('model not consistent: transient.isgroundingline should be a scalar (1 or 0)');
    317         end
    318 end
    319 %}}}
    320 %STEADYSTATE{{{1
    321 if md.private.solution==SteadystateSolutionEnum,
    322 
    323         %NDT
    324         if md.timestepping.time_step~=0,
    325                 message(['model not consistent: for a steadystate computation, timestepping.time_step must be zero.']);
    326         end
    327 
    328         %eps:
    329         if isnan(md.diagnostic.reltol),
    330                 message(['model not consistent: for a steadystate computation, diagnostic.reltol (relative convergence criterion) must be defined!']);
    331         end
    332 end
    333 %}}}
    334 %GROUNDINGLINEMIGRATION2D{{{1
    335 if md.private.solution==GroundinglineMigration2dSolutionEnum,
    336         if strcmpi(md.cluster.name,'none'),
    337                 message(['model not consistent: ' md.private.solution ' is only implemented in parallel mode !'])
    338         end
    339 
    340         if md.timestepping.time_step<=0,
    341                 message('model not consistent: field timestepping.time_step must be positive for a transient run')
    342         end
    343 
    344         if (md.mesh.dimension~=2),
    345                 message(['model not consistent: for a ' md.private.solution ' computation, the grounding line module is only implemented in 2d !'])
    346         end
    347 
    348         if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0),
    349                 message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between  0 and 1']);
    350         end
    351 end
    352 %}}}
    353 
    354 %Now check all analyses called for a given solution
    355 %ANALYSESCHECKS {{{1
    356 [analyses,numanalyses]=AnalysisConfiguration(md.private.solution);
    357 for i=1:numanalyses,
    358         checkforanalysis(md,analyses(i));
    359 end
    360 %}}}
    361 
    362         if modelconsistency==false, error(['model not consistent']); end
    363 end
    364 
    365 %analysis checks
    366 %checkforanalysis {{{1
    367 function checkforanalysis(md,analysis_enum)
    368 
    369         switch(analysis_enum),
    370                 case DiagnosticHorizAnalysisEnum,
    371                         % {{{2
    372                         %SINGULAR
    373                         if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2),
    374                                 message(['model not consistent: model ' md.miscellaneous.name ' is not well posed (singular). You need at least one node with fixed velocity!'])
    375                         end
    376 
    377                         %ROTATED SPC
    378                         %CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES
    379                         if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6),
    380                                 message(['model not consistent: model ' md.miscellaneous.name ' has problem with rotated spc. Each line of diagnostic.referential should contain either only NaN values or no NaN values']);
    381                         end
    382                         %CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL
    383                         if any(sum(isnan(md.diagnostic.referential),2)==0),
    384                                 pos=find(sum(isnan(md.diagnostic.referential),2)==0);
    385                                 if any(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))),
    386                                         dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))
    387                                         message(['model not consistent: model ' md.miscellaneous.name ' has problem with rotated spc. Vectors in diagnostic.referential (colums 1 to 3 and 4 to 6) must be orthogonal']);
    388                                 end
    389 
    390                         end
    391                         %CHECK THAT ROTATION IS IN THE (X,Y) PLANE FOR 2D MODELS
    392                         if md.mesh.dimension==2,
    393                                 pos=find(sum(isnan(md.diagnostic.referential),2)==0  & md.flowequation.vertex_equation==2);
    394                                 if any(md.diagnostic.referential(pos,3:5)~=0);
    395                                         message(['model not consistent: model ' md.miscellaneous.name ' has problem with rotated spc. The rotation should be in the (x,y) plane for 2D diagnostic models (nodeonmacayeal)']);
    396                                 end
    397                         end
    398 
    399                         %INITIAL VELOCITY
    400                         if ~isnan(md.initialization.vx) & ~isnan(md.initialization.vy),
    401                                 fields={'initialization.vx','initialization.vy'};
    402                                 checknan(md,fields);
    403                                 checksize(md,fields,[md.mesh.numberofvertices 1]);
    404                         end
    405 
    406                         %ELEMENTSTYPE
    407                         %Check the size of flowequation.element_equation
    408                         fields={'flowequation.element_equation'};
    409                         checksize(md,fields,[md.mesh.numberofelements 1]);
    410                         %Check the values of flowequation.element_equation
    411                         checkvalues(md,{'flowequation.element_equation'},[0:7]);
    412                         if (md.mesh.dimension==2),
    413                                 checkvalues(md,{'flowequation.element_equation'},[1 2]);
    414                         end
    415                         if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0),
    416                                 message(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
    417                         end
    418                        
    419                         %VERTICESTYPE
    420                         %Check the size of verticess_type
    421                         fields={'flowequation.vertex_equation'};
    422                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    423                         %Check the values of flowequation.vertex_equation
    424                         checkvalues(md,{'flowequation.vertex_equation'},[0:7]);
    425                         if (md.mesh.dimension==2),
    426                                 checkvalues(md,{'flowequation.vertex_equation'},[1 2]);
    427                         end
    428                         if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0),
    429                                 message(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
    430                         end
    431                         %}}}
    432                 case DiagnosticVertAnalysisEnum,
    433                         % {{{2
    434                         if md.mesh.dimension==3,
    435                                 % No checks for now
    436                         end
    437                         %}}}
    438                 case DiagnosticHutterAnalysisEnum,
    439                         % {{{2
    440                         %HUTTER ON ICESHELF WARNING
    441                         if any(md.flowequation.element_equation==1 & md.mask.elementonfloatingice),
    442                                 disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
    443                         end
    444                         %}}}
    445                 case BedSlopeAnalysisEnum,
    446                         % {{{2
    447                         % No checks for now
    448                         %}}}
    449                 case SurfaceSlopeAnalysisEnum,
    450                         % {{{2
    451                         % No checks for now
    452                         %}}}
    453                 case PrognosticAnalysisEnum,
    454                         % {{{2
    455                         %INITIAL VELOCITIES
    456                         fields={'initialization.vx','initialization.vy'};
    457                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    458                         checknan(md,fields);
    459 
    460                         %FORCINGS
    461                         fields={'surfaceforcings.mass_balance'};
    462                         checkforcing(md,fields);
    463 
    464                         %DIRICHLET
    465                         fields={'prognostic.spcthickness'};
    466                         checkforcing(md,fields);
    467 
    468                         %DIRICHLET IF THICKNESS <= 0
    469                         %Check the size of prognostic.spcthickness
    470                         if any(md.geometry.thickness<=0),
    471                                 pos=find(md.geometry.thickness<=0);
    472                                 if any(isnan(md.prognostic.spcthickness(pos))),
    473                                         message(['model not consistent: model ' md.miscellaneous.name ' has some nodes with 0 thickness']);
    474                                 end
    475                         end
    476 
    477                         %}}}
    478                 case HydrologyAnalysisEnum,
    479                         % {{{2
    480                         fields={'hydrology.spcwatercolumn'};
    481                         checkforcing(md,fields);
    482                         fields={'initialization.watercolumn'};
    483                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    484                         %}}}
    485                 case ThermalAnalysisEnum,
    486                         % {{{2
    487                         %EXTRUSION
    488                         if (md.mesh.dimension==2),
    489                                 if md.private.solution==TransientSolutionEnum,
    490                                         return;
    491                                 else
    492                                         message(['model not consistent: for a ' EnumToString(md.private.solution) ' computation, the model must be 3d, extrude it first!'])
    493                                 end
    494                         end
    495 
    496                         %CHECK THAT SPCTEMPERATURE IS AN APPROPRIATE FORCING
    497                         fields={'thermal.spctemperature'};
    498                         checkforcing(md,fields);
    499 
    500                         %CHECK THAT WE ARE NOT FULLY CONSTRAINED
    501                         if ~any(~isnan(md.thermal.spctemperature))
    502                                 message(['model not consistent: model ' md.miscellaneous.name ' is totally constrained for temperature, no need to solve!']);
    503                         end
    504 
    505                         %VELOCITIES AND PRESSURE
    506                         fields={'initialization.vx','initialization.vy','initialization.vz','initialization.pressure','basalforcings.geothermalflux'};
    507                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    508                         checknan(md,fields);
    509 
    510                         %THERMAL TRANSIENT
    511                         if md.timestepping.time_step~=0,
    512 
    513                                 %DT and NDT
    514                                 fields={'timestepping.time_step','timestepping.final_time'};
    515                                 checkgreaterstrict(md,fields,0);
    516 
    517                                 %INITIAL TEMPERATURE, MELTING AND ACCUMULATION
    518                                 fields={'initialization.temperature','basalforcings.melting_rate'};
    519                                 checksize(md,fields,[md.mesh.numberofvertices 1]);
    520                                 checknan(md,fields);
    521 
    522                                 %INITIAL TEMPERATURE
    523                                 fields={'initialization.temperature'};
    524                                 checkgreater(md,fields,0)
    525 
    526                                 %CHECK SPCTEMPERATURE that are not NaN are >0.
    527                                 if find(any(md.thermal.spctemperature(find(~isnan(md.thermal.spctemperature(1:md.mesh.numberofvertices,:))))<=0)),
    528                                         message(['model not consistent: model ' md.miscellaneous.name ' is constrained with negative or nil temperatures!']);
    529                                 end
    530 
    531                         end
    532                         %}}}
    533                 case EnthalpyAnalysisEnum,
    534                         % {{{2
    535                         %EXTRUSION
    536                         if (md.mesh.dimension==2),
    537                                 if md.private.solution==TransientSolutionEnum,
    538                                         return;
    539                                 else
    540                                         message(['model not consistent: for a ' EnumToString(md.private.solution) ' computation, the model must be 3d, extrude it first!'])
    541                                 end
    542                         end
    543 
    544                         %VELOCITIES AND PRESSURE
    545                         fields={'initialization.vx','initialization.vy','initialization.vz','initialization.pressure','basalforcings.geothermalflux'};
    546                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    547                         checknan(md,fields);
    548 
    549                         %THERMAL TRANSIENT
    550                         if md.timestepping.time_step~=0,
    551 
    552                                 %TIMESTEPPING
    553                                 fields={'timestepping.time_step','timestepping.final_time'};
    554                                 checkgreaterstrict(md,fields,0);
    555 
    556                                 %INITIAL TEMPERATURE, MELTING, ACCUMULATION AND WATERFRACTION
    557                                 fields={'initialization.temperature','basalforcings.melting_rate','initialization.waterfraction'};
    558                                 checksize(md,fields,[md.mesh.numberofvertices 1]);
    559                                 checknan(md,fields);
    560 
    561                                 %INITIAL TEMPERATURE
    562                                 fields={'initialization.temperature'};
    563                                 checkgreater(md,fields,0)
    564                         end
    565                         %}}}
    566                 case MeltingAnalysisEnum,
    567                         % {{{2
    568                         % No checks for now
    569                         %}}}
    570                 case BalancethicknessAnalysisEnum,
    571                         % {{{2
    572                         %VELOCITIES MELTING AND ACCUMULATION
    573                         fields={'initialization.vx','initialization.vy','basalforcings.melting_rate','balancethickness.thickening_rate'};
    574                         checksize(md,fields,[md.mesh.numberofvertices 1]);
    575                         checknan(md,fields);
    576 
    577                         %Triangle with zero velocity
    578                         if any(sum(abs(md.initialization.vx(md.elements)),2)==0 & sum(abs(md.initialization.vy(md.elements)),2)==0)
    579                                 message('model not consistent: at least one triangle has all its vertices with a zero velocity');
    580                         end
    581                         %}}}
    582                 case FlaimAnalysisEnum,
    583                         % {{{2
    584                         % No checks for now
    585                         %}}}
    586                 otherwise
    587                         disp(['WARNING: no check implemented for analysis ' EnumToString(analysis_enum) '!']);
    588                         message('stop')
    589         end
    590 
    591 end
    592 %}}}
    593 
    594 %checks additional functions
    595 %checklength {{{1
    596 function checklength(md,fields,fieldlength)
    597         %CHECKSIZE - check length of a field
    598         for i=1:length(fields),
    599                 eval(['field=md.' fields{i} ';']);
    600                 if length(field)~=fieldlength,
    601                         message(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
    602                 end
    603         end
    604 end
    605 %}}}
    606 %checksize {{{1
    607 function checksize(md,fields,fieldsize)
    608         %CHECKSIZE - check size of a field
    609         for i=1:length(fields),
    610                 eval(['field=md.' fields{i} ';']);
    611                 if isnan(fieldsize(1)),
    612                         if (size(field,2)~=fieldsize(2)),
    613                                 message(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
    614                         end
    615                 elseif isnan(fieldsize(2)),
    616                         if (size(field,1)~=fieldsize(1)),
    617                                 message(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
    618                         end
    619                 else
    620                         if ((size(field)~=fieldsize(1)) |  (size(field,2)~=fieldsize(2)))
    621                                 %WARNING: new version
    622                                 if strcmp(fields{i},'cm_min'),
    623                                         disp('!!! ');
    624                                         disp('!!! WARNING: cm_min must now be of size [md.mesh.numberofvertices x 1]. Update your parameter file as follows:');
    625                                         disp('!!! ');
    626                                         disp('!!! md.inversion.min_parameters=md.inversion.min_parameters*ones(md.mesh.numberofvertices,1);');
    627                                         disp('!!! ');
    628                                 end
    629                                 %WARNING: new version
    630                                 if strcmp(fields{i},'cm_max'),
    631                                         disp('!!! ');
    632                                         disp('!!! WARNING: cm_max must now be of size [md.mesh.numberofvertices x 1]. Update your parameter file as follows:');
    633                                         disp('!!! ');
    634                                         disp('!!! md.inversion.max_parameters=md.inversion.max_parameters*ones(md.mesh.numberofvertices,1);');
    635                                         disp('!!! ');
    636                                 end
    637                                 message(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
    638                         end
    639                 end
    640         end
    641 end
    642 %}}}
    643 %checknan {{{1
    644 function checknan(md,fields)
    645         %CHECKNAN - check nan values of a field
    646         for i=1:length(fields),
    647                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    648                 if any(isnan(field)),
    649                         message(['model not consistent: NaN values found in field ' fields{i}]);
    650                 end
    651         end
    652 end
    653 %}}}
    654 %checkreal{{{1
    655 function checkreal(md,fields)
    656         %CHECKREAL - check real values of a field
    657         for i=1:length(fields),
    658                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    659                 if any(~isreal(field)),
    660                         message(['model not consistent: complex values found in field ' fields{i}]);
    661                 end
    662         end
    663 end
    664 %}}}
    665 %checkgreaterstrict{{{1
    666 function checkgreaterstrict(md,fields,lowerbound)
    667         %CHECKGREATERSTRICT - check values of a field
    668         for i=1:length(fields),
    669                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    670                 if any(field<=lowerbound),
    671                         message(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
    672                 end
    673         end
    674 end
    675 %}}}
    676 %checkgreater{{{1
    677 function checkgreater(md,fields,lowerbound)
    678         %CHECKGREATER - check values of a field
    679         for i=1:length(fields),
    680                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    681                 if any(field<lowerbound),
    682                         message(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
    683                 end
    684         end
    685 end
    686 %}}}
    687 %checklessstrict{{{1
    688 function checklessstrict(md,fields,upperbound)
    689         %CHECKLESSSTRICT - check values of a field
    690         for i=1:length(fields),
    691                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    692                 if any(field>=upperbound),
    693                         message(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
    694                 end
    695         end
    696 end
    697 %}}}
    698 %checkless{{{1
    699 function checkless(md,fields,upperbound)
    700         %CHECKLESS - check values of a field
    701         for i=1:length(fields),
    702                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    703                 if any(field>upperbound),
    704                         message(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
    705                 end
    706         end
    707 end
    708 %}}}
    709 %checkvalues {{{1
    710 function checkvalues(md,fields,values)
    711         %CHECKVALUE - check that a field has specified values
    712         for i=1:length(fields),
    713                 eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']);
    714                 if any(~ismember(field,values)),
    715                         message(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
    716                 end
    717         end
    718 end
    719 %}}}
    720 %checkforcing {{{1
    721 function checkforcing(md,fields)
    722 
    723         for i=1:length(fields),
    724                 eval(['field=md.' fields{i} ';']);
    725                 if size(field,1)==md.mesh.numberofvertices,
    726                         if ~size(field,2)==1,
    727                                 message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' should have only one column as there are md.mesh.numberofvertices lines']);
    728                         end
    729                 elseif size(field,1)==md.mesh.numberofvertices+1
    730                         if any(field(end,:)~=sort(field(end,:))),
    731                                 message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' columns should be chronological']);
    732                         end
    733                         if any(field(end,1:end-1)==field(end,2:end)),
    734                                 message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' columns must not contain duplicate timesteps']);
    735                         end
    736                 else
    737                         message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' should have md.mesh.numberofvertices or md.mesh.numberofvertices+1 lines']);
    738                 end
    739         end
    740 end
    741 %}}}
    742 
    743 %error messages
    744 %modelconsistency{{{1
    745 function flag=modelconsistency(flag_in)
    746 
    747         persistent consistency;
    748 
    749         if nargin==1 & nargout==0,
    750                 %OK model is inconsistent, set flag as false
    751                 consistency=flag_in;
    752         elseif nargin==0 & nargout==1,
    753                 if isempty(consistency),
    754                         %modelinconsistent has never been called, model is consistent
    755                         consistency=true;
    756                 end
    757         else
    758                 message('Bad usage');
    759         end
    760 
    761         flag=consistency;
    762 end
    763 %}}}
    764 %message{{{1
    765 function message(string)
    766 
    767         disp(string);
    768         modelconsistency(false);
    769 end
    770 %}}}
     7%nothing for now
Note: See TracChangeset for help on using the changeset viewer.