The MATLAB code for First order upwind seems to be on the right track, here are some more things to consider:
- The time loop has an increment dT=dT+0.001; at the end, which doesn't make sense in this context since dT is your time step and should remain constant as defined at the beginning.
- The boundary conditions are not enforced in your code. You should set the edges of the computational domain to zero at each time step.
Here is a revised version of the code that takes care of these considerations:
% Define the domain and grid size
n = 320;
Lx = 4;
Ly = 4;
dx = Lx / n;
dy = Ly / n;
x = linspace(0, Lx, n+1);
y = linspace(0, Ly, n+1);
[X, Y] = meshgrid(x, y);
% Define the time domain
dt = 0.001;
t_final = 4;
Nt = round(t_final / dt);
% Define the initial condition
c = zeros(n+1, n+1);
for i = 2:n
for j = 2:n
if (x(i) > 0.9 && x(i) < 1.1) && (y(j) > 0.9 && y(j) < 1.1)
c(j, i) = sin(5 * (x(i) - 0.9)) * sin(5 * (y(j) - 0.9));
end
end
end
% Define velocities
u = 0.5;
v = 0.5;
% FVM implementation
for t = 1:Nt
c_new = c;
% Compute fluxes and update the scalar field c
for i = 2:n
for j = 2:n
% Compute fluxes using upwind scheme
flux_x = u * (c(j, i) - c(j, i-1)) / dx;
flux_y = v * (c(j, i) - c(j-1, i)) / dy;
% Update c based on net fluxes
c_new(j, i) = c(j, i) - dt * (flux_x + flux_y);
end
end
% Apply boundary conditions
c_new(1, :) = 0; % Bottom boundary
c_new(end, :) = 0; % Top boundary
c_new(:, 1) = 0; % Left boundary
c_new(:, end) = 0; % Right boundary
c = c_new;
% Plot at specific time steps
if mod(t, Nt / 4) == 0 || t == 1
figure;
surf(X, Y, c, 'EdgeColor', 'none');
title(sprintf('FVM with Upwind Scheme at t = %.3f', t * dt));
xlabel('x');
ylabel('y');
zlabel('c');
view(3);
colorbar;
pause(0.1); % To visualize the evolution
end
end
Hope this helps in resolving your query!
