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

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
Line 34: Line 34:
 
clear thkall uvelall vvelall
 
clear thkall uvelall vvelall
  
xn = repmat( [0:5e3:2e5] / 1e3, 41, 1 );        %% make a grid for plots (NOTE: you'll have to change this
+
xn = repmat( [0:10e3:2e5] / 1e3, 21, 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 !!! )
+
yn = repmat( [0:10e3:2e5]' / 1e3, 1, 21);        %% is you change the grid spacing in the config file !!! )
  
xs = repmat( [2.5e3:5e3:2e5-2.5e3] / 1e3, 40, 1 );
+
xs = repmat( [5e3:10e3:2e5-5e3] / 1e3, 20, 1 );
ys = repmat( [2.5e3:5e3:2e5-2.5e3]' / 1e3, 1, 40 );
+
ys = repmat( [5e3:10e3:2e5-5e3]' / 1e3, 1, 20 );
  
  

Revision as of 12:48, 7 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:10e3:2e5] / 1e3, 21, 1 );        %% make a grid for plots (NOTE: you'll have to change this
yn = repmat( [0:10e3:2e5]' / 1e3, 1, 21);        %% is you change the grid spacing in the config file !!! )
 
xs = repmat( [5e3:10e3:2e5-5e3] / 1e3, 20, 1 );
ys = repmat( [5e3:10e3:2e5-5e3]' / 1e3, 1, 20 );
 
 
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)' )