% diffusionOCFE.m
% This program runs a code to solve the diffusion problem
% 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 diffusion coefficient
iwhich = 3

% 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 = 8
% 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) = x(i);
end

while change>tol
iter=iter+1;

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

% set the matrices for the differential equation
for k=1:ne
for j=2:np-1
index = (np-1)*(k-1) + j;
ans = diffusivity(c(index),iwhich);
value = ans(1); derivative = ans(2);
index0 = (np-1)*(k-1);
sum = 0.;
for lll = 1:np
sum = sum + aoc(j,lll)*c(index0+lll);
end
for i=1:np
Aocfe(j,i,k) = value*boc(j,i) + derivative*2*sum*aoc(j,i);
end
Bocfe(index) = derivative*sum*sum;
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 (using the fact that Aocfe starts out as zero each iteration
Aocfe(1,1,1) = 1;
Bocfe(1) = 0;
Aocfe(np,np,ne) = 1;
Bocfe(nt) = 1;

% 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*')
xlabel('z')
ylabel('c')
c

flux1 = 0;
flux2 = 0;
index0 = (np-1)*(ne-1);
for i=1:np
flux1 = flux1+aoc(1,i)*c(i);
flux2 = flux2+aoc(np,i)*c(index0+i);
end
flux1=flux1/delx(1);
flux2=flux2/delx(ne);
ans = diffusivity(c(1),iwhich);
flux1=-ans(1)*flux1
ans = diffusivity(c(nt),iwhich);
flux2=-ans(1)*flux2

************************************************************************
% calculate the diffusivity as a function of concentration
% and its derivative with respect to concentration

function ans = diffusivity(c,iwhich)
if iwhich==1
% use for constant case
ans(1) = 1;
ans(2) = 0.;

elseif iwhich == 2
% use for linear case (problem is nonlinear)
ans(1) = 1+c; %protect against diff = 0
ans(2) = 1.;
else
% use for exponential case
ddd = 0.135*exp(2*c);
ans(1) = 1 + ddd;
ans(2) = 2*ddd;
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

************************************************************************
function planar(nx,kon)
% this subroutine computes the matrices for collocation without
% symmetry, using collocation points in table 4-3.
% input variable
% nx = number of collocation points, including two end points
%
% output variables, must be declared global
% aoc = matrix for first derivative, eq. 4-103
% boc = matrix for second derivative, eq. 4-103
% qinv = matrix for q inverse, eq. 4-101
% xoc = vector of collocations points, from table 4-3
% woc = vector of weights, eq. 4-106
%
global xoc woc aoc boc qinv
if(nx<3)
nx=3;
fprintf('\n ***************************************** \n')
fprintf('The number of collocation points must be between \n')
fprintf('3 and 7 (inclusive). nx has been changed to %g\n\n',nx)
end
if(nx>7)
nx=7;
fprintf('\n ***************************************** \n')
fprintf('The number of collocation points must be between \n')
fprintf('3 and 7 (inclusive). nx has been changed to %g\n\n',nx)
end
% now 3<=nx<=7
if kon==0
fprintf('The total number of collocation points is %g\n\n',nx)
end

xoc(1)=0.;
if (nx==3)
xoc(2)=0.5;
elseif (nx==4)
xoc(2)=0.21132486540519;
xoc(3)=1-xoc(2);
elseif (nx==5)
xoc(2)=0.11270166537926;
xoc(3)=0.5;
xoc(4)=1-xoc(2);
elseif (nx==6)
xoc(2)=0.06943184420297;
xoc(3)=0.33000947820757;
xoc(4)=1-xoc(3);
xoc(5)=1-xoc(2);
else
xoc(2)=0.04691007703067;
xoc(3)=0.23076534494716;
xoc(4)=0.5;
xoc(5)=1-xoc(3);
xoc(6)=1-xoc(2);
end
xoc(nx)=1.;
format long
for i=1:nx,
r(i,i) = 0.0;
aoc(i,i) = 0.0;
s(i) = 1.0;
boc(i,i) = 0.0 ;
for j=1:nx,
if (i~=j)
r(i,j) = 1.0/(xoc(i)-xoc(j)) ;
s(i) = s(i)*r(i,j);
end
end
for j=1:nx,
jx = nx-j+1;
if (jx<j), break, end
if (jx==j)
aoc(i,i) = aoc(i,i)+r(i,j);
end
if (jx>j)
aoc(i,i) = aoc(i,i)+r(i,j)+r(i,jx);
end
end
end
for i=1:nx
for j=1:nx
if (i~=j)
aoc(i,j) = s(j)*r(i,j)/s(i);
boc(i,j) = 2.0*aoc(i,j)*(aoc(i,i)-r(i,j));
boc(i,i) = boc(i,i)+r(i,j)*(aoc(i,i)-r(i,j));
end
end
end
for i=1:nx
qinv(1,i) = s(i);
k = 1;
woc(i) = 0.0 ;
for j=1:nx
if (j~=i)
l = k;
k = k+1;
qinv(k,i) = qinv(l,i);
while (l>1)
m = l-1;
qinv(l,i) = qinv(m,i)-xoc(j)*qinv(l,i);
l = m;
end
qinv(1,i) = -xoc(j)*qinv(1,i);
end
end
for j=1:nx
woc(i) = woc(i)+qinv(j,i)/j;
end
end

if kon == 0
fprintf('A matrix (aoc) \n')
disp(aoc)
fprintf('B matrix (boc) \n')
disp(boc)
fprintf('Q matrix (qinv) \n')
disp(qinv)
fprintf('Collocation points (xoc) \n')
disp(xoc)
fprintf('W matrix (woc) \n')
disp(woc)
end