% This file can be run to check that the current version of Pattyn model is giving 
% coherent results. This test deals with an iceshelf with icefront for a 3d 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/test11_penalties_iceshelf_allmodels: ';

% Create model;

tolerance=10^-13;


md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'all','');
md=parameterize(md,'Square.par');
md=extrude(md,10,3);
md=setelementstype(md,'pattyn','Pattyn.exp','fill','macayeal','stokes','Stokes.exp');

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

%Load old velocities and compare with the new model
load archive_velocity 
pos=find(archive_velocity-velocity);
error_vel=abs(norm((archive_velocity(pos)-velocity(pos))./archive_velocity(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 diagnostic model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from diagnostic model are comform with the archive version'));
end



%Compute solution with Ice model
md=solve(md,'thermalsteady','ice');
temperature=md.temperature;
melting_stokes=md.melting;

%Load old temperatures and compare with the new model
load archive_temperature
post=find(archive_temperature-temperature);
error_temp=abs(norm((archive_temperature(post)-temperature(post))./archive_temperature(post),2));

disp(sprintf('numerical difference between old and new version of Ice is : %d', error_temp));
if (error_temp>tolerance);
	logoutput(logstring,sprintf('%s\n','    ERROR. Results from thermal model differ from the archive version'));
else
	logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from thermal model are comform with the archive version'));
end


%Load old melting and compare with the new model
load archive_melting
post=find(archive_melting-melting_stokes);
error_melting=abs(norm((archive_melting(post)-melting_stokes(post))./archive_melting(post),2));

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

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