Changeset 908

Timestamp:

06/11/09 14:54:23 (16 years ago)

Author:

seroussi

Message:

simplified parameter file

File:

• issm/trunk/test/Validation/RoundIceSheet/Circ.par (modified) (1 diff)

issm/trunk/test/Validation/RoundIceSheet/Circ.par

@@ -1 @@
-
 %Ok, start defining model parameters here
 
-%material parameters
-        md.g=9.81;
-        md.rho_ice=910;
-        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;
+disp('      creating thickness');
+hmin=0.01;
+hmax=2756.7;
+radius=(sqrt((md.x).^2+(md.y).^2));
+radiusmax=max(radius);
+md.thickness=hmin*ones(size(md.x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8);
+md.firn_layer=10*ones(md.numberofgrids,1);
+md.bed=0*md.thickness;
+md.surface=md.bed+md.thickness;
 
-%Solution parameters
-        %parallelization 
-        client_server_mode='no';
-        md.cluster='none';
-        md.np=2;
-        md.time=1;
-        md.exclusive=0;
+disp('      creating drag');
+md.drag_type=2; %0 none 1 plastic 2 viscous
+md.drag=20*ones(md.numberofgrids,1); %q=1. %no drag is specified in the analytical solution
+%Take care of iceshelves: no basal drag
+pos=find(md.elementoniceshelf);
+md.drag(md.elements(pos,:))=0;
+md.p=ones(md.numberofelements,1);
+md.q=ones(md.numberofelements,1);
 
-        %statics
-        md.eps_rel=0.01;
-        md.eps_abs=10;
-        md.penalty_offset=3.5;
-        md.lowmem=1;
-        if md.numberofgrids<1000000,
-        md.sparsity=.001;
-        else
-        md.sparsity=.0001;
-        end
+disp('      creating temperatures');
+md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
 
-        %dynamics
-        md.dt=1*md.yts; %1 year
-        md.ndt=md.dt*10; 
-        md.artificial_diffusivity=1;
+disp('      creating flow law paramter');
+md.B=6.81*10^(7)*ones(md.numberofgrids,1); %to have the same B as the analytical solution 
+md.n=3*ones(md.numberofelements,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 accumulation rates');
+md.accumulation=0.3*ones(md.numberofgrids,1)/md.yts; %1m/a
 
-        disp('      creating thickness');
-        hmin=0.01;
-        hmax=2756.7;
-        radius=(sqrt((md.x).^2+(md.y).^2));
-        radiusmax=max(radius);
-        md.thickness=hmin*ones(size(md.x,1),1)+hmax*(4*((1/2)^(4/3)*ones(size(md.x,1),1)-((radius)./(2*radiusmax)).^(4/3))).^(3/8);
-        md.firn_layer=10*ones(md.numberofgrids,1);
-        md.bed=0*md.thickness;
-        md.surface=md.bed+md.thickness;
+disp('      creating velocities');
+constant=0.3;
+md.vx_obs=constant/2*md.x.*(md.thickness).^-1;
+md.vy_obs=constant/2*md.y.*(md.thickness).^-1;
+md.vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
 
-        disp('      creating drag');
-        md.drag_type=2; %0 none 1 plastic 2 viscous
-        md.drag=20*ones(md.numberofgrids,1); %q=1. %no drag is specified in the analytical solution
-        %Take care of iceshelves: no basal drag
-        pos=find(md.elementoniceshelf);
-        md.drag(md.elements(pos,:))=0;
-        md.p=ones(md.numberofelements,1);
-        md.q=ones(md.numberofelements,1);
+%Deal with boundary conditions:
+disp('      boundary conditions for diagnostic model: ');
+%Build gridonicefront, array of boundary grids belonging to the icefront:
+gridonicefront=double(md.gridonboundary);
 
-        disp('      creating temperatures');
-        md.observed_temperature=(273-20)*ones(md.numberofgrids,1);
+md.gridondirichlet_diag=zeros(md.numberofgrids,1);
+pos=find(md.gridonboundary & ~gridonicefront);md.gridondirichlet_diag(pos)=1;
+md.dirichletvalues_diag=zeros(md.numberofgrids,2);
 
-        disp('      creating flow law paramter');
-        md.B=6.81*10^(7)*ones(md.numberofgrids,1); %to have the same B as the analytical solution 
-        md.n=3*ones(md.numberofelements,1);
+radius=sqrt((md.x).*md.x+(md.y).*md.y);
+pos=find(radius==min(radius));md.gridondirichlet_diag(pos)=1;
+md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center
 
-        disp('      creating accumulation rates');
-        md.accumulation=0.3*ones(md.numberofgrids,1)/md.yts; %1m/a
+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('      creating velocities');
-        constant=0.3;
-        md.vx_obs=constant/2*md.x.*(md.thickness).^-1;
-        md.vy_obs=constant/2*md.y.*(md.thickness).^-1;
-
-        %md.vx_obs=zeros(md.numberofgrids,1);
-        %md.vy_obs=zeros(md.numberofgrids,1);
-        md.vel_obs=(sqrt((md.vx_obs).^2+(md.vy_obs).^2));
-        
-        %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',1);
-        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);
-        
-        radius=sqrt((md.x).*md.x+(md.y).*md.y);
-        pos=find(radius==min(radius));md.gridondirichlet_diag(pos)=1;
-        md.x(pos)=0; md.y(pos)=0; %the closest node to the center is changed to be exactly at the center
-
-        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 thermal model: ');
-        md.melting=zeros(md.numberofgrids,1);
-        md.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperatures
-        md.dirichletvalues_thermal=md.observed_temperature;
-        md.geothermalflux=zeros(md.numberofgrids,1); 
-        pos=find(md.elementonicesheet);
-        %md.geothermalflux(md.elements(pos,:))=50*10^-3; %50 mW/m^2
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-        ServerDisconnect;
-end   
+disp('      boundary conditions for thermal model: ');
+md.gridondirichlet_thermal=zeros(md.numberofgrids,1);