# Sods problem: Oscillations with WENO schemes

 Register Blogs Members List Search Today's Posts Mark Forums Read

 May 31, 2015, 17:57 Sods problem: Oscillations with WENO schemes #1 New Member   Join Date: May 2015 Posts: 1 Rep Power: 0 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. Code: ```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!