Index: /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par =================================================================== --- /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par (revision 14044) +++ /issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par (revision 14045) @@ -1,2 +1,6 @@ +% Parameters to change/Try +friction_coefficient = 10; % default [10] +Temp_change = 0; % default [0 K] + %Name and hemisphere md.miscellaneous.name='PIG'; @@ -49,5 +53,5 @@ disp(' Interpolating temperatures'); -md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15; +md.initialization.temperature=InterpFromGridToMesh(x1,y1,temp,md.mesh.x,md.mesh.y,0)+273.15+Temp_change; disp(' Interpolating surface mass balance'); @@ -56,5 +60,5 @@ disp(' Construct basal friction parameters'); -md.friction.coefficient=10*ones(md.mesh.numberofvertices,1); +md.friction.coefficient=friction_coefficient*ones(md.mesh.numberofvertices,1); pos=find(md.mask.vertexonfloatingice); md.friction.coefficient(pos)=0; %no friction applied on floating ice Index: /issm/workshop/2012/Talks/09_PIG_Application/runme.m =================================================================== --- /issm/workshop/2012/Talks/09_PIG_Application/runme.m (revision 14044) +++ /issm/workshop/2012/Talks/09_PIG_Application/runme.m (revision 14045) @@ -1,24 +1,29 @@ -steps=[1 2 3 4]; - -%SeaRISe dataset for Antarctica. Can be downloaded at: -%http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica - +% MacAyeal Steps +steps=[1 2 3 4 5]; + +% Blatter-Pattyn Steps +% steps=[6 7]; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% SeaRISe dataset for Antarctica. Can be downloaded at: +% http://websrv.cs.umt.edu/isis/index.php/Present_Day_Antarctica searise='./Data/Antarctica_5km_withshelves_v0.75.nc'; -%Rignot et al, 2011. NSIDC available velocities (at http://nsidc.org/data/nsidc-0484.html) +% Rignot et al, 2011. NSIDC available velocities (at http://nsidc.org/data/nsidc-0484.html) nsidc_vel='./Data/Antarctica_ice_velocity.nc'; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Create Mesh if any(steps==1) disp(' Step 1: Mesh creation'); - - %Generate initial uniform mesh (resolution = 10000 m) + % Generate an initial uniform mesh (resolution = 10000 m) md=bamg(model,'domain','Exp_Par/DomainOutline.exp','hmax',10000,'MaxCornerAngle',1); % Get velocities (Note: You can use ncdisp('file') to see an ncdump) xmin = ncreadatt(nsidc_vel,'/','xmin'); - xmin = strtrim(xmin); %this is a string, and we need to recover the double value - xmin = xmin(1:end-2); %get rid of the unit - xmin = str2num(xmin); %convert to double + xmin = strtrim(xmin); % this is a string, and we need to recover the double value + xmin = xmin(1:end-2); % get rid of the unit + xmin = str2num(xmin); % convert to double ymax = ncreadatt(nsidc_vel,'/','ymax'); ymax = strtrim(ymax); ymax = ymax(1:end-2); ymax = str2num(ymax); nx = ncreadatt(nsidc_vel,'/','nx'); @@ -30,19 +35,24 @@ y=(ymax-ny*spacing)+(0:1:ny)'*spacing; y=double(y); - %interpolate velocities onto coarse mesh + % Interpolate velocities onto coarse mesh vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0); vy_obs=InterpFromGridToMesh(x,y,flipud(vy'),md.mesh.x,md.mesh.y,0); vel_obs=sqrt(vx_obs.^2+vy_obs.^2); - %Adapt the mesh to minimize error in velocity interpolation + % Adapt the mesh to minimize error in velocity interpolation md=bamg(md,'hmax',400000,'hmin',5000,'gradation',1.7,'field',vel_obs,'err',8); - %convert x,y coordinates (Polar stereo) to lat/lon + % Convert x,y coordinates (Polar stereo) to lat/lon [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,+1,39,71); + % Save model save ./Models/Pig.Mesh_generation md; end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Masks + if any(steps==2) + disp(' Step 2: Mask (floating vs grounded)'); md = loadmodel('./Models/Pig.Mesh_generation'); @@ -63,21 +73,31 @@ md=setmask(md,elementoniceshelf,''); + % Save model save ./Models/Pig.Mask md; end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Parameterization + if any(steps==3) + disp(' Step 3: Parameterization'); md = loadmodel('./Models/Pig.Mask'); - md = parameterize(md,'./Exp_Par/Pig.par'); + + % Use a MacAyeal flow model md = setflowequation(md,'macayeal','all'); - + + % Save model save ./Models/Pig.Parameterization md; end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Control Method on Friction Coefficient if any(steps==4) disp(' Step 4: Control method friction'); md = loadmodel('./Models/Pig.Parameterization'); - %Control general + % Control general md.inversion.iscontrol=1; md.inversion.nsteps=20; @@ -86,5 +106,5 @@ md.verbose=verbose('solution',true,'control',true); - %Cost functions + % Cost functions md.inversion.cost_functions=[[103*ones(md.inversion.nsteps/2,1); 101*ones(md.inversion.nsteps/2,1)] 501*ones(md.inversion.nsteps,1)]; md.inversion.cost_functions_coefficients=ones(md.mesh.numberofvertices,2); @@ -92,5 +112,5 @@ md.inversion.cost_functions_coefficients(:,2)=8e-15; - %Controls + % Controls md.inversion.control_parameters={'FrictionCoefficient'}; md.inversion.gradient_scaling=50*ones(md.inversion.nsteps,1); @@ -98,10 +118,10 @@ md.inversion.max_parameters=200*ones(md.mesh.numberofvertices,1); - %Additional parameters + % Additional parameters md.diagnostic.restol=0.01; md.diagnostic.reltol=0.1; md.diagnostic.abstol=NaN; - %Go solve + % Solve md.solver=addoptions(md.solver,NoneAnalysisEnum,asmoptions); md.solver=addoptions(md.solver,DiagnosticVertAnalysisEnum,jacobiasmoptions); @@ -110,7 +130,71 @@ md=solve(md,DiagnosticSolutionEnum); - %Update model friction fields accordingly + % Update model friction fields accordingly md.friction.coefficient=md.results.DiagnosticSolution.FrictionCoefficient; + % Save model save ./Models/Pig.Control_drag md; end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Plot MacAyeal Results +if any(steps==5) + + disp(' Plotting MacAyeal Model Results') + md = loadmodel('./Models/Pig.Control_drag'); + + plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',... + 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,... + 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,... + 'FontSize#all',12,... + 'data',md.initialization.vel,'title','Observed velocity',... + 'data',md.results.DiagnosticSolution.Vel,'title','Modeled Velocity',... + 'data',md.geometry.bed,'title','Bed elevation',... + 'data',md.results.DiagnosticSolution.FrictionCoefficient,'title','Friction Coefficient',... + 'colorbar#all','on','colorbartitle#1-2','[m/yr]',... + 'caxis#1-2',([1.5,4000]),... + 'colorbartitle#3','[m]', 'log#all',10); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Use a Blatter-Pattyn Flow Model +if any(steps==6) + + % Load Model + md = loadmodel('./Models/Pig.Control_drag'); + md.inversion.iscontrol=0; + + disp(' Extruding mesh') + number_of_layers=3; + md=extrude(md, number_of_layers, 0.9); + + disp(' Using Batter-Pattyn Ice Flow Model') + md=setflowequation(md, 'pattyn', 'all'); + + % Solve + md=solve(md,DiagnosticSolutionEnum); + + % Save Model + save ./Models/Pig.Control_drag_pattyn md; +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Plot Blatter-Pattyn Results +if any(steps==7) + + disp(' Plotting Pattyn Model Results') + md = loadmodel('./Models/Pig.Control_drag_pattyn'); + + plotmodel(md,'nlines',2,'ncols',2,'unit#all','km','axis#all','equal',... + 'xlim#all',[min(md.mesh.x) max(md.mesh.x)]/10^3,... + 'ylim#all',[min(md.mesh.y) max(md.mesh.y)]/10^3,... + 'FontSize#all',12,... + 'data',md.initialization.vel,'title','Observed velocity',... + 'data',md.results.DiagnosticSolution.Vel,'title','Modeled Velocity',... + 'data',md.geometry.bed,'title','Bed elevation',... + 'data',md.friction.coefficient,'title','Friction Coefficient',... + 'caxis#1-2',([1.5,4000]),... + 'colorbar#all','on','view#all',2,... + 'colorbartitle#1-2','[m/yr]','colorbartitle#3','[m]',... + 'layer#1-3',(md.mesh.numberoflayers),'layer#4',1, 'log#all', 10); +end