source: issm/trunk/test/Verification/IceShelfIceFrontM3d_2/runme.m@ 487

function runme(varargin)
% This file can be run to check that the current version of issm is giving 
% coherent results. This test deals with an icesheet with icefront for a 3d 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 analysis_types to be tested
if nargin==1,
        packages=varargin{1};
else
        packages={'macayeal','ice','cielo_serial'};
end
solutions={'diagnostic','thermalsteady','thermaltransient','prognostic','transient'};

%Initialize log message for nightly runs.
testname='IceShelfIceFrontM3d_2';
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 solution requested
                if (~(strcmpi(package,'ice') | strcmpi(solution,'diagnostic')) | strcmpi(solution,'transient')),
                        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',100000);
                md=geography(md,'all','');
                md=parameterize(md,'Square.par');
                md=extrude(md,8,4);
                md=setelementstype(md,'macayeal','all');
                if strcmpi(package,'cielo_parallel'), md.cluster='wilkes'; end
                if md.numberofgrids==832
                        load Velocities; md.vx=0.8*vx; md.vy=0.8*vy; md.vz=0.8*vz; md.temperature=temperature-1; md.pressure=pressure;
                end

                %compute solution
                if (strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel')),
                        md=solve(md,'analysis_type',solution,'package','cielo');
                else
                        md=solve(md,'analysis_type',solution,'package',package);
                end

                %compute field to be checked
                if strcmpi(solution,'diagnostic'),
                        fields={'vy','vz'};
                elseif strcmpi(solution,'thermalsteady'),
                        fields={'temperature','melting'};
                elseif strcmpi(solution,'thermaltransient'),
                        fields={'thermaltransient_results(end).temperature','thermaltransient_results(end).melting'};
                elseif strcmpi(solution,'prognostic'),
                        fields={'new_thickness'};
                elseif strcmpi(solution,'transient'),
                        fields={'transient_results(end).vy','transient_results(end).vz','transient_results(end).temperature','transient_results(end).melting','transient_results(end).thickness'};
                end

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

                for k=1:length(fields),
                        field=fields{k};

                        %compare to archive
                        eval(['Archive=Archive' package solution '_field' num2str(k) ';']);
                        eval(['error_diff=abs(norm((Archive(find(Archive))-md.' field  '(find(Archive)))./Archive(find(Archive)),2));']);

                        %disp test result
                        if (error_diff>tolerance);
                                disp(sprintf(['\n\nERROR   (difference=%-7.2g > %g) --> test: %-25s solution: %-16s package: %-14s field: ' field  '.\n\n'],error_diff,tolerance,testname,solution,package));
                        else
                                disp(sprintf(['\n\nSUCCESS (difference=%-7.2g < %g) --> test: %-25s solution: %-16s package: %-14s field: ' field  '.\n\n'],error_diff,tolerance,testname,solution,package));
                        end

                end
        end
end