Changeset 22793

Timestamp:

05/17/18 12:19:20 (7 years ago)

Author:

adhikari

Message:

CHG: tutorial 1 completed

Location:

issm/trunk-jpl/examples/EsaWahr

Files:

• RoundDomain.exp (added)
• runme.m (modified) (8 diffs)
• wahr.m (added)

issm/trunk-jpl/examples/EsaWahr/runme.m

@@ -1 +1 @@
 
 clear all;
-steps=[4];
+steps=[1:6]; % [1:6]; 
 
-if any(steps==1)
-        disp('   Step 1: Mesh creation');
+if any(steps==0)
+        disp('   Step 0: Mesh creation');
 
-        %Generate initial uniform mesh 
+        % Generate initial uniform mesh 
         md=roundmesh(model,100000,10000);  % Domain radius and element size (meters) 
 
         save ./Models/EsaWahr.Mesh md;
         
+        plotmodel(md,'data','mesh');
+        %export_fig('Fig0.pdf'); 
+
+end
+
+if any(steps==1)
+        disp('   Step 1: Anisotropic mesh creation');
+
+        % Generate initial uniform mesh 
+        md=roundmesh(model,100000,1000);  % Domain radius and element size (meters) 
+        %plotmodel(md,'data','mesh'); 
+
+        % Define a fake field to have anisotropic meshing 
+        disc_radius=20*1000; % km => m 
+        rad_dist=sqrt(md.mesh.x.^2+md.mesh.y.^2);       % radial distance [m] 
+        field = abs(rad_dist-disc_radius); 
+        %plotmodel(md,'data',field); 
+
+        md = bamg(md,'field',field,'err',50,'hmax',10000); % error in field in m 
+
+        save ./Models/EsaWahr.Mesh md;
+        
         plotmodel (md,'data','mesh');
+        %export_fig('Fig1.pdf'); 
 
 end
@@ -20 +43 @@
 
         % primer 
-        nv=md.mesh.numberofvertices; 
-        ne=md.mesh.numberofelements; 
         rho_w_i=md.materials.rho_freshwater/md.materials.rho_ice; 
 
@@ -30 +51 @@
 
         % Define a 20-km radius disc and unload it by 1 meter water height equivalent uniformly  
-        md.esa.deltathickness = zeros(ne,1); 
+        md.esa.deltathickness = zeros(md.mesh.numberofelements,1); 
         disc_radius=20; % km 
         rad_dist=sqrt(x_cent.^2+y_cent.^2)/1000;        % radial distance in km 
@@ -38 +59 @@
         
         plotmodel (md,'data',md.esa.deltathickness,'title','Ice height equivalent [m]');
+        %export_fig('Fig2.pdf'); 
 
 end
@@ -67 +89 @@
         % Miscellaneous: 
         md.miscellaneous.name='EsaWahr';
+        %% IGNORE BUT DO NOT DELETE %% 
         
         save ./Models/EsaWahr.Parameterization md; 
@@ -86 +109 @@
         md=solve(md,'Esa');
 
-        save ./Models/EsaWahr.EsaSolutions md; 
+        save ./Models/EsaWahr.Solution md; 
 
 end
 
-return; 
+if any(steps==5)
+        disp('   Step 5: Plot solutions');
+        md = loadmodel('./Models/EsaWahr.Solution');
 
-% high-res mesh...
-% plot step. 
-% analytical step. 
-
-        
+        % m => mm 
         vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
-        horz_n = md.results.EsaSolution.EsaNmotion*1000; % [mm] 
-        horz_e = md.results.EsaSolution.EsaEmotion*1000; % [mm] 
+        horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 
+        horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 
         horz = sqrt(horz_n.^2+horz_e.^2); 
 
-        % grid data from the center 
-        xi=[0:100:60000]; % 100 m resolution 
-        yi=zeros(1,length(xi)); 
-        vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
-        horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
-
-        return; 
-
-        % figure S1a {{{ 
+        % plot 
         set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
         figure('Position', [100, 100, 800, 600]);
@@ -116 +129 @@
                 'xTickLabel#all',[],'yTickLabel#all',[],...
                 'colormap#all','jet','colorbarcornerposition#all','south',...
-                'expdisp#all','../Exp/RoundDomain.exp',...
+                'expdisp#all','./RoundDomain.exp',...
                 'gridded#all',1,...
                 'axispos#1',[0.02 0.505 0.47 0.47],...
@@ -132 +145 @@
                 'axispos',[0.505 0.02 0.47 0.47],...
                 'colorbarpos',[0.53,0.065,0.18,0.02],'colorbartitle#4','East-west [mm]'); 
-        % }}} 
-        %export_fig('./Fig/Fig_S1a.pdf'); 
 
-        % figure S1b {{{
-        set(0,'DefaultAxesFontSize',24,'DefaultAxesLineWidth',1,'DefaultTextFontSize',24,'DefaultLineMarkerSize',6); 
+        %export_fig('Fig5.pdf'); 
+
+end
+
+if any(steps==6)
+        disp('   Step 6: Compare results against Wahr semi-analytic solutions');
+        md = loadmodel('./Models/EsaWahr.Solution');
+
+        % numerical solutions 
+        % m => mm 
+        vert = md.results.EsaSolution.EsaUmotion*1000; % [mm] 
+        horz_n = md.results.EsaSolution.EsaYmotion*1000; % [mm] 
+        horz_e = md.results.EsaSolution.EsaXmotion*1000; % [mm] 
+        horz = sqrt(horz_n.^2+horz_e.^2); 
+        % grid data from the center 
+        xi=[0:500:100000]; % 500 m resolution 
+        yi=zeros(1,length(xi)); 
+        vert_track=griddata(md.mesh.x,md.mesh.y,vert,xi,yi,'linear'); 
+        horz_track=griddata(md.mesh.x,md.mesh.y,horz,xi,yi,'linear'); 
+
+        % semi-analytic solution (Wahr et al., 2013, JGR, Figure 1) 
+        disc_radius = 20*1000; % 20 km 
+        [vert_wahr, horz_wahr] = wahr(disc_radius,xi,md.esa.love_h,md.esa.love_l);
+
+        % plot 
+        set(0,'DefaultAxesFontSize',16,'DefaultAxesLineWidth',1,'DefaultTextFontSize',16,'DefaultLineMarkerSize',6); 
         figure1=figure('Position', [100, 100, 700, 400]);
         ylabel_1={'0',' ','1','','2','','3',''}; 
         axes1 = axes('Layer','top','Position',[0.1 0.15 0.8 0.8],...
-                'XTick',[0:10:60],'xlim',[0 60],...
+                'XTick',[0:10:100],'xlim',[0 100],...
                 'ylim',[0 3.5],'Ytick',[0:0.5:3.5],'yticklabel',ylabel_1); 
                 box(axes1,'on'); hold(axes1,'all'); grid on; 
                 xlabel(axes1,'Radial distance [km]'); 
                 ylabel(axes1,'Displacement [mm]');
-                plot([20 20],[0 3.5],'-.k','linewidth',2,'parent',axes1); 
-                h1=plot(xi/1000,vert_track*917/1000,'color',[0.8500 0.3250 0.0980],'linewidth',3,'parent',axes1); 
-                h2=plot(xi/1000,horz_track*917/1000,'b','linewidth',3,'parent',axes1); 
+                plot([20 20],[0 3.5],'-k','linewidth',2,'parent',axes1); 
+                % analytic soln 
+                h3=plot(xi/1000,vert_wahr,'-r','linewidth',5,'parent',axes1); 
+                h4=plot(xi/1000,horz_wahr,'-m','linewidth',5,'parent',axes1); 
+                % ISSM 
+                h1=plot(xi/1000,vert_track*917/1000,'-b','linewidth',3,'parent',axes1); 
+                h2=plot(xi/1000,horz_track*917/1000,'-c','linewidth',3,'parent',axes1); 
                 % first box 
                 ag1 = gca;
-                leg1a = legend(ag1,[h1,h2],'Vertical (upward)','Horizontal (away from disc)',1);
-                set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',24); 
+                leg1a = legend(ag1,[h3,h1,h4,h2],'Vertical (Wahr)','Vertical (ISSM)','Horizontal (Wahr)','Horizontal (ISSM)',1);
+                set(leg1a,'location','east','Orientation','Vertical','Box','Off','FontSize',16); 
                 % 
         set(gcf,'color','w');
-        % }}} 
-        %export_fig('./Fig/Fig_S1b.pdf'); 
+        
+        %export_fig('Fig6.pdf'); 
 
+end