Numerical Bug in HLLC Solver for Sod Shock Tube Problem

Hi all,

I was playing around with implementing Toro's HLLC algorithm in MATLAB. I think there is something silly I'm overlooking because I'm not getting the right answer when modeling Sod's shock tube. My pressure approaches zero where the initial discontinuity is located.

In the attached code, I'm using Test 1 from Chapter 10 of his book (i.e., Riemann Solvers and Numerical Methods for Fluid Dynamics). Please note I'm using the FVM and using ghost cells at the boundaries. Also note for my conservation of energy equation, I'm using rho*E as my conserved variable, not E as he uses in his book.

Any help would be much appreciated and I'd be happy to share the code with the community once it's working.

Thanks,

C

clear; clc; close all;
%% Parameters
nx = 50;                 % Number of physical cells
gamma = 1.4;              % Ratio of specific heats
CFL = 0.3;                % CFL number
t_final = 0.2;            % Final simulation time
diaphragm = 0.3;          % Initial discontinuity position
x = linspace(0, 1, nx);   % Physical cell centers
dx = x(2) - x(1);
% Include ghost cells (2 extra)
U = zeros(nx + 2, 3);     % [rho, rho*u, rho*E]
x_g = linspace(-dx, 1+dx, nx + 2);  % Extended x with ghost cell placeholders
%% Initial Conditions (Sod problem)
rhoL = 1.0; uL = 0.75; pL = 1.0;
rhoR = 0.125; uR = 0.0; pR = 0.1;
% Encode energy
E = @(rho, u, p) p ./ (rho .* (gamma - 1)) + 0.5 * u.^2;
% Assign initial values
for i = 1:nx+2
    if x_g(i) < diaphragm
        U(i,:) = [rhoL, rhoL*uL, rhoL*E(rhoL, uL, pL)];
    else
        U(i,:) = [rhoR, rhoR*uR, rhoR*E(rhoR, uR, pR)];
    end
end
%% Time integration loop
t = 0;
while t < t_final
    
    [rho, u, p, ~] = decode(U, gamma);
    a = sqrt(gamma * p ./ rho);
    dt = CFL * dx / max(abs(u) + a);
    if t + dt > t_final
        dt = t_final - t;
    end
    dF = compute_dF(U, gamma);
    % Only update interior cells
    U(2:end-1,:) = U(2:end-1,:) - dt * dF / dx;
    U = correctBC(U);  % Apply boundary conditions
    if any(isnan(U(:)))
        error("NaN detected in solution at time t = %g", t);
    end
    t = t + dt;
end
%% Plot results (exclude ghost cells)
[rho, u, p, ~] = decode(U, gamma);
x_phys = x_g(2:end-1);  % Strip ghost cells
figure;
plot(x_phys, rho(2:end-1), 'b-', 'LineWidth', 2); hold on;
plot(x_phys, u(2:end-1),   'r-', 'LineWidth', 2);
plot(x_phys, p(2:end-1),   'g-', 'LineWidth', 2);
legend('\rho', 'u', 'p');
xlabel('x'); ylabel('Solution');
title('1D Euler Equations with HLLC');
grid on;

%% HLLC Flux and Supporting Functions
function dF = compute_dF(U, gamma)
nx = size(U, 1) - 2;  % Exclude ghost cells
[rho, u, p, E] = decode(U, gamma);
% Reconstruct primitive variables at interfaces
[rhoL, rhoR] = iPlusHalf(rho);
[uL, uR]     = iPlusHalf(u);
[pL, pR]     = iPlusHalf(p);
% Reconstruct conserved variables
UL = encode(rhoL, uL, pL, gamma);
UR = encode(rhoR, uR, pR, gamma);
% Compute fluxes at interfaces
F = zeros(nx + 1, 3);
for i = 1:nx+1
    F(i,:) = hllc_flux(UL(i,:), UR(i,:), gamma);
end
% Flux differences for physical cells
dF = zeros(nx, 3);
for i = 1:nx
    dF(i,:) = F(i+1,:) - F(i,:);
end
end
function F = hllc_flux(UL, UR, gamma)
a = @(p, rho) sqrt(gamma * p ./ rho);
[rhoL, uL, pL, EL] = decode(UL, gamma);
[rhoR, uR, pR, ER] = decode(UR, gamma);
pmin = min(pL,pR);
pmax = max(pL,pR);
qmax = pmax/ pmin;
quser = 2.0;
aL = a(pL, rhoL); aR = a(pR, rhoR);
rhoBar = 0.5 * (rhoL + rhoR);
aBar = 0.5 * (aL + aR);
ppv = 0.5 * (pL + pR) - 0.5 * (uR - uL) * rhoBar * aBar;
pstar = 0;
if (qmax < quser) && (pmin<=ppv) && (ppv <= pmax)
    pstar = ppv;
else
    if ppv < pmin
        % use the two-expansion approximation, TRRS
        zz = (gamma - 1.0) / (2.0 * gamma);
        pstar = ((aL + aR - (gamma-1.0)/2.0*(uR - uL)) / (aL/(pL^zz) + aR/(pR^zz)))^ (1.0 / zz);
    else
        % use the TSRS solution
        arg_gL = (2.0/(gamma+1)/rhoL)/((gamma-1.0)/(gamma+1.0)*pL + ppv);
        arg_gR = (2.0/(gamma+1)/rhoR)/((gamma-1.0)/(gamma+1.0)*pR + ppv);
        % Add checks to prevent negative arguments to sqrt
        arg_gL = max(arg_gL, 0);
        arg_gR = max(arg_gR, 0);
        gL = sqrt(arg_gL);
        gR = sqrt(arg_gR);
        % Add check to prevent division by zero
        denom = gL + gR;
        if denom == 0
            denom = 1e-8; % Small value to avoid division by zero
        end
        pstar = (gL * pL + gR * pR - (uR - uL))/(denom);
    end
end
qL = computeQ(gamma, pstar, pL);
qR = computeQ(gamma, pstar, pR);
SL = uL - aL * qL;
SR = uR + aR * qR;
denom = max(abs(rhoL*(SL - uL) - rhoR*(SR - uR)), 1e-8);
Sstar = (pR - pL + rhoL*uL*(SL - uL) - rhoR*uR*(SR - uR)) / denom;
FL = flux(rhoL, uL, pL, EL);
FR = flux(rhoR, uR, pR, ER);
UStarL = intermediate_state(SL, uL, Sstar, rhoL, EL, pL);
UStarR = intermediate_state(SR, uR, Sstar, rhoR, ER, pR);
F = flux_selection(SL, Sstar, SR, FL, FR, UStarL, UL, UStarR, UR);
end
function q = computeQ(gamma, pstar, pk)
if pstar <= pk
    q = 1;
else
    q = sqrt(1 + (gamma+1)/(2*gamma)*(pstar/pk - 1));
end
end
function F = flux(rho, u, p, E)
F = zeros(1, 3);
F(1) = rho * u;
F(2) = rho * u^2 + p;
F(3) = u * (rho*E + p);
end
function UStar = intermediate_state(S, u, Sstar, rho, E, p)
den = S - u;
UStar(1) = rho * (S - u) / (S - Sstar);
UStar(2) = UStar(1) * Sstar;
UStar(3) = UStar(1) * (E + (Sstar - u)*(Sstar + p / (rho * den)));
end
function F = flux_selection(SL, Sstar, SR, FL, FR, UStarL, UL, UStarR, UR)
if 0 <= SL
    F = FL;
elseif SL <= 0 && 0 <= Sstar
    F = FL + SL * (UStarL - UL);
elseif Sstar <= 0 && 0 <= SR
    F = FR + SR * (UStarR - UR);
else
    F = FR;
end
end
function [rho, u, p, E] = decode(U, gamma)
rho = U(:,1);
u   = U(:,2) ./ rho;
E   = U(:,3) ./ rho;
p = (rho .* E - 0.5 * rho .* u.^2) * (gamma - 1);
p = max(p, 1e-8);   % Prevent negative pressure
end
function U = encode(rho, u, p, gamma)
U(:,1) = rho;
U(:,2) = rho .* u;
U(:,3) = p ./ (gamma - 1) + 0.5 * rho .* u.^2;
end
function U = correctBC(U)
U(1,:)     = U(2,:);
U(end,:)   = U(end-1,:);
end
function [thetaL, thetaR] = iPlusHalf(theta)
n = length(theta) - 2;  % physical cells
thetaL = zeros(n + 1, 1);  % one per interface
thetaR = zeros(n + 1, 1);
for i = 1:n
    left  = theta(i);
    right = theta(i+1);
    slope = minmod(right - left, left - theta(i));
    thetaL(i) = left + 0.5 * slope;
    thetaR(i) = right - 0.5 * slope;
end
% Last interface fallback (between cell nx and nx+1)
thetaL(end) = theta(n);
thetaR(end) = theta(n+1);
end
function m = minmod(a, b)
m = 0.5 * (sign(a) + sign(b)) .* min(abs(a), abs(b));
end

[Torsten]

Just a sidenote:

If you are in need of an excellent solver for hyperbolic PDEs, I suggest CLAWPACK which is available under

https://www.clawpack.org/

Answer 1

[CHowdy]

Hi all,

I figured out the issue. The code works as intended now. I'm only using a forward Euler time integrator, but it could be modified for other explicit time integrators. Also, different reconstruction schemes could be used for the face values. If anybody wants to make it more modular and share for the community, that would be swell. Please see below:

% Sod Shocktube problem
clear; clc; close all;
%% Parameters
nx = 50;                 % Number of physical cells
gamma = 1.4;              % Ratio of specific heats
CFL = 0.3;                % CFL number
t_final = 0.2;            % Final simulation time
diaphragm = 0.3;          % Initial discontinuity position
x = linspace(0, 1, nx);   % Physical cell centers
dx = x(2) - x(1);
% Include ghost cells (2 extra)
U = zeros(nx + 2, 3);     % [rho, rho*u, rho*E]
x_g = linspace(-dx, 1+dx, nx + 2);  % Extended x with ghost cell placeholders
%% Initial Conditions (Sod problem)
rhoL = 1.0; uL = 0.75; pL = 1.0;
rhoR = 0.125; uR = 0.0; pR = 0.1;
% Encode energy
E = @(rho, u, p) p ./ (rho .* (gamma - 1)) + 0.5 * u.^2;

% Assign initial values
for i = 1:nx+2
    if x_g(i) < diaphragm
        U(i,:) = [rhoL, rhoL*uL, rhoL*E(rhoL, uL, pL)];
    else
        U(i,:) = [rhoR, rhoR*uR, rhoR*E(rhoR, uR, pR)];
    end
end

%% Time integration loop
t = 0;
while  t < t_final
    
    [rho, u, p, ~] = decode(U, gamma);
    % Compute speed of sound 
    a = sqrt(gamma * p ./ rho);
    dt = CFL * dx / max(abs(u) + a);
    if t + dt > t_final
        dt = t_final - t;
    end
    % Compute divF
    dFdx = compute_divF(U, gamma, dx);
    % Only update interior cells
    U(2:end-1,:) = U(2:end-1,:) - dt * dFdx;
    U = correctBC(U);  % Apply boundary conditions
    if any(isnan(U(:)))
        error("NaN detected in solution at time t = %g", t);
    end
    t = t + dt;
end

%% Plot results (exclude ghost cells)
[rho, u, p, ~] = decode(U, gamma);
x_phys = x_g(2:end-1);  % Strip ghost cells
figure;
plot(x_phys, rho(2:end-1), 'b-', 'LineWidth', 2); hold on;
plot(x_phys, u(2:end-1),   'r-', 'LineWidth', 2);
plot(x_phys, p(2:end-1),   'g-', 'LineWidth', 2);
legend('\rho', 'u', 'p');
xlabel('x'); ylabel('Solution');
title('1D Euler Equations with HLLC');
grid on;

%% HLLC Flux and Supporting Functions
function dFdx = compute_divF(U, gamma, dx)
nx = size(U, 1) - 2;  % Exclude ghost cells
[rho, u, p, ~] = decode(U, gamma);
% Reconstruct primitive variables at left and right side of faces
[rhoL, rhoR] = iPlusHalf(rho);
[uL, uR]     = iPlusHalf(u);
[pL, pR]     = iPlusHalf(p);
% Reconstruct conserved variables at interfaces
UL = encode(rhoL, uL, pL, gamma);
UR = encode(rhoR, uR, pR, gamma);
% Compute fluxes at interfaces using HLLC method
F = zeros(nx + 1, 3);
for i = 1:nx+1
    F(i,:) = hllc_flux(UL(i,:), UR(i,:), gamma);
end
% Compute divergence for physical cells
dFdx = zeros(nx, 3);
for i = 1:nx
    dFdx(i,:) = (F(i+1,:) - F(i,:))/dx;
end
end
function F = hllc_flux(UL, UR, gamma)
a = @(p, rho) sqrt(gamma * p ./ rho);
[rhoL, uL, pL, EL] = decode(UL, gamma);
[rhoR, uR, pR, ER] = decode(UR, gamma);
aL = a(pL, rhoL); aR = a(pR, rhoR);
rhoBar = 0.5 * (rhoL + rhoR);
aBar = 0.5 * (aL + aR);
ppv = 0.5 * (pL + pR) - 0.5 * (uR - uL) * rhoBar * aBar;
pstar = max(0, ppv);
qL = computeQ(gamma, pstar, pL);
qR = computeQ(gamma, pstar, pR);
SL = uL - aL * qL;
SR = uR + aR * qR;
denom = rhoL*(SL - uL) - rhoR*(SR - uR);
Sstar = (pR - pL + rhoL*uL*(SL - uL) - rhoR*uR*(SR - uR)) / denom;
FL = flux(rhoL, uL, pL, EL);
FR = flux(rhoR, uR, pR, ER);
UStarL = intermediate_state(SL, uL, Sstar, rhoL, EL, pL);
UStarR = intermediate_state(SR, uR, Sstar, rhoR, ER, pR);
F = flux_selection(SL, Sstar, SR, FL, FR, UStarL, UL, UStarR, UR);
end

function q = computeQ(gamma, pstar, pk)
if pstar <= pk
    q = 1;
else
    q = sqrt(1 + (gamma+1)/(2*gamma)*(pstar/pk - 1));
end
end

function F = flux(rho, u, p, E)
F = zeros(1, 3);
F(1) = rho * u;
F(2) = rho * u^2 + p;
F(3) = u * (rho*E + p);
end

function UStar = intermediate_state(S, u, Sstar, rho, E, p)
den = S - u;
UStar(1) = rho * (S - u) / (S - Sstar);
UStar(2) = UStar(1) * Sstar;
UStar(3) = UStar(1) * (E + (Sstar - u)*(Sstar + p / (rho * den)));
end

function F = flux_selection(SL, Sstar, SR, FL, FR, UStarL, UL, UStarR, UR)
if 0 <= SL
    F = FL;
elseif SL <= 0 && 0 <= Sstar
    F = FL + SL * (UStarL - UL);
elseif Sstar <= 0 && 0 <= SR
    F = FR + SR * (UStarR - UR);
else
    F = FR;
end
end

function [rho, u, p, E] = decode(U, gamma)
rho = U(:,1);
u   = U(:,2) ./ rho;
E   = U(:,3) ./ rho;
p = (rho .* E - 0.5 * rho .* u.^2) * (gamma - 1);
p = max(p, 1e-8);   % Prevent negative pressure
end

function U = encode(rho, u, p, gamma)
U(:,1) = rho;
U(:,2) = rho .* u;
U(:,3) = p ./ (gamma - 1) + 0.5 * rho .* u.^2;
end

function U = correctBC(U)
% Zero-gradient (extrapolation) BC
U(1,:)     = U(2,:);
U(end,:)   = U(end-1,:);
end

function [thetaL, thetaR] = iPlusHalf(theta)
n = length(theta) - 2;  % Number of physical cells
thetaL = zeros(n + 1, 1);  % Interface values on the left
thetaR = zeros(n + 1, 1);  % Interface values on the right
for i = 1:n
    left  = theta(i);
    right = theta(i+1);
    slope = 0; %minmod(right - left, left - theta(i+2));
    thetaL(i) = left + 0.5 * slope;
    thetaR(i) = right - 0.5 * slope;
end

% Last interface fallback (between cell nx and nx+1)
thetaL(end) = theta(end-1);
thetaR(end) = theta(end);
end
% function m = minmod(a, b)
%     m = 0.5 * (sign(a) + sign(b)) .* min(abs(a), abs(b));
% end