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^-12;

%check usage
if nargin~=1,
	help ismodelselfconsistent
	error('ismodelselfconsistent error message: wrong usage');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   TRANSIENT   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%check analysis
if md.analysis_type==TransientAnalysisEnum,
	if md.dt<=0,
		error('model not consistent: field dt must be positive for a transient run')
	end

	%recursive call to ismodelselfconsistent
	analysis=[DiagnosticAnalysisEnum PrognosticAnalysisEnum ThermalAnalysisEnum];
	for i=1:length(analysis),
		md.analysis_type=analysis;
		ismodelselfconsistent(md);
	end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   COMMON CHECKS   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%COUNTER
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
if isempty(md.name),
	error(['model is not correctly configured: missing name!']);
end

%ELEMENTSTYPE
if size(md.elements_type,1)~=md.numberofelements | size(md.elements_type,2)~=2,
	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,1)==HutterFormulationEnum) + (md.elements_type(:,1)==MacAyealFormulationEnum)  + (md.elements_type(:,1)==PattynFormulationEnum)))
	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
end
if any(ones(md.numberofelements,1)-((md.elements_type(:,2)==StokesFormulationEnum) + (md.elements_type(:,2)==NoneFormulationEnum)))
	error(['model not consistent: types of elements have not been set properly, run setelementstype first'])
end
if strcmpi(md.type,'2d'),
	if (ismember(PattynFormulationEnum,md.elements_type(:,1)) |  ismember(StokesFormulationEnum,md.elements_type(:,2))),
		error(['model not consistent: for a 2d model, only MacAyeal''s and Hutter''s elements are allowed']);
	end
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
if (md.analysis_type==DiagnosticAnalysisEnum & any(ismember(MacAyealFormulationEnum,md.elements_type(:,1)) & ismember(PattynFormulationEnum,md.elements_type(:,1))))
	error(['model not consistent: coupling MacAyeal/Pattyn not implemented yet']);
end
if (md.isstokes & md.analysis_type==TransientAnalysisEnum());
	error(['model not consistent: Stokes transient not implemented yet']);
end

%ICEFRONT
if strcmpi(md.type,'2d'),
	fields={'pressureload'};
	checksize(md,fields,[NaN 3]);
elseif strcmpi(md.type,'3d'),
	fields={'pressureload'};
	checksize(md,fields,[NaN 5]);
end

%NO NAN
fields={'numberofelements','numberofgrids','x','y','z','drag','drag_type','p','q',...
	'rho_ice','rho_water','B','elementoniceshelf','surface','thickness','bed','g','lowmem','sparsity','nsteps','maxiter',...
	'tolx','np','eps_res','exclusive','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checknan(md,fields);

%FIELDS >= 0 
fields={'numberofelements','numberofgrids','elements','drag','drag_type','p','q',...
	'rho_ice','rho_water','B','elementoniceshelf','thickness','g','eps_res','eps_rel','eps_abs','nsteps','maxiter','tolx','exclusive',...
	'sparsity','lowmem','n','gridonbed','gridonsurface','elementonbed','elementonsurface','deltaH','DeltaH','timeacc','timedec'};
checkgreater(md,fields,0);

%FIELDS > 0
fields={'numberofelements','numberofgrids','elements','drag_type','p',...
	'rho_ice','rho_water','B','thickness','g','eps_res','eps_rel','eps_abs','maxiter','tolx',...
	'sparsity','deltaH','DeltaH','timeacc','timedec'};
checkgreaterstrict(md,fields,0);

%SIZE NUMBEROFELEMENTS
fields={'elements','p','q','elementoniceshelf','n','elementonbed'};
checklength(md,fields,md.numberofelements);

%SIZE NUMBEROFGRIDS
fields={'x','y','z','B','drag','spcvelocity','melting','accumulation','surface','thickness','bed','gridonbed','gridonsurface'};
checklength(md,fields,md.numberofgrids);

%THICKNESS = SURFACE - BED
if any((md.thickness-md.surface+md.bed)>tolerance),
	error(['model not consistent: model ' md.name ' violates the equality thickness=surface-bed!']);
end

%RIFTS
if md.numrifts,
	if ~strcmpi(md.type,'2d'),
		error(['model not consistent: models with rifts are only supported in 2d for now!']);
	end
end
if ~isstruct(md.rifts),
	if ~isnan(md.rifts),
		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
	end
end

%ARTIFICIAL DIFFUSIVITY
if ~ismember(md.artificial_diffusivity,[0 1]),
	error('model not consistent: artificial_diffusivity should be a scalar (1 or 0)');
end

%PARAMETEROUTPUT
if md.numoutput~=length(md.parameteroutput),
	error('model not consistent: numoutput should be the same size as parameteroutput');
end

%CONNECTIVITY
if strcmpi(md.type,'2d'),
	if md.connectivity<9, 
		error('model not consistent: connectivity should be at least 9 for 2d models');
	end
end
if strcmpi(md.type,'3d'),
	if md.connectivity<24, 
		error('model not consistent: connectivity should be at least 24 for 3d models');
	end
end

%LOWMEM = 0 or 1
if ((md.lowmem ~= 1) & (md.lowmem~=0)),
	error(['model not consistent: model ' md.name ' lowmem field should be 0 or 1']);
end

%PARALLEL
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  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%QMU
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 md.npart~=md.numberofgrids,
			error(['model not consistent: user supplied partition should have same size as md.npart']);
		end
	end
	if md.eps_rel>10^-3,
		error(['model not consistent: for qmu analysis, eps_rel should be least than 10^-5, 10^-15 being a better value']);
	end
end

%DIAGNOSTIC
if md.analysis_type==DiagnosticAnalysisEnum,

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

	%HUTTER ON ICESHELF WARNING
	if any(md.elements_type(:,1)==HutterFormulationEnum & 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

end

%PROGNOSTIC
if md.analysis_type==PrognosticAnalysisEnum,

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

end

%STEADYSTATE
if md.analysis_type==SteadystateAnalysisEnum,

	%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 strcmpi(md.type,'2d'),
		error(['model not consistent: for a steadystate computation, model needs to be 3d']);
	end
end

%THERMAL STEADY AND THERMAL TRANSIENT
if md.analysis_type==ThermalAnalysisEnum,

	%EXTRUSION
	if strcmp(md.type,'2d'),
		error(['model not consistent: for a ' md.analysis_type ' computation, the model must be 3d, extrude it first!'])
	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==ThermalAnalysisEnum & md.dt~=0

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

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

end

%CONTROL
if md.control_analysis,

	%CONTROL TYPE
	if ~ischar(md.control_type),
		error('model not consistent: control_type should be a string');
	end

	%LENGTH CONTROL FIELDS
	fields={'maxiter','optscal','fit','cm_jump'};
	checklength(md,fields,md.nsteps);

	%FIT
	if sum((double(md.fit==1) + double(md.fit==0) + double(md.fit==2))==1)~=md.nsteps
		error('model not consistent: wrong fits: fit should be a vector of size nsteps holding 0, 1 and 2 only')
	end

	%OBSERVED VELOCITIES
	fields={'vx_obs','vy_obs'};
	checksize(md,fields,[md.numberofgrids 1]);
	checknan(md,fields);

	%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
end

%QMU
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 %end function

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CHECK FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function checklength(md,fields,fieldlength)
	%CHECKSIZE - check length of a field
	for i=1:length(fields),
		if length(eval(['md.' fields{i}]))~=fieldlength,
			error(['model not consistent: field ' fields{i} ' length should be ' num2str(fieldlength)]);
		end
	end
end

function checksize(md,fields,fieldsize)
	%CHECKSIZE - check size of a field
	for i=1:length(fields),
		if isnan(fieldsize(1)),
			if (size(eval(['md.' fields{i}]),2)~=fieldsize(2)),
				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(2)) ' columns']);
			end
		elseif isnan(fieldsize(2)),
			if (size(eval(['md.' fields{i}]),1)~=fieldsize(1)),
				error(['model not consistent: field ' fields{i} ' should have ' num2str(fieldsize(1)) ' rows']);
			end
		else
			if ((size(eval(['md.' fields{i}]),1)~=fieldsize(1)) |  (size(eval(['md.' fields{i}]),2)~=fieldsize(2)))
				error(['model not consistent: field ' fields{i} ' size should be ' num2str(fieldsize(1)) ' x ' num2str(fieldsize(2))]);
			end
		end
	end
end

function checknan(md,fields)
	%CHECKNAN - check nan values of a field
	for i=1:length(fields),
		if any(isnan(eval(['md.' fields{i}]))),
			error(['model not consistent: NaN values in field ' fields{i}]);
		end
	end
end

function checkreal(md,fields)
	%CHECKREAL - check real values of a field
	for i=1:length(fields),
		if any(eval(['~isreal(md.' fields{i} ')'])),
			error(['model not consistent: complex values in field ' fields{i}]);
		end
	end
end

function checkgreaterstrict(md,fields,lowerbound)
	%CHECKGREATERSTRICT - check values of a field
	for i=1:length(fields),
		if any(eval(['md.' fields{i} '<=' num2str(lowerbound) ])),
			error(['model not consistent: field ' fields{i} ' should have values stricly above ' num2str(lowerbound)]);
		end
	end
end

function checkgreater(md,fields,lowerbound)
	%CHECKGREATER - check values of a field
	for i=1:length(fields),
		if any(eval(['md.' fields{i} '<' num2str(lowerbound) ])),
			error(['model not consistent: field ' fields{i} ' should have values above ' num2str(lowerbound)]);
		end
	end
end

function checklessstrict(md,fields,upperbound)
	%CHECKLESSSTRICT - check values of a field
	for i=1:length(fields),
		if any(eval(['md.' fields{i} '>=' num2str(upperbound) ])),
			error(['model not consistent: field ' fields{i} ' should have values stricly below ' num2str(upperbound)]);
		end
	end
end

function checkless(md,fields,upperbound)
	%CHECKLESS - check values of a field
	for i=1:length(fields),
		if any(eval(['md.' fields{i} '>' num2str(upperbound) ])),
			error(['model not consistent: field ' fields{i} ' should have values below ' num2str(upperbound)]);
		end
	end
end