% Lee Addams
% HW 6: FINITE ELEMENT MODEL
% FOR TWO-COMPONENT LINEAR KINETICS using
% Crank Nicholsen Implicit Model for
% Advection Dispersion Eqn in 1-D
% GENERALIZED FOR ADSORBTION, DECAY, etc.

% NOTE: the node k=0 cannot exist in MATLAB;
%	so the first node is k=1 (B.C) and the last
% 	is k=51 (B.C)

% identify constants
clear variables
DT=.05;		% Time step
DX=1;		% Node interval
Por=.3;		% Porosity effective
Por2=1.2;	% Porosity effective **2nd component

Cb1=0.0;	% Boundary cond at x=50 m and initial cond at all x *Component 1**
Cf1=100.0;	% Influx Conc. at x=0 t>0  [mg/L]		    *Component 1** 
Cb2=0.0;	% Boundary cond at x=50 m and initial cond at all x *Component 2**
Cf2=200.0;	% Influx Conc. at x=0 t>0  [mg/L]		    *Component 2** 	
DH=2.5;		% hydrodynamic dispersion coeff [m2/d]
Q=1.5;		% GW *flux* (not lin velocity) [m/d]
DH2=0.0;	% hydrodynamic dispersion coeff [m2/d] ** 2nd Component
Q2=0.0;		% GW *flux* (not lin velocity) [m/d]    ** 2nd Component

L=50		% Limit of modeled flow [m]
Tmax=5		% Max time
%% DECAY=0.05	% Decay Constant K (linear kinetics)  **NOT USED!
Kf=0.2		% Forward rate constant [1/day]
Kb=0.1		% Backward rate constant [day^-1]

% initialize the relevant variables
% generally, among a,b,c,d,e,f, etc., these parameters will
% be capitalized if they change with inverting a matrix, and
% will be lowercase if they are constant and only dependent on the 
% variables of the problem given above...

N=L/DX+1		% Number of Nodes
Cold=ones(N,1)*Cb1;	% Initialize for t=0
Cold(1)=Cf1;
Vold=ones(N,1)*Cb2;	% Initialize for t=0
Vold(1)=Cf2;
% Define Integral values specific to the Galerkin-Chapeau Basis Function
% (=1 at node k, goes to zero linearly at node k-1, k+1)

Up=DX/6;Us=DX/6;
Up1=-.5;Us1=.5;
Up2=-1/DX;Us2=-1/DX;
Tvar=2*DX/3;Tvar1=0;Tvar2=2/DX;

Pp=(Kf*Up)+(DH*Up2)+(Q*Up1);
Pt=(Kf*Tvar)+(DH*Tvar2)+(Q*Tvar1);
Ps=(Kf*Us)+(DH*Us2)+(Q*Us1);

Rp=(Kb*Up)+(DH2*Up2)+(Q2*Up1);
Rt=(Kb*Tvar)+(DH2*Tvar2)+(Q2*Tvar1);
Rs=(Kb*Us)+(DH2*Us2)+(Q2*Us1);


a1=(2*Por*Up)+(DT*Pp);	% (IM parameter)   a1,b1,c1 correspond to the 1st component 
b1=2*Por*Tvar+(DT*Pt);  % (IM parameter)   , first eqn
c1=(2*Por*Us)+(DT*Ps);	% (IM parameter)

a2=-1*DT*Kb*Up;		% (IM parameter)   a2,b2,c2 correspond to the 2nd component	
b2=-1*DT*Kb*Tvar;	% (IM parameter)   , 1st eqn	
c2=-1*DT*Kb*Us;		% (IM parameter)

a3=-1*DT*Kf*Up; 	% (IM parameter)   a3,b3,c3 correspond to the 2nd component	
b3=-1*DT*Kf*Tvar;	% (IM parameter)   , 2nd eqn	
c3=-1*DT*Kf*Us;		% (IM parameter)

a4=(DT*Rp)+(2*Por2*Up);	 % (IM parameter)   a4,b4,c4 correspond to the 1st component 
b4=(DT*Rt)+2*Por2*Tvar;  % (IM parameter)   , 2nd eqn
c4=(DT*Rs)+(2*Por2*Us);	 % (IM parameter)


% Now each eqn has 6 unknowns, the three nodes of "C"
% and the three nodes of "V"... the a1..c2 are the coefficients
% This is a bi-tridiagonal matrix... two eqns each with 6 unknowns.
% 2(n-1) algebraic eqns and 2(n-1) unknowns...
% Use von Rosenbergs conventions for naming...


J=[0:DX:L]   	% Vector describing position of N=51 nodes

% Time step loop--marching through time at DT time steps

NPRINT=0;

% 	Open output file and write the initial parameters
%	Boundary, Initial Conditions; Time step; Node Interval;
%	Dispersion Coefficient; GW flux
fid= fopen('FEM2comp.out','w+')
fprintf(fid,'\nFINITE ELEMENT APPROXIMATION\n');
fprintf(fid,'Two Component Linear Kinetics\n');
fprintf(fid,'\tLeft Boundary Condition Cf1=%.3f [mg/L]\n',Cf1);
fprintf(fid,'\tLeft Boundary Condition Cf2=%.3f [mg/L](2nd component)\n',Cf2);
fprintf(fid,'\tRight Boundary Condition Cb1=%.3f [mg/L]\n',Cb1);
fprintf(fid,'\tRight Boundary Condition Cb2=%.3f [mg/L] (2nd component)\n',Cb2);

fprintf(fid,'\tInitial Condition at x=0  : C(0)=%.3f [mg/L]\n',Cf1);
fprintf(fid,'\tInitial Condition at x>0  : C(i)=%.3f [mg/L]\n',Cb1);
fprintf(fid,'\tGW flux Q=%.3f [m/d]\n ',Q);
fprintf(fid,'\tGW flux Q2=%.3f [m/d]\n ',Q2);
fprintf(fid,'\tForward Reaction Rate Constant Kf=%.3f [1/day]\n ',Kf);
fprintf(fid,'\tReverse Reaction Rate Constant Kb=%.3f [1/day]\n ',Kb);
fprintf(fid,'\tHydrodynamic Dispersion DH=%.3f [m/d]\n ',DH);
fprintf(fid,'\tHydrodynamic Dispersion DH2=%.3f [m/d]\n ',DH2);
fprintf(fid,'\tEffective Porosity Component1 Por=%.3f\n ',Por);
fprintf(fid,'\tEffective Porosity Component2 Por2=%.3f\n ',Por2);
fprintf(fid,'\tTime Step DT=%.3f [days]\n',DT);
fprintf(fid,'\tNode Interval DX=%.3f [m]\n',DX);
fprintf(fid,'\tTotal Nodes N=%1g\n',N);
fprintf(fid,'\n----------------------------\n Intermediate Parameters:')
fprintf(fid,'\nPp=\t%.4f\tRp=%.4f',Pp,Rp);
fprintf(fid,'\nPt=\t%.4f\tRt=%.4f',Pt,Rt);
fprintf(fid,'\nPs=\t%.4f\tRs=%.4f',Ps,Rs);
fprintf(fid,'\na1=\t%.4f\ta2=\t%.4f\ta3=\t%.4f\ta4=\t%.4f',a1,a2,a3,a4);
fprintf(fid,'\nb1=\t%.4f\tb2=\t%.4f\tb3=\t%.4f\tb4=\t%.4f',b1,b2,b3,b4);
fprintf(fid,'\nc1=\t%.4f\tc2=\t%.4f\tc3=\t%.4f\tc4=\t%.4f',c1,c2,c3,c4);

% TIME LOOP by DT up to Tmax;

for i=DT:DT:Tmax    % 

	Cnew(1)=Cf1; Cnew(N)=Cb1; 	% Boundary conditions for component 1
	Vnew(1)=Cf2; Vnew(N)=Cb2; 	% Boundary conditions for component 2
	NPRINT=NPRINT+1;
	
	%  Loop to calculate D,E,F (forward)
	% "D" is R.H.S that before was Dcurly
	%  These variables are all simplifications of the 
	%  matrix we will wish to invert for each time step
	
	%  D1(2) does not fit the "general" case and must be calculated
	%  separately
	
	D1=(4*Por*Tvar-b1)*Cold(2)+(4*Por*Us-c1)*Cold(3)-b2*Vold(2)-c2*Vold(3)-2*DT*Pp*Cnew(1)+2*DT*Kb*Up*Vnew(1);
	D2=(4*Por2*Tvar-b4)*Vold(2)+(4*Por2*Us-c4)*Vold(3)-b3*Cold(2)-c3*Cold(3)-2*DT*Rp*Vnew(1)+2*DT*Kf*Up*Cnew(1);
			% The Vnew and Cnew terms contain the right B.C.
			
		
		bt1=b1;bt2=b2;bt3=b3;bt4=b4;dt1=D1;dt2=D2;DEN=bt1*bt4-bt2*bt3;	
		E1(2)=[bt4*c1-bt2*c3]/DEN;
		E2(2)=[bt4*c2-bt2*c4]/DEN;
		E3(2)=[bt1*c3-bt3*c1]/DEN;	
		E4(2)=[bt1*c4-bt3*c2]/DEN;
		F1(2)=[bt4*dt1-bt2*dt2]/DEN;F2(2)=[bt1*dt2-bt3*dt1]/DEN;
		
	for k=3:1:N-2	% Get the Dcurly from node 3 to node N-1
			% Get the E,F from node 3 to node N-1
		D1=(4*Por*Up-a1)*Cold(k-1)+(4*Por*Tvar-b1)*Cold(k)+(4*Por*Us-c1)*Cold(k+1)-a2*Vold(k-1)-b2*Vold(k)-c2*Vold(k+1);
		D2=(4*Por2*Up-a4)*Vold(k-1)+(4*Por2*Tvar-b4)*Vold(k)+(4*Por2*Us-c4)*Vold(k+1)-a3*Cold(k-1)-b3*Cold(k)-c3*Cold(k+1);
			bt1=b1-a1*E1(k-1)-a2*E3(k-1);
			bt2=b2-a1*E2(k-1)-a2*E4(k-1);	
			bt3=b3-a3*E1(k-1)-a4*E3(k-1);
			bt4=b4-a3*E2(k-1)-a4*E4(k-1);
			dt1=D1-a1*F1(k-1)-a2*F2(k-1);
			dt2=D2-a3*F1(k-1)-a4*F2(k-1);
			DEN=bt1*bt4-bt2*bt3;
		E1(k)=[bt4*c1-bt2*c3]/DEN;
		E2(k)=[bt4*c2-bt2*c4]/DEN;
		E3(k)=[bt1*c3-bt3*c1]/DEN;	
		E4(k)=[bt1*c4-bt3*c2]/DEN;
		F1(k)=[bt4*dt1-bt2*dt2]/DEN;F2(k)=[bt1*dt2-bt3*dt1]/DEN;
		
	end
	
	%  Loop to calculate all the Cj's, Vj's (reverse)
	%  "Inversion" algorithm for the tridiagonal matrix
	D1=(4*Por*Up-a1)*Cold(N-2)+(4*Por*Tvar-b1)*Cold(N-1)-a2*Vold(N-2)-b2*Vold(N-1)-2*DT*Ps*Cnew(N)+2*DT*Kb*Us*Vnew(N);
	D2=(4*Por2*Up-a4)*Vold(N-2)+(4*Por2*Tvar-b4)*Vold(N-1)-a3*Cold(N-2)-b3*Cold(N-1)-2*DT*Rs*Vnew(N)+2*DT*Kf*Us*Cnew(N);
	E1(N-1)=0;E2(N-1)=0;E3(N-1)=0;E4(N-1)=0;
			bt1=b1-a1*E1(N-2)-a2*E3(N-2);
			bt2=b2-a1*E2(N-2)-a2*E4(N-2);	
			bt3=b3-a3*E1(N-2)-a4*E3(N-2);
			bt4=b4-a3*E2(N-2)-a4*E4(N-2);
			dt1=D1-a1*F1(N-2)-a2*F2(N-2);
			dt2=D2-a3*F1(N-2)-a4*F2(N-2);
			DEN=bt1*bt4-bt2*bt3;
	
	F1(N-1)=[bt4*dt1-bt2*dt2]/DEN;
	F2(N-1)=[bt1*dt2-bt3*dt1]/DEN;

		Cnew(N-1)=F1(N-1);	%  Get the N-1 node
		Vnew(N-1)=F2(N-1);	
	for j=N-1:-1:2	
		Cnew(j)=F1(j)-E1(j)*Cnew(j+1)-E2(j)*Vnew(j+1);	% Start going backward...
		Vnew(j)=F2(j)-E3(j)*Cnew(j+1)-E4(j)*Vnew(j+1);
	end

	Cold=Cnew;Vold=Vnew; 

	if NPRINT==20;
		fprintf(fid,'\nTime step=\t%.3f days\n\n',i)
		fprintf(fid,'<< C >>\n')
		fprintf(fid,'%.6f\t',Cnew)
		fprintf(fid,'\n<< V >>\n')
		fprintf(fid,'%.6f\t',Vnew)
		fprintf(fid,'\n\n')
		NPRINT=0;
	end
	
%	Collect Snapshot of days 1,3,5;
	
	if i==1
		C1=Cnew;
	elseif i==3
		C3=Cnew;
	elseif i==5
		C5=Cnew;
		V5=Vnew;
	end
end
i

% 	Plot Analytical results for comparison
figure(1)
clf;
%	Crank-Nicholsen numerical Soln indicated with "+"

plot(J,C5,J,V5);
hold on;
xlabel('distance [m]')
ylabel('Concentration [mg/L]')
title('HW5- Two Component Linear Kinetics');
gtext('Day 5,Component 1');gtext('Day 5, Component 2');
fclose('all');
end

 	



