% 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. md=model; md=mesh(md,'../Exp/Square.exp',100000); md=geography(md,'',''); md=parameterize(md,'../Par/SquareThermal.par'); md=extrude(md,30,1); %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%) md=setelementstype(md,'Pattyn','all'); %Thermal boundary conditions pos1=find(md.elementonbed); md.spctemperature(md.elements(pos1,1:3),1)=1; md.spctemperature(md.elements(pos1,1:3),2)=10; pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6),1)=1; md.spctemperature(md.elements(pos2,4:6),2)=0; md.vz=0.1*ones(md.numberofgrids,1); md.vel=sqrt( md.vx.^2+ md.vy.^2+ md.vz.^2); md.pressure=zeros(md.numberofgrids,1); %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,'analysis_type',ThermalSolutionEnum); %plot results comp_temp=zeros(md.numberofgrids,1); comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value; relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100; relative(find(comp_temp==md.observed_temperature))=0; plotmodel(md,'data',comp_temp,'title','modeled temperature','data','observed_temperature','view',3,'title','analytical temperature','view',3,'data',comp_temp-md.observed_temperature,'title','absolute error','view',3,'data',relative,'title','relative error [%]','view',3) %Fields and tolerances to track changes field_names ={'AdvectionTemperature'}; field_tolerances={1e-13}; field_values ={comp_temp};