Changeset 3510

Timestamp:

04/12/10 08:46:50 (15 years ago)

Author:

Mathieu Morlighem

Message:

updated Bumps tests

Location:

issm/trunk/examples/Bumps/Bump2_bed

Files:

• Square.par (modified) (3 diffs)
• runme.m (modified) (6 diffs)

issm/trunk/examples/Bumps/Bump2_bed/Square.par

@@ -1 +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.cluster=oshostname;
+        md.np=8;
         md.time=1;
-        md.exclusive=0;
 
         %statics
-        md.lowmem=1;
         md.eps_rel=0.01;
         md.eps_abs=10; %m/yr
-        md.penalty_offset=4;
-        md.penalty_melting=10^7;
-        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');
-        hmin=500;
-        hmax=500;
-        ymin=min(md.y);
-        ymax=max(md.y);
-        md.thickness=hmax+(hmin-hmax)*(md.y-ymin)/(ymax-ymin);
+        md.thickness=500*ones(md.numberofgrids,1);
         md.firn_layer=10*ones(md.numberofgrids,1);
-        md.bed=-di*md.thickness-50; %*10^(-4)*min(0,md.y-900000);;
+        md.bed=-md.rho_ice/md.rho_water*md.thickness+5;
 
         md.surface=md.bed+md.thickness;
 
         %Add bumps to the ice sheet :
-        pos=find(sqrt((md.x-500000).^2+(md.y-500000).^2)<100000); %for a cicular bump
-        md.bed(pos)=md.bed(pos)+150*cos(pi/200000*(sqrt((md.y(pos)-500000).^2+(md.x(pos)-500000).^2))); %for a circular bump
+        radius= 500000;
+        center=[500000 500000];
+        height=150;
+        pos=find(sqrt((md.x-center(1)).^2+(md.y-center(2)).^2)<radius); %for a cicular bump
+        md.bed(pos)=md.bed(pos)+height*cos(pi/(2*radius)*(sqrt((md.y(pos)-center(1)).^2+(md.x(pos)-center(2)).^2))); %for a circular bump
 
         %Change the thickness to follow the new bed
+        %height=-15;
+        %md.surface(pos)=md.surface(pos)+height*cos(pi/(2*radius)*(sqrt((md.y(pos)-center(1)).^2+(md.x(pos)-center(2)).^2))); %for a circular bump
         md.thickness=md.surface-md.bed;
 
         disp('      creating velocities');
-        md.vx_obs=zeros(md.numberofgrids,1);
-        md.vy_obs=zeros(md.numberofgrids,1);
+        md.vx_obs=   0*ones(md.numberofgrids,1);
+        md.vy_obs=100*ones(md.numberofgrids,1);
         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
         
@@ -73 +37 @@
         md.drag_type=2; %0 none 1 plastic 2 viscous
         md.drag=10*ones(md.numberofgrids,1); %q=1.
-        %Take care of iceshelves: no basal drag
         pos=find(md.elementoniceshelf);
         md.drag(md.elements(pos,:))=0;
@@ -86 +49 @@
         md.n=3*ones(md.numberofelements,1);
 
-        disp('      creating accumulation rates');
-        md.accumulation=ones(md.numberofgrids,1)/md.yts; %1m/a
-        md.melting=0*ones(md.numberofgrids,1)/md.yts; %1m/a
-
         %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=100*ones(md.numberofgrids,2)*[0,0;0,1];
-
-        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.gridondirichlet_thermal=ones(md.numberofgrids,1); %surface temperature
-        md.dirichletvalues_thermal=md.observed_temperature;
-        md.geothermalflux=zeros(md.numberofgrids,1); 
-
-
-% Some Cielo code, ignore.
-if strcmp(md.cluster,'yes')
-        ServerDisconnect;
-end   
+        md=SetIceShelfBC(md,'Front.exp');
+        md.verbose=10;

issm/trunk/examples/Bumps/Bump2_bed/runme.m

@@ -8 +8 @@
 % Create macayeal model;
 
-mdm=model;
-mdm=mesh(mdm,'DomainOutline.exp',150000);
+mdm=mesh(model,'DomainOutline.exp',150000);
 mdm=geography(mdm,'','');
 mdm=parameterize(mdm,'Square.par');
@@ -16 +15 @@
 
 %Compute solution with Ice model
-mdm=solve(mdm,'diagnostic','ice');
+mdm.cluster=oshostname;
+mdm=solve(mdm,'analysis_type','diagnostic');
 
 save modelmacayeal mdm
@@ -22 +22 @@
 % Create pattyn model;
 
-mdp=model;
-mdp=mesh(mdp,'DomainOutline.exp',150000);
+mdp=mesh(model,'DomainOutline.exp',150000);
 mdp=geography(mdp,'','');
 mdp=parameterize(mdp,'Square.par');
@@ -30 +29 @@
 
 %Compute solution with Ice model
-mdp=solve(mdp,'diagnostic','ice');
+mdp=solve(mdp,'analysis_type','diagnostic');
 
 save modelpattyn mdp
@@ -36 +35 @@
 % Create stokes model;
 
-mds=model;
-mds=mesh(mds,'DomainOutline.exp',150000);
+mds=mesh(model,'DomainOutline.exp',150000);
 mds=geography(mds,'','');
 mds=parameterize(mds,'Square.par');
@@ -44 +42 @@
 
 %Compute solution with Ice model
-mds=solve(mds,'diagnostic','ice');
+mds=solve(mds,'analysis_type','diagnostic');
 
 save modelstokes mds