Can you help me write this question in finite volume method? This is my first time writing FVM. I wrote using FDM below :(

%Upwind scheme
n=320;
X=n+1;
Y=n+1;
u=0.5;
v=0.5;
Lx=4;
Ly=4;
dX=Lx/n;
dY=Ly/n;
c(:,:,1)=zeros(Y,X);
c(:,:,2)=zeros(Y,X);
for a=(2+0.9/dX):1:(1.1/dX)
    for b=(2+0.9/dY):1:(1.1/dY)
        x=(a-1)*dX;
        y=(b-1)*dY;
        c(a,b,1)=sin(5*pi*(x-0.9))*sin(5*pi*(y-0.9));
    end
end


dT=0.001;
for T=1:1:4
    for n=2:1:(round(T/dT)+1)
        for j=2:1:(Lx/dX)
            for i=2:1:(Ly/dY)
                c(i,j,2)=c(i,j,1)+(v*dT/dX)*(c(i-1,j,1)-c(i,j,1))+(u*dT/dY)*(c(i,j-1,1)-c(i,j,1));
                
            end
        end
        c(:,:,1)=c(:,:,2);
    end
    dT=dT+0.001;
    fig = figure();
    mesh(0:dX:Lx,0:dY:Ly,c(:,:,2));
    xlabel('X');
    ylabel('Y');
    zlabel('c');
    grid on;
    axis([0 4 0 4 -1 1]);
end

Answer 1

[Brahmadev]

Hi @Shengwei Bie,

The MATLAB code for First order upwind seems to be on the right track, here are some more things to consider:

1. 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.
2. 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!