% ps07_3.m is the mfile answer to problem set 6 question 3 created by
% T. KENNA 19981202 1155. It does a first upwind differencing method 
% calculation of pure advection of a step function in a linear domain.
%
% ***NOTE*** You need a file called show.m to run this file.
% 01:01:15 Modif'd by DMG to get rid of those annoying clc's
% 02:11:01 Modif'd by DMG to fix a minor textual error.
% 02:11:08 Modif'd by DMG to remove clc's and the use of SHOW.M
% 06:11:09 Modif'd by DMG to correct the name of the file in the header
% Modif'd 2008:11:13 DMG to change the name of the file in the header

clear all
close all

global G
global g
global n

% PART 1
 
%clc
disp('We first start by desigining our 1-D system; to make it 25,000 km');
disp('long requires that with a dx of 50,000 m, we need an nx of 500.  In');
disp('order to insure that we choose a stable time step, we must return to');
disp('our Von Neumann stability analysis the criteria needed for u > 0 is');
disp('that c+2*d <= 0. Since there is no explicit diffusivity, K and also');
disp('d are both equal to zero.  With the variables selected we can now run');
disp('the model out for 10 years with u = .01 m/s. At the end of the tenth');
disp('year all model parameters will be printed on the graph as well as the');
disp('numeric and theoretical diffusivity which are calculated according to');
disp('the equations given in the notes and assignment. So let''s go.');
disp('<return>');
pause

tyr=24*3600*365.25; 
dx=50000; 
dt=200000; 
u=.01;
nt=10*tyr/dt;
nx=500;
K=0;
d=K*dt/dx^2;
c=u*dt/dx;

%			FUDM weights
w0=1-u*dt/dx;
wm=u*dt/dx;
%			Initialize arrays
g=zeros(nx,1);
g(nx/2)=1; 
gold=g; 
x=dx*[1:nx]/1.e6;
m=zeros(10,2);
j=1;
oldnow=0;

% Start the 10 year run

for i=1:nt
	now=ceil(i*dt/tyr);
	if now ~= oldnow
		subplot(3,1,1)
			p=plot(x,g,'g');
			set(p,'LineWidth',2)
			grid on
			axis([0 25 0 .1]);
			xlabel('Distance')
			ylabel('Concentration')
		subplot(3,1,2)
			p=plot(now,round(sum(g)-1),'r*');
			hold on
			xlabel('Year')
			ylabel('Mass Production')
			grid on
		subplot(3,1,3)
			p=plot(now,1/max(g)^2,'b*');
			hold on
			grid on
			xlabel('Year')
			ylabel('1/C_m_a_x^2')
		m(j,1)=(now);
		m(j,2)=1/max(g)^2;
		j=j+1;
%		pause(1)
		oldnow=now;
	end
	g(1)=g(now);
	g(2:nx)=gold(2:nx)*w0+gold(1:nx-1)*wm;
	gold=g;
end

% Plot the data, calculate the slope, theoretical and numeric diffusivities

subplot(3,1,1)
title('FUDM Diffusion only')
text(1,.085,sprintf('Length of system = %.0f km',dx*nx/1000))
text(1,.07,sprintf('dx = %.0f m',dx))
text(1,.055,sprintf('dt = %.0f s',dt))
text(1,.04,sprintf('u = %.2f m/s',u))
text(1,.025,sprintf('Number of time steps = %.0f',nt))
text(1,.010,sprintf('Number of space steps = %.0f',nx))
text(12,.085,sprintf('c + 2*d = %.2f',c+2*d))

subplot(3,1,3)
[a]=linfit(m(:,1),m(:,2),1);
y=a(2).*m(:,1)+a(1);
plot(m(:,1),y,'b--')
slope=a(2)/tyr;
Ka=(slope*(dx^2))/pi^2;
Kt=u*dx/2*(1-c);
text(1.4,300,sprintf('K_a = %.1f\n',Ka))
text(1.4,225,sprintf('K_t = %.1f\n',Kt))
text(1.4,175,sprintf('Slope = %.3e\n',slope))
Apparent_vs_Theoretical(1,1)=Ka;
Apparent_vs_Theoretical(1,2)=Kt;
pause

%clc
disp(' ')
disp('For this particular system, we observe a K_a of 305 and a K_t of 240.');
disp('Before discussing the results, let''s do all the model runs. Next is');
disp('to double the velocity from u = .01 to u = .02. <return>');
pause

% PART 2

figure(2);clf
tyr=24*3600*365.25; 
dx=50000; 
dt=200000; 
u=.02;
nt=10*tyr/dt;
nx=500;
K=0;
d=K*dt/dx^2;
c=u*dt/dx;

%			FUDM weights
w0=1-u*dt/dx;
wm=u*dt/dx;
%			Initialize arrays
g=zeros(nx,1);
g(nx/2)=1; 
gold=g; 
x=dx*[1:nx]/1.e6;
m=zeros(10,2);
j=1;
oldnow=0;

% Start the 10 year run

for i=1:nt
	now=ceil(i*dt/tyr);
	if now ~= oldnow
		subplot(3,1,1)
			p=plot(x,g,'g');
			set(p,'LineWidth',2)
			grid on
			axis([0 25 0 .1]);
			xlabel('Distance')
			ylabel('Concentration')
		subplot(3,1,2)
			p=plot(now,round(sum(g)-1),'r*');
			hold on
			xlabel('Year')
			ylabel('Mass Production')
			grid on
		subplot(3,1,3)
			p=plot(now,1/max(g)^2,'b*');
			hold on
			grid on
			xlabel('Year')
			ylabel('1/C_m_a_x^2')
		m(j,1)=(now);
		m(j,2)=1/max(g)^2;
		j=j+1;
%		pause(1)
		oldnow=now;
	end
	g(1)=g(now);
	g(2:nx)=gold(2:nx)*w0+gold(1:nx-1)*wm;
	gold=g;
end

% Plot the data, calculate the slope, theoretical and numeric diffusivities

subplot(3,1,1)
title('FUDM Diffusion only: Velocity is Doubled')
text(1,.085,sprintf('Length of system = %.0f km',dx*nx/1000))
text(1,.07,sprintf('dx = %.0f m',dx))
text(1,.055,sprintf('dt = %.0f s',dt))
text(1,.04,sprintf('u = %.2f m/s',u))
text(1,.025,sprintf('Number of time steps = %.0f',nt))
text(1,.010,sprintf('Number of space steps = %.0f',nx))
text(12,.085,sprintf('c + 2*d = %.2f',c+2*d))

subplot(3,1,3)
[a]=linfit(m(:,1),m(:,2),1);
y=a(2).*m(:,1)+a(1);
plot(m(:,1),y,'b--')
slope=a(2)/tyr;
Ka=(slope*(dx^2))/pi^2;
Kt=u*dx/2*(1-c);
text(1.4,725,sprintf('K_a = %.1f\n',Ka))
text(1.4,575,sprintf('K_t = %.1f\n',Kt))
text(1.4,450,sprintf('Slope = %.3e\n',slope))
Apparent_vs_Theoretical(2,1)=Ka;
Apparent_vs_Theoretical(2,2)=Kt;

%clc
disp(' ')
disp('This time K_a=586 and K_t=460. In order to run the next system');
disp('(double the space step), and keep the system length at 25,000 km');
disp('I changed dx to 100000 and had to change nx to 250; u was changed');
disp('back to .01 m/s. <return>');
pause
 
% PART 3

figure(3);clf
tyr=24*3600*365.25; 
dx=100000; 
dt=200000; 
u=.01;
nt=10*tyr/dt;
nx=250;
K=0;
d=K*dt/dx^2;
c=u*dt/dx;

%			FUDM weights
w0=1-u*dt/dx;
wm=u*dt/dx;
%			Initialize arrays
g=zeros(nx,1);
g(nx/2)=1; 
gold=g; 
x=dx*[1:nx]/1.e6;
m=zeros(10,2);
j=1;
oldnow=0;

% Start the 10 year run

for i=1:nt
	now=ceil(i*dt/tyr);
	if now ~= oldnow
		subplot(3,1,1)
			p=plot(x,g,'g');
			set(p,'LineWidth',2)
			grid on
			axis([0 25 0 .1]);
			xlabel('Distance')
			ylabel('Concentration')
		subplot(3,1,2)
			p=plot(now,round(sum(g)-1),'r*');
			hold on
			xlabel('Year')
			ylabel('Mass Production')
			grid on
		subplot(3,1,3)
			p=plot(now,1/max(g)^2,'b*');
			hold on
			grid on
			xlabel('Year')
			ylabel('1/C_m_a_x^2')
		m(j,1)=(now);
		m(j,2)=1/max(g)^2;
		j=j+1;
%		pause(1)
		oldnow=now;
	end
	g(1)=g(now);
	g(2:nx)=gold(2:nx)*w0+gold(1:nx-1)*wm;
	gold=g;
end

% Plot the data, calculate the slope, theoretical and numeric diffusivities

subplot(3,1,1)
title('FUDM Diffusion only: Space step is Doubled')
text(1,.085,sprintf('Length of system = %.0f km',dx*nx/1000))
text(1,.07,sprintf('dx = %.0f m',dx))
text(1,.055,sprintf('dt = %.0f s',dt))
text(1,.04,sprintf('u = %.2f m/s',u))
text(1,.025,sprintf('Number of time steps = %.0f',nt))
text(1,.010,sprintf('Number of space steps = %.0f',nx))
text(12,.085,sprintf('c + 2*d = %.2f',c+2*d))

subplot(3,1,3)
[a]=linfit(m(:,1),m(:,2),1);
y=a(2).*m(:,1)+a(1);
plot(m(:,1),y,'b--')
slope=a(2)/tyr;
Ka=(slope*(dx^2))/pi^2;
Kt=u*dx/2*(1-c);
text(1.4,150,sprintf('K_a = %.1f\n',Ka))
text(1.4,115,sprintf('K_t = %.1f\n',Kt))
text(1.4,90,sprintf('Slope = %.3e\n',slope))
Apparent_vs_Theoretical(3,1)=Ka;
Apparent_vs_Theoretical(3,2)=Kt;
pause

%clc
disp(' ')
Apparent_vs_Theoretical
disp('Above are the apparent vs theoretical diffusivities for all three');
disp('in each case, the apparent diffusivity was larger than the');
disp('theoretical diffusivity. If we examine rows 1 and 2 we can observe');
disp('the effect of doubling the velocity; this causes a near doubling of');
disp('the numeric diffusivity. A closer look at the equation for numeric');
disp('diffusivity (Kn=u*(dx-c*dt)/2 over the velocities encountered (say');
disp('-.02 to +.02) shows that the relationship between Kn and u is nearly');
disp('linear. In actuality, over a larger range of u the relationship is');
disp('parabolic. If we examine rows 1 and 3, we can observe that the');
disp('effect of doubling the space step yields a numeric diffusivity that');
disp('is slightly more than twice the starting value. If we hold u');
disp('constant at .01 and let dx range from 50,000 to 100,000 again');
disp('observe a near linear trend. (this trend appears to stay linear over');
disp('a wide range of dx values. <return>');
pause

%clc
disp(' ')
disp('If we examine the notes, we observe that when solving for the');
disp('theoretical numeric diffusivity we do a Taylor series expansion for');
disp('both the space and time derivatives. In both expansions, terms');
disp('higher than the quadratic are ignored which would tend to');
disp('underestimate the true value of the series. I would also add that');
disp('the series for the space step alternates between positive and');
disp('negative while the time step series is positive only. As a result,');
disp('the effect of the polynomial degree one decides to use will have a');
disp('slightly different effect on time and space. I suspect that if we');
disp('did our expansions to degree 5 instead of degree 2, our theoretical');
disp('diffusivities would approach our apparent diffusivities. The values');
disp('calculated with the slope method are probably more accurate;');
disp('however, depending on the range of diffusivities, the effect of the');
disp('difference between the two methods is small. <return>');
pause

%clc
disp(' ')
disp('End of problem 3');