[Ironiker]

Hey everyone!

When applying my WENO code for the Sods Shock Tube problem, I get some minor oscillations in the solutions:


I am using a central compact WENO scheme with SSP time stepping, which works well in the scalar case (Transport, Burgers, etc.). Initial data and numerical constants are taken from the paper
http://epubs.siam.org/doi/abs/10.1137/S1064827599359461

Thus the mistake must be in the numerical process of cell updating/ calculation the right hand side. I am just copying the necessary code fragments to maintain clarity.

- U is a matrix evaluated rho, m and E at every cell midpoint.
- The function WENO gives back the reconstructed value at the left side and at the right side for each cell for u1, u2, and u3.
- To calculate solutions at cell boundaries I am shifting the solution to the right, calculate the numerical flux and then shift that to the left due to conservation.
- F and the Jacobian matrix are written down explicitly and checked twice.

U = [rho; m; E]; 

%%%%% Calculating righthandside for time stepping %%%%%

function du = rhs(u)
vL2,vR2 = WENO(U)
z = f_num(shmatrixperiodl(vR2),vL2);
y = shmatrixperiodr(z);
du = -(y-z);
end

%%%%% Shift function routines %%%%%

function v = shmatrixperiodl(U)
v = [U(:,1) U(:,1:(end-1))];
end


function v = shmatrixperiodr(U)
v = [U(:,2:end) U(:,end)];
end

%%%%% Calculating numerical flux %%%%%

function Z = f_num(A,B)
% Numerical flux function
% Here: local lambdax-Friedrichs flux
[m,n]=size(A);
C=sparse(m,n);
for i=1:n
C(:,i) = max(abs(eig(Jacobmatrix(A(:,i)))), abs(eig(Jacobmatrix(B(:,i)))) );
end
Z = (F(A)+F(B)-C.*(B-A))./2;

%%%%% Flux Function %%%%%

function Y = F(U)
global gamma
Y = [U(2,:); 
     0.5*(3-gamma)*U(2,:).^2./U(1,:)+(gamma-1)*U(3,:); 
     gamma*U(3,:).*U(2,:)./U(1,:)+0.5*(1-gamma)*U(2,:).^3./(U(1,:).^2)]; 
end


%%%%% Jacobian %%%%%%

function Y = Jacobmatrix(U)
global gamma
Y = [ 0                                      1                                    0 
    -0.5*(3-gamma)*U(2).^2./(U(1).^2)      (3-gamma)*U(2)./U(1)                  (gamma-1);
    -gamma*U(3).*U(2)./(U(1).^2)+(gamma-1)*U(2).^3./(U(1).^3)  gamma*U(3)./U(1)+1.5*(1-gamma)*U(2).^2/(U(1).^2)     gamma*U(2)./U(1)];
end

I would be glad for some helpful advices!