% FEm.PSE template, 1D unsteady conduction.

global X1;        % array for node coordinates

% uniform discretization, M=4
nnodes = 9;
X1 = linspace(1,2,nnodes);

% physical constants (input, dimensional)
k = 10;	                    % conductivity, Btu/(hr*ft*degF)
rho = 1;                    % density, ft^3/lbm
Cp = 1;                     % specific heat, Btu/(lbm*degF)
h = 20;			    % convection coefficient, Btu/(hr*ft^2*degF)
Tr = 1500;		    % convection reference temperature, degF
Tb = 306.8528;		    % diri temperature, degF

% initial condition (input, dimensional)
T_initial = Tb * ones(nnodes,1);

% physical constants (calculated)
Tmax = Tr;                  % max reference dimensional temperature
Tmin = Tb;                  % min reference dimensional temperature
pot_Tr = (Tr-Tmin)/(Tmax-Tmin); % convection reference temperature in
                                % potential form
L = X1(nnodes) - X1(1);     % length of the domain
Nu = h*L/k;                 % Nusselt number
alpha = rho*Cp/k;           % thermal diffusivity
t_ref = L^2 * alpha;        % non-dimensional time
r_inner = X1(1);            % inner radius
r_inner_nd = r_inner/L;     % non-dimesnional inner radius

% time loop parameters
delta_t = .0001;            % dimensional time step (input)
delta_t_nd = delta_t/t_ref; % non-dimenstional time step
theta = .5;
last_timestep = 10;

% transform initial condition to potential form
pot_T_initial = (T_initial - Tmin) ./ (Tmax-Tmin);

% load the FE matrix library
load femlib;

% The required matrices are constant, and need be computed only
% once, outside the time loop
MASS = asjac1D([],[],[X1],1,'A3000L',[]);
DIFF = asjac1D([],[],[X1],-1,'A3011L',[]);
% make a matrix to point to the convection node
H_pointer = zeros(nnodes,nnodes);
H_pointer(1,1) = 1;
HBC = r_inner_nd * Nu * H_pointer;

% convective contribution to b at node 1
b = zeros(nnodes,1);
b(1) = r_inner_nd * Nu * pot_Tr;

% make a matrix to store the solutions at each time step
Q_history = zeros(nnodes,last_timestep+1);

% initialize Qn with the non-dimensional temperature profile
Qn = pot_T_initial;

Q_history(:,1) = Qn;     % store the initial condition

% time loop
for step = 1:last_timestep
  SYS = MASS + theta*delta_t_nd*(DIFF+HBC);
  RESn = (DIFF + HBC)*Qn - b;
  RHS = -delta_t_nd * RESn;
  
  % enforce the dirichlet bc on the last node by making the coefficient
  % in the nnodes,nnodes position very large
  SYS(nnodes,nnodes) = SYS(nnodes,nnodes) * 1E15;
    
  % solve for deltaQ
  deltaQ = SYS\RHS;
  
  % update the solution, but don't use deltaQ on the diri node
  deltaQ(nnodes) = 0;
  Qn = Qn + deltaQ;
  
  % store the solution for this time step
  Q_history(:,step+1) = Qn;
end

% convert the non-dimensional Q_history to dimensional form
T_history = Q_history*(Tmax-Tmin) + Tmin;

% To make a movie of the time history of this solution, run the
% axi_uns_movie script or append its contents to the bottom of this file