% For Use in GES 233 / CEE 260C 
% FEM2comp 
% FINITE ELEMENT MODEL
% FOR TWO-COMPONENT LINEAR KINETICS using
% Crank Nicholson Implicit Method for
% Advection Dispersion Reaction Eqns in 1-D
% (Addams, 1998 modified by Gorelick, 1999)
% NOTE: node k=0 cannot exist in MATLAB;
%	so the first node is k=1 (B.C) and the last
% 	is k=N (B.C)

clear variables

% Concentrations for equation 1 are C's
% Concentrations for equation 2 are V's

% PARAMETERS IN FEM2comp MODEL--CHANGE THESE,SAVE and Run

% System Length and Time and space and time increments
L=60;		% Limit of modeled flow [m]
Tmax=5.0;		% Max time
DT=0.01;		% Time step
DX=1;		   % Node interval

% Effective Porosity
Por1=0.3;		% Porosity effective equation1
Por2=1.2;	% Porosity effective equation2

% Specific discharge 
Q1=10;		% GW *flux* (not lin velocity) [m/d] equation 1
Q2=0.0;		% GW *flux* (not lin velocity) [m/d] equation 2

%Dispersion Parameters
alfa1=2.0;  % dispersivity for equation 1
alfa2=0.0;  % dispersivity for equation 2
DH1=Q1*alfa1;	% hydrodynamic dispersion coeff [m2/d] for equation 1
DH2=Q2*alfa2;	% hydrodynamic dispersion coeff [m2/d] for equation 2 

%Kinetic Rate Constants
Kf=800;		% Forward rate constant [1/day]
Kb=400;		% Backward rate constant [1/day]

%Boundary and Initial Conditions
Cb1=0.0;	% Boundary cond at x=L m and initial cond at all x for equation 1
Cf1=100.0;	% Influx Conc. at x=0 t>0  [mg/L] for equation 1
Vb2=0.0;	% Boundary cond at x=L m and initial cond at all x for equation 2
Vf2=200.0;	% Influx Conc. at x=0 t>0  [mg/L] for equation 2  	

% initialize variables
% generally, among a,b,c,d,e,f, etc., these parameters will
% be capitalized if they change in 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)*Vb2;	% Initialize for t=0
Vold(1)=Vf2;

% 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)+(DH1*Up2)+(Q1*Up1);
Pt=(Kf*Tvar)+(DH1*Tvar2)+(Q1*Tvar1);
Ps=(Kf*Us)+(DH1*Us2)+(Q1*Us1);

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


a1=(2*Por1*Up)+(DT*Pp);	% (IM parameter)   a1,b1,c1 correspond to the 1st component 
b1=2*Por1*Tvar+(DT*Pt);  % (IM parameter)   , first eqn
c1=(2*Por1*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 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 Vf2=%.3f [mg/L](2nd component)\n',Vf2);
fprintf(fid,'\tRight Boundary Condition Cb1=%.3f [mg/L]\n',Cb1);
fprintf(fid,'\tRight Boundary Condition Vb2=%.3f [mg/L] (2nd component)\n',Vb2);

fprintf(fid,'\tInitial Condition at x=0  : Cf1(0)=%.3f [mg/L]\n',Cf1);
fprintf(fid,'\tInitial Condition at x>0  : Cb1(i)=%.3f [mg/L]\n',Cb1);
fprintf(fid,'\tInitial Condition at x=0  : Vf2(0)=%.3f [mg/L]\n',Vf2);
fprintf(fid,'\tInitial Condition at x>0  : Vb2(i)=%.3f [mg/L]\n',Vb2);
fprintf(fid,'\tGW flux Q1=%.3f [m/d]\n ',Q1);
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 Coef. DH1=%.3f [m/d]\n ',DH1);
fprintf(fid,'\tHydrodynamic Dispersion Coef. DH2=%.3f [m/d]\n ',DH2);
fprintf(fid,'\tEffective Porosity Component1 Por1=%.3f\n ',Por1);
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+DT    % 

	Cnew(1)=Cf1; Cnew(N)=Cb1; 	% Boundary conditions for component 1
	Vnew(1)=Vf2; Vnew(N)=Vb2; 	% 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*Por1*Tvar-b1)*Cold(2)+(4*Por1*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*Por1*Up-a1)*Cold(k-1)+(4*Por1*Tvar-b1)*Cold(k)+(4*Por1*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*Por1*Up-a1)*Cold(N-2)+(4*Por1*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; See FEM2comp.out
%  Will write the Day 5 C's to screen as a column of values	
if i==1
%      C1=Cnew;
%      V1=Vnew;
%	elseif i==3
%      C3=Cnew;
%      V3=Vnew;
   elseif i==5
      disp('Vs first and Cs second')
      disp('Vs Day 5')
%     V5=Vnew'
      V5=Vnew
      disp('Cs Day 5')
      C5=Cnew
%     C5=Cnew'
%     cvonc=[C5',V5'];
% Plot graph at day 5
clf;
plot(J,C5,J,V5);
hold on;
xlabel('distance [m]')
ylabel('Concentration [mg/L]')
title('Two Component Linear Kinetics');
gtext('Day 5, Component 1');gtext('Day 5, Component 2');
fclose('all');
	end
end



