
% This file can be run to check that the conduction is correctly modeled.
% There is no velocity (no advection) the only thermal boundary conditions are an imposed temperature
% at the lower and upper surface. The result must be a linear temperature from the upper to the lower
% surface. if it is not the case, something is thermal modeling has been changed...
% 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,'Shelf.exp','');
md=parameterize(md,'Square.par');
md=extrude(md,11,1);
md=setelementstype(md,'Pattyn','all');
thermalboundarycondition;

%analytical results
%d2T/dz2=0 T(bed)=10 T(surface)=0  => T=10/thickness*z-10*surface/thickness
%each layer of the 3d mesh must have a constant value
Tlayer=10:-1:0;
for layer=1:md.numlayers
	md.observed_temperature((layer-1)*md.numberofgrids2d+1:layer*md.numberofgrids2d)=Tlayer(layer);
end

%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)
