%The aim of this program is to compare a model with an analytical solution given in MacAyeal EISMINT : Lessons in Ice-Sheet Modeling
%The constants used in the file circ.par are the ones used in the analytical solution.
%This solution is valid for a circular contour

%You can adjust the parameters numberofpoints and density to refine the mesh and check the precision of your solution. You can change the number of layers of the domain with the parameter layernumber.

%To begin with the numerical model
md=model;

numberofgrids=20;
density=190000;
layernumber=4;

%First we have to create the domain outline 
contour=CreateDomainOutline(numberofgrids,750000);

expgen('DomainOutline.exp',contour,1);
md=mesh(md,'DomainOutline.exp',density)

%move the closest node to the center
[mini pos]=min(md.x.^2+md.y.^2);
md.x(pos)=0;
md.y(pos)=0;

%Then we have to create the geography and the parameters of the mesh
md=geography(md,'',''); %We can not test iceshelves nor ice rises with this analytical solution

%We now create a file Front.exp needed to parameterize but the analytical solution can't support an icefront so this file is empty
front=CreateDomainOutline(numberofgrids,1000000);

expgen('Front.exp',front,1);

md=parameterize(md,'Circ.par');

%Calculation of the analytical 2d velocity field
constant=0.3;
vx_obs=constant/2*md.x.*(md.thickness).^-1;
vy_obs=constant/2*md.y.*(md.thickness).^-1;
vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));

%We extrude the model to have a 3d model
md=extrude(md,layernumber,3);
md=setelementstype(md,'Pattyn','all');

%Spc the grids on the bed
pos=find(md.gridonbed);
md.gridondirichlet_diag(pos)=1;

%Now we can solve the problem 
md=solve(md,'diagnostic','ice');

%Calculate the depth averaged velocity field (2d):
vel=zeros(md.numberofgrids2d,1);

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

%Plot of the velocity from the exact and calculated solutions
figure(1)
subplot(2,2,1)
p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel,'FaceColor','interp','EdgeColor','none');
title('Modelled velocity','FontSize',14,'FontWeight','bold')
colorbar; 
caxis([0 200]);
   
subplot(2,2,2)
p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
vel_obs,'FaceColor','interp','EdgeColor','none');
title('Analytical velocity','FontSize',14,'FontWeight','bold')
colorbar; 
caxis([0 200]);

subplot(2,2,3)
hold on;
plot(sqrt((md.x(1:md.numberofgrids2d)).^2+(md.y(1:md.numberofgrids2d)).^2),vel,'r.');
plot(sqrt((md.x2d).^2+(md.y2d).^2),vel_obs,'b.');
title('Analytical vs calculated velocity','FontSize',14,'FontWeight','bold');
xlabel('distance to the center of the icesheet [m]','FontSize',14,'FontWeight','bold');
ylabel('velocity [m/yr]','FontSize',14,'FontWeight','bold');
legend('calculated velocity','exact velocity');
axis([0 750000 0 200]);
hold off;

subplot(2,2,4)
p=patch('Faces',md.elements2d,'Vertices',[md.x2d md.y2d],'FaceVertexCData',...
abs(vel-vel_obs)./vel_obs*100,'FaceColor','interp','EdgeColor','none');
title('Relative misfit [%]','FontSize',14,'FontWeight','bold')
colorbar;
caxis([0 100]);

