CFD Basics: Code Vectorization

     This post is about comparing 2 codes to solve the 2D Laplace equation using finite difference method. A sample code mentioned under "Code 01" uses nested loops. We must, wherever possible avoid nested loops. The solution to the code is shown in Fig. 1. The second code uses vectorization instead of the nested loops. The vectorized version is mentioned under  "Code 02".

Code 01

clear

clc

close all

%% Parameters

Lx = 1; % Length of the domain in the x-direction

Ly = 1; % Length of the domain in the y-direction

Nx = 201; % Number of grid points in the x-direction

Ny = 201; % Number of grid points in the y-direction

dx = Lx / (Nx - 1); % Grid spacing in the x-direction

dy = Ly / (Ny - 1); % Grid spacing in the y-direction

%% Initialize temperature matrix

T = zeros(Nx, Ny);

T(:, 1) = 100;

T(:, Nx) = 0;

T(1, :) = 25;

T(Ny, :) = 50;

%% Gauss-Seidel iteration

max_iter = 50000; % Maximum number of iterations

tolerance = 1e-10; % Convergence tolerance

error = inf; % Initialize error

iter = 0; % Iteration counter

while error > tolerance && iter < max_iter

T_old = T;

% Solve Laplace equation using Gauss-Seidel iterations

for i = 2:Nx-1

for j = 2:Ny-1

T(i, j) = ((T(i+1, j) + T(i-1, j))*dy^2 + (T(i, j+1) + T(i, j-1))*dx^2) / (2*(dx^2 + dy^2));

end

end

% Compute error

error = max(abs(T(:) - T_old(:)));

% Increment iteration counter

iter = iter + 1;

end

%% plotting

[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);

contourf(X, Y, T')

axis equal

colormap jet

colorbar

clim([0 100])

title('Temperature Distribution Non Vec')

xlabel('x')

ylabel('y')

zlabel('Temperature (T)')

colorbar

Code 02

clear

clc

close all

%% Parameters

Lx = 1; % Length of the domain in the x-direction

Ly = 1; % Length of the domain in the y-direction

Nx = 201; % Number of grid points in the x-direction

Ny = 201; % Number of grid points in the y-direction

dx = Lx / (Nx - 1); % Grid spacing in the x-direction

dy = Ly / (Ny - 1); % Grid spacing in the y-direction

i = 2:Nx-1;

j = 2:Ny-1;

%% Initialize temperature matrix

T = zeros(Nx, Ny);

T(:, 1) = 100;

T(:, Nx) = 0;

T(1, :) = 25;

T(Ny, :) = 50;

%% Gauss-Seidel iteration

max_iter = 50000; % Maximum number of iterations

tolerance = 1e-10; % Convergence tolerance

error = inf; % Initialize error

iter = 0; % Iteration counter

while error > tolerance && iter < max_iter

T_old = T;

% Solve Laplace equation using Gauss-Seidel iterations

T(i, j) = ((T(i+1, j) + T(i-1, j))*dy^2 + (T(i, j+1) + T(i, j-1))*dx^2) / (2*(dx^2 + dy^2));

% Compute error

error = max(abs(T(:) - T_old(:)));

% Increment iteration counter

iter = iter + 1;

end

%% plotting

[X, Y] = meshgrid(0:dx:Lx, 0:dy:Ly);

contourf(X, Y, T')

axis equal

colormap jet

colorbar

clim([0 100])

title('Temperature Distribution Vec')

xlabel('x')

ylabel('y')

zlabel('Temperature (T)')

colorbar

Result

Results say that vectorized code is ~1.5x faster than nested looped code for 200x200 matrix. Simulation results are now presented.

Fig. 1, Vectorized VS Nested Loops

Thank you for reading! I hope you learned something new! If you like this blog and want to hire me as your PhD student, please get in touch!

1D Laplace's equation using Finite Difference Method

This post is about FDM for Laplace Equation with various boundary conditions.

MATLAB Code (1D, Dirichlet Boundary Conditions)

%% initialize the workspace, clear the command window

clear; clc

%% finite difference 1D laplace dirichlet boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, u(L) = 4 %% Ax=b%%

N = 4 ; %number of grid points

L = 1; %length of domain

dx = L/(N-1); %element size

%% initialize variables %%

l = linspace(0,L,N); %independent

u=zeros(1,N); %dependent

%% boundary conditions %%

u(1)=10;

u(end)=4;

%% b vector %%

b=zeros(N-2,1);

b(1) = b(1) - u(1);

b(end) = b(end) - u(end);

%% A matrix

A = -2*eye(N-2,N-2);

for i=1:N-2

    if i<N-2

        A(i,i+1)=1;

    end

    if i>1

        A(i,i-1)=1;

    end

end

%% solve for unknowns %%

x = A\b;

%% fill the u vector with unknowns %%

u(2:end-1) = x;

%% plot the results %%

hold on; grid on; box on, grid minor

set(gca,'FontSize',40)

set(gca, 'FontName', 'Times New Roman')

ylabel('u','FontSize',44)

xlabel('l','FontSize',44)

plot(l,u,'-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)

MATLAB Code (1D, Mixed Boundary Conditions)

%% initialize the workspace, clear the command window

clear; clc

%% finite difference 1D laplace mixed boundary conditions %% d2u/dx2 = 0 %% u(o) = 10, du/dx(L) = 4 %% Ax=b%%

N = 5; %number of grid points

L = 1; %length of domain

dx = L/(N-1); %element size

a = 4;

%% initialize variables %%

l = linspace(0,L,N); %independent

u=zeros(1,N); %dependent

%% dirichlet boundary condition %%

u(1) = 10;

%% b vector %%

b=zeros(N-1,1);

b(1) = b(1) - u(1);

b(end) = b(end) + dx*a; %neumann boundary condition added to b vector

%% A matrix

A = -2*eye(N-1,N-1);

for i=1:N-1

    if i<N-1

        A(i,i+1)=1;

    end

    if i>1

        A(i,i-1)=1;

    end

end

A(N-1,N-1) = -1; %neumann boundary condition added to A matrix

%% solve for unknowns %%

x = A\b;

% fill the u vector with unknowns %%

u(2:end) = x;

%% plot the results %%

hold on; grid on; box on, grid minor

set(gca,'FontSize',40)

set(gca, 'FontName', 'Times New Roman')

ylabel('u','FontSize',44)

xlabel('l','FontSize',44)

plot(l,u','-o','color',[0 0 0],'LineWidth',2,'MarkerSize',20)