%
% Solution of unsteady 1D diffusion equation withinin the domain [0,1] using
% explicit second order finite difference method.
% Dirichlet and Neumann boundary conditions are implemented.
%
% This is a part of teaching materials provided for the course
% "Computational Heat & Fluid Flow (ME 605)" taught at IIT Goa
% during the winter semester of 2018
%
% Author: Y Sudhakar, IIT Goa
% email: sudhakar@iitgoa.ac.in
%
% Tested with Matlab R2018b on 16 January 2020
%
% Assignment for the students
% 1. check solution at large time steps -- recall condition of stability
% 2. Implement Robin boundary condition
% 3. Implement one-sided 2nd order approximation for Neumann BC
% 4. Write the iteration number and error-norm in a file
% 5. Introduce a source term and solve the problem
%

clc; clear all; close all;

% --------------- Set the following input parameters ------------------
num_mesh = 11;                  % number of mesh points
del_t    = 0.00001;             % size of the time step
alpha    = 1.0;                 % Thermal diffusivity
bc       = 'Dirichlet';         % Boundary condition at the right boundary; set either 'Dirichlet', 'Neumann' or 'Robin'
% ---------------------------------------------------------------------

del_x    = 1.0/(num_mesh-1.0);       % mesh size (Delta_x)
xmesh    = zeros(num_mesh,1);        % Mesh points
u_old    = zeros(num_mesh,1);        % solution at the previous time step
u_new    = zeros(num_mesh,1);        % Solution at the new time step
steady_term = zeros(num_mesh,1);     % RHS of the discretized governing equation (approaches zero atsteady state)

%del_t = del_x*del_x*0.25/alpha; %-- to demonstrate stability
del_t = del_x*del_x*0.52/alpha; %-- to demonstrate stability

% Compute the location of mesh points
for i=1:num_mesh
    xmesh(i) = (i-1) * del_x;
end

% Solve the discretized equation until steady state is achieved
norm = 100.0;   % measure of difference between values at two consequetive time steps
iter = 0;
converged = false;
delx_sq = del_x^2;
while converged==false
    % the discretized equation
    for i=2:num_mesh-1
        u_new(i) = u_old(i) + alpha * del_t / delx_sq * (u_old(i-1)+u_old(i+1)-2.0*u_old(i));
    end
    u_new = apply_bc(u_new,bc,num_mesh);
    
    for i=2:num_mesh-1
        steady_term(i) = abs(alpha*(u_new(i+1)-2.0*u_new(i)+u_new(i-1))/(del_x^2.0));
    end
    
    norm = max(steady_term); % compute the difference between two solutions
    
    iter = iter + 1;
    if norm <0.00001 || iter == 100 % to demonstrate stability
    %if norm <0.00001
        converged = true;
    end
    
    % Uncomment the following lines to uncommented for stability
    if rem(iter,25) == 0
        str = strcat(num2str(iter),'iter');
        plot(xmesh,u_new,'x-','DisplayName',str,'LineWidth',5);
        hold on;
    end
    legend('FontSize',36)
    
    % Soution not converged -- copy current value into old
    for i=1:num_mesh
        u_old(i) = u_new(i);
    end
end

plot(xmesh,u_new,'x-');
xlabel('x');
ylabel('Temperature');

%----------------------------------------------------------------
% This function applies boundary conditions on the given array
function u = apply_bc(u,bc,num_mesh)
    u(1) = 100.0;
    if strcmp(bc,'Dirichlet')
        u(num_mesh) = 50.0;
    elseif strcmp(bc,'Neumann')
        u(num_mesh) = u(num_mesh-1);
    else
        msg = 'Set correct boundary condition at the right side boundary';
        error(msg);
    end
end
