% Task 1, Computational Fluid Dynamics of turbulent flow
%
% Outline of 2D stationary heat conduction code - a help to get
% you started with Matlab.
% 
% By Håkan Nilsson (hani@tfd.chalmerse.se), 1999-01-20

% Set parameters
max_e= 0.01;   % Maximum error - Check if this is small enough!!!
k    = 1;      % Conduction coeff - vary this to see how it affects physics
               % and numerics.
s    = -1.5;   % Source Term
L    = 1.;     % Length of area
H    = 0.5;    % Heigth of area

% Grid parameters
% The number of nodes in each direction are set by N_i and N_j.
% If the boundary nodes are put on the boundaries, the number of
% grid points becomes N_im1=N_i-1 and N_jm1=N_j-1
% General advice: Think about how you number your nodes with
% i,j / rows,colums etc. Mixing this up is one of the most common
% problems. Matlab uses Origo for the i-j plane in the upper left
% corner of the matrix with i going vertical and j going horizontal.
% In CFD, the Origo for the i-j plane is usually placed in the
% lower left of the grid with i going horizontal and j going vertical.
N_i = 10;     % Cells in i (usually equivalent to x) direction
N_j = 10;     % Cells in j (usually equivalent to y) direction
N_im1=N_i-1;   % Grid points in i direction
N_jm1=N_j-1;   % Grid points in j direction

% Create the grid. You should refine the grid in regions where you
% expect large temperature gradients!
for j = 1:N_jm1
  for i = 1:N_im1
    x_grid(i,j)=   ;
    y_grid(i,j)=   ;
  end
end

% Calculate the node positions, in the center of the control volumes
for j = 1:N_j
  for i = 1:N_i
    x_node(i,j)=   ;
    y_node(i,j)=   ;
  end
end

% Calculate constant geometric properties, for example delta_x(i,j)
% and delta_y(i,j)

% Initialize things.
T    = ones(N_i,N_j);  % Note that T is a property of the nodes,
                       % therefore N_i and N_j is used.

% Prepare boundary conditions, for example
T(1,:)=10;

% Main iteration loop

error=1;
while (error > max_e)
  Told=T;

  % Depending on the case the following should 
  % be done here or inside/outside the loop:
  %
  % Fix boundary conditions ...
  % Calculate a_P, a_E, a_W, a_S, a_N ...
  % Calculate b ...
  % Loop over all interior nodes and update T - One Gauss-Seidel sweep

  error=norm(T-Told)
end

% Plot your results. See 'help mesh' in matlab

