source: issm/trunk/test/Validation/MacAyealVsPattyn/test1_iceshelf/runme.m@ 4871

% This file can be run to compare MacAyeal and Pattyn's model ie a 2d and a 3d model.
% This test deals with a  mesh includind an iceshelf ans an icesheet. The geometry 
% is square. Just run this file in Matlab, with a properly setup Ice code. 

md=model;
md=mesh(md,'DomainOutline.exp',50000);
md=geography(md,'Shelf.exp','');
md=parameterize(md,'Square.par');
md=setelementstype(md,'Macayeal','all');

%Compute solution for a 2d model
md=solve(md,'analysis_type',DiagnosticSolutionEnum);
vel_2d=zeros(md.numberofgrids,1);
vel_2d(md.results.DiagnosticSolution.Vel.index)=md.results.DiagnosticSolution.Vel.value;

%Compute solution for a 3d model
md=extrude(md,5,3);
md=setelementstype(md,'Pattyn','Pattyn.exp','fill','Macayeal');
md=solve(md,'analysis_type',DiagnosticSolutionEnum);

%Calculate the average velocity on each grid
vel_3d=zeros(md.numberofgrids2d,1);
vel=vel_3d;
vel(md.results.DiagnosticSolution.Vel.index)=md.results.DiagnosticSolution.Vel.value;
grid_vel=0;

for i=1:md.numberofgrids2d
        for j=1:(md.numlayers-1)
                grid_vel=grid_vel+1/(2*(md.numlayers-1))*(vel(i+j*md.numberofgrids2d,1)+vel(i+(j-1)*md.numberofgrids2d,1));
        end
        vel_3d(i,1)=grid_vel;
        grid_vel=0;
end

vel_diff=(vel_2d-vel_3d)./vel_2d;
vel_diff(find(vel_2d==vel_3d))=0;

%Plot of the velocity from the macayeal and pattyn model 
figure(1)
subplot(2,2,1)
p1=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_2d,'FaceColor','interp','EdgeColor','none');
title('MacAyeal model [m/yr]','FontSize',14,'FontWeight','bold')
colorbar; 
   
subplot(2,2,2)
p2=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_3d,'FaceColor','interp','EdgeColor','none');
title('Pattyn model [m/yr]','FontSize',14,'FontWeight','bold')
colorbar; 

subplot(2,2,3)
p3=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_diff*100,'FaceColor','interp','EdgeColor','none');
title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
colorbar;

subplot(2,2,4)
p3=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_2d-vel_3d,'FaceColor','interp','EdgeColor','none');
title('Absolute misfit [m/yr]','FontSize',14,'FontWeight','bold')
colorbar;