1
$\begingroup$

I'm trying to solve the advection equation $$m_t+(\alpha m)_x=0$$ with $m(0,\cdot)=m_0$ numerically using the first order Godunov scheme. Hence I write $$m_j^{i+1}=m_j^i-\frac{dt}{dx}(m_{j+1/2}^i\alpha_{j+1/2}^i-m_{j-1/2}^i\alpha_{j-1/2}^i)$$ where $m_{j+1/2}^i=m_{j+1}^i$ if $\alpha_{j+1/2}^i<0$ and $m_j^i$ otherwise.

Here's a Matlab code with $\alpha=1/2$ everywhere:

%% parameters
N=99; M=99; 

%Allocate memory
x=linspace(0,1,M+1); 
mk=zeros(N+1,M+1); a_ph=ones(N,M)*0.5; 

%% Set initial conditions
sigma=0.08;
mk(1,:)=1/sqrt(2*pi*sigma^2)*exp(-(x-0.5).^2/(2*sigma^2));

q=dt/h;

for i=1:N
    ai=ak_ph(i,2:M); aim=a_ph(i,1:M-1); %local var
    B=diag([0,ai.*(ai>=0),0]) + diag([0,ai.*(ai<0)],1) ...
       - diag([0,aim.*(aim<0),0]) - diag([aim.*(aim>=0),0],-1);
    mk(i+1,:)=(eye(M+1)-q*B)*mk(i,:)';
end

[X1,X2]=meshgrid(x1,x2);
surf(X1,X2,mk','FaceColor','interp','EdgeColor','none')

When I compute the solution this is what I get:

enter image description here

I expect a wave travelling at a uniform velocity but the amplitude gets smaller in the meantime which is strange. There are also ripples at the peak of the wave. Does anyone have a suggestion?

$\endgroup$
  • $\begingroup$ Welcome to SciComp.SE! It seems that you are using two different accounts, which is why you cannot edit your own posts (which you should be able to do, independent of your reputation). You might want to register and then contact the SE moderators and ask them to merge the accounts. $\endgroup$ – Christian Clason Aug 19 '15 at 17:09
  • 2
    $\begingroup$ Can you plot some snapshots of your solution in time (say 3) instead of doing the surface plot? I assume that the oscillation at the peak of the wave is simply a plotting artifact but it's hard to see if the shape of your wave is changing the way it should. $\endgroup$ – Kyle Mandli Aug 20 '15 at 18:01
  • $\begingroup$ Could be numerical diffusion? Is the peak also getting broader? $\endgroup$ – boyfarrell Aug 20 '15 at 20:18
  • $\begingroup$ This method is diffusive, which is why you see the amplitude decreasing. Whatever reference you used to learn about the method surely also discusses this effect. $\endgroup$ – David Ketcheson Aug 21 '15 at 17:44
  • 2
    $\begingroup$ @DavidKetcheson: I think the $x$-axis of the plot is increasing from right-to-left, which is the mirror image of the typical orientation seen in plots, hence why the wave is moving to the left. $\endgroup$ – Geoff Oxberry Aug 21 '15 at 22:27
2
$\begingroup$

Godunov's scheme is a tricky scheme that chooses the stencil based on the direction of the propagation of the wave. Because of the hyperbolic character of the equation, we know in advection problems that the information should go from left to right if the wave speed is positive and vice versa for negative wave speed. Godunov's scheme described this trick with the mathematical expression that you have written that in your expression like "$m_{j+1/2}^i=m_{j+1}^i$ if $\alpha_{j+1/2}^i<0$ and $m_j^i$ otherwise". If the sign is reversed it will lead to unphysical oscillations!

Why this scheme is diffusive and how to get rid of that.

For simplicity I'm assuming wave speed $\alpha$ =1; Then the exact solution of this problem, can be written as $$m^{n+1}_i=m^n_{i-1}$$

For positive $\alpha$ the finite difference equation can be written as $$m^{n+1}_i=m^n_i-CFL*(m^n_i-m^n_{i-1})$$ or $$m^{n+1}_I=CFL*m^n_{i-1}+m^n_i(1-CFL)$$ here CFL=Courant–Friedrichs–Lewy number =$\alpha \frac{\delta t}{\delta x}$ If we do Von-Newman stability analysis on that problem, it will restrict our $CFL$ number less than one.

Still our discretized result is not matching with the exact result. Let put $CFL$ =1, then the result will match with exact result. the second term $m^n_i(1-CFL)$ is called numerical viscosity, that is the term that diffuses the result. You can notice that that term increases with reducing the $CFL$ number. That means artificial viscosity is increasing, will decreases the signal amplitude (you have observed this in your code) .

Please try to change the $CFL$ number and analyze your result you will understand it clearly. If you want to see the effect of grid spacing then you should try some higher order TVD schemes or ENO with different initial conditions.

Please note that analysis I discussed here is only valid for forward Euler time discretization and may not be valid for other time discretizations. For that we should do phase velocity plot and numerical diffusion coefficient plot for that scheme. As a rule of thumb we can say odd order schemes are diffusive and even order schemes are dispersive. In most of the good numerical schemes this diffusion and dispersion can be controlled by changing CFL or cell Reynolds number and grid size.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.