Homec4science
Diffusion Computational physics (master) Browse

compphys/exo4/heat_problem.m
09de8b211167master

Last Change
Raw File

heat_problem.m
View Options

global N
N = 25
L = 10.0
global dx
dx = L / (N+1);
A = zeros(N^2);
b = zeros(N^2,1);
sd_cg = true
j_precondition = true
global T1
global T2
T1 = 2
T2 = 1
% bind values to the hamiltonian
for i = 1:N
for j = 1:N
% diagonal
A(grid_k(i,j), grid_k(i,j)) = 4;
% off diagonal, letting boundary unchanged
if i > 1
A(grid_k(i,j), grid_k(i-1,j)) = -1;
end
if j > 1
A(grid_k(i,j), grid_k(i,j-1)) = -1;
end
if j < N
A(grid_k(i,j), grid_k(i,j+1)) = -1;
end
if i < N
A(grid_k(i,j), grid_k(i+1,j)) = -1;
end
end
end
% boundary conditions
% for each x distretization (j): find x boundary
for j = 1:N
% determine b
[Tb, nb] = bars(1,j);
b(grid_k(1,j)) = Tb; % i = 1
A(grid_k(1,j), grid_k(1,j)) = contacts(i,j) + nb;
[Tb, nb] = bars(N,j);
b(grid_k(N,j)) = Tb; % i = N
A(grid_k(N,j), grid_k(N,j)) = contacts(i,j) + nb;
end
% for each y distretization (i): find x boundary
for i = 1:N
% determine b
[Tb, nb] = bars(i,1);
b(grid_k(i,1)) = Tb; % j = 1
A(grid_k(i,1), grid_k(i,1)) = contacts(i,j) + nb;
[Tb, nb] = bars(i,N);
b(grid_k(i,N)) = Tb; % j = N
A(grid_k(i,N), grid_k(i,N)) = contacts(i,j) + nb;
end
% Jacobi preconditioning
if j_precondition
fprintf("Old conditioning: %g\n", cond(A))
M = diag(1 ./ diag(A));
A = M * A;
b = M * b;
fprintf("New conditioning: %g\n", cond(A))
end
% solve
X = A \ b;
if sd_cg
[X_sd, err_sd] = solve_SD(A, b, 100000);
fprintf("Number of iterations for SD: %d\n", length(err_sd))
[X_cg, err_cg] = solve_CG(A, b, 1000);
fprintf("Number of iterations for CG: %d\n", length(err_cg))
fprintf("|X - X_SD| = %g\n", norm(X_sd - X))
fprintf("|X - X_CG| = %g\n", norm(X_cg - X))
end
T = reshape(X, N, N)';
[Sx, Sy] = eval_stream(T);
% plot solution
figure
hold on
contourf(x(1:N), y(1:N), T, 15, 'LineStyle','None')
quiver(x(1:N), y(1:N), Sx, Sy,'k')
xlabel('x [m]')
ylabel('y [m]')
hold off
% tools
% hypothesis, i,j are boundary indexes
function [val, nb] = bars(i,j)
global T1
global T2
nb = 0;
val = 0;
if is_Bar1(i,j)
val = val + T1;
nb = nb + 1;
end
if is_Bar2(i,j)
val = val + T2;
nb = nb + 1;
end
end
function nb = contacts(i,j)
global N
nb = (i > 1) + (j > 1) + (i < N) + (j < N);
end
function flag = is_Bar1(i,j)
global N
flag = (j == 1) && (i <= N/2);
end
function flag = is_Bar2(i,j)
global N
flag = (j == N) && (i >= N/2);
end
function k = grid_k(i, j)
global N
k = (i-1) * N + j;
end
function pos = x(j)
global dx
pos = (j - 0.5) * dx;
end
function pos = y(i)
global dx
pos = (i - 0.5) * dx;
end
% T : NxN reshaped resultant temperature grid
% Sx : Nx(N+1)
% Sy : (N+1)xN
function [Sx, Sy] = eval_stream(T)
global N
global dx
global T1
global T2
Sx = zeros(N,N);
Sy = zeros(N,N);
% internal
% for each cell compute left and bottom current
for i = 2:N-1
for j = 2:N-1
Sx(i,j) = -(T(i,j+1) - T(i,j-1))/(2*dx);
Sy(i,j) = -(T(i+1,j) - T(i-1,j))/(2*dx);
end
end
% boundaries (! particular case)
for i = 1:N
% bottom and upper are already null
% left
Sx(i,1) = -(T(i,1) - T1)/(2*dx) * is_Bar1(i,1);
% right
Sx(i,N) = -(T2 - T(i,N))/(2*dx) * is_Bar2(i,N);
end
% angles
end
c4science · Help