blown up solution for linear advection in upwind method with finite difference

I am going to solve this advection equation regarding the flow simulation of an energy tower $$$$Y_t+ v(t) Y_x =0$$$$ with the following boundary conditions which depend on the direction of the flow \begin{align} Y(x,0) = Y_0 \end{align} And \begin{align} Y(0,t) = Y_{bot} \text{ if }v(t)>0 \\ Y(1,t) = Y_{top} \text{ if }v(t)<0 \end{align} And for the simulations, I have written a Matlab code in the upwind method with different inputs to $$Y_0, Y_{bot}$$ and $$Y_{top}$$. But almost all the solution seems to blow up in time. Can someone tell me what has gone wrong here?

t_end = 10;
n_x = 10;

y0 = 1;  % Inputs for the initial condition Y(x,0)
y_bot = 0; %  Inputs for the boundary condition Y(0,t)
y_top = 1; %  Inputs for the boundary condition Y(1,t)

dx = 1/n_x;

v0 =1; % initial velocity
CFL = 1;
dt = CFL*(dx)/((abs(v0)));

n_t = round(t_end/dt)
x = linspace(0,1,n_x);  % space grid
t = linspace(0,t_end,n_t); % time grid

Y = zeros(n_x, n_t);

Y(:,1)= y0;  % initial condition Y_0
Y(1,:)=y_bot; % boundary condition at the bottom of the tower
% Y(n_x,:)= y_top; % boundary condition at the top of the tower

v = sort(randi([-20 50],1,n_t),'descend'); % generating a random array of velocities in the decending order
v(1) = v0; % make the initial velocity as v0

for n=2:n_x
for j=2:n_x-1

if v(n) >= 0
Y(1,n) = y_bot;
else
if v(n_x) < 0
Y(n_x,n) = y_top;
end
end

if v(n) >= 0
Y(j,n) = Y(j,n-1)-v(n)*dt/dx*(Y(j,n-1)-Y(j-1,n-1));
else
if v(n)< 0
Y(j,n) = Y(j,n-1)-v(n)*dt/dx*(Y(j+1,n-1)-Y(j,n-1));
end
end
end
end


• The line if v(n_x) < 0 is meaningless, always true. Even corrected to if v(n) < 0 it is always true if it is reached. Commented Nov 21, 2022 at 7:03

You are using n_x at at least two places as boundary for the time index.

Loops are slow, so avoid them where not necessary. The time loop can be condensed to

for n=2:n_t
if v(n) >= 0
Y(1,n) = y_bot;
Y(2:n_x,n) = Y(2:n_x,n-1)-v(n)*dt/dx*(Y(2:n_x,n-1)-Y(1:n_x-1,n-1));
else
Y(n_x,n) = y_top;
Y(1:n_x-1,n) = Y(1:n_x-1,n-1)-v(n)*dt/dx*(Y(2:n_x,n-1)-Y(1:n_x-1,n-1));
end
end


Your CFL-dt calculation needs to take the highest velocity into account. You compute it for v0=1, but the actual highest values in magnitude are 50 and -20. Being sneaky and setting v(1)=v0 does not change the general problem, the second largest velocity is still too large.

• Thank you for the feedback. But to take the highest velocity by taking max(v) to calculate time step $dt$, according to the CFL. how can I define $v = sort(randi([-20 50],1,n_t)$, since $n_t$ depends on $dt$?
– TMW
Commented Nov 21, 2022 at 9:14
• You take the theoretical maximum 50. If you want truth in advertizing, linearly stretch the random array after construction to exactly fit the boundaries. Commented Nov 21, 2022 at 9:25
• Thank you. This is a question on a bit improved step, I want to know if my $v(t)$ is also developing inside the loop and how I can get the max velocity.? This happens when $v(t)$ also comes out as a solution of an ODE w.r.t. time with some initial velocity $v_0$
– TMW
Commented Nov 21, 2022 at 10:21