[TurbJet]

Hello,

I am solving a nonlinear hyperbolic equation using simple upwind scheme with discontinuity in the solution. But when a peak occurred, the solution blows up, just as the one in the attached figure, the right-hand side peak keeps increasing until inf eventually. When I tried Lax-Friedrichs, it all went well: the peak reaches 100 and stop increasing, which is expected.

Upwind is much more diffusive than LF, so I think this blow-up should not happened. But I really don't know why this is happening. Is it upwind invalid for this situation?

Any thoughts?

1.jpg

Thx.

[FMDenaro]

Provided that you are fulfilling the stability constraint, You have for sure a bug in the code.

[TurbJet]

Originally Posted by FMDenaro

Provided that you are fulfilling the stability constraint, You have for sure a bug in the code.

With fixed mesh size, I tried smaller time steps, but they did not change anything. So I think it has nothing wrong with my CFL condition.


On the other hand, I checked and checked my code, I don't see any bugs.

The equation I am solving is

,

where, and c are constant, but varies with x.

Anyway, denoting , than I write this code for upwind:



where unew stands for the u for next time level,

I don't know where it goes wrong.

Codes for the upwind (for MATLAB):

for i = 2:N+2
unew(i) = u(i) - b* (v0(i)*(u(i)-u(i)^2/c) - v0(i-1)*(u(i-1)-u(i-1)^2/c));
end

[FMDenaro]

The upwind criterion is applied considering the characteristic direction, therefore you have to consider both the signs of the advective velocity. Are you sure that your problem has always a positive advection component?

[TurbJet]

Originally Posted by FMDenaro

The upwind criterion is applied considering the characteristic direction, therefore you have to consider both the signs of the advective velocity. Are you sure that your problem has always a positive advection component?

Yes, the convective velocity will always be positive.

[TurbJet]

Originally Posted by FMDenaro

The upwind criterion is applied considering the characteristic direction, therefore you have to consider both the signs of the advective velocity. Are you sure that your problem has always a positive advection component?

I'll just paste my entire code here. It's not a long code; if you have time, please take a look.

clearvars; clc; close all

L = 50; % length of computation domain
N = 100; % set up mesh interval
dx = L/N; % mesh size

% coefficients
a = 0;
c = 100;

t = 0; % starting time level
T = 1; % ending time level
dt = 0; % time step size

iter = 0; % iteration counter
% set up time level for saving the data for plotting
tc = 1; % time advancing counter
timeLevel = [0:0.01:1]; % time levels for saving solution
tl = length(timeLevel);

% set up spatial coordinate
x = -dx : dx : 50+dx;
for i = 1:N+3
% find index for x = 0, 10, 30, 40
if (x(i) == 0)
i0 = i;
elseif (abs(x(i) - 10.0) < 1e-6)
i10 = i;
elseif (abs(x(i) - 30.0) < 1e-6)
i30 = i;
elseif (abs(x(i) - 40.0) < 1e-6)
i40 = i;
end
end

% allocate arrays for u(x,t), v_0(x), unew(x,t)
u = zeros(N+3, 1);
unew = zeros(N+3, 1);
v0 = zeros(N+3, 1);

% initializing u(x,t)
for i = 1:N+3
if ((i0 < i) && (i < i10))
u(i) = 50;
end
end

% initializing v_0(x)
for i = 1:N+3
if ((i30 < i) && (i < i40))
v0(i) = 10;
else
v0(i) = 100;
end
end

uarr = zeros(N+3, tl);
uarr(:, 1) = u; tc = tc + 1;


while (t < T)

% time advancing
iter = iter + 1;

a = max(v0.*(1-2*u/c));

if (t < timeLevel(tc) && t + dt > timeLevel(tc))
dt = timeLevel(tc) - t;
else
% dt = dx/(a + 1e-8); % plus a small number in case v = 0
dt = 0.0005;
end
t = t + dt;

% create shorthand
b = dt/dx;

for i = 2:N+2
unew(i) = u(i) - b* (v0(i)*(u(i)-u(i)^2/c) - v0(i-1)*(u(i-1)-u(i-1)^2/c));
end

% apply periodic BC
unew(1) = unew(N+2); unew(N+3) = unew(2);
% update solution
u = unew;

% save solution at given time level
if (t == timeLevel(tc))
uarr(:, tc) = u;
tc = tc + 1;
end


end

[FMDenaro]

Ok, let me some time, I never solved numerically this equation, I will have a look. However, I suggest to write the code in terms of the FV formulation, computing first the fluxes and then performing the update.

[FMDenaro]

I insert the plot into the code to see the evolution step-by-step. I see the left branch becoming a shock wave and the right branch becoming a rarefaction wave. The simulation run until a small wiggle is produced at the basis of the rarefaction wave, near the value u=0.
Therefore, I assume that your solution produces some very small negative value that will be amplified owing to the assumed positive upwind. That is a numerical instability. I suggest to translate the initial solution adding a positive constant (also simply 1) to avoid having small negative velocity. Otherwise, check the negative event and set the value to zero.
I suggest also to replicate the case illustrated in the book of Leveque a page 206, you could also check for the source with the Godunov scheme.

[TurbJet]

Originally Posted by FMDenaro

I insert the plot into the code to see the evolution step-by-step. I see the left branch becoming a shock wave and the right branch becoming a rarefaction wave. The simulation run until a small wiggle is produced at the basis of the rarefaction wave, near the value u=0.
Therefore, I assume that your solution produces some very small negative value that will be amplified owing to the assumed positive upwind. That is a numerical instability. I suggest to translate the initial solution adding a positive constant (also simply 1) to avoid having small negative velocity. Otherwise, check the negative event and set the value to zero.
I suggest also to replicate the case illustrated in the book of Leveque a page 206, you could also check for the source with the Godunov scheme.

First of all, thank you for your time! But I am not sure if I see this wiggle you mentioned. This is what I got when the solution start to blow up, I plot both velocity v and density u.
2.jpg1.jpg
I highlighted the negative velocity and some wiggle(I suppose). Are these what you mentioned?

[FMDenaro]

Originally Posted by TurbJet

First of all, thank you for your time! But I am not sure if I see this wiggle you mentioned. This is what I got when the solution start to blow up, I plot both velocity v and density u.
Attachment 61938Attachment 61939
I highlighted the negative velocity and some wiggle(I suppose). Are these what you mentioned?

Just plot u at each time step and you will see creating a small wiggles at the expansion front (about x=30) that will be amplified. The shock is well resolved.
Thus, check again your code for some bug.

[TurbJet]

Originally Posted by FMDenaro

Just plot u at each time step and you will see creating a small wiggles at the expansion front (about x=30) that will be amplified. The shock is well resolved.
Thus, check again your code for some bug.

You mean this one?
1.jpg2.jpg

If you do mean this one, I think this is what expected to happen, since the wave is kind of "congested" in this section.

[FMDenaro]

Yes, it is.
Are you sure about the way you set the periodic BC.s?

[TurbJet]

Originally Posted by FMDenaro

Yes, it is

If you mean this, I don't think this is the problem because just as I said, it's expected.

[FMDenaro]

Originally Posted by TurbJet

If you mean this, I don't think this is the problem because just as I said, it's expected.

But that is implied by the periodic BC.s you set?
I don't know the physics of the problem but you should see a coalescence in the characteristic lines to generate a discontinuity

[TurbJet]

Originally Posted by FMDenaro

But that is implied by the periodic BC.s you set?
I don't know the physics of the problem but you should see a coalescence in the characteristic lines to generate a discontinuity

I think you're right. But with time goes by, the wave will congested round x=30, as the follows which I computed with Lax-Friedrichs scheme at different time level.

1.jpg2.jpg3.jpg

The peak (or as you said, wiggle) at the most-right will gets thicker and thicker, eventually entire wave will reform into a higher, thinner square wave, comparing to the initial solution. The maximum peak should stop at around u=100, which will never generate negative velocity.

But with upwind, I don't see "congestion" at x=30, so for this peak, instead of getting thicker, it gets higher and higher. This is what I am wondering about.

[FMDenaro]

Originally Posted by TurbJet

I think you're right. But with time goes by, the wave will congested round x=30, as the follows which I computed with Lax-Friedrichs scheme at different time level.

Attachment 61942Attachment 61943Attachment 61944

The peak (or as you said, wiggle) at the most-right will gets thicker and thicker, eventually entire wave will reform into a higher, thinner square wave, comparing to the initial solution. The maximum peak should stop at around u=100, which will never generate negative velocity.

But with upwind, I don't see "congestion" at x=30, so for this peak, instead of getting thicker, it gets higher and higher. This is what I am wondering about.

Ok, the problem depends on the function v(x), it should be better described in the x,t plane by looking at the characteristic curves.
As I wrote before, I never solved this specific problem that is slightly different from the example in the Leveque book.
The rarefaction front should encounter a low or vanishing velocity in such a way to generate compression.
However, I see that continuing the run, after some time the u plot shows negative values.
Could you write me the initial condition u(x,0) and the function v(x)?

[TurbJet]

Originally Posted by FMDenaro

Ok, the problem depends on the function v(x), it should be better described in the x,t plane by looking at the characteristic curves.
As I wrote before, I never solved this specific problem that is slightly different from the example in the Leveque book.
The rarefaction front should encounter a low or vanishing velocity in such a way to generate compression.
However, I see that continuing the run, after some time the u plot shows negative values.
Could you write me the initial condition u(x,0) and the function v(x)?

Sure.

[TurbJet]

Originally Posted by FMDenaro

Ok, the problem depends on the function v(x), it should be better described in the x,t plane by looking at the characteristic curves.
As I wrote before, I never solved this specific problem that is slightly different from the example in the Leveque book.
The rarefaction front should encounter a low or vanishing velocity in such a way to generate compression.
However, I see that continuing the run, after some time the u plot shows negative values.
Could you write me the initial condition u(x,0) and the function v(x)?

I guess you are right for the "low velocity" and "compression" part. The front edge does suffer a low velocity.

[FMDenaro]

Originally Posted by TurbJet

Sure.

Well, I see that v0 produces a variation for x>=30 but the initial condition is upward. That makes me suppose that the starting of the compression region is downwind with respect to the used stencil [i-1,i]. Somehow, that does not respect the upwind criterion. I am not sure if this is the problem, just check the sign of the characteristic. This problem has differences from the standard non-linear Burgers equation.

[FMDenaro]

I wrote few lines about your problem, I suggest to read 16.4.1 and 16.4.2 in the book of Leveque. You will see that the definition of the upwind criterion is quite more complex than you supposed.