| 1 | %Which steps to be performed
|
|---|
| 2 | steps=[1] ;
|
|---|
| 3 |
|
|---|
| 4 | %Run Steps
|
|---|
| 5 |
|
|---|
| 6 | %Mesh Generation #1
|
|---|
| 7 | if 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;
|
|---|
| 74 | end
|
|---|
| 75 |
|
|---|
| 76 | %Masks #2
|
|---|
| 77 | if any(steps==2)
|
|---|
| 78 |
|
|---|
| 79 | md = loadmodel('./Models/PIG.Mesh_generation');
|
|---|
| 80 |
|
|---|
| 81 | % Load SeaRISe dataset for Antarctica
|
|---|
| 82 | % http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica
|
|---|
| 83 | searise='../Data/Antarctica_5km_withshelves_v0.75.nc';
|
|---|
| 84 |
|
|---|
| 85 | %read thickness mask from SeaRISE
|
|---|
| 86 | x1=double(ncread(searise,'x1'));
|
|---|
| 87 | y1=double(ncread(searise,'y1'));
|
|---|
| 88 | thkmask=double(ncread(searise,'thkmask'));
|
|---|
| 89 |
|
|---|
| 90 | %interpolate onto our mesh vertices
|
|---|
| 91 | groundedice=double(InterpFromGridToMesh(x1,y1,thkmask',md.mesh.x,md.mesh.y,0));
|
|---|
| 92 | groundedice(groundedice<=0)=-1;
|
|---|
| 93 | clear thkmask;
|
|---|
| 94 |
|
|---|
| 95 | %fill in the md.mask structure
|
|---|
| 96 | md.mask.groundedice_levelset=groundedice; %ice is grounded for mask equal one
|
|---|
| 97 | md.mask.ice_levelset=-1*ones(md.mesh.numberofvertices,1);%ice is present when negatvie
|
|---|
| 98 |
|
|---|
| 99 | %ploting
|
|---|
| 100 | plotmodel(md,'data',md.mask.groundedice_levelset,'title','grounded/floating','data',md.mask.ice_levelset,'title','ice/no-ice')
|
|---|
| 101 |
|
|---|
| 102 | % Save model
|
|---|
| 103 | save ./Models/PIG.SetMask md;
|
|---|
| 104 | end
|
|---|
| 105 |
|
|---|
| 106 | %Parameterization #3
|
|---|
| 107 | if any(steps==3)
|
|---|
| 108 |
|
|---|
| 109 | md = loadmodel('./Models/PIG.SetMask');
|
|---|
| 110 | md = parameterize(md,'./Pig.par');
|
|---|
| 111 |
|
|---|
| 112 | % Use a MacAyeal flow model
|
|---|
| 113 | md = setflowequation(md,'SSA','all');
|
|---|
| 114 |
|
|---|
| 115 | % Save model
|
|---|
| 116 | save ./Models/PIG.Parameterization md;
|
|---|
| 117 | end
|
|---|
| 118 |
|
|---|
| 119 | %Control Method #4
|
|---|
| 120 | if any(steps==4)
|
|---|
| 121 |
|
|---|
| 122 | md = loadmodel('./Models/PIG.Parameterization');
|
|---|
| 123 |
|
|---|
| 124 | % Control general
|
|---|
| 125 | md.inversion.iscontrol=1;
|
|---|
| 126 | md.inversion.maxsteps=20;
|
|---|
| 127 | md.inversion.maxiter=40;
|
|---|
| 128 | md.inversion.dxmin=0.1;
|
|---|
| 129 | md.inversion.gttol=1.0e-4;
|
|---|
| 130 | md.verbose=verbose('solution',true,'control',true);
|
|---|
| 131 |
|
|---|
| 132 | % Cost functions
|
|---|
| 133 | md.inversion.cost_functions=[101 103 501];
|
|---|
| 134 | md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,3);
|
|---|
| 135 | md.inversion.cost_functions_coefficients(:,1)=1;
|
|---|
| 136 | md.inversion.cost_functions_coefficients(:,2)=1;
|
|---|
| 137 | md.inversion.cost_functions_coefficients(:,3)=8e-15;
|
|---|
| 138 |
|
|---|
| 139 | % Controls
|
|---|
| 140 | md.inversion.control_parameters={'FrictionCoefficient'};
|
|---|
| 141 | md.inversion.min_parameters=1*ones(md.mesh.numberofvertices,1);
|
|---|
| 142 | md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1);
|
|---|
| 143 |
|
|---|
| 144 | % Additional parameters
|
|---|
| 145 | md.stressbalance.restol=0.01;
|
|---|
| 146 | md.stressbalance.reltol=0.1;
|
|---|
| 147 | md.stressbalance.abstol=NaN;
|
|---|
| 148 |
|
|---|
| 149 | % Solve
|
|---|
| 150 | md.toolkits=toolkits;
|
|---|
| 151 | md.cluster=generic('name',oshostname,'np',2);
|
|---|
| 152 | md=solve(md,StressbalanceSolutionEnum);
|
|---|
| 153 |
|
|---|
| 154 | % Update model friction fields accordingly
|
|---|
| 155 | md.friction.coefficient=md.results.StressbalanceSolution.FrictionCoefficient;
|
|---|
| 156 |
|
|---|
| 157 | plotmodel(md,'data',md.friction.coefficient)
|
|---|
| 158 |
|
|---|
| 159 | % Save model
|
|---|
| 160 | save ./Models/PIG.Control_drag md;
|
|---|
| 161 | end
|
|---|
| 162 |
|
|---|
| 163 | %Plot #5
|
|---|
| 164 | if any(steps==5)
|
|---|
| 165 |
|
|---|
| 166 | md = loadmodel('./Models/PIG.Control_drag');
|
|---|
| 167 |
|
|---|
| 168 | plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',...
|
|---|
| 169 | 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,...
|
|---|
| 170 | 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,...
|
|---|
| 171 | 'FontSize#all',12,...
|
|---|
| 172 | 'data',md.initialization.vel,'title','Observed velocity',...
|
|---|
| 173 | 'data',md.results.StressbalanceSolution.Vel,'title','Modeled Velocity',...
|
|---|
| 174 | 'data',md.geometry.base,'title','Bed elevation',...
|
|---|
| 175 | 'data',md.results.StressbalanceSolution.FrictionCoefficient,'title','Friction Coefficient',...
|
|---|
| 176 | 'colorbar#all','on','colorbartitle#1-2','[m/yr]',...
|
|---|
| 177 | 'caxis#1-2',([1.5,4000]),...
|
|---|
| 178 | 'colorbartitle#3','[m]', 'log#1-2',10);
|
|---|
| 179 | end
|
|---|
| 180 |
|
|---|
| 181 | %HO #6
|
|---|
| 182 | if any(steps==6)
|
|---|
| 183 |
|
|---|
| 184 | % Load Model
|
|---|
| 185 |
|
|---|
| 186 | % Disable inversion
|
|---|
| 187 |
|
|---|
| 188 | % Extrude Mesh
|
|---|
| 189 |
|
|---|
| 190 | % Set Flowequation
|
|---|
| 191 |
|
|---|
| 192 | % Solve
|
|---|
| 193 |
|
|---|
| 194 | % Save Model
|
|---|
| 195 |
|
|---|
| 196 | end
|
|---|
| 197 |
|
|---|
| 198 | %Plot #7
|
|---|
| 199 | if any(steps==7)
|
|---|
| 200 |
|
|---|
| 201 | mdHO = loadmodel('./Models/PIG.ModelHO');
|
|---|
| 202 | mdSSA = loadmodel('./Models/PIG.Control_drag');
|
|---|
| 203 |
|
|---|
| 204 | basal=find(mdHO.mesh.vertexonbase);
|
|---|
| 205 | surf=find(mdHO.mesh.vertexonsurface);
|
|---|
| 206 |
|
|---|
| 207 | plotmodel(mdHO,'nlines',3,'ncols',2,'unit#all','km','axis#all','equal',...
|
|---|
| 208 | 'xlim#all',[min(mdHO.mesh.x) max(mdHO.mesh.x)]/10^3,...
|
|---|
| 209 | 'ylim#all',[min(mdHO.mesh.y) max(mdHO.mesh.y)]/10^3,...
|
|---|
| 210 | 'FontSize#all',12,...
|
|---|
| 211 | 'data',mdHO.initialization.vel,'title','Observed velocity',...
|
|---|
| 212 | 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.initialization.vel(surf)),'title','(HO-observed) velocities',...
|
|---|
| 213 | 'data',mdSSA.results.StressbalanceSolution.Vel,'title','Modeled SSA Velocity',...
|
|---|
| 214 | 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdSSA.results.StressbalanceSolution.Vel),'title','(HO-SSA) velocities',...
|
|---|
| 215 | 'data',mdHO.results.StressbalanceSolution.Vel,'title','Modeled HO surface Velocities',...
|
|---|
| 216 | 'data',(mdHO.results.StressbalanceSolution.Vel(surf)-mdHO.results.StressbalanceSolution.Vel(basal)),'title','(HOsurf-HO base) velocities',...
|
|---|
| 217 | 'caxis#1',([1.5,4000]),'caxis#3',([1.5,4000]),'caxis#5',([1.5,4000]),...
|
|---|
| 218 | 'colorbar#all','on','view#all',2,...
|
|---|
| 219 | 'colorbartitle#all','[m/yr]',...
|
|---|
| 220 | 'layer#5',1, 'log#1', 10,'log#3', 10,'log#5', 10);
|
|---|
| 221 | end
|
|---|
![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
