source: issm/trunk/test/Verification/test1_icesheet_noicefront/runme.m@ 16

% 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 icesheet without 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/test1_icesheet_noicefront: ';

% Create model;

tolerance=10^-13;


md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'','');
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 solution with Hutter elements in Ice model
md=setelementstype(md,'hutter','all');
md=solve(md,'diagnostic','ice');
vel_hutter=md.vel;


%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_acc(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 with Hutter elements
load archive_hutter
pos=find(archive_hutter-vel_hutter);
error_vel=abs(norm((archive_hutter(pos)-vel_hutter(pos))./archive_hutter(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 with Hutter elements differ from the archive version'));
else
        logoutput(logstring,sprintf('%s\n','    SUCCESS. Results from Ice model with Hutter elements are comform with the archive version'));
end