Changeset 515

Timestamp:

05/20/09 08:34:42 (16 years ago)

Author:

Mathieu Morlighem

Message:

fixed test geothermal flux

Location:

issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux

Files:

• Square.par (modified) (4 diffs)
• runme.m (modified) (1 diff)

issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/Square.par

@@ -1 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-        md.g=9.8;
-        md.rho_ice=917;
-        md.rho_water=1023;
-        di=md.rho_ice/md.rho_water;
-        md.yts=365*24*3600;
-        md.heatcapacity=2009;
-        md.thermalconductivity=2.2; %W/mK
-        md.beta=9.8*10^-8;
-
-%Solution parameters
         %parallelization 
         md.cluster='none';
-        md.np=2;
-        md.time=1;
-        md.exclusive=0;
-
-        %statics
-        md.eps_rel=0.01;
-        md.eps_abs=10;
-        %md.penalty_horiz=10^25; %penalty parameters for rifts
-        %md.penalty_vert=10^15; %penalty parameters for rifts
-        md.lowmem=1;
-        if md.numberofgrids<1000000,
-        md.sparsity=.001;
-        else
-        md.sparsity=.0001;
-        end
-
-        %dynamics
-        md.dt=1*md.yts; %1 year
-        md.ndt=md.dt*10; 
-        md.artificial_diffusivity=1;
-
-        %control
-        md.control_type={'drag'}; %'drag', 'B'
-        md.nsteps=5;
-        md.tolx=10^-4;
-        md.maxiter=20;
-        md.optscal=10;
-        md.fit='logarithmic'; %'absolute','relative','logarithmic'
-        md.meanvel=1000/md.yts; %1000 meters/year
-        md.epsvel=eps;
-
 
         disp('      creating thickness');
-        hmax=1000;
-        md.thickness=hmax*ones(md.numberofgrids,1)
+        h=1000;
+        md.thickness=h*ones(md.numberofgrids,1);
         md.firn_layer=10*ones(md.numberofgrids,1);
-        md.surface=zeros(md.numberofgrids,1);
-        md.bed=-md.thickness;
-
-%       hmin=700;
-%       hmax=1000;
-%       ymin=min(md.y);
-%       ymax=max(md.y);
-%       md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
-%       md.firn_layer=10*ones(md.numberofgrids,1);
-%       md.surface=zeros(md.numberofgrids,1);
-%       md.bed=-md.thickness;
+        md.bed=zeros(md.numberofgrids,1);
+        md.surface=md.thickness;
 
         disp('      creating velocities');
@@ -83 +31 @@
         md.B=paterson(md.observed_temperature);
         md.n=3*ones(md.numberofelements,1);
-
+        
         disp('      creating accumulation rates');
         md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
@@ -90 +38 @@
         %Deal with boundary conditions:
         
-        disp('      boundary conditions for diagnostic model: ');
-        %Build gridonicefront, array of boundary grids belonging to the icefront:
-        gridinsideicefront=ArgusContourToMesh(md.elements,md.x,md.y,expread('Front.exp',1),'node',2);
-        gridonicefront=double(md.gridonboundary & gridinsideicefront);
-
-        md.gridondirichlet_diag=zeros(md.numberofgrids,1);
-        pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
-        md.dirichletvalues_diag=zeros(md.numberofgrids,2);
-
-        pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-        md.segmentonneumann_diag=md.segments(pos,:);
-        md.neumannvalues_diag=NaN*ones(length(md.segmentonneumann_diag),1); %dynamic boundary conditions (water pressure)
-
-        disp('      boundary conditions for prognostic model');
-        md.gridondirichlet_prog=zeros(md.numberofgrids,1);
-        md.dirichletvalues_prog=zeros(md.numberofgrids,1);
-        pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-        md.segmentonneumann_prog=md.segments(pos,:);
-        md.neumannvalues_prog=zeros(size(md.segmentonneumann_prog,1),1);
-        md.neumannvalues_prog(:)=NaN; %free radiation
-        
-        pos=find(gridonicefront(md.segments(:,1)) | gridonicefront(md.segments(:,2)));
-        md.segmentonneumann_prog2=md.segments(pos,:);
-        md.neumannvalues_prog2=zeros(size(md.segmentonneumann_prog2,1),1);
-        md.neumannvalues_prog2(:)=NaN; %free radiation
+        disp('      boundary conditions for diagnostic model');
+        md=SetIceShelfBC(md,'Front.exp');
         
         disp('      boundary conditions for thermal model');
@@ -120 +45 @@
         md.dirichletvalues_thermal=md.observed_temperature;
         md.geothermalflux=zeros(md.numberofgrids,1); 
-        pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %1 mW/m^2
-
-
-        
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-        ServerDisconnect;
-end   
+        pos=find(md.elementonicesheet);md.geothermalflux(md.elements(pos,:))=0.1; %100 mW/m^2

issm/trunk/test/Validation/ThermalTests/Simplegeothermalflux/runme.m

@@ -19 +19 @@
 %the result is linear with depth and is equal to 0 on the upper surface (See BC)
 %d2T/dz2=0  -k*dT/dz(bed)=G  T(surface)=0  => T=-G/k*(z-surface)
-md.observed_temperature=-0.1/md.thermalconductivity*md.z; %G=0.1 W/m2
+md.observed_temperature=-0.1/md.thermalconductivity*(md.z-md.surface); %G=0.1 W/m2
 
 %modeled  results
-md=solve(md,'thermalsteady');
+md=solve(md,'analysis_type','thermal','sub_analysis_type','steady','package','cielo');
 
 %plot results