% SS_CSTR_2_rxn_scan_Q.m
% This MATLAB m-file is modified from SS_CSTR_2_rxn.m to
% "scan" the concentrations of each species at steady
% state as a function of the volumetric flow rate through
% the reactor.
%
% K. Beers
% MIT ChE
% 6/26/2002
%
% ----------------------------------------------------------------------
% SS_CSTR_2_rxn.m
%
% This MATLAB m-file demonstrates the use of the optimization
% toolkit function fsolve to solve sets of nonlinear algebraic
% equations. This program calculates the steady-state
% concentrations of four species in an isothermal CSTR with
% two reactions :
% A + B ---> C
% B + C ---> D
% The MATLAB optimization toolkit function fsolve uses a trust-region
% algorithm to minimize the 2-norm of the function vector. Since this
% scalar function has global minima at the solutions with f(x) = 0, this
% method will find a solution unless the Jacobian becomes singular, in
% which case local minima may exist that are not solutions.
%
% The user must provide at least an initial guess to the solution, and
% the name of a function that takes as input the value of x and returns
% the vector f(x). An option exists for providing a routine to calculate
% the Jacobian analytically, but this example does not do this.
%
% For large, sparse systems, the approximation of the Jacobian can be
% made much more efficient if the fsolve routine is told which elements
% of the Jacobian may be non-zero. This can be done by providing a
% sparse, "structure" matrix that contains non-zero values (for example 1's)
% only at those positions in which the Jacobian is allowed to have non-zero
% elements. Since this example problem involves only four unknowns, the code
% that would provide such a matrix (here for a tri-diagonal problem) is
% commented out.
%
% K. Beers
% MIT ChE
% 6/26/2002


function iflag_main = SS_CSTR_2_rxn();

iflag_main = 0;

% We now set the values of the system parameters, which we
% store in the data stucture Param.

% the volumetric flow rate through the reactor (base case)
Param.Q = 1;
% the total reactor volume
Param.V = 10;
% the rate constant of the first reaction
Param.k1 = 1;
% the rate constant of the second reaction
Param.k2 = 1;
% the inlet concentration of A
Param.c_A_in = 1;
% the inlet concentration of B
Param.c_B_in = 2;
% the inlet concentration of C
Param.c_C_in = 0;
% the inlet concentration of D
Param.c_D_in = 0;


% Next, the options for fsolve are set using the command
% optimset. Here we only decrease the tolerance from the
% default value and use the small-scale algorithm. The
% default is to use a "large-scale" algorithm that is less
% memory intensive for very large systems but is less robust.
% Since this example system has only four unknowns, we turn the
% 'LargeScale' option off.
options = optimset('TolFun',1e-8,'LargeScale','off');

% With these options selected, we have fsolve estimate the
% Jacobian numerically. For large problems with sparse
% Jacobians, we can improve the efficiency by supplying
% a sparse matrix SJac that has non-zero values only at
% those locations where we expect the Jacobian is also
% non-zero. The following comment lines demonstrate the
% code for doing this for a tri-diagonal Jacobian.
%
% SJac = spalloc(Nvar,Nvar,3*Nvar);
% for i=1:Nvar
% SJac(i,i-1) = 1;
% SJac(i,i) = 1;
% SJac(i,i+1) = 1;
% end
% options = optimset('JacobPattern',SJac);
%
% Alternatively, if we write our function routine to return the
% Jacobian as well as the function value (as we do below), then
% the following flag tells fsolve to use this Jacobian instead
% of using numerical approximation. This reduces the number
% of function evaluations during the calculation, but requires
% more programming effort.
% options = optimset(options,'Jacobian','on');
%
% Another option is to provide a function that returns the
% product of multiplying an input vector by the Jacobian.
% This is controlled by the JacobMult option, and helps avoid
% memory problems associated with storing the Jacobian matrix
% in memory. This is primarily used for large problems.
%
% We now turn off the display of the solver.
options = optimset(options,'Display','off');



% We now allocate memory to store the trial values
% of the volumetric flow rate Q and the output
% values of the concentrations.
Q_vect = logspace(-3,3,100)';
results = zeros(length(Q_vect),4);

% Now, we run a FOR loop over each calcualtion. We
% provide as the initial guess the inlet concentration.
% Since we expect this to be near the solution for
% high values of the flow rate (low residence times),
% we start the calculations at the upper end of Q and
% move downwards.

% initial guess
x_guess = [Param.c_A_in; Param.c_B_in; ...
Param.c_C_in; Param.c_D_in];

for iQ = length(Q_vect):-1:1

% set value of volumetric flow rate
Param.Q = Q_vect(iQ);

% We now call fsolve to perform the calculation. Note that we
% also pass along the values of the system parameter structure
% Param to the routine through the argument list.

[x,fval,exitflag,output,Jacob] = ...
fsolve(@SS_CSTR_2_rxn_calc_f,x_guess,options,Param);

if(exitflag ~= 1)
disp(Param.Q);
disp(['fsolve finished with exitflag = ', int2str(exitflag)]);
end

% store results in array
results(iQ,:) = x';

% update x_guess to results of last simulation
x_guess = x;

end

% Now, plot the results in a master figure.
figure;
% [A] vs. Q
subplot(2,2,1);
semilogx(Q_vect,results(:,1));
xlabel('Q');
ylabel('[A]');
axis tight;
% [B] vs. Q
subplot(2,2,2);
semilogx(Q_vect,results(:,2));
xlabel('Q');
ylabel('[B]');
axis tight;
% [C] vs. Q
subplot(2,2,3);
semilogx(Q_vect,results(:,3));
xlabel('Q');
ylabel('[C]');
axis tight;
% [D] vs. Q
subplot(2,2,4);
semilogx(Q_vect,results(:,4));
xlabel('Q');
ylabel('[D]');
axis tight;
% add master title using gtext
gtext('CSTR, A + B --> C, B + C --> D');

iflag_main = 1;

return;


% ======================================================
% This function evaluates the values of the equations for the
% current estimate of the solution. If we also calculated
% the Jacobian analytically, we would return Jac - the Jacobian
% matrix - as a second return argument.

function [f] = SS_CSTR_2_rxn_calc_f(x,Param);

% First, we unstack the vector x into more
% meaningful names.
c_A = x(1);
c_B = x(2);
c_C = x(3);
c_D = x(4);

% We now evalue the values of the function vector f(x)
% for each of the species' mass balances.

f = zeros(4,1);

% mass balance on A
f(1) = (Param.Q/Param.V) * (Param.c_A_in - c_A) ...
- Param.k1*c_A*c_B;

% mass balance on B
f(2) = (Param.Q/Param.V) * (Param.c_B_in - c_B) ...
- Param.k1*c_A*c_B - Param.k2*c_B*c_C;

% mass balance on C
f(3) = (Param.Q/Param.V) * (Param.c_C_in - c_C) ...
+ Param.k1*c_A*c_B - Param.k2*c_B*c_C;

% mass balance on D
f(4) = (Param.Q/Param.V) * (Param.c_D_in - c_D) ...
+ Param.k2*c_B*c_C;


return;