% axialOCFE.m
% This program runs a code to solve the axial dispersion
% equations using the method of orthogonal collocation on finite elements.

global xoc woc aoc boc qinv np ne nt delx Aocfe Bocfe

% set linear to 1 if problem is linear - this avoids one unnecessary iteration
linear = 2
% set the type of reaction rate
iwhich = 2
%set nx, Peclet (Pe), and Damkohler (Da) before calling
Pe = 100.
Da = 1.
peinv = 1/Pe;

% np is the total number of collocation points in one finite element (z direction)
np = 4
% ne is the number of elements in the z-direction
ne = 4
% This can be changed to produce a variable grid
for k=1:ne
delx(k) = 1/ne;
end
delx(1)

% nt is the total number of points
nt = (np-1)*ne + 1
% number of interior points is np-2 per element
nint=(np-2)*ne;

% create the collocation matrices for the fluid
planar(np,0)

% set up the z values for plotting purposes
num=0;
xstart = 0.;
for k=1:ne
for j=1:np-1
num=num+1;
x(num)=xstart+delx(k)*xoc(j);
end
xstart = xstart+delx(k);
end
x(nt) = 1.

% iterate
tol = 1.e-12;
change = 1.;
iter=0;
% initial guess
for i=1:nt
c(i) = 1-x(i);
end

while change>tol
iter=iter+1;

% set the matrices to zero
Bocfe = zeros(nt);
Aocfe = zeros(np,np,ne);

% set the matrices for the differential equation
for k=1:ne
for j=2:np-1
for i=1:np
Aocfe(j,i,k) = boc(j,i)*peinv/(delx(k)*delx(k)) - aoc(j,i)/delx(k);
end
index = (np-1)*(k-1) + j;
ans = rate(c(index),iwhich);
raterxn = ans(1);
drate = ans(2);
Aocfe(j,j,k) = Aocfe(j,j,k) - Da*drate;
Bocfe(index) = Da*(raterxn - drate*c(index));
end
end

% set the continuity conditions
for k=1:ne-1
for i=1:np
Aocfe(np,i,k) = aoc(np,i)/delx(k);
end
end
for k=2:ne
for i=1:np
Aocfe(1,i,k) = -aoc(1,i)/delx(k);
end
end
for k=2:ne
term = Aocfe(np,np,k-1)+Aocfe(1,1,k);
Aocfe(np,np,k-1) = term;
Aocfe(1,1,k) = term; % this is redundant,since it isn't used
Bocfe((np-1)*(k-1)+1) = 0.;
end
% set the boundary conditions
for i=1:np
Aocfe(1,i,1) = -peinv*aoc(1,i)/delx(1);
end
Aocfe(1,1,1) = Aocfe(1,1,1) + 1;
Bocfe(1) = 1;
for i=1:np
Aocfe(np,i,ne) = -aoc(np,i);
end
Bocfe(nt) = 0;

% Solve one iteration.
OCFElud(np,ne);
OCFEfas(np,ne,nt);
c11 = Bocfe;

% Calculate the criterion to stop the iterations.
sum=0.;
for i=1:nt
sum = sum + abs(c11(i) - c(i));
c(i) = c11(i);
end
if linear==1,break,end
change = sum
end

% Plot the solution.
%plot(x,c,'-r*')
% Find the coefficients
%dd1 = qinv*c';
%xp = 0:0.02:1; %plotting points
%cinterp = planar_interp(nx,dd1,xp,0,1);
%plot(xoc,c,'*b',xp,cinterp,'-b')
%xlabel('z')
%ylabel('c')
c

************************************************************
% rate.m
% calculate the rate and derivative for a given concentration

function ans = rate(c,iwhich)

if iwhich==1
% use for linear reaction
ans(1) = c;
ans(2) = 1.;
elseif iwhich == 2
% use for second order reaction
ans(1) = c*c;
ans(2) = 2.*c;
else
% use for non-isothermal reaction
%ggg = 18.; bbb = 0.3;
ggg = 30; bbb = 0.4;
t = 1 + bbb - bbb*c;
eee = exp(ggg*(1 - 1/t));
ans(1) = c*eee;
ans(2) = eee*(1 - bbb*ggg*c/(t*t));
end

************************************************************
% This subroutine does the LU decomposition for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% and ne is the number of elements.

function OCFElud(np,ne)
global Aocfe

n1 = np - 1;
for l=1:ne
for k=1:n1
k1=k+1;
for i=k1:np
s = Aocfe(i,k,l)/Aocfe(k,k,l);
Aocfe(i,k,l) = s;
for j=k1:np
Aocfe(i,j,l) = Aocfe(i,j,l) - s*Aocfe(k,j,l);
end
end
end
if l<ne
Aocfe(1,1,l+1) = Aocfe(np,np,l);
end
end

************************************************************
% This subroutine does the fore and aft sweep for the
% orthogonal collocation method on finite elements.
% This version is for a single unknown, giving a block
% diagonal matrix with a single entry of overlap.
% This subroutine is in the Appendix
% of "Nonlinear Analysis in Chemical Engineering".
% It assumes diagonal pivoting is OK.
% The matrix Aocfe is stored as Aocfe(np,np,ne), where
% np is the number of collocation points in one element
% ne is the number of elements.
% nt is the total number of unknowns, nt = (np-1)*ne + 1
% To use OCFEfas, you must first run OCFElud.
% The right-hand side is in Bocfe(nt)

function OCFEfas(np,ne,nt)
global Aocfe Bocfe

np1=np;
% forward sweep
for l=1:ne
for i=2:np1
i2=(l-1)*(np-1)+i;
s = 0.;
i1=i-1;
for j=1:i1
j2=i2-i+j;
s=s+Aocfe(i,j,l)*Bocfe(j2);
end
Bocfe(i2)=Bocfe(i2)-s;
end
end
% back substitution
for l1=1:ne
l=ne-l1+1;
if (l==ne)
Bocfe(nt)=Bocfe(nt)/Aocfe(np,np,ne);
end
n1 = np-1;
for k=1:n1
i=n1+1-k;
i2=(l-1)*(np-1)+i;
m=i+1;
n2=n1+1;
s=0.;
for j=m:n2
j2=(l-1)*(np-1)+j;
s=s+Aocfe(i,j,l)*Bocfe(j2);
end
Bocfe(i2)=(Bocfe(i2) - s)/Aocfe(i,i,l);
end
end