Changeset 4986
- Timestamp:
- 08/04/10 15:29:53 (15 years ago)
- 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 2 2 3 3 %Run steps 4 steps=[ 3];4 steps=[1 2 3]; 5 5 6 6 %steps available: … … 13 13 modeldatapath='./../../Data/'; 14 14 %clustername=oshostname(); 15 clustername=' cosmos';15 clustername='wilkes'; 16 16 type='2d'; 17 17 … … 24 24 cluster.name='larsen'; cluster.np=8; cluster.queue='long'; cluster.time=12*60; 25 25 else 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; 27 28 end%}}} 28 29 … … 30 31 t1=clock; 31 32 32 %Mesh generation {{{133 %Mesh generation {{{1 33 34 num=1; 34 35 if ismember(num,steps), 35 36 message(rs,num) 37 38 %Mesh parameters 39 resolution=15000; 36 40 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); 43 44 44 %Initil mesh creation45 disp(' Interpolating fields');46 md=model;47 md=bamg(md,'domain',domain,'hmax',hmin,'splitcorner',0,'KeepVertices',0);48 49 %Mesh adaptation50 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 velocities60 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 77 45 savemodel(rs,num,md); 78 46 end %}}} 79 %Parameterization {{{147 %Parameterization {{{1 80 48 num=num+1; 81 49 if ismember(num,steps), … … 83 51 md=loadmodel(rs,'Mesh generation'); 84 52 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); 89 56 gridonicesheet=zeros(md.numberofgrids,1); 90 gridonicesheet(find(land_type> 2))=1;57 gridonicesheet(find(land_type>0))=1; 91 58 elementonicesheet=zeros(md.numberofelements,1); 92 59 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'); 97 64 md=setelementstype(md,'macayeal','all'); 98 md=setcluster(md,cluster);99 65 100 66 savemodel(rs,num,md); 101 end %}}}102 %Control Method drag 1{{{167 end %}}} 68 %Control method drag 1 {{{ 103 69 num=num+1; 104 70 if ismember(num,steps), … … 106 72 md=loadmodel(rs,'Parameterization'); 107 73 108 md.eps_res=0.01;109 md.eps_rel=0.1; 74 disp('Compute 2d drag coefficient'); 75 110 76 md.eps_abs=NaN; 111 77 md.eps_cm=NaN; 112 78 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); 128 90 md.cm_jump=.99*ones(md.nsteps,1); 129 91 md.control_type='drag'; … … 131 93 md.cm_min=1; 132 94 md.cm_max=200; 133 md.cm_noisedmp=10^-15; %13 is too much134 95 135 96 %spc all 0 velocities. … … 137 98 npos=length(pos); 138 99 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;141 100 142 101 md=setcluster(md,cluster); 143 md.control_analysis=1;144 102 md=solve(md,'analysis_type','diagnostic'); 145 103 md.drag=md.results.diagnostic.parameter; 146 104 147 105 savemodel(rs,num,md); 148 end %}}}106 end %}}} 149 107 150 108 %timing
Note:
See TracChangeset
for help on using the changeset viewer.