source: issm/trunk-jpl/examples/Pig/runme.m@ 18202

Last change on this file since 18202 was 18202, checked in by schlegel, 11 years ago

CHG: update runmes for tutorials
  • Property svn:executable set to *
File size: 6.7 KB
Line 
1%Which steps to be performed
2steps=[1:5] ;
3
4%Run Steps
5
6% {{{ Mesh Generation #1
7if any(steps==1)
8
9 md.miscellaneous.name='PIG.Mesh_generation';
10
11 %Mesh parameters
12 domain =['./DomainOutline.exp'];
13 hmax=40000; % maximum element size of the final mesh
14 hmin=5000; % minimum element size of the final mesh
15 hinit=10000; % element size for the initial mesh
16 gradation=1.7; % maximum size ratio between two neighboring elements
17 err=8; % maximum error between interpolated and control field
18
19 % Generate an initial uniform mesh (resolution = hinit m)
20 md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);
21
22 %ploting
23 plotmodel(md,'data','mesh')
24
25 % Load Velocities
26 % http://nsidc.org/data/nsidc-0484.html
27 nsidc_vel='../Data/Antarctica_ice_velocity.nc';
28
29 % Get necessary data to build up the velocity grid
30 xmin = ncreadatt(nsidc_vel,'/','xmin');
31 xmin = strtrim(xmin); % this is a string, and we need to recover the double value
32 xmin = xmin(1:end-2); % get rid of the unit
33 xmin = str2num(xmin); % convert to double
34
35 ymax = ncreadatt(nsidc_vel,'/','ymax');
36 ymax = strtrim(ymax);
37 ymax = ymax(1:end-2);
38 ymax = str2num(ymax);
39
40 nx = ncreadatt(nsidc_vel,'/','nx');
41 ny = ncreadatt(nsidc_vel,'/','ny');
42
43 spacing = ncreadatt(nsidc_vel,'/','spacing');
44 spacing = strtrim(spacing);
45 spacing = spacing(1:end-2);
46 spacing = str2num(spacing);
47
48 % Get velocities (Note: You can use ncdisp('file') to see an ncdump)
49 vx = double(ncread(nsidc_vel,'vx'));
50 vy = double(ncread(nsidc_vel,'vy'));
51
52 x=xmin+(0:1:nx)'*spacing;
53 x=double(x);
54 y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
55 y=double(y);
56
57 % Interpolate velocities onto coarse mesh
58 vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
59 vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
60 vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
61 clear vx vy x y;
62
63 % Adapt the mesh to minimize error in velocity interpolation
64 md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
65
66 %ploting
67 plotmodel(md,'data','mesh')
68
69 % Convert x,y coordinates (Polar stereo) to lat/lon
70 [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
71
72 % Save model
73 save ./Models/PIG.Mesh_generation md;
74end
75% }}}
76
77% {{{ Masks #2
78if any(steps==2)
79
80 md = loadmodel('./Models/PIG.Mesh_generation');
81
82 % Load SeaRISe dataset for Antarctica
83 % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
84 searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
85
86 %read thickness mask from SeaRISE
87 x1=double(ncread(searise,'x1'));
88 y1=double(ncread(searise,'y1'));
89 thkmask=double(ncread(searise,'thkmask'));
90
91 %interpolate onto our mesh vertices
92 groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
93 groundedice(groundedice<=0)=-1;
94 clear thkmask;
95
96 %fill in the md.mask structure
97 md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one
98 md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie
99
100 %ploting
101 plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
102
103 % Save model
104 save ./Models/PIG.SetMask md;
105end
106% }}}
107
108% {{{ Parameterization #3
109if any(steps==3)
110
111 md = loadmodel('./Models/PIG.SetMask');
112 md = parameterize(md,'./Pig.par');
113
114 % Use a MacAyeal flow model
115 md = setflowequation(md,'SSA','all');
116
117 % Save model
118 save ./Models/PIG.Parameterization md;
119end
120% }}}
121
122% {{{ Control Method #4
123if any(steps==4)
124
125 md = loadmodel('./Models/PIG.Parameterization');
126
127 % Control general
128 md.inversion.iscontrol=1;
129 md.inversion.maxsteps=20;
130 md.inversion.maxiter=40;
131 md.inversion.dxmin=0.1;
132 md.inversion.gttol=1.0e-4;
133 md.verbose=verbose('solution',true,'control',true);
134
135 % Cost functions
136 md.inversion.cost_functions=[101 103 501];
137 md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
138 md.inversion.cost_functions_coefficients(:,1)=1;
139 md.inversion.cost_functions_coefficients(:,2)=1;
140 md.inversion.cost_functions_coefficients(:,3)=8e-15;
141
142 % Controls
143 md.inversion.control_parameters={'FrictionCoefficient'};
144 md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
145 md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
146
147 % Additional parameters
148 md.stressbalance.restol=0.01;
149 md.stressbalance.reltol=0.1;
150 md.stressbalance.abstol=NaN;
151
152 % Solve
153 md.toolkits=toolkits;
154 md.cluster=generic('name',oshostname,'np',2);
155 md=solve(md,StressbalanceSolutionEnum);
156
157 % Update model friction fields accordingly
158 md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
159
160 plotmodel(md,'data',md.friction.coefficient)
161
162 % Save model
163 save ./Models/PIG.Control_drag md;
164end
165% }}}
166
167% {{{ Plot #5
168if any(steps==5)
169
170 md = loadmodel('./Models/PIG.Control_drag');
171
172 plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
173 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
174 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
175 'FontSize#all',12,...
176 'data',md.initialization.vel,'title','Observed velocity',...
177 'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
178 'data',md.geometry.base,'title','Bed elevation',...
179 'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
180 'colorbar#all','on','colorbartitle#1-2','[m/yr]',...
181 'caxis#1-2',([1.5,4000]),...
182 'colorbartitle#3','[m]', 'log#1-2',10);
183end
184% }}}
185
186% {{{ HO #6
187if any(steps==6)
188
189 % Load Model
190
191 % Disable inversion
192
193 %Extrude Mesh
194
195 % Set Flowequation
196
197 % Solve
198
199 % Save Model
200
201end
202% }}}
203
204% {{{ Plot #7
205if any(steps==7)
206
207 mdHO = loadmodel('./Models/PIG.ModelHO');
208 mdSSA = loadmodel('./Models/PIG.Control_drag');
209
210 basal=find(mdHO.mesh.vertexonbase);
211 surf=find(mdHO.mesh.vertexonsurface);
212
213 plotmodel(md,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',...
214 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
215 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
216 'FontSize#all',12,...
217 'data',mdHO.initialization.vel,'title','Observed velocity',...
218 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
219 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
220 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
221 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
222 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
223 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
224 'colorbar#all','on','view#all',2,...
225 'colorbartitle#all','[m/yr]',...
226 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
227end
228% }}}
Note: See TracBrowser for help on using the repository browser.