Changeset 14045

Timestamp:

11/28/12 15:38:57 (12 years ago)

Author:

Eric.Larour

Message:

CHG: last minute changes

Location:

issm/workshop/2012/Talks/09_PIG_Application

Files:

• Exp_Par/Pig.par (modified) (3 diffs)
• runme.m (modified) (7 diffs)

issm/workshop/2012/Talks/09_PIG_Application/Exp_Par/Pig.par ¶

@@ +1 @@
+% Parameters to change/Try
+friction_coefficient = 10; % default [10]
+Temp_change          =  0;  % default [0 K]
+
 %Name and hemisphere
 md.miscellaneous.name='PIG';
@@ -49 +53 @@
 
 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 +60 @@
 
 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

issm/workshop/2012/Talks/09_PIG_Application/runme.m ¶

@@ -1 @@
-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 +35 @@
         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 +73 @@
         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 +106 @@
         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 +112 @@
         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 +118 @@
         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 +130 @@
         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