Changeset 4986

Timestamp:

08/04/10 15:29:53 (15 years ago)

Author:

Mathieu Morlighem

Message:

updated runme

File:

• issm/trunk/scripts/runme.m (modified) (8 diffs)

issm/trunk/scripts/runme.m

@@ -1 @@
-%AGU 09: second attempt (after Moca09) at  control methods on entire Antarctica.
+%SaeRISE Antarctica
 
 %Run steps
-steps=[3];
+steps=[1 2 3];
 
 %steps available: 
@@ -13 +13 @@
 modeldatapath='./../../Data/';
 %clustername=oshostname();
-clustername='cosmos';
+clustername='wilkes';
 type='2d';
 
@@ -24 +24 @@
         cluster.name='larsen';   cluster.np=8;    cluster.queue='long';  cluster.time=12*60;
 else
-        cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
+        cluster.name=oshostname; cluster.np=8;   cluster.queue='long';  cluster.time=12*60;
+        %cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
 end%}}}
 
@@ -30 +31 @@
 t1=clock;
 
-%Mesh generation{{{1
+%Mesh generation {{{1
 num=1;
 if ismember(num,steps),
         message(rs,num)
+        
+        %Mesh parameters
+        resolution=15000;
 
-        %Mesh parameters
-        domain =['DomainOutline.exp'];
-        hmin=3000;
-        hmax=15*10^3;
-        gradation=2;
-        err=[4 40];
+        %Initial mesh
+        md=model;
+        md=mesh(md,'DomainOutline.exp',resolution);
 
-        %Initil mesh creation
-        disp('   Interpolating fields'); 
-        md=model;
-        md=bamg(md,'domain',domain,'hmax',hmin,'splitcorner',0,'KeepVertices',0);
-
-        %Mesh adaptation
-        disp('   Interpolating fields'); 
-        vxpath   =[modeldatapath '/surfvelx.mat'];
-        vypath   =[modeldatapath '/surfvely.mat'];
-        md.vx_obs=InterpFromFile(md.x,md.y,vxpath,0);
-        md.vy_obs=InterpFromFile(md.x,md.y,vypath,0);
-        md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-        thicknesspath=[modeldatapath '/thk.mat'];
-        md.thickness=InterpFromFile(md.x,md.y,thicknesspath,0);
-
-        %Complete with balanced velocities
-        surfacepath  =[modeldatapath '/usrf.mat'];
-        balancedvelpath=[modeldatapath '/balvelmag.mat'];
-        md.surface=InterpFromFile(md.x,md.y,surfacepath,0);
-        [md.surface_slopex md.surface_slopey]=slope(md);
-        md.surface_slopex=averaging(md,md.surface_slopex,0);
-        md.surface_slopey=averaging(md,md.surface_slopey,0);
-        balancedvelocities=InterpFromFile(md.x,md.y,balancedvelpath,0);
-        velx=balancedvelocities.*md.surface_slopex./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
-        vely=balancedvelocities.*md.surface_slopey./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
-        pos=find(md.vel_obs==0);
-        md.vx_obs(pos)=velx(pos);
-        md.vy_obs(pos)=vely(pos);
-        md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
-
-        disp('   Mesh adaptation'); 
-        md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',[md.vel_obs md.thickness],'err',[4 40],'splitcorner',0,'KeepVertices',0);
-        
         savemodel(rs,num,md);
 end %}}}
-%Parameterization{{{1
+%Parameterization {{{1
 num=num+1;
 if ismember(num,steps),
@@ -83 +51 @@
         md=loadmodel(rs,'Mesh generation');
 
-        %{
-        landcoverfile=[modeldatapath '/landcover.mat'];
-        load(landcoverfile);
-        land_type=InterpFromGridToMesh(x1,y1,landcover,md.x,md.y,0);
+        %Find ice sheets and ice shelves
+        landcoveringpath=[modeldatapath '/thkmask.mat'];
+        land_type=InterpFromFile(md.x,md.y,landcoveringpath,0);
         gridonicesheet=zeros(md.numberofgrids,1);
-        gridonicesheet(find(land_type>2))=1;
+        gridonicesheet(find(land_type>0))=1;
         elementonicesheet=zeros(md.numberofelements,1);
         elementonicesheet(find(sum(gridonicesheet(md.elements(:,:)),2)==3))=1;
-        md=geography(md,double(~elementonicesheet),'');
-        %}
-        md=geography(md,'','');
-        md=parameterize(md,'Greenland.par');
+        elementoniceshelf=1-elementonicesheet;
+        md=geography(md,elementoniceshelf,'');
+
+        md=parameterize(md,'Antarctica.par');
         md=setelementstype(md,'macayeal','all');
-        md=setcluster(md,cluster);
         
         savemodel(rs,num,md);
-end%}}}
-%Control Method drag 1{{{1
+end %}}}
+%Control method drag 1 {{{
 num=num+1;
 if ismember(num,steps),
@@ -106 +72 @@
         md=loadmodel(rs,'Parameterization');
 
-        md.eps_res=0.01;
-        md.eps_rel=0.1;
+        disp('Compute 2d drag coefficient');
+
         md.eps_abs=NaN;
         md.eps_cm=NaN;
         md.verbose=0;
-
-        %50  (2-0) cmmax=120: misfit(md)  = 8.6417
-        %100 (2-3) cmmax=120: misfit(md)  = 7.4280
-        %100 (2-3) cmmax=200: misfit(md)  = 7.4146
-        %100 (4-3) cmmax=200: misfit(md)  = 8.6725
-        %100 (2-0) cmmax=200: misfit(md)  = 8.4770
-        %100 (2-3) cmmax=200: misfit(md)  = 7.3538
-        md=parametercontroldrag(md,'nsteps',100);
-        %fit=[0;0;2;2];
-        %md.fit=repmat(fit,md.nsteps,1);
-        md.fit=[];
-        md.fit(1:floor(md.nsteps/2),1)=2;
-        md.fit(floor(md.nsteps/2):md.nsteps,1)=3;
-        md.optscal(1:floor(md.nsteps/2))=500;
-        md.optscal(floor(md.nsteps/2):md.nsteps)=100;
+        
+        md.drag(:)=30;
+        md.control_analysis=1;
+        md.nsteps=10;
+        md.maxiter=10*ones(md.nsteps,1);
+        fit=[0;0;2;0];
+        md.fit=repmat(fit,md.nsteps,1);
+        md.fit(1:30)=2;
+        md.fit(31:60)=0;
+        md.fit(md.nsteps+1:end)=[];
+        md.optscal=30*ones(md.nsteps,1);
         md.cm_jump=.99*ones(md.nsteps,1);
         md.control_type='drag';
@@ -131 +93 @@
         md.cm_min=1;
         md.cm_max=200;
-        md.cm_noisedmp=10^-15; %13 is too much
 
         %spc all 0 velocities.
@@ -137 +98 @@
         npos=length(pos);
         md.spcvelocity(pos,:)=[ones(npos,1) ones(npos,1) ones(npos,1) zeros(npos,1) zeros(npos,1) zeros(npos,1)];
-        in=ContourToMesh(md.elements,md.x,md.y,expread('trackscontours.exp',1),'node',1);
-        md.weights(find(in))=0.1;
 
         md=setcluster(md,cluster);
-        md.control_analysis=1;
         md=solve(md,'analysis_type','diagnostic');
         md.drag=md.results.diagnostic.parameter;
 
         savemodel(rs,num,md);
-end%}}}
+end %}}}
 
 %timing