[zen3gr]

Hello everybody.

I'm trying to learn some basic theory on CFD. Currently I'm studying the MacCormak scheme, I've read on wikipedia that it should give very accurate results in case of non linear PDEs so I'm trying it for the 1D Burger's Inviscid Equation. In the picture attached, the variables of matlab used in the code below are presented with the corresponding mathematical equation. In the second picture, you can see a graph of the solution with two methods: with MacCormak method and a simple forward in time, forward in space with CFL condition; it's clear that the former is pretty bad in contrast with the second method that seems pretty correct.

Matlab Codes:

--MacCormak

L=2; %length of domain
nx = 41; %x nodes
dx = L./(nx-1);
nt =70; %nt is the number of timesteps we want to calculate
dt = 0.025; %dt is the amount of time each timestep covers (delta t)
%IC
u = ones(1,nx); %matrix of 1 row and nx columns with all values = 1
u(1,0.5/dx : 1/dx+1)=2;
%initializing matrices
partial_t=ones(1,nx);
u_predicted=partial_t;
partial_t_dt=u_predicted;
partial_ave=partial_t_dt;

for t=1:1:nt
u_p=u; % u_p matrix is where u values from previous time step are stored and used for partial_t and u_predicted

for n=1:1:nx-1 %predictor
partial_t(n)=-u_p(n)*((u_p(n+1)-u_p(n))/(dx));
u_predicted(n)=u_p(n)+partial_t(n)*dt;

end

for n=2:nx %corrector

partial_t_dt(n)=-u_predicted(n)*((u_predicted(n)-u_predicted(n-1))/(dx));
partial_ave(n)=0.5*(partial_t(n)+partial_t_dt(n));
u(n)=u_p(n)+partial_ave(n)*dt;
end


plot(linspace(0,L,nx),u) %%Plot the results
ylim([0 3]);
pause(.1)

end

--CFL - Simple Scheme

L=2;
nx = 41;
dx = L./(nx-1);
nt =70;
%CFL
sigma=.5;
dt=sigma*dx
%IC
u = ones(1,nx);
u(1,round(.5/dx):round(1 / dx + 1))=2;
for t=1:1:nt
u_p=u;
for n=2:1:nx
u(n)=u_p(n)-u_p(n)*dt/dx*(u_p(n)-u_p(n-1));
end
fprintf('t=%d -> u_prediction=%.2f \n', t, u(20));
plot(linspace(0,L,nx),u) %%Plot the results
ylim([0 3]);
%dokimh gia na fainetai h timh ths U mono stous komvous

pause(.1)
% hold on
end

Have you any thoughts about the bad performance of the MacCormak Method? I think the implementation is right so I've ran out of ideas about what maybe is wrong.

Thanks in advance for your time

[FMDenaro]

Originally Posted by zen3gr

Hello everybody.

I'm trying to learn some basic theory on CFD. Currently I'm studying the MacCormak scheme, I've read on wikipedia that it should give very accurate results in case of non linear PDEs so I'm trying it for the 1D Burger's Inviscid Equation. In the picture attached, the variables of matlab used in the code below are presented with the corresponding mathematical equation. In the second picture, you can see a graph of the solution with two methods: with MacCormak method and a simple forward in time, forward in space with CFL condition; it's clear that the former is pretty bad in contrast with the second method that seems pretty correct.

Matlab Codes:

--MacCormak


--CFL - Simple Scheme


Have you any thoughts about the bad performance of the MacCormak Method? I think the implementation is right so I've ran out of ideas about what maybe is wrong.

Thanks in advance for your time

And you really want to learn the basic CFD reading wikepedia??

I strongly suggest to read non fundamental textbooks, you can find a lot of basic reading.
For example, read this text and try to repeat the proposed test first on the linear wave equation then on the non linear Burgers equation.
https://www.google.com/url?sa=t&rct=...5CwoG7Gt8UmeMX

[zen3gr]

Originally Posted by FMDenaro

And you really want to learn the basic CFD reading wikepedia??

I strongly suggest to read non fundamental textbooks, you can find a lot of basic reading.
For example, read this text and try to repeat the proposed test first on the linear wave equation then on the non linear Burgers equation.
https://www.google.com/url?sa=t&rct=...5CwoG7Gt8UmeMX

Thanks for the reply.

Actually I'm reading the book by Anderson (Computational Fluid Dynamics, which covers, the very basic things required for a solid base) in combination with some video lectures for that purpose. I only quoted Wikipedia because it states that MacCormak should give satisfying results and that's not the case with my attempt.

I would like a comment on the actual code. I'm gonna study the eBook you offered to get some more insight on the numerical schemes and try to troubleshoot it on my own.

Thanks anyway, really appreciate your time!

[duri]

Originally Posted by zen3gr

Hello everybody.

I'm trying to learn some basic theory on CFD. Currently I'm studying the MacCormak scheme, I've read on wikipedia that it should give very accurate results in case of non linear PDEs so I'm trying it for the 1D Burger's Inviscid Equation. In the picture attached, the variables of matlab used in the code below are presented with the corresponding mathematical equation. In the second picture, you can see a graph of the solution with two methods: with MacCormak method and a simple forward in time, forward in space with CFL condition; it's clear that the former is pretty bad in contrast with the second method that seems pretty correct.

Matlab Codes:

--MacCormak


--CFL - Simple Scheme


Have you any thoughts about the bad performance of the MacCormak Method? I think the implementation is right so I've ran out of ideas about what maybe is wrong.

Thanks in advance for your time

One of the method is wrong. Not sure which one because initial condition and time step is not given, I guess most likely simple is wrong. Plot the analytical solution and check the results.
As initial step try to validate constant wave speed (Linear convection equation) and check that the initial condition have traveled exact distance on given time using x=ct relation.

[zen3gr]

Originally Posted by duri

One of the method is wrong. Not sure which one because initial condition and time step is not given, I guess most likely simple is wrong. Plot the analytical solution and check the results.
As initial step try to validate constant wave speed (Linear convection equation) and check that the initial condition have traveled exact distance on given time using x=ct relation.

Thanks for your input.
I think the simple/CFL code captures the essence of phenomenon pretty correct.
Time step is 0.025 for both cases as is the I.C. (u is everywhere 1 except some points that its value is two("hat" function), see attached graph).

I've solved the linear convection equation before proceeding to the nonlinear case; I'm not sure how the bold text will help me with the implementation of maccormak.

Thanks again for your time!

[zen3gr]

I've tried to solve 1d Linear convection equation with the maccormak scheme; it turns out I should have done it before attempting to solve the nonlinear one but I thought it would be too easy, that's not the case unfortunately.

In the pictures attached you can see the solution for two schemes: maccormak and a simple one with plain forward relations for time and space derivatives as well as the initial conditions and a interesting graph I found at Anderson's book about CFD and Heat Transfer about mathematical dispersion. I suppose that this kind of behaviour is typical for the second order MacCormak scheme (the wiggles of the solution compared to the simple scheme one).

My question is: what can be done about it? Is something wrong with my code?

Matlab Code:

clear all
clc
L=2;
nx = 41;
dx = L./(nx-1);
nt = 70; %nt is the number of timesteps we want to calculate
dt = .025; %dt is the amount of time each timestep covers (delta t)
c=1;

%IC
u = ones(1,nx);
u(1,0.5/dx : 1/dx+1)=2;
u_cfl=u;
for t=1:1:15
u_p=u;
u_cfl_p=u_cfl;
for n=2:1:nx-1
%---MacCormak---%
partial_t(n)=-(c/dx)*(u_p(n+1)-u_p(n));
u_bar(n)=u_p(n)+partial_t(n)*dt;
derivative_bar(n)=-(c/dx)*(u_bar(n)-u_bar(n-1));
u(n)=u_p(n)+0.5*(partial_t(n)+derivative_bar(n))*d t;
%---Simple FTFS scheme---%
u_cfl(n)=u_cfl_p(n)-c*dt/dx*(u_cfl_p(n)-u_cfl_p(n-1));
end
plot(linspace(0,2,nx),u,linspace(0,2,nx),u_cfl,'r' )
ylim([0 3]);
legend('MacCormak','Simple FTFS');
pause(.1)
end

EDIT: For dx=dt => C (from CFL condition) =1 the solution using MacCormak Scheme looks pretty satisfying, except the part "behind the movement of the wave". (4th picture attached)


Thanks in advance.

[tas38]

Originally Posted by zen3gr

I've tried to solve 1d Linear convection equation with the maccormak scheme; it turns out I should have done it before attempting to solve the nonlinear one but I thought it would be too easy, that's not the case unfortunately.

In the pictures attached you can see the solution for two schemes: maccormak and a simple one with plain forward relations for time and space derivatives as well as the initial conditions and a interesting graph I found at Anderson's book about CFD and Heat Transfer about mathematical dispersion. I suppose that this kind of behaviour is typical for the second order MacCormak scheme (the wiggles of the solution compared to the simple scheme one).

My question is: what can be done about it? Is something wrong with my code?

Matlab Code:


EDIT: For dx=dt => C (from CFL condition) =1 the solution using MacCormak Scheme looks pretty satisfying, except the part "behind the movement of the wave". (4th picture attached)


Thanks in advance.

I have added an alternative implementation of MacCormack to your code below. It exhibits the dispersive (but non-diffusive) errors which you would expect. I think your implementation may have some subtle errors.

clear all
clc
L=2;
nx = 41;
dx = L./(nx-1);
nt = 70; %nt is the number of timesteps we want to calculate
dt = .025; %dt is the amount of time each timestep covers (delta t)
c=1;

%IC
u1 = ones(1,nx);
u_pred = zeros(1,nx);
u_pred(1) = c;
u1(1,0.5/dx : 1/dx+1)=2;
u_cfl=u1;
u2 = u1;
for t=1:1:15
u_p1=u1;
u_p2=u2;
u_cfl_p=u_cfl;
for n=2:1:nx-1
%---MacCormak---%
partial_t(n)=-(c/dx)*(u_p1(n+1)-u_p1(n));
u_pred(n) = u2(n)-c*dt/dx*(u2(n+1)-u2(n));
u_bar(n)=u_p1(n)+partial_t(n)*dt;
derivative_bar(n)=-(c/dx)*(u_bar(n)-u_bar(n-1));
u1(n)=u_p1(n)+0.5*(partial_t(n)+derivative_bar(n))*dt;
u2(n)=0.5*(u2(n)+u_pred(n)) - 0.5*c*dt/dx*(u_pred(n)-u_pred(n-1));
%---Simple FTFS scheme---%
u_cfl(n)=u_cfl_p(n)-c*dt/dx*(u_cfl_p(n)-u_cfl_p(n-1));
end
plot(linspace(0,2,nx),u1,linspace(0,2,nx),u_cfl,'r',linspace(0,2,nx),u2,'o')
ylim([0 3]);
legend('MacCormack','Simple FTFS','alt-MacCormack');
pause(.1)
end

[FMDenaro]

Originally Posted by zen3gr

I've tried to solve 1d Linear convection equation with the maccormak scheme; it turns out I should have done it before attempting to solve the nonlinear one but I thought it would be too easy, that's not the case unfortunately.

In the pictures attached you can see the solution for two schemes: maccormak and a simple one with plain forward relations for time and space derivatives as well as the initial conditions and a interesting graph I found at Anderson's book about CFD and Heat Transfer about mathematical dispersion. I suppose that this kind of behaviour is typical for the second order MacCormak scheme (the wiggles of the solution compared to the simple scheme one).

My question is: what can be done about it? Is something wrong with my code?

Matlab Code:


EDIT: For dx=dt => C (from CFL condition) =1 the solution using MacCormak Scheme looks pretty satisfying, except the part "behind the movement of the wave". (4th picture attached)


Thanks in advance.

You have a bug in the code. Even if Godunov theorem states that a linear second order scheme is not monotone, if you work on a smooth initial condition the numerical oscillations are disregardable.

[zen3gr]

Originally Posted by tas38

I have added an alternative implementation of MacCormack to your code below. It exhibits the dispersive (but non-diffusive) errors which you would expect. I think your implementation may have some subtle errors.

Thank you very much, your answer has been very illustrating. I guess that line "u_pred(1)=c" makes all the difference. I've adjusted it to my implementation and the results are pretty good. I'm going to try to solve the MacCormak scheme for a few nodes by hand and figure out why the term is equal to c. Any guidance would be much appreciated.

clear all
clc
L=2;
nx = 41;
dx = L./(nx-1);
nt = 70; %nt is the number of timesteps we want to calculate
dt = .025; %dt is the amount of time each timestep covers (delta t)
c=1;

%IC
u1 = ones(1,nx);
u_pred = zeros(1,nx);
u_pred(1) = c;
u_bar(1)=c;
u1(1,0.5/dx : 1/dx+1)=2;
u_cfl=u1;
u2 = u1;
for t=1:1:15
u_p1=u1;
u_p2=u2;
u_cfl_p=u_cfl;
for n=2:1:nx-1

%---ANDERSON_MacCormak---%
        
        u_bar(n)=u_p(n)-(c/dx)*(u_p(n+1)-u_p(n))*dt;  
            derivative_bar(n)=-(c/dx)*(u_bar(n)-u_bar(n-1));
            u(n)=u_p(n)+0.5*(-(c/dx)*(u_p(n+1)-u_p(n))+derivative_bar(n))*dt;
%---MacCormak---%
u_pred(n) = u2(n)-c*dt/dx*(u2(n+1)-u2(n));
u2(n)=0.5*(u2(n)+u_pred(n)) - 0.5*c*dt/dx*(u_pred(n)-u_pred(n-1));
%---Simple FTFS scheme---%
u_cfl(n)=u_cfl_p(n)-c*dt/dx*(u_cfl_p(n)-u_cfl_p(n-1));
end
plot(linspace(0,2,nx),u1,linspace(0,2,nx),u_cfl,'r',linspace(0,2,nx),u2,'o')
ylim([0 3]);
legend('MacCormack','Simple FTFS','alt-MacCormack');
pause(.1)
end

What I would like to ask now, before going further with the MacCormak method, trying to solve the nonlinear equation, is if it's wrong to use values from previous time steps. I mean in the %---ANDERSON_MacCormak---% section I've used values from time step n to evaluate the while on your the formula, tas38, the is evaluated using values from the same time level. The results,, for both cases seem to be identical. However, I'm thinking that wouldn't be a "prediction" if it was supposed to be calculated from previous time steps. So some clarification about that would help me greatly.

Thanks!

[FMDenaro]

Originally Posted by zen3gr

Thank you very much, your answer has been very illustrating. I guess that line "u_pred(1)=c" makes all the difference. I've adjusted it to my implementation and the results are pretty good. I'm going to try to solve the MacCormak scheme for a few nodes by hand and figure out why the term is equal to c. Any guidance would be much appreciated.

clear all
clc
L=2;
nx = 41;
dx = L./(nx-1);
nt = 70; %nt is the number of timesteps we want to calculate
dt = .025; %dt is the amount of time each timestep covers (delta t)
c=1;

%IC
u1 = ones(1,nx);
u_pred = zeros(1,nx);
u_pred(1) = c;
u_bar(1)=c;
u1(1,0.5/dx : 1/dx+1)=2;
u_cfl=u1;
u2 = u1;
for t=1:1:15
u_p1=u1;
u_p2=u2;
u_cfl_p=u_cfl;
for n=2:1:nx-1

%---ANDERSON_MacCormak---%
        
        u_bar(n)=u_p(n)-(c/dx)*(u_p(n+1)-u_p(n))*dt;  
            derivative_bar(n)=-(c/dx)*(u_bar(n)-u_bar(n-1));
            u(n)=u_p(n)+0.5*(-(c/dx)*(u_p(n+1)-u_p(n))+derivative_bar(n))*dt;
%---MacCormak---%
u_pred(n) = u2(n)-c*dt/dx*(u2(n+1)-u2(n));
u2(n)=0.5*(u2(n)+u_pred(n)) - 0.5*c*dt/dx*(u_pred(n)-u_pred(n-1));
%---Simple FTFS scheme---%
u_cfl(n)=u_cfl_p(n)-c*dt/dx*(u_cfl_p(n)-u_cfl_p(n-1));
end
plot(linspace(0,2,nx),u1,linspace(0,2,nx),u_cfl,'r',linspace(0,2,nx),u2,'o')
ylim([0 3]);
legend('MacCormack','Simple FTFS','alt-MacCormack');
pause(.1)
end

What I would like to ask now, before going further with the MacCormak method, trying to solve the nonlinear equation, is if it's wrong to use values from previous time steps. I mean in the %---ANDERSON_MacCormak---% section I've used values from time step n to evaluate the while on your the formula, tas38, the is evaluated using values from the same time level. The results,, for both cases seem to be identical. However, I'm thinking that wouldn't be a "prediction" if it was supposed to be calculated from previous time steps. So some clarification about that would help me greatly.

Thanks!

Be aware that on the liner wave equation you should get exactly the Lax-Wendroff scheme. Try that a compare the results.