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

Last change on this file was 25986, checked in by jdquinn, 4 years ago

CHG: Cleanup

  • Property svn:executable set to *
File size: 5.7 KB
RevLine 
[21129]1steps=[1];
[18198]2
[25986]3if any(steps==1) %Mesh Generation #1
[18198]4 %Mesh parameters
5 domain =['./DomainOutline.exp'];
[25915]6 hinit=10000; % element size for the initial mesh
7 hmax=40000; % maximum element size of the final mesh
8 hmin=5000; % minimum element size of the final mesh
9 gradation=1.7; % maximum size ratio between two neighboring elements
10 err=8; % maximum error between interpolated and control field
[18198]11
12 % Generate an initial uniform mesh (resolution = hinit m)
[20733]13 md=bamg(model,'domain',domain,'hmax',hinit);
[18198]14
15 % Load Velocities
[25915]16 nsidc_vel='../Data/Antarctica_ice_velocity.nc';
[18198]17
18 % Get necessary data to build up the velocity grid
[25915]19 xmin = ncreadatt(nsidc_vel,'/','xmin');
20 ymax = ncreadatt(nsidc_vel,'/','ymax');
21 spacing = ncreadatt(nsidc_vel,'/','spacing');
22 nx = double(ncreadatt(nsidc_vel,'/','nx'));
23 ny = double(ncreadatt(nsidc_vel,'/','ny'));
24 vx = double(ncread(nsidc_vel,'vx'));
25 vy = double(ncread(nsidc_vel,'vy'));
[20733]26
27 % Read coordinates
[25915]28 xmin = strtrim(xmin);
29 xmin = str2num(xmin(1:end-2));
30 ymax = strtrim(ymax);
31 ymax = str2num(ymax(1:end-2));
[18198]32 spacing = strtrim(spacing);
[25915]33 spacing = str2num(spacing(1:end-2));
[20733]34
35 % Build the coordinates
36 x=xmin+(0:1:nx)'*spacing;
37 y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
[25915]38
[18198]39 % Interpolate velocities onto coarse mesh
40 vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
41 vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0);
42 vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
43 clear vx vy x y;
44
45 % Adapt the mesh to minimize error in velocity interpolation
46 md=bamg(md,'hmax',hmax,'hmin',hmin,'gradation',gradation,'field',vel_obs,'err',err);
[25915]47
[18198]48 %ploting
49 plotmodel(md,'data','mesh')
50
51 % Save model
[20754]52 save ./Models/PIG_Mesh_generation md;
[25986]53end
[18198]54
[25986]55if any(steps==2) %Masks #2
[25915]56 md = loadmodel('./Models/PIG_Mesh_generation');
[18198]57
[25915]58 % Load SeaRISe dataset for Antarctica
[18198]59 % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
60 searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
[25915]61
[18198]62 %read thickness mask from SeaRISE
63 x1=double(ncread(searise,'x1'));
64 y1=double(ncread(searise,'y1'));
65 thkmask=double(ncread(searise,'thkmask'));
[25915]66
[18198]67 %interpolate onto our mesh vertices
68 groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
69 groundedice(groundedice<=0)=-1;
70 clear thkmask;
71
72 %fill in the md.mask structure
[24864]73 md.mask.ocean_levelset=groundedice; %ice is grounded for mask equal one
[18198]74 md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie
75
76 %ploting
[24864]77 plotmodel(md,'data',md.mask.ocean_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
[25915]78
[18198]79 % Save model
[20754]80 save ./Models/PIG_SetMask md;
[25986]81end
[18198]82
[25986]83if any(steps==3) %Parameterization #3
[20754]84 md = loadmodel('./Models/PIG_SetMask');
[18198]85 md = parameterize(md,'./Pig.par');
86
87 % Use a MacAyeal flow model
88 md = setflowequation(md,'SSA','all');
[25915]89
[18198]90 % Save model
[20754]91 save ./Models/PIG_Parameterization md;
[25986]92end
[18198]93
[25986]94if any(steps==4) %Control Method #4
[20754]95 md = loadmodel('./Models/PIG_Parameterization');
[18198]96
97 % Control general
98 md.inversion.iscontrol=1;
[21607]99 md.inversion.maxsteps=20;
[18198]100 md.inversion.maxiter=40;
101 md.inversion.dxmin=0.1;
102 md.inversion.gttol=1.0e-4;
[21607]103 md.verbose=verbose('control',true);
[18198]104
105 % Cost functions
106 md.inversion.cost_functions=[101 103 501];
107 md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
108 md.inversion.cost_functions_coefficients(:,1)=1;
109 md.inversion.cost_functions_coefficients(:,2)=1;
110 md.inversion.cost_functions_coefficients(:,3)=8e-15;
111
112 % Controls
113 md.inversion.control_parameters={'FrictionCoefficient'};
114 md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
115 md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
116
117 % Additional parameters
118 md.stressbalance.restol=0.01;
119 md.stressbalance.reltol=0.1;
120 md.stressbalance.abstol=NaN;
121
122 % Solve
123 md.toolkits=toolkits;
124 md.cluster=generic('name',oshostname,'np',2);
[21057]125 md=solve(md,'Stressbalance');
[18198]126
127 % Update model friction fields accordingly
128 md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
129
130 plotmodel(md,'data',md.friction.coefficient)
131
132 % Save model
[20754]133 save ./Models/PIG_Control_drag md;
[25986]134end
[18198]135
[25986]136if any(steps==5) %Plot #5
[20754]137 md = loadmodel('./Models/PIG_Control_drag');
[18198]138
[21607]139 plotmodel(md,...
[18198]140 'data',md.initialization.vel,'title','Observed velocity',...
141 'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
142 'data',md.geometry.base,'title','Bed elevation',...
143 'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
[21607]144 'colorbar#all','on','colorbartitle#1-2','(m/yr)',...
[18198]145 'caxis#1-2',([1.5,4000]),...
[21607]146 'colorbartitle#3','(m)', 'log#1-2',10);
[25986]147end
[18198]148
[25986]149if any(steps==6) %Higher-Order #6
[18198]150 % Load Model
[20951]151
[18269]152 % Disable inversion
[20951]153
[18260]154 % Extrude Mesh
[20951]155
[18269]156 % Set Flowequation
[20951]157
[18198]158 % Solve
[20951]159
[18198]160 % Save Model
161
[25915]162end % step 6 end
[18198]163
[25986]164if any(steps==7) %Plot #7
[20754]165 mdHO = loadmodel('./Models/PIG_ModelHO');
166 mdSSA = loadmodel('./Models/PIG_Control_drag');
[18198]167
168 basal=find(mdHO.mesh.vertexonbase);
169 surf=find(mdHO.mesh.vertexonsurface);
170
[20795]171 plotmodel(mdHO,'nlines',3,'ncols',2,'axis#all','equal',...
[25915]172 'data',mdHO.initialization.vel,'title','Observed velocity',...
173 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
174 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
175 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
176 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
177 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
178 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
179 'colorbar#all','on','view#all',2,...
180 'colorbartitle#all','(m/yr)',...
181 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
[25986]182end
Note: See TracBrowser for help on using the repository browser.