source: issm/trunk/test/NightlyRun/test1302.m@ 9045

Last change on this file since 9045 was 9045, checked in by seroussi, 14 years ago

fixed spctemperature

  • Property svn:executable set to *
File size: 2.2 KB
RevLine 
[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]4printingflag=false;
[5046]5
6md=model;
7md=mesh(md,'../Exp/Square.exp',100000);
8md=geography(md,'','');
9md=parameterize(md,'../Par/SquareThermal.par');
10md=extrude(md,30,1); %NB: the more one extrudes, the better (10-> relative~0.35%, 20->0.1%, 30->0.05%)
11md=setelementstype(md,'Pattyn','all');
12
13%Thermal boundary conditions
[9045]14pos1=find(md.elementonbed); md.spctemperature(md.elements(pos1,1:3))=10;
15pos2=find(md.elementonsurface); md.spctemperature(md.elements(pos2,4:6))=0;
[8304]16md.vz=0.1*ones(md.numberofnodes,1);
[5046]17md.vel=sqrt( md.vx.^2+ md.vy.^2+ md.vz.^2);
[8304]18md.pressure=zeros(md.numberofnodes,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
22alpha=0.1/md.yts*md.rho_ice*md.heatcapacity/md.thermalconductivity; %alpha=w rho_ice c /k and w=0.1m/an
23A=10/(exp(alpha*(-1000))-1); %A=T(bed)/(exp(alpha*bed)-1) with bed=-1000 T(bed)=10
24B=-A;
25md.observed_temperature=A*exp(alpha*md.z)+B;
26
27%modeled results
[8630]28md.cluster=generic('name',oshostname(),'np',2);
[8295]29md=solve(md,ThermalSolutionEnum);
[5046]30
31%plot results
[8304]32comp_temp=zeros(md.numberofnodes,1);
[5046]33comp_temp(md.results.ThermalSolution.Temperature.index)=md.results.ThermalSolution.Temperature.value;
34relative=abs((comp_temp-md.observed_temperature)./md.observed_temperature)*100;
35relative(find(comp_temp==md.observed_temperature))=0;
[5271]36plotmodel(md,'data',comp_temp,'title','Modeled temperature [K]','data','observed_temperature','view',3,...
37 'title','Analytical temperature [K]','view',3,'data',comp_temp-md.observed_temperature,...
38 'title','Absolute error [K]','view',3,'data',relative,'title','Relative error [%]','view',3,...
[5422]39 'figposition','mathieu','FontSize#all',20)
40if 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]44end
[5098]45
46%Fields and tolerances to track changes
47field_names ={'AdvectionTemperature'};
48field_tolerances={1e-13};
49field_values ={comp_temp};
Note: See TracBrowser for help on using the repository browser.