% sor2d: Solves model elliptic problem using SOR ...

% Discretization parameters ...
level = 6;
n = 2^level + 1;
h = 1.0 / (n - 1);

% Define optimal SOR parameter ...
omega = 2.0 / (1.0 + sin(pi * h))

% Vectors of x , y coordinates for plotting ...
[x, y] = meshgrid(linspace(0.0, 1.0, n), linspace(0.0, 1.0, n));

% Frequencies for exact solution ...
omega_x = pi;
omega_y = pi;

[uxct f] = calc_exact(omega_x, omega_y, n, h);

% Set boundary conditions and initialize solution to random values between
% 0 and 1 ...
u = zeros(n,n);
u(2:n-1,2:n-1) = rand(n-2,n-2);

% Number of relaxation sweeps ...
nsweep = 1000;

% Set to 1 to enable .avi creation
avienable = 1;
avisecs = 20;
aviframerate = 25;
aviframes = avisecs * aviframerate;
avifilename = 'sor.avi';

% AVI initialization and output ...
if avienable
   avifreq = round(nsweep / aviframes);
   aviobj = VideoWriter(avifilename);
   open(aviobj);

   uerr = uxct - u;
   surf(x, y, u);
   title('Sweep 0');

   framecount = 5 * aviframerate
   for iframe = 1 : framecount
      writeVideo(aviobj, getframe(gcf));
   end
end

% For each relaxation sweep ...
for isweep = 1 : nsweep
   rnorm = 0.0;
   for i = 2 : n - 1
      for j = 2 : n - 1
         % Update unknown and accumulate residual ...
         rij = (u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - ...
                4*u(i,j)) / h^2 - f(i,j);
         rnorm = rnorm + rij^2;
         % Compute GS value ...
         uij_gs = u(i,j) + 0.25 * h^2 * rij;
         % Overrelaxation step ...
         u(i,j) = omega * uij_gs + (1.0 - omega) * u(i,j);
      end
   end
   % Compute residual and error norms and output ...
   rnorm = sqrt(rnorm / n^2);
   uerr = uxct - u;
   fprintf('sweep: %d  |res|: %g  |uerr|: %g\n', isweep, rnorm, ...
           norm(uerr)/n);

   % AVI output ...
   if mod(isweep, avifreq) == 0
      surf(x, y, u);
      titlestr = sprintf('Sweep %g\n', isweep);
      title(titlestr);
      writeVideo(aviobj, getframe(gcf));
   end
end

% Close avi file ...
if avienable
   close(aviobj);
   fprintf('Created video file: %s\n', avifilename);
end

