% This file can be run to see the response to a bump on the bed of an ice sheet. The bump is not 
% transmitted to the surface of the ice. The results of the three models are saved in order to be 
% compared. The geometry is square and the ice sheet is flat. The velocity of the icesheet is 
% 100 m/yr in one particular direction on every border. Just run this file in Matlab, i
% with a properly setup Ice code. 
% Take care to have the same density and number of layers to compare the results easily.

% Create macayeal model;

mdm=model;
mdm=mesh(mdm,'DomainOutline.exp',150000);
mdm=geography(mdm,'','');
mdm=parameterize(mdm,'Square.par');
mdm=extrude(mdm,6,3);
mdm=setelementstype(mdm,'macayeal','all');

%Compute solution with Ice model
mdm=solve(mdm,'diagnostic','ice');

save modelmacayeal mdm

% Create pattyn model;

mdp=model;
mdp=mesh(mdp,'DomainOutline.exp',150000);
mdp=geography(mdp,'','');
mdp=parameterize(mdp,'Square.par');
mdp=extrude(mdp,6,3);
mdp=setelementstype(mdp,'pattyn','all');

%Compute solution with Ice model
mdp=solve(mdp,'diagnostic','ice');

save modelpattyn mdp

% Create stokes model;

mds=model;
mds=mesh(mds,'DomainOutline.exp',150000);
mds=geography(mds,'','');
mds=parameterize(mds,'Square.par');
mds=extrude(mds,6,3);
mds=setelementstype(mds,'pattyn','all','stokes','all');

%Compute solution with Ice model
mds=solve(mds,'diagnostic','ice');

save modelstokes mds
