Poisson Solve Gauss-Seidel and SOR (MATLAB) 01

PoissonSolve-SOR_01.m

% Same as 00, but we explictlly calcuclate the RHS from the R,G Components. 
 
clear all;
format long; 

% Square unit domain
N = 64;
dx = 1 / (N-1);
img_path = "C:\Users\Niall\Pictures\rock.png";
img_form = "png";

% Use 64^2 Image as input "vector field" to derive RHS from.
im = imread(img_path, img_form);
% R,G --> U,V components
im_x = im(:, :, 1);
im_y = im(:, :, 2);
% Cast
im_x = cast(im_x, "double");
im_y = cast(im_y, "double");
% Normalize
im_x = im_x / 255;
im_y = im_y / 255;
%im = im * (0.5 * dx);

% Unlike _01 we explictlly RHS calc divergence (using R,G components)
% This better shows what a "real" pressure solution looks like. 
rhs = zeros(N, N);
for i=2:N-1
    for j=2:N-1
        xx = im_x(i, j+1) - im_x(i, j-1);
        yy = im_y(i+1, j) - im_y(i-1, j);
        rhs(i, j) = (0.5 * dx) * (xx + yy); 
    end
end

n_iter      = 512; 
resid_k     = zeros(n_iter, 1); 
resid_k_l10 = zeros(n_iter, 1); 
% Presure and RHS vectors (Assume they are 2D Grids, not "real" matrices)
pn = zeros(N, N);
%b = im;
b = rhs; 
% Spatial Grid vectors (per dim) 
x = 0:dx:1;
y = 0:dx:1;

% Gauss-Seidel + SOR.
for it = 1:n_iter
    for i=2:N-1
        for j=2:N-1
            pn_k = (1/4) * (pn(i+1, j) + pn(i-1, j) + pn(i, j+1) + pn(i, j-1) + b(i,j));
            pn(i, j) = (1-1.8) * pn(i, j) + 1.8 * pn_k; % SOR \omega = 1.8
            %pn(i, j) = pn_k; % GS 
        end
    end

    % Calc Residual (it)
    r = zeros(N, N);
    for i=2:N-1
        for j=2:N-1
            r(i,j) = (b(i, j) + pn(i+1, j) + pn(i-1, j) + pn(i, j+1) + pn(i, j-1) - 4*pn(i, j));
        end
    end
    resid_k(it)     = norm(r, "inf");
    resid_k_l10(it) = log10(resid_k(it));
end

% Plot Solution (Pressure) Surface
figure(1) 
subplot(3,1,1) 
surf(x, y, pn)
title('Solution Surf');
xlabel('X');
ylabel('Y');
zlabel('Solution');

% Plot Residual Norm Log10
subplot(3,1,2); 
plot((1:n_iter), resid_k_l10);
title('Residual Plot Log');
xlabel('Iteration K');
ylabel('Residual Norm log10');

% Plot Residual Norm 
subplot(3,1,3); 
plot((1:n_iter), resid_k);
title('Residual Plot');
xlabel('Iteration K');
ylabel('Residual Norm');