source: issm/trunk/test/Verification/IceSheetNoIceFrontP3d_18/updatearchive.m@ 451

% This file can be run to update the velocity archives  of the test1.
% 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 Ice code. 

% The archive files will be saved in this directory but will not commited to ice1. 
% Just commit the result if you want to.

%packages and solutions to be tested
packages={'macayeal','ice','cielo_serial','cielo_parallel'};
solutions={'diagnostic','thermalsteady','thermaltransient','prognostic','transient'};

%go through the solutions requested
testname='IceSheetNoIceFrontP3d_18';
for i=1:length(packages),
        package=packages{i};

        for j=1:length(solutions),
                solution=solutions{j};

                %check package
                if ~(strcmpi(package,'macayeal') | strcmpi(package,'ice') | strcmpi(package,'cielo_serial') | strcmpi(package,'cielo_parallel'));
                        error(['package: ' package  ' in test: ' testname  ' not supported yet']);
                        %check solution
                elseif ~(strcmpi(solution,'diagnostic') | strcmpi(solution,'thermalsteady') | strcmpi(solution,'thermaltransient') |...
                                strcmpi(solution,'prognostic') | strcmpi(solution,'transient'));
                        error(['solution: ' solution  ' in test: ' testname  ' not supported yet']);
                end
                %check solution requested
                if (~(strcmpi(package,'ice') | (strcmpi(package,'cielo_serial') & 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,'','');
                md=parameterize(md,'Square.par');
                md=extrude(md,8,4);
                md=setelementstype(md,'pattyn','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,solution,'cielo');
                else
                        md=solve(md,solution,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

                %save new archive
                for k=1:length(fields),
                        field=fields{k};
                        eval(['Archive' package solution '_field' num2str(k) '=md.' field  ';']);
                end
                eval(['save Archive' package solution ' Archive' package solution '_field*']);
        end
end