Difference between revisions of "Script for comparing output from blob experiments"

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
(New page: <source lang=matlab> %% for spreading blob of ice experiments, compare thickness evolution %% profiles from 1st-order upwinding scheme w/ those from a more %% sophisticated scheme (increme...)
 
Line 20: Line 20:
 
vvelall = netcdf.getVar(ncid,vvelid);
 
vvelall = netcdf.getVar(ncid,vvelid);
  
thk0 = thkall(:,:,2);          %% get initial fields for which vels are nonzero
+
thk0 = thkall(:,:,1);          %% get initial fields for which vels are nonzero
 
uvel0 = uvelall(:,:,1,2);
 
uvel0 = uvelall(:,:,1,2);
 
vvel0 = vvelall(:,:,1,2);
 
vvel0 = vvelall(:,:,1,2);
  
thk1 = thkall(:,:,end);        %% get same fields from end of run
+
thk1 = thkall(:,:,50);        %% get same fields from middle of run
uvel1 = uvelall(:,:,1,end);
+
uvel1 = uvelall(:,:,1,50);
vvel1 = vvelall(:,:,1,end);
+
vvel1 = vvelall(:,:,1,50);
 +
 
 +
thk2 = thkall(:,:,end);        %% get same fields from end of run
 +
uvel2 = uvelall(:,:,1,end);
 +
vvel2 = vvelall(:,:,1,end);
  
 
clear thkall uvelall vvelall
 
clear thkall uvelall vvelall
  
xn = repmat( [0:5e3:2e5] / 1e3, 41, 1 );        %% make a grid for plots
+
xn = repmat( [0:5e3:2e5] / 1e3, 41, 1 );        %% make a grid for plots (NOTE: you'll have to change this
yn = repmat( [0:5e3:2e5]' / 1e3, 1, 41);
+
yn = repmat( [0:5e3:2e5]' / 1e3, 1, 41);       %% is you change the grid spacing in the config file !!! )
  
 
xs = repmat( [2.5e3:5e3:2e5-2.5e3] / 1e3, 40, 1 );
 
xs = repmat( [2.5e3:5e3:2e5-2.5e3] / 1e3, 40, 1 );
Line 39: Line 43:
 
figure(1), clf, hold on                                    %% plot INITIAL fields
 
figure(1), clf, hold on                                    %% plot INITIAL fields
 
subplot(2,2,1), hold on
 
subplot(2,2,1), hold on
contourf( xn, yn, thk0, 100 ), shading flat, colorbar( 'horiz' )
+
contourf( xn, yn, thk0, 100 ), shading flat, colorbar
 
contour( xn, yn, thk0, [0:100:1000], 'k-' )
 
contour( xn, yn, thk0, [0:100:1000], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init thickness (m)' )
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init thickness (m)' )
 
subplot(2,2,2), hold on
 
subplot(2,2,2), hold on
contourf( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), 100 ), shading flat, colorbar('horiz')
+
contourf( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), 100 ), shading flat, colorbar
 
contour( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), [0:10:100], 'k-' )
 
contour( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), [0:10:100], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
Line 51: Line 55:
 
figure(1), hold on                                          %% plot FINAL fields
 
figure(1), hold on                                          %% plot FINAL fields
 
subplot(2,2,3), hold on
 
subplot(2,2,3), hold on
contourf( xn, yn, thk1, 100 ), shading flat, colorbar( 'horiz' )
+
contourf( xn, yn, thk2, 100 ), shading flat, colorbar
contour( xn, yn, thk1, [0:100:1000], 'k-' )
+
contour( xn, yn, thk2, [0:100:1000], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
 
subplot(2,2,4), hold on
 
subplot(2,2,4), hold on
contourf( xs, ys, sqrt( uvel1.^2 + vvel1.^2 ), 100 ), shading flat, colorbar('horiz')
+
contourf( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel1.^2 + vvel1.^2 ), [0:10:100], 'k-' )
+
contour( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), [0:10:100], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
Line 64: Line 68:
 
figure(3), clf, hold on
 
figure(3), clf, hold on
 
plot( xn(20,:), thk0(20,:), 'b-' )
 
plot( xn(20,:), thk0(20,:), 'b-' )
plot( xn(20,:), thk1(20,:), 'r-' )
+
plot( xn(20,:), thk1(20,:), 'g-' )
 +
plot( xn(20,:), thk2(20,:), 'r-' )
 +
legend( 'initial(fo)', 'middle(fo)', 'final(fo)' )
 
xlabel( 'x (km)' ), ylabel( 'z (m)' ), title( 'elevation (m)' ), box on
 
xlabel( 'x (km)' ), ylabel( 'z (m)' ), title( 'elevation (m)' ), box on
 
 
  
  
Line 84: Line 88:
 
vvelall = netcdf.getVar(ncid,vvelid);
 
vvelall = netcdf.getVar(ncid,vvelid);
  
thk0 = thkall(:,:,1);
+
thk0 = thkall(:,:,1);         %% get initial fields for which vels are nonzero
 
uvel0 = uvelall(:,:,1,2);
 
uvel0 = uvelall(:,:,1,2);
 
vvel0 = vvelall(:,:,1,2);
 
vvel0 = vvelall(:,:,1,2);
  
thk1 = thkall(:,:,end);
+
thk1 = thkall(:,:,50);       %% get same fields from middle of run
uvel1 = uvelall(:,:,1,end);
+
uvel1 = uvelall(:,:,1,50);
vvel1 = vvelall(:,:,1,end);
+
vvel1 = vvelall(:,:,1,50);
  
clear thkall uvelall vvelall
+
thk2 = thkall(:,:,end);        %% get same fields from end of run
 +
uvel2 = uvelall(:,:,1,end);
 +
vvel2 = vvelall(:,:,1,end);
  
 +
clear thkall uvelall vvelall
  
 
figure(2), clf, hold on                                            %% plot INITIAL fields
 
figure(2), clf, hold on                                            %% plot INITIAL fields
Line 109: Line 116:
 
figure(2), hold on                                                    %% plot FINAL fields
 
figure(2), hold on                                                    %% plot FINAL fields
 
subplot(2,2,3), hold on
 
subplot(2,2,3), hold on
contourf( xn, yn, thk1, 100 ), shading flat, colorbar
+
contourf( xn, yn, thk2, 100 ), shading flat, colorbar
contour( xn, yn, thk1, [0:100:1000], 'k-' )
+
contour( xn, yn, thk2, [0:100:1000], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
 
subplot(2,2,4), hold on
 
subplot(2,2,4), hold on
contourf( xs, ys, sqrt( uvel1.^2 + vvel1.^2 ), 100 ), shading flat, colorbar
+
contourf( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel1.^2 + vvel1.^2 ), [0:10:100], 'k-' )
+
contour( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), [0:10:100], 'k-' )
 
axis xy, axis equal, axis tight, box on
 
axis xy, axis equal, axis tight, box on
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
 
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
Line 122: Line 129:
 
%% plot initial and final thickness (profile)
 
%% plot initial and final thickness (profile)
 
figure(3), hold on
 
figure(3), hold on
plot( xn(20,:), thk0(20,:), 'b--' )
+
plot( xn(20,:), thk0(20,:), 'b.--' )
plot( xn(20,:), thk1(20,:), 'r--' )
+
plot( xn(20,:), thk1(20,:), 'g.--' )
 
+
plot( xn(20,:), thk2(20,:), 'r.--' )
legend( 'initial(fo)', 'final(fo)', 'initial(remap)', 'final(remap)' )
+
legend( 'initial(fo)', 'middle(fo)', 'final(fo)', 'initial(remap)', 'middle(remap)', 'final(remap)' )
 
</source>
 
</source>

Revision as of 15:18, 6 August 2009

%% for spreading blob of ice experiments, compare thickness evolution
%% profiles from 1st-order upwinding scheme w/ those from a more
%% sophisticated scheme (incremental remapping)
 
 
%% first, data from fo scheme
clear all; close all
 
ncid = netcdf.open( 'blob.out.fo.nc', 'NC_NOWRITE' );
 
%% get field ids
thkid = netcdf.inqVarID( ncid, 'thk' );
uvelid = netcdf.inqVarID( ncid, 'uvelhom' );
vvelid = netcdf.inqVarID( ncid, 'vvelhom' );
 
%% get fields
thkall = netcdf.getVar(ncid,thkid);
uvelall = netcdf.getVar(ncid,uvelid);
vvelall = netcdf.getVar(ncid,vvelid);
 
thk0 = thkall(:,:,1);          %% get initial fields for which vels are nonzero
uvel0 = uvelall(:,:,1,2);
vvel0 = vvelall(:,:,1,2);
 
thk1 = thkall(:,:,50);        %% get same fields from middle of run
uvel1 = uvelall(:,:,1,50);
vvel1 = vvelall(:,:,1,50);
 
thk2 = thkall(:,:,end);        %% get same fields from end of run
uvel2 = uvelall(:,:,1,end);
vvel2 = vvelall(:,:,1,end);
 
clear thkall uvelall vvelall
 
xn = repmat( [0:5e3:2e5] / 1e3, 41, 1 );        %% make a grid for plots (NOTE: you'll have to change this
yn = repmat( [0:5e3:2e5]' / 1e3, 1, 41);        %% is you change the grid spacing in the config file !!! )
 
xs = repmat( [2.5e3:5e3:2e5-2.5e3] / 1e3, 40, 1 );
ys = repmat( [2.5e3:5e3:2e5-2.5e3]' / 1e3, 1, 40 );
 
 
figure(1), clf, hold on                                     %% plot INITIAL fields
subplot(2,2,1), hold on
contourf( xn, yn, thk0, 100 ), shading flat, colorbar
contour( xn, yn, thk0, [0:100:1000], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init thickness (m)' )
subplot(2,2,2), hold on
contourf( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), [0:10:100], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init speed (m/a)' )
 
figure(1), hold on                                          %% plot FINAL fields
subplot(2,2,3), hold on
contourf( xn, yn, thk2, 100 ), shading flat, colorbar
contour( xn, yn, thk2, [0:100:1000], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
subplot(2,2,4), hold on
contourf( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), [0:10:100], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
 
%% plot initial and final thickness (profile)
figure(3), clf, hold on
plot( xn(20,:), thk0(20,:), 'b-' )
plot( xn(20,:), thk1(20,:), 'g-' )
plot( xn(20,:), thk2(20,:), 'r-' )
legend( 'initial(fo)', 'middle(fo)', 'final(fo)' )
xlabel( 'x (km)' ), ylabel( 'z (m)' ), title( 'elevation (m)' ), box on
 
 
 
%% second, data from remapping scheme
ncid = netcdf.open( 'blob.out.remap.nc', 'NC_NOWRITE' );
 
%% get field ids
thkid = netcdf.inqVarID( ncid, 'thk' );
uvelid = netcdf.inqVarID( ncid, 'uvelhom' );
vvelid = netcdf.inqVarID( ncid, 'vvelhom' );
 
%% get fields
thkall = netcdf.getVar(ncid,thkid);
uvelall = netcdf.getVar(ncid,uvelid);
vvelall = netcdf.getVar(ncid,vvelid);
 
thk0 = thkall(:,:,1);          %% get initial fields for which vels are nonzero
uvel0 = uvelall(:,:,1,2);
vvel0 = vvelall(:,:,1,2);
 
thk1 = thkall(:,:,50);        %% get same fields from middle of run
uvel1 = uvelall(:,:,1,50);
vvel1 = vvelall(:,:,1,50);
 
thk2 = thkall(:,:,end);        %% get same fields from end of run
uvel2 = uvelall(:,:,1,end);
vvel2 = vvelall(:,:,1,end);
 
clear thkall uvelall vvelall
 
figure(2), clf, hold on                                             %% plot INITIAL fields
subplot(2,2,1), hold on
contourf( xn, yn, thk0, 100 ), shading flat, colorbar
contour( xn, yn, thk0, [0:100:1000], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init thickness (m)' )
subplot(2,2,2), hold on
contourf( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel0.^2 + vvel0.^2 ), [0:10:100], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'init speed (m/a)' )
 
figure(2), hold on                                                    %% plot FINAL fields
subplot(2,2,3), hold on
contourf( xn, yn, thk2, 100 ), shading flat, colorbar
contour( xn, yn, thk2, [0:100:1000], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final thickness (m)' )
subplot(2,2,4), hold on
contourf( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), 100 ), shading flat, colorbar
contour( xs, ys, sqrt( uvel2.^2 + vvel2.^2 ), [0:10:100], 'k-' )
axis xy, axis equal, axis tight, box on
xlabel( 'x (km)' ), ylabel( 'y (km)' ), title( 'final speed (m/a)' )
 
 
%% plot initial and final thickness (profile)
figure(3), hold on
plot( xn(20,:), thk0(20,:), 'b.--' )
plot( xn(20,:), thk1(20,:), 'g.--' )
plot( xn(20,:), thk2(20,:), 'r.--' )
legend( 'initial(fo)', 'middle(fo)', 'final(fo)', 'initial(remap)', 'middle(remap)', 'final(remap)' )