Changeset 20733

Timestamp:

06/14/16 11:42:36 (9 years ago)

Author:

seroussi

Message:

CHG: finished update pig test for the workshop

Location:

issm/trunk-jpl/examples/Pig

Files:

• CheatSheet.m (modified) (2 diffs)
• Pig.par (modified) (6 diffs)
• runme.m (modified) (8 diffs)

issm/trunk-jpl/examples/Pig/CheatSheet.m ¶

@@ -1 @@
-if any(steps==6)
-        % Load Model
+if step==6
         md = loadmodel('./Models/PIG.Control_drag');
         md.inversion.iscontrol=0;
@@ -6 +5 @@
         disp('   Extruding mesh')
         number_of_layers=3;
-        md=extrude(md,number_of_layers,0.9);
+        md=extrude(md,number_of_layers,1);
 
         disp('   Using HO Ice Flow Model')
         md=setflowequation(md, 'HO', 'all');
 
-        % Solve
         md=solve(md,StressbalanceSolutionEnum);
 
-        % Save Model
         save ./Models/PIG.ModelHO md;
 end

issm/trunk-jpl/examples/Pig/Pig.par ¶

@@ -17 +17 @@
 smb   = ncread(ncdata,'presprcp')';
 gflux = ncread(ncdata,'bheatflx_fox')';
-
-disp('   Loading velocities data from NetCDF');
-nsidc_vel='../Data/Antarctica_ice_velocity.nc';
-
-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
-
-ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-ymax = strtrim(ymax);  
-ymax = ymax(1:end-2);  
-ymax = str2num(ymax); 
-
-nx = ncreadatt(nsidc_vel,'/','nx');
-ny = ncreadatt(nsidc_vel,'/','ny');
-
-spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-spacing = strtrim(spacing);  
-spacing = spacing(1:end-2);  
-spacing = str2num(spacing); 
-
-velx = double(ncread(nsidc_vel,'vx'));
-vely = double(ncread(nsidc_vel,'vy'));
-
-x2=xmin+(0:1:nx)'*spacing; 
-x2=double(x2);
-
-y2=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-y2=double(y2);
 
 %Geometry
@@ -70 +40 @@
 
 %Set min thickness to 1 meter
-pos0=find(md.geometry.thickness<=0);
+pos0=find(md.geometry.thickness<=1);
 md.geometry.thickness(pos0)=1;
 md.geometry.surface=md.geometry.thickness+md.geometry.base;
@@ -79 +49 @@
 clear temp;
 
+disp('   Loading velocities data from NetCDF');
+nsidc_vel='../Data/Antarctica_ice_velocity.nc';
+xmin    = ncreadatt(nsidc_vel,'/','xmin');
+ymax    = ncreadatt(nsidc_vel,'/','ymax');
+spacing = ncreadatt(nsidc_vel,'/','spacing');
+nx      = double(ncreadatt(nsidc_vel,'/','nx'));
+ny      = double(ncreadatt(nsidc_vel,'/','ny'));
+velx    = double(ncread(nsidc_vel,'vx'));
+vely    = double(ncread(nsidc_vel,'vy'));
+% Read coordinates
+xmin = strtrim(xmin);  
+xmin = str2num(xmin(1:end-2)); 
+ymax = strtrim(ymax);  
+ymax = str2num(ymax(1:end-2));  
+spacing = strtrim(spacing);
+spacing = str2num(spacing(1:end-2));  
+% Build the coordinates
+x2=xmin+(0:1:nx)'*spacing;
+y2=(ymax-ny*spacing)+(0:1:ny)'*spacing;
+
 disp('   Set observed velocities')
-vx_obs=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0);
-vy_obs=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0);
+md.initialization.vx=InterpFromGridToMesh(x2,y2,flipud(velx'),md.mesh.x,md.mesh.y,0);
+md.initialization.vy=InterpFromGridToMesh(x2,y2,flipud(vely'),md.mesh.x,md.mesh.y,0);
+md.initialization.vz=zeros(md.mesh.numberofvertices,1);
+md.initialization.vel=sqrt(md.initialization.vx.^2+md.initialization.vy.^2);
 clear velx vely;
-
-vel_obs=sqrt(vx_obs.^2+vy_obs.^2);
-md.initialization.vx=vx_obs;
-md.initialization.vy=vy_obs;
-md.initialization.vz=zeros(md.mesh.numberofvertices,1);
-md.initialization.vel=vel_obs;
 
 disp('   Set Pressure');
@@ -99 +85 @@
 %Forcings
 disp('   Interpolating surface mass balance');
-md.smb.mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
-md.smb.mass_balance=md.smb.mass_balance*md.materials.rho_water/md.materials.rho_ice;
+mass_balance=InterpFromGridToMesh(x1,y1,smb,md.mesh.x,md.mesh.y,0);
+md.smb.mass_balance=mass_balance*md.materials.rho_water/md.materials.rho_ice;
 clear smb;
 
@@ -118 +104 @@
 
 md.inversion=m1qn3inversion();
-md.inversion.vx_obs=vx_obs;
-md.inversion.vy_obs=vy_obs;
-md.inversion.vel_obs=vel_obs;
+md.inversion.vx_obs=md.initialization.vx;
+md.inversion.vy_obs=md.initialization.vy;
+md.inversion.vel_obs=md.initialization.vel;
 
 disp('   Set boundary conditions');
@@ -126 +112 @@
 md.basalforcings.floatingice_melting_rate = zeros(md.mesh.numberofvertices,1);
 md.basalforcings.groundedice_melting_rate = zeros(md.mesh.numberofvertices,1);
-md.thermal.spctemperature     = [md.initialization.temperature;1]; %impose observed temperature on surface
-md.masstransport.spcthickness    = NaN*ones(md.mesh.numberofvertices,1);
+md.thermal.spctemperature                 = md.initialization.temperature;
+md.masstransport.spcthickness             = NaN*ones(md.mesh.numberofvertices,1);

issm/trunk-jpl/examples/Pig/runme.m ¶

@@ -1 @@
-%Which steps to be performed
-steps=[1] ;
+step=[7];
 
-%Run Steps
-
-%Mesh Generation #1
-if any(steps==1)
-
-        md.miscellaneous.name='PIG.Mesh_generation';
-
+if step==1   %Mesh Generation #1
         %Mesh parameters
         domain =['./DomainOutline.exp'];
+        hinit=10000;   % element size for the initial mesh
         hmax=40000;    % maximum element size of the final mesh
         hmin=5000;     % minimum element size of the final mesh
-        hinit=10000;   % element size for the initial mesh
         gradation=1.7; % maximum size ratio between two neighboring elements
         err=8;         % maximum error between interpolated and control field
 
         % Generate an initial uniform mesh (resolution = hinit m)
-        md=bamg(model,'domain',domain,'hmax',hinit,'MaxCornerAngle',1);
+        md=bamg(model,'domain',domain,'hmax',hinit);
 
         %ploting
@@ -24 +17 @@
 
         % Load Velocities
-        % http://nsidc.org/data/nsidc-0484.html
         nsidc_vel='../Data/Antarctica_ice_velocity.nc';         
 
         % Get necessary data to build up the velocity grid
-        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    = ncreadatt(nsidc_vel,'/','xmin');
+        ymax    = ncreadatt(nsidc_vel,'/','ymax');
+        spacing = ncreadatt(nsidc_vel,'/','spacing');
+        nx      = double(ncreadatt(nsidc_vel,'/','nx'));
+        ny      = double(ncreadatt(nsidc_vel,'/','ny'));
+        vx      = double(ncread(nsidc_vel,'vx'));
+        vy      = double(ncread(nsidc_vel,'vy'));
+
+        % Read coordinates
+        xmin = strtrim(xmin);  
+        xmin = str2num(xmin(1:end-2)); 
+        ymax = strtrim(ymax);  
+        ymax = str2num(ymax(1:end-2));  
+        spacing = strtrim(spacing);
+        spacing = str2num(spacing(1:end-2));  
+
+        % Build the coordinates
+        x=xmin+(0:1:nx)'*spacing;
+        y=(ymax-ny*spacing)+(0:1:ny)'*spacing;
         
-        ymax = ncreadatt(nsidc_vel,'/','ymax'); 
-        ymax = strtrim(ymax);  
-        ymax = ymax(1:end-2);  
-        ymax = str2num(ymax); 
-        
-        nx = ncreadatt(nsidc_vel,'/','nx');
-        ny = ncreadatt(nsidc_vel,'/','ny');
-        
-        spacing = ncreadatt(nsidc_vel,'/','spacing'); 
-        spacing = strtrim(spacing);
-        spacing = spacing(1:end-2);  
-        spacing = str2num(spacing); 
-        
-        % Get velocities (Note: You can use ncdisp('file') to see an ncdump)
-        vx = double(ncread(nsidc_vel,'vx'));
-        vy = double(ncread(nsidc_vel,'vy'));
-        
-        x=xmin+(0:1:nx)'*spacing; 
-        x=double(x);
-        y=(ymax-ny*spacing)+(0:1:ny)'*spacing; 
-        y=double(y);
-
         % Interpolate velocities onto coarse mesh
         vx_obs=InterpFromGridToMesh(x,y,flipud(vx'),md.mesh.x,md.mesh.y,0);
@@ -67 +52 @@
         plotmodel(md,'data','mesh')
 
-        % Convert x,y coordinates (Polar stereo) to lat/lon
-        [md.mesh.lat,md.mesh.long]=xy2ll(md.mesh.x,md.mesh.y,-1);
-        
         % Save model
         save ./Models/PIG.Mesh_generation md;
 end
 
-%Masks #2
-if any(steps==2) 
+if step==2  %Masks #2
 
         md = loadmodel('./Models/PIG.Mesh_generation'); 
@@ -104 +85 @@
 end
 
-%Parameterization #3
-if any(steps==3) 
+if step==3  %Parameterization #3
 
         md = loadmodel('./Models/PIG.SetMask');
@@ -117 +97 @@
 end
 
-%Control Method #4
-if any(steps==4)
+if step==4  %Control Method #4
 
         md = loadmodel('./Models/PIG.Parameterization');
@@ -161 +140 @@
 end
 
-%Plot #5
-if any(steps==5)
+if step==5 %Plot #5
 
         md = loadmodel('./Models/PIG.Control_drag');
@@ -179 +157 @@
 end
 
-%HO #6
-if any(steps==6)
+if step==6  %Higher-Order #6
 
         % Load Model
@@ -196 +173 @@
 end
 
-%Plot #7
-if any(steps==7)
+if step==7  %Plot #7
 
         mdHO = loadmodel('./Models/PIG.ModelHO');