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

% 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/test10_iceshelf_icefront_stokes: ';

% 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','all','stokes','all');

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

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



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

%Load old velocities and compare with the new model
load archive_temperature
post=find(archive_temperature-temp_stokes);
error_temp=abs(norm((archive_temperature(post)-temp_stokes(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 velocities 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