% This file can be run to check that the current version of macayeal and ice models are giving 
% coherent results. This test deals with an iceshelf with icefront for a 2d model. The geometry 
% is square. Just run this file in Matlab, with a properly setup Ice 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^-14. 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.

%Initialize log message for nightly runs.
logstring='Verification/test3_iceshelf_icefront: ';

% Create model;

tolerance=10^-13;


md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'macayeal','all');

%Compute solution with MacAyeal's model
md=solve(md,'diagnostic','macayeal');
vel_macayeal=md.vel;

%Compute solution with Ice model
md.acceleration=0;
md=solve(md,'diagnostic','ice');
vel_ice=md.vel;

%Compute solution with accelerated Ice model
md.acceleration=1;
md=solve(md,'diagnostic','ice');
vel_acc=md.vel;

%Compute prognostic with Ice model
md=solve(md,'prognostic','ice');
thickness=md.new_thickness;



%Load old velocities and compare with the new model

%Macayeal
load archive_macayeal
pos=find(archive_macayeal-vel_macayeal);
error_vel=abs(norm((archive_macayeal(pos)-vel_macayeal(pos))./archive_macayeal(pos),2));

disp(sprintf('numerical difference between old and new version of MacAyeal is : %d', error_vel));

if (error_vel>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR.  Results from MacAyeal model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from MacAyeal model are comform with the archive version'));
end

%Ice
load archive_ice
pos=find(archive_ice-vel_ice);
error_vel=abs(norm((archive_ice(pos)-vel_ice(pos))./archive_ice(pos),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_vel));

if (error_vel>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from Ice model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from Ice model are comform with the archive version'));
end

%Accelerated Ice
load archive_acc
pos=find(archive_acc-vel_acc);
error_vel=abs(norm((archive_acc(pos)-vel_acc(pos))./archive_ice(pos),2));

disp(sprintf('numerical difference between old and new version of Accelerated Ice is : %d', error_vel));

if (error_vel>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from Accelerated Ice model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from Accelerated Ice model are comform with the archive version'));
end

%Ice prognostic
load archive_thickness
pos=find(archive_thickness-thickness);
error_thickness=abs(norm((archive_thickness(pos)-thickness(pos))./archive_thickness(pos),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_thickness));

if (error_thickness>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from prognostic model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from prognostic model are comform with the archive version'));
end


