% FFEA template file - Peclet problem

clear all;
clf reset;
format long e;      % double precision exponent output
profile on;

global X1;          % array for node coordinates
load femlib;        % load FEA libarary

% constant data
              % non-dimensional peclet number
              % inverse of Peclet number

% spatial data
				% left boundary
				% right boundary
    			% number of nodes
                % uniform grid
%X1 = ratlist(XL,XR,nnodes-1,0.75); % non-uniform grid

% temporal data
					% implicitness factor
					% max number of Newton iterations
    			    % Newton convergence criteria
					% initial time (s)
					% final time (s)
					% time step size
            	    % number of time steps

% boundary conditions
			% left dirichlet BC - non-dimensional
			% right dirichlet BC - non-dimensional

%initial condition
Q_ini = X1';			% linear distribution
Q_ini(1,1) = QL;		% make sure BC's are exactly enforced
Q_ini(end,1) = QR;

Q_n = Q_ini;                    % Initialize Qn with IC data
Q_np1 = Q_ini;					% Initialize Qnp1 with guess of IC data
DELTA_Q = Q_np1 - Q_n;			% Form temporal change in Q
delta_Q = pi*ones(nnodes,1);	% Initialize delta Q
num_its = 1;                    % Initialize iteration counter

Q_big(:,1) = Q_ini;             % Initialize matrix for solution storage
X_big(:,1) = X1';               % Initialize matrix for x storage
T_big(:,1) = zeros(nnodes,1);   % Initialize matrix for time storage

%**************************************************************************
%OUTER TIME STEPPING LOOP
for t_loop=1:nsteps
    %INNER ITERATION LOOP   
    while ((max(abs(delta_Q))>newt_crit)&(num_its <= max_its))
        
        % Form temporal MASS contribution to total residual FQp using DELTA_Q
        
        
        % Form spatial residual contributions to total residual FQp
        % {RES}_np1 term
        
        % {RES)_n term
        
        
        % Form total residual term FQp
        
        
        %*******************************************************************
        % Form temporal MASS contribution to jacobian JACp
        
        % Form spatial residual contributions to jacobian JACp
        
        % Form total jacobiaion term JACp
        
        
        
        %*******************************************************************
        % Modify JACp for dirichlet data
        
        
        % Modify FQp for Dirichlet data
        
        
        % Solve for incremental change in Q_np1
        %delta_Q = JACp \ -FQp;
        delta_Q = PINV(JACp)*(-FQp);
        
        % Update Q_np1 with incremental change
        
        
        % Recover new value of DELTA_Q
        
        
        % echo stuff
        status = [t_loop max(abs(delta_Q)) max(abs(DELTA_Q))]
        
        % increment numits counter
        num_its = num_its + 1;
    end
    
    % Moving to new time station, the current Q now becomes the previous Q
    Q_big(:,t_loop+1) = Q_np1;
    X_big(:,t_loop+1) = X1';
    T_big(:,t_loop+1) = t_loop*DELTA_T*ones(nnodes,1);
    Q_n = Q_np1;
    DELTA_Q = Q_np1 - Q_n;			% Form temporal change in Q
    delta_Q = pi*ones(nnodes,1);	% Initialize delta Q
    num_its = 1;
    
end

% calculate exact solution
x_exact = XL:0.01:XR;
q_exact = (1-exp(Pe.*x_exact))./(1-exp(Pe));


figure(1)
plot(x_exact,q_exact);
hold on
plot(X1,Q_big(:,end),'ro-')
legend('Exact','k=1',2)
hold off
title(['Steady Peclet Problem - Pe=' num2str(Pe)])
xlabel('x')
ylabel('q')

% pause
% figure(2)
% surf(T_big,X_big,Q_big)
% view(40,22)
% title(['Time Evolution of Peclet Problem - Pe=' num2str(Pe) '  \Theta=' num2str(theta) '  \Delta t=' num2str(DELTA_T)])
% xlabel('t')
% ylabel('x')
% zlabel('q')
