Index: /issm/trunk/src/m/model/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 9725) +++ /issm/trunk/src/m/model/ismodelselfconsistent.m (revision 9726) @@ -5,766 +5,3 @@ % ismodelselfconsistent(md), -%tolerance we use in litmus checks for the consistency of the model -tolerance=10^-9; -%check usage {{{1 -if nargin~=1, - help ismodelselfconsistent - message('ismodelselfconsistent message message: wrong usage'); -end -%}}} - -%initialize consistency as true -modelconsistency(true); - -% Flaim doesn't need the entire model, so take an early return -%FLAIM {{{1 -if (md.private.solution == FlaimSolutionEnum), - if ~exist(md.flaim.tracks,'file') - message(['model not consistent: flaim.tracks file ''' md.flaim.tracks ''' must exist.']); - end - % probably going to need some checks on flaim.flightreqs here - if (numel(md.flaim.criterion) ~= md.mesh.numberofvertices) && (numel(md.flaim.criterion) ~= md.mesh.numberofelements) - 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.']); - end - return; -end -%}}} - -% Common checks -%VERBOSE{{{1 -if ~strcmp('verbose',class(md.verbose)) - disp('WARNING: md.verbose must now be an object of class ''verbose'''); - disp(' To update your model use the following command:'); - disp(' '); - disp(' md.verbose=verbose;'); - disp(' '); - message(['field verbose should be of class ''verbose'' ']); -end -%}}} -%NAME{{{1 -if isempty(md.miscellaneous.name), - message(['model is not correctly configured: missing name!']); -end -%}}} -%ELEMENTS{{{1 -fields={'elements'}; -if (md.mesh.dimension==2), - checksize(md,fields,[md.mesh.numberofelements 3]); -else - checksize(md,fields,[md.mesh.numberofelements 6]); -end -if any(~ismember(1:md.mesh.numberofvertices,sort(unique(md.elements(:))))); - message('orphan nodes have been found. Check the mesh'); -end -%}}} - -%DG {{{1 -if md.prognostic.stabilization==3; - %CHECK THE COMPATIBILITY OF THE EDGES - fields={'edges'}; - checksize(md,fields,[NaN 4]); - fields={'edges(:,1:3)'}; - checknan(md,fields); - checkgreaterstrict(md,fields,0); -end -%}}} -%PRESSURELOAD{{{1 -if (md.mesh.dimension==2), - fields={'diagnostic.icefront'}; - checksize(md,fields,[NaN 4]); -elseif md.mesh.dimension==3, - fields={'diagnostic.icefront'}; - checksize(md,fields,[NaN 6]); -else - message('dim should be either 2 3'); -end -checkvalues(md,{'diagnostic.icefront(:,end)'},[0 1 2]); -%}}} -%NO NAN {{{1 -fields={'mesh.numberofelements','mesh.numberofvertices','x','y','z','friction.coefficient','friction.p','friction.q',... - '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',... - 'diagnostic.restol','diagnostic.maxiter','materials.rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface','elementconnectivity'}; -checknan(md,fields); -%}}}} -%FIELDS >= 0 {{{1 -fields={'mesh.numberofelements','mesh.numberofvertices','elements','friction.coefficient','friction.p','friction.q',... - 'materials.rho_ice','materials.rho_water','materials.rheology_B','mask.elementonfloatingice','geometry.thickness','constants.g','diagnostic.restol','diagnostic.maxiter','diagnostic.reltol',... - 'diagnostic.abstol','inversion.nsteps','inversion.maxiter_per_step',... - 'settings.lowmem','materials.rheology_n','nodeonbed','nodeonsurface','elementonbed','elementonsurface'}; -checkgreater(md,fields,0); -%}}} -%FIELDS > 0 {{{1 -fields={'mesh.numberofelements','mesh.numberofvertices','elements','friction.p',... - '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'}; -checkgreaterstrict(md,fields,0); -%}}} -%SIZE NUMBEROFELEMENTS {{{1 -fields={'friction.p','friction.q','mask.elementonfloatingice','materials.rheology_n','elementonbed'}; -checksize(md,fields,[md.mesh.numberofelements 1]); -%}}} -%SIZE NUMBEROFNODES {{{1 -fields={'x','y','z','materials.rheology_B','friction.coefficient','basalforcings.melting_rate','geometry.surface','geometry.thickness','geometry.bed','nodeonbed','nodeonsurface'}; -checksize(md,fields,[md.mesh.numberofvertices 1]); -%}}} -%OTHER SIZES {{{1 -fields={'diagnostic.referential'}; -checksize(md,fields,[md.mesh.numberofvertices 6]); -if ~isempty(md.diagnostic.requested_outputs), - if(size(md.diagnostic.requested_outputs,2)~=1), - message(['model ' md.miscellaneous.name ' requested outputs should be a column vector']); - end -end -%}}} -%THICKNESS = SURFACE - BED {{{1 -if any((md.geometry.thickness-md.geometry.surface+md.geometry.bed)>tolerance), - message(['model not consistent: model ' md.miscellaneous.name ' violates the equality thickness=surface-bed!']); -end -%GROUNDING LINE MIGRATION {{{1 -if ~ismember({md.groundingline.migration},{'None' 'AgressiveMigration' 'SoftMigration'}), - message(['model not consistent: model ' md.miscellaneous.name ' groundingline.migration field should be ''None'' ''AgressiveMigration'' or ''SoftMigration''']); -end -if ~strcmp(md.groundingline.migration,'None'), - if (md.mesh.dimension==3 | strcmpi(md.cluster.name,'none')), - message(['model ' md.miscellaneous.name ' requesting grounding line migration, but grounding line module only implemented for 2d models and parallel runs!']); - end - if isnan(md.geometry.bathymetry), - message(['model not consistent: model ' md.miscellaneous.name ' requesting grounding line migration, but bathymetry is absent!']); - end - pos=find(md.mask.vertexongroundedice); - if any(md.geometry.bed(pos)-md.geometry.bathymetry(pos)), - message(['model not consistent: model ' md.miscellaneous.name ' bathymetry not equal to bed on grounded ice !']); - end - pos=find(md.mask.vertexonfloatingice); - if any(md.geometry.bathymetry(pos)-md.geometry.bed(pos)>tolerance), - message(['model not consistent: model ' md.miscellaneous.name ' bathymetry superior to bed on floating ice !']); - end - -end - -%}}} -%RIFTS{{{1 -if md.rifts.numrifts, - if ~(md.mesh.dimension==2), - message(['model not consistent: models with rifts are only supported in 2d for now!']); - end - if ~isstruct(md.rifts.riftstruct), - message(['model not consistent: md.rifts should be a structure!']); - end - if ~isempty(find(md.mesh.segmentmarkers>=2)), - %We have segments with rift markers, but no rift structure! - message(['model not consistent: model ' md.miscellaneous.name ' should be processed for rifts (run meshprocessrifts)!']); - end - %Check that rifts are filled with proper material - checkvalues(md,{'rifts.riftstruct.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]); -else - if ~isnans(md.rifts.riftstruct), - message(['model not consistent: md.rifts.riftstruct shoud be NaN since md.rifts.numrifts is 0!']); - end -end -%}}} -%FLAGS (0 or 1){{{1 -if ~ismember(md.prognostic.stabilization,[0 1 3]), - message('model not consistent: prognostic.stabilization should be a scalar (0 or 1 or 3)'); -end -if ~ismember(md.thermal.stabilization,[0 1 2]), - message('model not consistent: thermal.stabilization should be a scalar (0 or 1 or 2)'); -end -if ~ismember(md.balancethickness.stabilization,[0 1 3]), - message('model not consistent: balancethickness.stabilization should be a scalar (0 or 1 or 3)'); -end -if ~ismember(md.settings.lowmem,[0 1]), - message(['model not consistent: model ' md.miscellaneous.name ' settings.lowmem field should be 0 or 1']); -end -if ~ismember(md.timestepping.time_adapt,[0 1]), - message(['model not consistent: model ' md.miscellaneous.name ' time_adapt field should be 0 or 1']); -end -if ~ismember(md.prognostic.hydrostatic_adjustment,{'Absolute' 'Incremental'}), - message(['model not consistent: model ' md.miscellaneous.name ' prognostic.hydrostatic_adjustment field should be AbsoluteEnum or IncrementalEnum']); -end -if ~ismember({md.materials.rheology_law},{'None' 'Paterson' 'Arrhenius'}), - message(['model not consistent: model ' md.miscellaneous.name ' rheology_law field should be ''none'' ''paterson'' or ''arrhenius''']); -end -%}}} -%CONNECTIVITY {{{1 -if (md.mesh.dimension==2), - if md.mesh.average_vertex_connectivity<9, - message('model not consistent: connectivity should be at least 9 for 2d models'); - end -end -if md.mesh.dimension==3, - if md.mesh.average_vertex_connectivity<24, - message('model not consistent: connectivity should be at least 24 for 3d models'); - end -end -%}}} -%PARALLEL{{{1 -IsConsistent(md.cluster); -%}}} -%CONTROL{{{1 -if md.inversion.iscontrol, - - %CONTROL TYPE - num_controls=numel(md.inversion.control_parameters); - num_costfunc=size(md.inversion.cost_functions,2); - if ~iscell(md.inversion.control_parameters) - message(['model not consistent: model ' md.miscellaneous.name ' inversion.control_parameters field should be a cell of strings']); - end - if ~ismember(md.inversion.control_parameters,{'BalancethicknessThickeningRate' 'FrictionCoefficient' 'MaterialsRheologyBbar' 'Vx' 'Vy'}); - message(['model not consistent: model ' md.miscellaneous.name ' inversion.control_parameters field should be ''BalancethicknessThickeningRate'' ''FrictionCoefficient'' ''MaterialsRheologyBbar'' ''Vx'' ''Vy''']); - end - - %LENGTH CONTROL FIELDS - fields={'inversion.maxiter_per_step','inversion.step_threshold'}; - checksize(md,fields,[md.inversion.nsteps 1]); - fields={'inversion.cost_functions'}; - checksize(md,fields,[md.inversion.nsteps num_costfunc]); - fields={'inversion.gradient_scaling'}; - checksize(md,fields,[md.inversion.nsteps num_controls]); - fields={'inversion.min_parameters','inversion.max_parameters'}; - checksize(md,fields,[md.mesh.numberofvertices num_controls]); - - %RESPONSES - checkvalues(md,{'inversion.cost_functions'},[101:105 201 501:503]); - - %WEIGHTS - fields={'inversion.cost_functions_coefficients'}; - checksize(md,fields,[md.mesh.numberofvertices num_costfunc]); - checkgreater(md,fields,0); - - %OBSERVED VELOCITIES - if md.private.solution==BalancethicknessSolutionEnum - fields={'inversion.thickness_obs'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - else - fields={'inversion.vx_obs','inversion.vy_obs'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - end - - %DIRICHLET IF THICKNESS <= 0 - if any(md.geometry.thickness<=0), - pos=find(md.geometry.thickness<=0); - if any(isnan(md.balancethickness.spcthickness(pos))), - message(['model not consistent: model ' md.miscellaneous.name ' has some nodes with 0 thickness']); - end - end -end -%}}} -%QMU {{{1 -if md.qmu.isdakota, - if md.qmu.params.evaluation_concurrency~=1, - message(['model not consistent: concurrency should be set to 1 when running dakota in library mode']); - end - if ~isempty(md.qmu.partition), - if numel(md.qmu.partition)~=md.mesh.numberofvertices, - message(['model not consistent: user supplied partition for qmu analysis should have size md.mesh.numberofvertices x 1 ']); - end - if find(md.qmu.partition)>=md.mesh.numberofvertices, - message(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']); - end - if min(md.qmu.partition)~=0, - message(['model not consistent: partition vector not indexed from 0 on']); - end - if max(md.qmu.partition)>=md.mesh.numberofvertices, - message(['model not consistent: partition vector cannot have maximum index larger than number of nodes']); - end - if ~isempty(find(md.qmu.partition<0)), - message(['model not consistent: partition vector cannot have values less than 0']); - end - if ~isempty(find(md.qmu.partition>=md.qmu.numberofpartitions)), - message(['model not consistent: partition vector cannot have values more than md.qmu.numberofpartitions-1']); - end - if max(md.qmu.partition)>=md.qmu.numberofpartitions, - message(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']); - end - end -end - -if strcmpi(md.private.solution,'qmu'), - if ~strcmpi(md.cluster.name,'none'), - if md.settings.waitonlock==0, - message(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']); - end - end -end -%}}} - -%Solution specific check -%TRANSIENT {{{1 -if (md.private.solution==TransientSolutionEnum), - - if md.timestepping.time_step<=0, - message('model not consistent: field timestepping.time_step must be positive for a transient run') - end - if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0), - message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between 0 and 1']); - end - if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0), - message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between 0 and 1']); - end - if ~ismember(md.transient.isdiagnostic,[0 1]), - message('model not consistent: transient.isdiagnostic should be a scalar (1 or 0)'); - end - if ~ismember(md.transient.isprognostic,[0 1]), - message('model not consistent: transient.isprognostic should be a scalar (1 or 0)'); - end - if ~ismember(md.transient.isthermal,[0 1]), - message('model not consistent: transient.isthermal should be a scalar (1 or 0)'); - end - if ~ismember(md.transient.isgroundingline,[0 1]), - message('model not consistent: transient.isgroundingline should be a scalar (1 or 0)'); - end -end -%}}} -%STEADYSTATE{{{1 -if md.private.solution==SteadystateSolutionEnum, - - %NDT - if md.timestepping.time_step~=0, - message(['model not consistent: for a steadystate computation, timestepping.time_step must be zero.']); - end - - %eps: - if isnan(md.diagnostic.reltol), - message(['model not consistent: for a steadystate computation, diagnostic.reltol (relative convergence criterion) must be defined!']); - end -end -%}}} -%GROUNDINGLINEMIGRATION2D{{{1 -if md.private.solution==GroundinglineMigration2dSolutionEnum, - if strcmpi(md.cluster.name,'none'), - message(['model not consistent: ' md.private.solution ' is only implemented in parallel mode !']) - end - - if md.timestepping.time_step<=0, - message('model not consistent: field timestepping.time_step must be positive for a transient run') - end - - if (md.mesh.dimension~=2), - message(['model not consistent: for a ' md.private.solution ' computation, the grounding line module is only implemented in 2d !']) - end - - if(md.timestepping.cfl_coefficient>1 | md.timestepping.cfl_coefficient<0), - message(['model not consistent: model ' md.miscellaneous.name ' cfl_coefficient field should between 0 and 1']); - end -end -%}}} - -%Now check all analyses called for a given solution -%ANALYSESCHECKS {{{1 -[analyses,numanalyses]=AnalysisConfiguration(md.private.solution); -for i=1:numanalyses, - checkforanalysis(md,analyses(i)); -end -%}}} - - if modelconsistency==false, error(['model not consistent']); end -end - -%analysis checks -%checkforanalysis {{{1 -function checkforanalysis(md,analysis_enum) - - switch(analysis_enum), - case DiagnosticHorizAnalysisEnum, - % {{{2 - %SINGULAR - if ~any((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy))==2), - message(['model not consistent: model ' md.miscellaneous.name ' is not well posed (singular). You need at least one node with fixed velocity!']) - end - - %ROTATED SPC - %CHECK THAT EACH LINES CONTAINS ONLY NAN VALUES OR NO NAN VALUES - if any(sum(isnan(md.diagnostic.referential),2)~=0 & sum(isnan(md.diagnostic.referential),2)~=6), - 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']); - end - %CHECK THAT THE TWO VECTORS PROVIDED ARE ORTHOGONAL - if any(sum(isnan(md.diagnostic.referential),2)==0), - pos=find(sum(isnan(md.diagnostic.referential),2)==0); - if any(dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6))), - dot(md.diagnostic.referential(pos,1:3),md.diagnostic.referential(pos,4:6)) - 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']); - end - - end - %CHECK THAT ROTATION IS IN THE (X,Y) PLANE FOR 2D MODELS - if md.mesh.dimension==2, - pos=find(sum(isnan(md.diagnostic.referential),2)==0 & md.flowequation.vertex_equation==2); - if any(md.diagnostic.referential(pos,3:5)~=0); - 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)']); - end - end - - %INITIAL VELOCITY - if ~isnan(md.initialization.vx) & ~isnan(md.initialization.vy), - fields={'initialization.vx','initialization.vy'}; - checknan(md,fields); - checksize(md,fields,[md.mesh.numberofvertices 1]); - end - - %ELEMENTSTYPE - %Check the size of flowequation.element_equation - fields={'flowequation.element_equation'}; - checksize(md,fields,[md.mesh.numberofelements 1]); - %Check the values of flowequation.element_equation - checkvalues(md,{'flowequation.element_equation'},[0:7]); - if (md.mesh.dimension==2), - checkvalues(md,{'flowequation.element_equation'},[1 2]); - end - if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0), - message(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']); - end - - %VERTICESTYPE - %Check the size of verticess_type - fields={'flowequation.vertex_equation'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - %Check the values of flowequation.vertex_equation - checkvalues(md,{'flowequation.vertex_equation'},[0:7]); - if (md.mesh.dimension==2), - checkvalues(md,{'flowequation.vertex_equation'},[1 2]); - end - if (md.flowequation.ismacayealpattyn==0 && md.flowequation.ishutter==0 && md.flowequation.isstokes==0), - message(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']); - end - %}}} - case DiagnosticVertAnalysisEnum, - % {{{2 - if md.mesh.dimension==3, - % No checks for now - end - %}}} - case DiagnosticHutterAnalysisEnum, - % {{{2 - %HUTTER ON ICESHELF WARNING - if any(md.flowequation.element_equation==1 & md.mask.elementonfloatingice), - disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n')); - end - %}}} - case BedSlopeAnalysisEnum, - % {{{2 - % No checks for now - %}}} - case SurfaceSlopeAnalysisEnum, - % {{{2 - % No checks for now - %}}} - case PrognosticAnalysisEnum, - % {{{2 - %INITIAL VELOCITIES - fields={'initialization.vx','initialization.vy'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %FORCINGS - fields={'surfaceforcings.mass_balance'}; - checkforcing(md,fields); - - %DIRICHLET - fields={'prognostic.spcthickness'}; - checkforcing(md,fields); - - %DIRICHLET IF THICKNESS <= 0 - %Check the size of prognostic.spcthickness - if any(md.geometry.thickness<=0), - pos=find(md.geometry.thickness<=0); - if any(isnan(md.prognostic.spcthickness(pos))), - message(['model not consistent: model ' md.miscellaneous.name ' has some nodes with 0 thickness']); - end - end - - %}}} - case HydrologyAnalysisEnum, - % {{{2 - fields={'hydrology.spcwatercolumn'}; - checkforcing(md,fields); - fields={'initialization.watercolumn'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - %}}} - case ThermalAnalysisEnum, - % {{{2 - %EXTRUSION - if (md.mesh.dimension==2), - if md.private.solution==TransientSolutionEnum, - return; - else - message(['model not consistent: for a ' EnumToString(md.private.solution) ' computation, the model must be 3d, extrude it first!']) - end - end - - %CHECK THAT SPCTEMPERATURE IS AN APPROPRIATE FORCING - fields={'thermal.spctemperature'}; - checkforcing(md,fields); - - %CHECK THAT WE ARE NOT FULLY CONSTRAINED - if ~any(~isnan(md.thermal.spctemperature)) - message(['model not consistent: model ' md.miscellaneous.name ' is totally constrained for temperature, no need to solve!']); - end - - %VELOCITIES AND PRESSURE - fields={'initialization.vx','initialization.vy','initialization.vz','initialization.pressure','basalforcings.geothermalflux'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %THERMAL TRANSIENT - if md.timestepping.time_step~=0, - - %DT and NDT - fields={'timestepping.time_step','timestepping.final_time'}; - checkgreaterstrict(md,fields,0); - - %INITIAL TEMPERATURE, MELTING AND ACCUMULATION - fields={'initialization.temperature','basalforcings.melting_rate'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %INITIAL TEMPERATURE - fields={'initialization.temperature'}; - checkgreater(md,fields,0) - - %CHECK SPCTEMPERATURE that are not NaN are >0. - if find(any(md.thermal.spctemperature(find(~isnan(md.thermal.spctemperature(1:md.mesh.numberofvertices,:))))<=0)), - message(['model not consistent: model ' md.miscellaneous.name ' is constrained with negative or nil temperatures!']); - end - - end - %}}} - case EnthalpyAnalysisEnum, - % {{{2 - %EXTRUSION - if (md.mesh.dimension==2), - if md.private.solution==TransientSolutionEnum, - return; - else - message(['model not consistent: for a ' EnumToString(md.private.solution) ' computation, the model must be 3d, extrude it first!']) - end - end - - %VELOCITIES AND PRESSURE - fields={'initialization.vx','initialization.vy','initialization.vz','initialization.pressure','basalforcings.geothermalflux'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %THERMAL TRANSIENT - if md.timestepping.time_step~=0, - - %TIMESTEPPING - fields={'timestepping.time_step','timestepping.final_time'}; - checkgreaterstrict(md,fields,0); - - %INITIAL TEMPERATURE, MELTING, ACCUMULATION AND WATERFRACTION - fields={'initialization.temperature','basalforcings.melting_rate','initialization.waterfraction'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %INITIAL TEMPERATURE - fields={'initialization.temperature'}; - checkgreater(md,fields,0) - end - %}}} - case MeltingAnalysisEnum, - % {{{2 - % No checks for now - %}}} - case BalancethicknessAnalysisEnum, - % {{{2 - %VELOCITIES MELTING AND ACCUMULATION - fields={'initialization.vx','initialization.vy','basalforcings.melting_rate','balancethickness.thickening_rate'}; - checksize(md,fields,[md.mesh.numberofvertices 1]); - checknan(md,fields); - - %Triangle with zero velocity - if any(sum(abs(md.initialization.vx(md.elements)),2)==0 & sum(abs(md.initialization.vy(md.elements)),2)==0) - message('model not consistent: at least one triangle has all its vertices with a zero velocity'); - end - %}}} - case FlaimAnalysisEnum, - % {{{2 - % No checks for now - %}}} - otherwise - disp(['WARNING: no check implemented for analysis ' EnumToString(analysis_enum) '!']); - message('stop') - end - -end -%}}} - -%checks additional functions -%checklength {{{1 -function checklength(md,fields,fieldlength) - %CHECKSIZE - check length of a field - for i=1:length(fields), - eval(['field=md.' fields{i} ';']); - if length(field)~=fieldlength, - message(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]); - end - end -end -%}}} -%checksize {{{1 -function checksize(md,fields,fieldsize) - %CHECKSIZE - check size of a field - for i=1:length(fields), - eval(['field=md.' fields{i} ';']); - if isnan(fieldsize(1)), - if (size(field,2)~=fieldsize(2)), - message(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']); - end - elseif isnan(fieldsize(2)), - if (size(field,1)~=fieldsize(1)), - message(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']); - end - else - if ((size(field)~=fieldsize(1)) | (size(field,2)~=fieldsize(2))) - %WARNING: new version - if strcmp(fields{i},'cm_min'), - disp('!!! '); - disp('!!! WARNING: cm_min must now be of size [md.mesh.numberofvertices x 1]. Update your parameter file as follows:'); - disp('!!! '); - disp('!!! md.inversion.min_parameters=md.inversion.min_parameters*ones(md.mesh.numberofvertices,1);'); - disp('!!! '); - end - %WARNING: new version - if strcmp(fields{i},'cm_max'), - disp('!!! '); - disp('!!! WARNING: cm_max must now be of size [md.mesh.numberofvertices x 1]. Update your parameter file as follows:'); - disp('!!! '); - disp('!!! md.inversion.max_parameters=md.inversion.max_parameters*ones(md.mesh.numberofvertices,1);'); - disp('!!! '); - end - message(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]); - end - end - end -end -%}}} -%checknan {{{1 -function checknan(md,fields) - %CHECKNAN - check nan values of a field - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(isnan(field)), - message(['model not consistent: NaN values found in field ' fields{i}]); - end - end -end -%}}} -%checkreal{{{1 -function checkreal(md,fields) - %CHECKREAL - check real values of a field - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(~isreal(field)), - message(['model not consistent: complex values found in field ' fields{i}]); - end - end -end -%}}} -%checkgreaterstrict{{{1 -function checkgreaterstrict(md,fields,lowerbound) - %CHECKGREATERSTRICT - check values of a field - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(field<=lowerbound), - message(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]); - end - end -end -%}}} -%checkgreater{{{1 -function checkgreater(md,fields,lowerbound) - %CHECKGREATER - check values of a field - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(field=upperbound), - message(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]); - end - end -end -%}}} -%checkless{{{1 -function checkless(md,fields,upperbound) - %CHECKLESS - check values of a field - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(field>upperbound), - message(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]); - end - end -end -%}}} -%checkvalues {{{1 -function checkvalues(md,fields,values) - %CHECKVALUE - check that a field has specified values - for i=1:length(fields), - eval(['field=reshape(md.' fields{i} ',[prod(size(md.' fields{i} ')) 1]);']); - if any(~ismember(field,values)), - message(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]); - end - end -end -%}}} -%checkforcing {{{1 -function checkforcing(md,fields) - - for i=1:length(fields), - eval(['field=md.' fields{i} ';']); - if size(field,1)==md.mesh.numberofvertices, - if ~size(field,2)==1, - message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' should have only one column as there are md.mesh.numberofvertices lines']); - end - elseif size(field,1)==md.mesh.numberofvertices+1 - if any(field(end,:)~=sort(field(end,:))), - message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' columns should be chronological']); - end - if any(field(end,1:end-1)==field(end,2:end)), - message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' columns must not contain duplicate timesteps']); - end - else - message(['model not consistent: model ' md.miscellaneous.name ': ' fields{i} ' should have md.mesh.numberofvertices or md.mesh.numberofvertices+1 lines']); - end - end -end -%}}} - -%error messages -%modelconsistency{{{1 -function flag=modelconsistency(flag_in) - - persistent consistency; - - if nargin==1 & nargout==0, - %OK model is inconsistent, set flag as false - consistency=flag_in; - elseif nargin==0 & nargout==1, - if isempty(consistency), - %modelinconsistent has never been called, model is consistent - consistency=true; - end - else - message('Bad usage'); - end - - flag=consistency; -end -%}}} -%message{{{1 -function message(string) - - disp(string); - modelconsistency(false); -end -%}}} +%nothing for now