Changeset 4986


Ignore:
Timestamp:
08/04/10 15:29:53 (15 years ago)
Author:
Mathieu Morlighem
Message:

updated runme

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/scripts/runme.m

    r3206 r4986  
    1 %AGU 09: second attempt (after Moca09) at  control methods on entire Antarctica.
     1%SaeRISE Antarctica
    22
    33%Run steps
    4 steps=[3];
     4steps=[1 2 3];
    55
    66%steps available:
     
    1313modeldatapath='./../../Data/';
    1414%clustername=oshostname();
    15 clustername='cosmos';
     15clustername='wilkes';
    1616type='2d';
    1717
     
    2424        cluster.name='larsen';   cluster.np=8;    cluster.queue='long';  cluster.time=12*60;
    2525else
    26         cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
     26        cluster.name=oshostname; cluster.np=8;   cluster.queue='long';  cluster.time=12*60;
     27        %cluster.name=oshostname; cluster.np=16;   cluster.queue='long';  cluster.time=12*60;
    2728end%}}}
    2829
     
    3031t1=clock;
    3132
    32 %Mesh generation{{{1
     33%Mesh generation {{{1
    3334num=1;
    3435if ismember(num,steps),
    3536        message(rs,num)
     37       
     38        %Mesh parameters
     39        resolution=15000;
    3640
    37         %Mesh parameters
    38         domain =['DomainOutline.exp'];
    39         hmin=3000;
    40         hmax=15*10^3;
    41         gradation=2;
    42         err=[4 40];
     41        %Initial mesh
     42        md=model;
     43        md=mesh(md,'DomainOutline.exp',resolution);
    4344
    44         %Initil mesh creation
    45         disp('   Interpolating fields');
    46         md=model;
    47         md=bamg(md,'domain',domain,'hmax',hmin,'splitcorner',0,'KeepVertices',0);
    48 
    49         %Mesh adaptation
    50         disp('   Interpolating fields');
    51         vxpath   =[modeldatapath '/surfvelx.mat'];
    52         vypath   =[modeldatapath '/surfvely.mat'];
    53         md.vx_obs=InterpFromFile(md.x,md.y,vxpath,0);
    54         md.vy_obs=InterpFromFile(md.x,md.y,vypath,0);
    55         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    56         thicknesspath=[modeldatapath '/thk.mat'];
    57         md.thickness=InterpFromFile(md.x,md.y,thicknesspath,0);
    58 
    59         %Complete with balanced velocities
    60         surfacepath  =[modeldatapath '/usrf.mat'];
    61         balancedvelpath=[modeldatapath '/balvelmag.mat'];
    62         md.surface=InterpFromFile(md.x,md.y,surfacepath,0);
    63         [md.surface_slopex md.surface_slopey]=slope(md);
    64         md.surface_slopex=averaging(md,md.surface_slopex,0);
    65         md.surface_slopey=averaging(md,md.surface_slopey,0);
    66         balancedvelocities=InterpFromFile(md.x,md.y,balancedvelpath,0);
    67         velx=balancedvelocities.*md.surface_slopex./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
    68         vely=balancedvelocities.*md.surface_slopey./sqrt(md.surface_slopex.^2+md.surface_slopey.^2);
    69         pos=find(md.vel_obs==0);
    70         md.vx_obs(pos)=velx(pos);
    71         md.vy_obs(pos)=vely(pos);
    72         md.vel_obs=sqrt(md.vx_obs.^2+md.vy_obs.^2);
    73 
    74         disp('   Mesh adaptation');
    75         md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',[md.vel_obs md.thickness],'err',[4 40],'splitcorner',0,'KeepVertices',0);
    76        
    7745        savemodel(rs,num,md);
    7846end %}}}
    79 %Parameterization{{{1
     47%Parameterization {{{1
    8048num=num+1;
    8149if ismember(num,steps),
     
    8351        md=loadmodel(rs,'Mesh generation');
    8452
    85         %{
    86         landcoverfile=[modeldatapath '/landcover.mat'];
    87         load(landcoverfile);
    88         land_type=InterpFromGridToMesh(x1,y1,landcover,md.x,md.y,0);
     53        %Find ice sheets and ice shelves
     54        landcoveringpath=[modeldatapath '/thkmask.mat'];
     55        land_type=InterpFromFile(md.x,md.y,landcoveringpath,0);
    8956        gridonicesheet=zeros(md.numberofgrids,1);
    90         gridonicesheet(find(land_type>2))=1;
     57        gridonicesheet(find(land_type>0))=1;
    9158        elementonicesheet=zeros(md.numberofelements,1);
    9259        elementonicesheet(find(sum(gridonicesheet(md.elements(:,:)),2)==3))=1;
    93         md=geography(md,double(~elementonicesheet),'');
    94         %}
    95         md=geography(md,'','');
    96         md=parameterize(md,'Greenland.par');
     60        elementoniceshelf=1-elementonicesheet;
     61        md=geography(md,elementoniceshelf,'');
     62
     63        md=parameterize(md,'Antarctica.par');
    9764        md=setelementstype(md,'macayeal','all');
    98         md=setcluster(md,cluster);
    9965       
    10066        savemodel(rs,num,md);
    101 end%}}}
    102 %Control Method drag 1{{{1
     67end %}}}
     68%Control method drag 1 {{{
    10369num=num+1;
    10470if ismember(num,steps),
     
    10672        md=loadmodel(rs,'Parameterization');
    10773
    108         md.eps_res=0.01;
    109         md.eps_rel=0.1;
     74        disp('Compute 2d drag coefficient');
     75
    11076        md.eps_abs=NaN;
    11177        md.eps_cm=NaN;
    11278        md.verbose=0;
    113 
    114         %50  (2-0) cmmax=120: misfit(md)  = 8.6417
    115         %100 (2-3) cmmax=120: misfit(md)  = 7.4280
    116         %100 (2-3) cmmax=200: misfit(md)  = 7.4146
    117         %100 (4-3) cmmax=200: misfit(md)  = 8.6725
    118         %100 (2-0) cmmax=200: misfit(md)  = 8.4770
    119         %100 (2-3) cmmax=200: misfit(md)  = 7.3538
    120         md=parametercontroldrag(md,'nsteps',100);
    121         %fit=[0;0;2;2];
    122         %md.fit=repmat(fit,md.nsteps,1);
    123         md.fit=[];
    124         md.fit(1:floor(md.nsteps/2),1)=2;
    125         md.fit(floor(md.nsteps/2):md.nsteps,1)=3;
    126         md.optscal(1:floor(md.nsteps/2))=500;
    127         md.optscal(floor(md.nsteps/2):md.nsteps)=100;
     79       
     80        md.drag(:)=30;
     81        md.control_analysis=1;
     82        md.nsteps=10;
     83        md.maxiter=10*ones(md.nsteps,1);
     84        fit=[0;0;2;0];
     85        md.fit=repmat(fit,md.nsteps,1);
     86        md.fit(1:30)=2;
     87        md.fit(31:60)=0;
     88        md.fit(md.nsteps+1:end)=[];
     89        md.optscal=30*ones(md.nsteps,1);
    12890        md.cm_jump=.99*ones(md.nsteps,1);
    12991        md.control_type='drag';
     
    13193        md.cm_min=1;
    13294        md.cm_max=200;
    133         md.cm_noisedmp=10^-15; %13 is too much
    13495
    13596        %spc all 0 velocities.
     
    13798        npos=length(pos);
    13899        md.spcvelocity(pos,:)=[ones(npos,1) ones(npos,1) ones(npos,1) zeros(npos,1) zeros(npos,1) zeros(npos,1)];
    139         in=ContourToMesh(md.elements,md.x,md.y,expread('trackscontours.exp',1),'node',1);
    140         md.weights(find(in))=0.1;
    141100
    142101        md=setcluster(md,cluster);
    143         md.control_analysis=1;
    144102        md=solve(md,'analysis_type','diagnostic');
    145103        md.drag=md.results.diagnostic.parameter;
    146104
    147105        savemodel(rs,num,md);
    148 end%}}}
     106end %}}}
    149107
    150108%timing
Note: See TracChangeset for help on using the changeset viewer.