
% This file can be run to check that the advection-diffusion  is correctly modeled.
% There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature
% at upper surface and an impose flux at its base.
% Just run this file in Matlab, with a properly setup Ice code. 

% Create model;

md=model;
md=mesh(md,'DomainOutline.exp',100000);
md=geography(md,'','');
md=parameterize(md,'Square.par');
md=extrude(md,20,1);   %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
md=setelementstype(md,'Pattyn','all'); 
thermalboundarycondition;

%analytical results
%d2T/dz2-w*rho_ice*c/k*dT/dz=0   T(surface)=0  T(bed)=10   => T=A exp(alpha z)+B
alpha=0.1/md.yts*md.rho_ice*md.heatcapacity/md.thermalconductivity;   %alpha=w rho_ice c /k  and w=0.1m/an
A=10/(exp(alpha*(-1000))-1);    %A=T(bed)/(exp(alpha*bed)-1)  with bed=-1000 T(bed)=10
B=-A;
md.observed_temperature=A*exp(alpha*md.z)+B;

%modeled  results
md=solve(md,'thermalsteady');

%plot results
relative=abs((md.temperature-md.observed_temperature)./md.observed_temperature)*100;
relative(find(md.temperature==md.observed_temperature))=0;
plotmodel(md,'data','temperature','title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',md.temperature-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3)
