[21129] | 1 | steps=[1];
|
---|
[18198] | 2 |
|
---|
[25986] | 3 | if 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] | 53 | end
|
---|
[18198] | 54 |
|
---|
[25986] | 55 | if 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] | 81 | end
|
---|
[18198] | 82 |
|
---|
[25986] | 83 | if 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] | 92 | end
|
---|
[18198] | 93 |
|
---|
[25986] | 94 | if 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] | 134 | end
|
---|
[18198] | 135 |
|
---|
[25986] | 136 | if 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] | 147 | end
|
---|
[18198] | 148 |
|
---|
[25986] | 149 | if 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] | 162 | end % step 6 end
|
---|
[18198] | 163 |
|
---|
[25986] | 164 | if 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] | 182 | end
|
---|