% filename run_radialFD.m
% This is the rhs of the PDEs for a reactor with radial dispersion.
% They are solved with the finite difference method in cylindrical geometry.

global alphap betap gammap Pem Peh Damk Biot tsurr ngrid delr delr2 rad

ngrid=9
% set parameters
alphap = 1.
betap = 0.6666667
gammap = 20.
Pem = 1
Peh = 1
Damk = 0.3
Biot = 1
tsurr = 0.92

% set the time interval for integration
zspan = [0 1]
options=odeset('RelTol',1.e-12,'AbsTol',1.e-6)

for i=1:ngrid
yi(1,i) = 1.;
yi(2,i) = 1.;
rad(i) = (i-1)/(ngrid-1);
end
delr = 1/(ngrid-1);
delr2=delr*delr;

% move variables to one long list
num = 0;
for nvar = 1:2
for k=1:ngrid
num=num+1;
y0(num) = yi(nvar,k);
end
end
y0

[z,solution] = ode15s('radialrhsFD',zspan, y0,options);

% move the variables into variable names
num = 0;
for nvar = 1:2
for k=1:ngrid
num=num+1;
ytemp(:,nvar,k) = solution(:,num);
end
end
for k=1:ngrid
conc(:,k) = ytemp(:,1,k);
temp(:,k) = ytemp(:,2,k);
end
nn = size(conc);
for kkk = 1:nn(1)
for k=1:ngrid
c(k) = conc(kkk,k);
t(k) = temp(kkk,k);
end

% calculate the average conc and temp
% must use an odd number of points
for i=2:2:ngrid-1
w(i) = 4.;
w(i+1) = 2.;
end
w(1) = 1;
w(ngrid) = 1;
sum1 = 0.;
sum2 = 0.;
for i=1:ngrid
sum1 = sum1 + w(i)*c(i)*rad(i);
sum2 = sum2 + w(i)*t(i)*rad(i);
end
cavg(kkk) = 2.*(delr/3.)*sum1;
tavg(kkk) = 2.*(delr/3.)*sum2;
end
plot(z,cavg,'-r',z,tavg,'-b')
xlabel('axial location')
legend('average concentration','average temperature')
title('Problem One')


************************************************
% filename radialrhsFD.m
% This is the rhs of the PDEs for a reactor with radial dispersion.
% They are solved with the finite difference method in cylindrical geometry.

function ydot=radialrhs(z,y)
global alphap betap gammap Pem Peh Damk Biot tsurr ngrid delr delr2 rad

% move y into variable names
num = 0;
for nvar = 1:2
for k=1:ngrid
num=num+1;
ytemp(nvar,k) = y(num);
end
end
c = ytemp(1,:);
t = ytemp(2,:);

% now compute the right-hand sides, for the time derivatives

% j = 1
ydot(1)=4.*(c(2)-c(1))/delr2;
ydot(1+ngrid)=4.*(t(2)-t(1))/delr2;
% j = ngrid
ydot(ngrid)=2.*(-c(ngrid)+c(ngrid-1))/delr2;
ydot(2*ngrid)=2.*(-t(ngrid)+t(ngrid-1)-Biot*delr*(t(ngrid)-tsurr))/delr2;
ydot(2*ngrid)=ydot(2*ngrid)-Biot*(t(ngrid)-tsurr);
% j = 2 to ngrid-1
for j=2:ngrid-1
ydot(j)=(c(j+1)-2*c(j)+c(j-1))/delr2 + (c(j+1)-c(j-1))/(2.*rad(j)*delr);
ydot(j+ngrid)=(t(j+1)-2*t(j)+t(j-1))/delr2 + (t(j+1)-t(j-1))/(2.*rad(j)*delr);
end

for j=1:ngrid
ydot(j) = (alphap/Pem)*ydot(j);
ydot(j+ngrid) = (alphap/Peh)*ydot(j+ngrid);
end

% next do the reaction terms
for j=1:ngrid
rate1 = Damk*c(j)*exp(gammap*(1.-1./t(j)));
ydot(j) = ydot(j) - rate1;
ydot(j+ngrid) = ydot(j+ngrid) + betap*rate1;
end

ydot = ydot';
************************************************

% plot_radialFD.m
% plots results of the reactor calculation
% plot versus radius, for a fixed ngrid, at several z.

global alphap betap gammap Pem Peh Damk Biot tsurr ngrid delr delr2 rad

ngrid=9
% set parameters
alphap = 1.
betap = 0.6666667
gammap = 20.
Pem = 1
Peh = 1
Damk = 0.3
Biot = 1
tsurr = 0.92

% run for six z
dataz = [0.2 0.4 0.5 0.6 0.7 1.0]
% set the time interval for integration
zspan = [0 dataz(1)]

for i=1:ngrid
yi(1,i) = 1.;
yi(2,i) = 1.;
rad(i) = (i-1)/(ngrid-1);
end
delr = 1/(ngrid-1);
delr2=delr*delr;

% move variables to one long list
num = 0;
for nvar = 1:2
for k=1:ngrid
num=num+1;
y0(num) = yi(nvar,k);
end
end

for kkk=1:6
[z,solution] = ode15s('radialrhsFD',zspan, y0);

% move the variables into variable names
num = 0;
for nvar = 1:2
for k=1:ngrid
num=num+1;
ytemp(:,nvar,k) = solution(:,num);
end
end
nn = size(ytemp)
for i=1:ngrid
c(kkk,i) = ytemp(nn(1),1,i);
t(kkk,i) = ytemp(nn(1),2,i);
end
% set parameters for next interval
if kkk==6, break, end
zspan = [dataz(kkk) dataz(kkk+1)]
for i=1:2*ngrid
y0(i) = solution(nn(1),i);
end
clear ytemp solution conc temp
y0
pause
end

% plot the solutions
plot(rad,c(1,:),rad,c(2,:),rad,c(3,:),rad,c(4,:),rad,c(5,:),rad,c(6,:))
%plot(rad,t(1,:),rad,t(2,:),rad,t(3,:),rad,t(4,:),rad,t(5,:),rad,t(6,:))
xlabel('radial location')
legend('z = 0.2','z = 0.4','z = 0.5','z = 0.6','z = 0.7','z = 1.0')
title('Problem One')