% This file can be run to check that the current version of issm is giving 
% coherent results. This test deals with an icesheet without icefront for a 2d model. The geometry 
% is square. Just run this file in Matlab, with a properly setup ISSM code. 
% The results of this test will indicate if there is a difference between current computations 
% and archived results.

% Errors  between archived results and the current version will get flagged if they are not within 
% a certain tolerance. The current tolerance is 10^-12. If you have good reasons to believe this 
% tolerance should be lowered (for example, if you are running single precision compilers?), feel 
% free to tweak the tolerance variable.

%packages and solutions to be tested
packages={'macayeal','ice','cielo_serial','cielo_parallel'};
solutions={'diagnostic','prognostic'};

%Initialize log message for nightly runs.
testname='IceSheetNoIceFrontM2d_16';
tolerance=10^-12;

%go through the solutions requested
for i=1:length(packages),
	package=packages{i};

	for j=1:length(solutions),
		solution=solutions{j};

		%check package
		if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
			error(['package: ' package  ' in test: ' testname  ' not supported yet']);
		%check solution
		elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'prognostic'));
			error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
		end
		%check solution requested
		if (~strcmpi(package,'ice') & strcmpi(solution,'prognostic')),
			disp(sprintf(['\nsolution: ' solution  ', with package: ' package  ', in test: ' testname  ', not supported yet.\n']));
			continue
		end

		%initialize model
		md=model;
		md=mesh(md,'DomainOutline.exp',50000);
		md=geography(md,'','');
		md=parameterize(md,'Square.par');
		md=setelementstype(md,'macayeal','all');
		if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
		if md.numberofgrids==388
			load Velocities; md.vx=0.5*vx; md.vy=0.5*vy;
		end

		%compute solution
		if strcmpi(package,'cielo_parallel') & strcmpi(solution,'diagnostic'),
			md=solve(md,'diagnostic_horiz','cielo');
		elseif strcmpi(package,'cielo_serial'),
			eval(['md=solve(md,''' solution ''',''cielo'');']);
		else
			eval(['md=solve(md,''' solution ''',''' package ''');']);
		end

		%compute field to be checked
		if strcmpi(solution,'diagnostic'),
			fieldtest=md.vel;
		elseif strcmpi(solution,'prognostic'),
			fieldtest=md.new_thickness';
		end

		%load archive
		eval(['load Archive' package solution ]);

		%compare to archive
		eval(['pos=find(Archive' package solution ');']);
		eval(['error_diff=abs(norm((Archive' package solution '(pos)-fieldtest(pos))./Archive' package solution '(pos),2));']);

		%disp test result
		if (error_diff>tolerance);
			disp(['test: ' testname  ', solution: ' solution  ', package: ' package ' -> ERROR (difference=' num2str(error_diff) ' > ' num2str(tolerance) ')'])
		else
			disp(['test: ' testname  ', solution: ' solution  ', package: ' package ' -> SUCCESS (difference=' num2str(error_diff) ' < ' num2str(tolerance) ')'])
		end
	end
end