[5046] | 1 | % This file can be run to check that the advection-diffusion is correctly modeled.
|
---|
| 2 | % There is u=v=0 and w=cst everywhere the only thermal boundary conditions are an imposed temperature
|
---|
| 3 | % at upper surface and an impose flux at its base.
|
---|
[5606] | 4 | printingflag=false;
|
---|
[5046] | 5 |
|
---|
| 6 | md=model;
|
---|
[9703] | 7 | md=setmesh(md,'../Exp/Square.exp',100000);
|
---|
[9641] | 8 | md=setmask(md,'','');
|
---|
[5046] | 9 | md=parameterize(md,'../Par/SquareThermal.par');
|
---|
| 10 | md=extrude(md,30,1); %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
|
---|
[9664] | 11 | md=setflowequation(md,'Pattyn','all');
|
---|
[5046] | 12 |
|
---|
| 13 | %Thermal boundary conditions
|
---|
[9632] | 14 | pos1=find(md.elementonbed); md.thermal.spctemperature(md.elements(pos1,1:3))=10;
|
---|
| 15 | pos2=find(md.elementonsurface); md.thermal.spctemperature(md.elements(pos2,4:6))=0;
|
---|
[9725] | 16 | md.initialization.vz=0.1*ones(md.mesh.numberofvertices,1);
|
---|
[9684] | 17 | md.initialization.vel=sqrt( md.initialization.vx.^2+ md.initialization.vy.^2+ md.initialization.vz.^2);
|
---|
[9725] | 18 | md.initialization.pressure=zeros(md.mesh.numberofvertices,1);
|
---|
[5046] | 19 |
|
---|
| 20 | %analytical results
|
---|
| 21 | %d2T/dz2-w*rho_ice*c/k*dT/dz=0 T(surface)=0 T(bed)=10 => T=A exp(alpha z)+B
|
---|
[9636] | 22 | alpha=0.1/md.constants.yts*md.materials.rho_ice*md.heatcapacity/md.thermalconductivity; %alpha=w rho_ice c /k and w=0.1m/an
|
---|
[5046] | 23 | A=10/(exp(alpha*(-1000))-1); %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10
|
---|
| 24 | B=-A;
|
---|
[9684] | 25 | md.initialization.temperature=A*exp(alpha*md.z)+B;
|
---|
[5046] | 26 |
|
---|
| 27 | %modeled results
|
---|
[8630] | 28 | md.cluster=generic('name',oshostname(),'np',2);
|
---|
[8295] | 29 | md=solve(md,ThermalSolutionEnum);
|
---|
[5046] | 30 |
|
---|
| 31 | %plot results
|
---|
[9725] | 32 | comp_temp=zeros(md.mesh.numberofvertices,1);
|
---|
[5046] | 33 | comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
|
---|
[9684] | 34 | relative=abs((comp_temp-md.initialization.temperature)./md.initialization.temperature)*100;
|
---|
| 35 | relative(find(comp_temp==md.initialization.temperature))=0;
|
---|
[9433] | 36 | plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data','temperature','view',3,...
|
---|
[9684] | 37 | 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.initialization.temperature,...
|
---|
[5271] | 38 | 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
|
---|
[5422] | 39 | 'figposition','mathieu','FontSize#all',20)
|
---|
| 40 | if printingflag,
|
---|
| 41 | set(gcf,'Color','w')
|
---|
| 42 | printmodel('thermaladvection','png','margin','on','marginsize',25,'frame','off','resolution',0.7,'hardcopy','off');
|
---|
[6088] | 43 | system(['mv thermaladvection.png ' ISSM_DIR '/website/doc_pdf/validation/Images/EISMINT ']);
|
---|
[5422] | 44 | end
|
---|
[5098] | 45 |
|
---|
| 46 | %Fields and tolerances to track changes
|
---|
| 47 | field_names ={'AdvectionTemperature'};
|
---|
| 48 | field_tolerances={1e-13};
|
---|
| 49 | field_values ={comp_temp};
|
---|