function ismodelselfconsistent(md),
%ISMODELSELFCONSISTENT - check that model forms a closed form solvable problem.
%
%   Usage:
%      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
	error('ismodelselfconsistent error message: wrong usage');
end
%}}}

%recursive call for TRANSIENT {{{1
if (md.analysis_type==Transient2DSolutionEnum | md.analysis_type==Transient3DSolutionEnum),
	if md.dt<=0,
		error('model not consistent: field dt must be positive for a transient run')
	end

	%recursive call to ismodelselfconsistent
	if (md.dim==2),
		analysis=[DiagnosticSolutionEnum PrognosticSolutionEnum];
	else
		analysis=[DiagnosticSolutionEnum PrognosticSolutionEnum ThermalSolutionEnum];
	end

	for i=1:length(analysis),
		md.analysis_type=analysis(i);
		ismodelselfconsistent(md);
	end
end
%}}}

% Common checks
%COUNTER {{{1
if md.counter<3,
	error(['model ' md.name ' is not correctly configured. You forgot one step in the following sequence (mesh, geography, parameterize,setelementstype)!']);
end
%}}}
%NAME{{{1
if isempty(md.name),
	error(['model is not correctly configured: missing name!']);
end
%}}}
%ELEMENTS{{{1
fields={'elements'};
if (md.dim==2),
	checksize(md,fields,[md.numberofelements 3]);
else
	checksize(md,fields,[md.numberofelements 6]);
end
%}}}
%ELEMENTSTYPE{{{1
%Check the size of elements_type
fields={'elements_type'};
checksize(md,fields,[md.numberofelements 1]);
%Check the values of elements_type
checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() PattynStokesApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
if (md.dim==2),
	checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]);
end
if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
	error(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
end
%}}}
%VERTICESTYPE{{{1
%Check the size of verticess_type
fields={'vertices_type'};
checksize(md,fields,[md.numberofgrids 1]);
%Check the values of vertices_type
checkvalues(md,{'vertices_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
if (md.dim==2),
	checkvalues(md,{'vertices_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]);
end
if (md.ismacayealpattyn==0 && md.ishutter==0 && md.isstokes==0),
	error(['model not consistent: no elements type set for this model. at least one of ismacayealpattyn, ishutter and isstokes need to be =1']);
end
%}}}
%DG {{{1
if md.prognostic_DG==1;
	%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.dim==2),
	fields={'pressureload'};
	checksize(md,fields,[NaN 4]);
elseif md.dim==3,
	fields={'pressureload'};
	checksize(md,fields,[NaN 6]);
else
	error('dim should be either 2 3');
end
checkvalues(md,{'pressureload(:,end)'},[WaterEnum() AirEnum()]);
%}}}
%NO NAN {{{1
fields={'numberofelements','numberofgrids','x','y','z','drag_coefficient','drag_type','drag_p','drag_q',...
	'rho_ice','rho_water','rheology_B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
	'tolx','np','eps_res','max_nonlinear_iterations','exclusive','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checknan(md,fields);
%}}}}
%FIELDS >= 0 {{{1
fields={'numberofelements','numberofgrids','elements','drag_coefficient','drag_type','drag_p','drag_q',...
	'rho_ice','rho_water','rheology_B','elementoniceshelf','thickness','g','eps_res','max_nonlinear_iterations','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
	'sparsity','lowmem','rheology_n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checkgreater(md,fields,0);
%}}}
%FIELDS > 0 {{{1
fields={'numberofelements','numberofgrids','elements','drag_type','drag_p',...
	'rho_ice','rho_water','rheology_B','thickness','g','max_nonlinear_iterations','eps_res','eps_rel','eps_abs','maxiter','tolx',...
	'sparsity','deltaH','DeltaH','timeacc','timedec'};
checkgreaterstrict(md,fields,0);
%}}}
%SIZE NUMBEROFELEMENTS {{{1
fields={'drag_p','drag_q','elementoniceshelf','rheology_n','elementonbed'};
checksize(md,fields,[md.numberofelements 1]);
%}}}
%SIZE NUMBEROFGRIDS {{{1
fields={'x','y','z','rheology_B','drag_coefficient','melting_rate','accumulation_rate','surface','thickness','bed','gridonbed','gridonsurface'};
checksize(md,fields,[md.numberofgrids 1]);
%}}}
%OTHER SIZES {{{1
fields={'spcvelocity'};
checksize(md,fields,[md.numberofgrids 6]);
%}}}
%THICKNESS = SURFACE - BED {{{1
if any((md.thickness-md.surface+md.bed)>tolerance),
	error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
end
%}}}
%RIFTS{{{1
if md.numrifts,
	if ~(md.dim==2),
		error(['model not consistent: models with rifts are only supported in 2d for now!']);
	end
	if ~isstruct(md.rifts),
		error(['model not consistent: md.rifts should be a structure!']);
	end
	if ~isempty(find(md.segmentmarkers>=2)),
		%We have segments with rift markers, but no rift structure!
		error(['model not consistent: model ' md.name ' should be processed for rifts (run meshprocessrifts)!']);
	end
	%Check that rifts are filled with proper material
	checkvalues(md,{'rifts.fill'},[WaterEnum() AirEnum() IceEnum() MelangeEnum()]);
else
	if ~isnans(md.rifts),
		error(['model not consistent: md.rifts shoud be NaN since md.numrifts is 0!']);
	end
end
%}}}
%FLAGS (0 or 1){{{1
if ~ismember(md.artificial_diffusivity,[0 1]),
	error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
end
if ~ismember(md.prognostic_DG,[0 1]),
	error('model not consistent: prognostic_DG should be a scalar (1 or 0)');
end
if ~ismember(md.lowmem,[0 1]),
	error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
end
%}}}
%SCALAR
if ~isscalar(md.control_type),
	error('model not consistent: md.control_type should be a scalar');
end
%}}}
%PARAMETEROUTPUT {{{1
if md.numoutput~=length(md.parameteroutput),
	error('model not consistent: numoutput should be the same size as parameteroutput');
end
%}}}
%CONNECTIVITY {{{1
if (md.dim==2),
	if md.connectivity<9, 
		error('model not consistent: connectivity should be at least 9 for 2d models');
	end
end
if md.dim==3,
	if md.connectivity<24, 
		error('model not consistent: connectivity should be at least 24 for 3d models');
	end
end
%}}}
%PARALLEL{{{1
if ~strcmpi(md.cluster,'none'),

	%NAN VALUES
	fields={'time','np'};
	checknan(md,fields);

	%FIELD > 0
	fields={'time','np'};
	checkgreaterstrict(md,fields,0);

end
%}}}

% Solution checks
%DIAGNOSTIC{{{1
if md.analysis_type==DiagnosticSolutionEnum,

	%HUTTER ON ICESHELF WARNING
	if any(md.elements_type==HutterApproximationEnum & md.elementoniceshelf),
		disp(sprintf('\n !!! Warning: Hutter''s model is not consistent on ice shelves !!!\n'));
	end

	%SINGULAR
	if ~any(sum(md.spcvelocity(:,1:2),2)==2),
		error(['model not consistent: model ' md.name ' is not well posed (singular). You need at least one grid with fixed velocity!'])
	end

	%DIRICHLET IF THICKNESS <= 0
	if any(md.thickness<=0),
		pos=find(md.thickness<=0);
		if any(find(md.spcthickness(pos,1)==0)),
			error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
		end
	end

	%INITIAL VELOCITY
	if length(md.vx)==md.numberofgrids & length(md.vy)==md.numberofgrids,
		fields={'vx','vy'};
		checknan(md,fields);
	end
end
%}}}
%PROGNOSTIC{{{1
if md.analysis_type==PrognosticSolutionEnum,

	%INITIAL VELOCITIES
	fields={'vx','vy'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

	%CHECK THAT WE ARE NOT FULLY CONSTRAINED
	if (md.dim==2),
		if isempty(find(~md.spcthickness(:,1))),
			error(['model not consistent: model ' md.name ' is totally constrained for prognostic, no need to solve!']);
		end
	end

end
%}}}
%STEADYSTATE{{{1
if md.analysis_type==SteadystateSolutionEnum,

	%NDT
	if md.dt~=0,
		error(['model not consistent: for a steadystate computation, dt must be zero.']);
	end

	%PRESSURE
	if isnans(md.pressure),
		error(['model not consistent: for a steadystate computation, the model must have an initial pressure, even lithostatic will do.']);
	end

	%eps: 
	if isnan(md.eps_rel),
		error(['model not consistent: for a steadystate computation, eps_rel (relative convergence criterion) must be defined!']);
	end

	%dim: 
	if (md.dim==2),
		error(['model not consistent: for a steadystate computation, model needs to be 3d']);
	end
end
%}}}
%THERMAL {{{1
%THERMAL STEADY AND THERMAL TRANSIENT
if md.analysis_type==ThermalSolutionEnum,

	%EXTRUSION
	if (md.dim==2),
		error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
	end

	%CHECK THAT WE ARE NOT FULLY CONSTRAINED
	if isempty(find(~md.spctemperature(:,1))),
		error(['model not consistent: model ' md.name ' is totally constrained for temperature, no need to solve!']);
	end

	%VELOCITIES AND PRESSURE
	fields={'vx','vy','vz','pressure'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

end

%THERMAL TRANSIENT
if md.analysis_type==ThermalSolutionEnum & md.dt~=0,

	%DT and NDT
	fields={'dt','ndt'};
	checkgreaterstrict(md,fields,0);

	%INITIAL TEMPERATURE, MELTING AND ACCUMULATION
	fields={'temperature','accumulation_rate','melting_rate'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

	%INITIAL TEMPERATURE
	fields={'temperature','spctemperature(:,2)','observed_temperature'};
	checkgreater(md,fields,0)

end
%}}}
%BALANCEDTHICKNESS{{{1
if md.analysis_type==BalancedthicknessSolutionEnum

	%VELOCITIES MELTING AND ACCUMULATION
	fields={'vx','vy','accumulation_rate','melting_rate','dhdt'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

	%SPC				 
	if ~md.prognostic_DG,
		if any(md.spcthickness(find(md.gridonboundary))~=1),		 
			error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcthickness']);			 
		end 
	end

	%Triangle with zero velocity
	if any(sum(abs(md.vx(md.elements)),2)==0 & sum(abs(md.vy(md.elements)),2)==0)
		error('model not consistent: at least one triangle has all its vertices with a zero velocity');
	end
end
%}}}
%BALANCEDVELOCITIES{{{1
if md.analysis_type==BalancedvelocitiesSolutionEnum

	%VELOCITIES MELTING AND ACCUMULATION
	fields={'vx','vy','accumulation_rate','melting_rate'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

	%SPC
	if any(md.spcvelocity(find(md.gridonboundary),[1:2])~=1),
		error(['model not consistent: model ' md.name ' should have all the nodes on boundary constrained in field spcvelocity']);
	end
end
%}}}
%CONTROL{{{1
if md.control_analysis,

	%CONTROL TYPE
	if ~isnumeric(md.control_type),
		error('model not consistent: control_type should be an enum');
	end
	checkvalues(md,{'control_type'},[DhDtEnum DragCoefficientEnum RheologyBbarEnum]);

	%LENGTH CONTROL FIELDS
	fields={'maxiter','optscal','cm_responses','cm_jump'};
	checksize(md,fields,[md.nsteps 1]);

	%RESPONSES
	checkvalues(md,{'cm_responses'},[SurfaceAbsVelMisfitEnum SurfaceRelVelMisfitEnum SurfaceLogVelMisfitEnum SurfaceLogVxVyMisfitEnum SurfaceAverageVelMisfitEnum ThicknessAbsMisfitEnum]);

	%WEIGHTS
	fields={'weights'};
	checksize(md,fields,[md.numberofgrids 1]);
	checkgreater(md,fields,0);

	%OBSERVED VELOCITIES
	if md.analysis_type==BalancedthicknessSolutionEnum
		fields={'thickness_obs'};
		checksize(md,fields,[md.numberofgrids 1]);
		checknan(md,fields);
	else
		fields={'vx_obs','vy_obs'};
		checksize(md,fields,[md.numberofgrids 1]);
		checknan(md,fields);
	end

	%DIRICHLET IF THICKNESS <= 0
	if any(md.thickness<=0),
		pos=find(md.thickness<=0);
		if any(find(md.spcthickness(pos,1)==0)),
			error(['model not consistent: model ' md.name ' has some grids with 0 thickness']);
		end
	end

	%parameters
	fields={'cm_noisedmp'};
	checknan(md,fields);
end
%}}}
%QMU {{{1
if md.qmu_analysis,
	if md.qmu_params.evaluation_concurrency~=1,
		error(['model not consistent: concurrency should be set to 1 when running dakota in library mode']);
	end
	if ~isempty(md.part),
		if numel(md.part)~=md.numberofgrids,
			error(['model not consistent: user supplied partition for qmu analysis should have size md.numberofgrids x 1 ']);
		end
		if find(md.part)>=md.numberofgrids,
			error(['model not consistent: user supplied partition should be indexed from 0 (c-convention)']);
		end
		if min(md.part)~=0,
			error(['model not consistent: partition vector not indexed from 0 on']);
		end
		if max(md.part)>=md.numberofgrids,
			error(['model not consistent: partition vector cannot have maximum index larger than number of grids']);
		end
		if ~isempty(find(md.part<0)),
			error(['model not consistent: partition vector cannot have values less than 0']);
		end
		if ~isempty(find(md.part>=md.npart)),
			error(['model not consistent: partition vector cannot have values more than md.npart-1']);
		end
		if max(md.part)>=md.npart,
			error(['model not consistent: for qmu analysis, partitioning vector cannot go over npart, number of partition areas']);
		end
	end
	if ~md.qmu_relax,
		if md.eps_rel>1.1*10^-5,
			error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
		end
	end
end

if strcmpi(md.analysis_type,'qmu'),
	if ~strcmpi(md.cluster,'none'),
		if md.waitonlock==0,
			error(['model is not correctly configured: waitonlock should be activated when running qmu in parallel mode!']);
		end
	end
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,
			error(['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)),
				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
			end
		elseif isnan(fieldsize(2)),
			if (size(field,1)~=fieldsize(1)),
				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
			end
		else
			if ((size(field)~=fieldsize(1)) |  (size(field,2)~=fieldsize(2)))
				error(['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)),
			error(['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)),
			error(['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),
			error(['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<lowerbound),
			error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
		end
	end
end
%}}}
%checklessstrict{{{1
function checklessstrict(md,fields,upperbound)
	%CHECKLESSSTRICT - 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),
			error(['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),
			error(['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)),
			error(['model not consistent: field ' fields{i} ' should have values in ' num2str(values)]);
		end
	end
end
%}}}
