clc;
n = 40; % number of elements
finalT = 0.2; % final time
boundary = 1; % switch between Dirichlet conditions (1) and Absorbing (0)
k = 3; % polynomial degree
gamma = 1.4; % specific heat ratio
CFL = 0.1; % Courant number for stability
flux = 1; % flux type, 0 = LF, 1 = HDG
nc = k+1; % number of quadrature points
left = -5; % left end of the domain
right = 5; % right end of the domain
plot_accurate = 1; % plot resolution
%% Quadrature setup
[xg, wg] = get_gauss_quadrature(nc);
xunit = get_gauss_lobatto_quadrature(k+1);
%% Mesh creation
kp1 = length(xunit);
y = zeros(2*n, 1);
h = zeros(n, 1);
for e = 1:n
y(2*e-1) = left + (right-left)*(e-1)/n;
y(2*e) = left + (right-left)*e/n;
h(e) = y(2*e) - y(2*e-1);
end
dx = (right-left)/n;
% Node points
x = zeros(kp1, n);
for e = 1:n
for i = 1:kp1
x(i,e) = y(2*e-1) + (y(2*e)-y(2*e-1))*(0.5 + 0.5*xunit(i));
end
end
% Time step calculation
a_max = sqrt(gamma*1.0/1.0); % max wave speed estimate
dt = CFL*dx/a_max;
if dt <= 0
dt = finalT/100;
end
num_steps = ceil(finalT/dt);
dt = finalT/num_steps;
%% Initial conditions - Sod shock tube
rho = zeros(kp1, n);
u = zeros(kp1, n);
p = zeros(kp1, n);
for e = 1:n
for j = 1:kp1
x_pos = x(j,e);
if x_pos < 0.0
rho(j,e) = 1.0;
u(j,e) = 0.0;
p(j,e) = 1.0;
else
rho(j,e) = 0.125;
u(j,e) = 0.0;
p(j,e) = 0.1;
end
end
end
% Conservation variables
E = p./(gamma-1) + 0.5*rho.*u.^2;
U = zeros(3, kp1, n);
U(1,:,:) = rho;
U(2,:,:) = rho.*u;
U(3,:,:) = E;
%% Basis functions and mass matrix
[psi, dpsi] = evaluate_lagrange_basis(xunit, xg);
Me = psi * diag(wg) * psi';
M1inv = zeros(kp1, kp1, n);
for e = 1:n
M1inv(:,:,e) = inv(0.5*dx*Me);
end
%% Time stepping
t = 0;
for step = 1:num_steps
% Runge-Kutta 3
k1 = rhs_corrected(U, dx, gamma, xunit, psi, dpsi, wg, M1inv);
U1 = U + dt*k1;
k2 = rhs_corrected(U1, dx, gamma, xunit, psi, dpsi, wg, M1inv);
U2 = 0.75*U + 0.25*U1 + 0.25*dt*k2;
k3 = rhs_corrected(U2, dx, gamma, xunit, psi, dpsi, wg, M1inv);
U = (1/3)*U + (2/3)*U2 + (2/3)*dt*k3;
t = t + dt;
end
%% 可视化部分修正
figure;
% 统一使用[]自动计算维度
x_plot = reshape(x, [], 1);
rho_plot = reshape(U(1,:,:), [], 1);
% 修正速度计算,先重塑再除法
rho_reshaped = reshape(U(1,:,:), [], 1);
u_reshaped = reshape(U(2,:,:), [], 1);
u_plot = u_reshaped ./ rho_reshaped;
% 修正压力计算
E_reshaped = reshape(U(3,:,:), [], 1);
p_plot = reshape((gamma-1)*(E_reshaped - 0.5*u_reshaped.^2./rho_reshaped), [], 1);
subplot(3,1,1);
plot(x_plot, rho_plot, 'b.', 'MarkerSize', 10);
title('Density Distribution');
xlabel('x'); ylabel('\rho');
grid on;
subplot(3,1,2);
plot(x_plot, u_plot, 'r.', 'MarkerSize', 10);
title('Velocity Distribution');
xlabel('x'); ylabel('u');
grid on;
subplot(3,1,3);
plot(x_plot, p_plot, 'g.', 'MarkerSize', 10);
title('Pressure Distribution');
xlabel('x'); ylabel('p');
grid on;
%% Corrected right-hand side function
function rhs = rhs_corrected(U, dx, gamma, xunit, psi, dpsi, weights, M1inv)
[~, kp1, n] = size(U);
rhs = zeros(size(U));
J = dx/2; % Jacobian
% Compute Euler fluxes at solution points
F = zeros(3, kp1, n);
for e = 1:n
for i = 1:kp1
rho = U(1,i,e);
u_val = U(2,i,e)/rho;
E_val = U(3,i,e);
p_val = (gamma-1)*(E_val - 0.5*rho*u_val^2);
F(1,i,e) = rho*u_val;
F(2,i,e) = rho*u_val^2 + p_val;
F(3,i,e) = u_val*(E_val + p_val);
end
end
% Element integral term: -∫dφ/dx·F dx
for e = 1:n
% Interpolate flux to quadrature points
F_q = zeros(3, length(weights));
for q = 1:length(weights)
F_interp = [0; 0; 0];
for i = 1:kp1
F_interp = F_interp + psi(i,q)*F(:,i,e);
end
F_q(:,q) = F_interp; % Corrected assignment
end
% Compute integral for each component
for comp = 1:3
integral = zeros(kp1,1);
for i = 1:kp1
for q = 1:length(weights)
integral(i) = integral(i) + weights(q)*dpsi(i,q)*F_q(comp,q);
end
end
rhs(comp,:,e) = -J * integral;
end
% Apply inverse mass matrix
for comp = 1:3
rhs(comp,:,e) = (M1inv(:,:,e) * rhs(comp,:,e)')';
end
end
% Boundary flux terms
for e = 1:n
% Left boundary
if e == 1
% Left boundary condition
rho_L = 1.0;
u_L = 0.0;
p_L = 1.0;
E_L = p_L/(gamma-1) + 0.5*rho_L*u_L^2;
U_L = [rho_L; rho_L*u_L; E_L];
% Current element left state
rho_R = U(1,1,e);
u_R = U(2,1,e)/rho_R;
E_R = U(3,1,e);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
U_R = [rho_R; rho_R*u_R; E_R];
% Numerical flux
F_num = roe_flux(U_L, U_R, gamma);
% Apply flux (using temporary variables)
Minv = M1inv(:,:,e);
boundary_vector = Minv(:,1);
contribution = F_num * boundary_vector';
rhs(:,:,e) = rhs(:,:,e) - contribution;
else
% Internal left boundary
U_L = [U(1,end,e-1); U(2,end,e-1); U(3,end,e-1)];
U_R = [U(1,1,e); U(2,1,e); U(3,1,e)];
F_num = roe_flux(U_L, U_R, gamma);
% Apply to left element
Minv_left = M1inv(:,:,e-1);
boundary_vector_left = Minv_left(:,end);
contribution_left = F_num * boundary_vector_left';
rhs(:,:,e-1) = rhs(:,:,e-1) + contribution_left;
% Apply to current element
Minv_current = M1inv(:,:,e);
boundary_vector_current = Minv_current(:,1);
contribution_current = F_num * boundary_vector_current';
rhs(:,:,e) = rhs(:,:,e) - contribution_current;
end
% Right boundary
if e == n
% Right boundary condition
rho_R = 0.125;
u_R = 0.0;
p_R = 0.1;
E_R = p_R/(gamma-1) + 0.5*rho_R*u_R^2;
U_R = [rho_R; rho_R*u_R; E_R];
% Current element right state
rho_L = U(1,end,e);
u_L = U(2,end,e)/rho_L;
E_L = U(3,end,e);
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
U_L = [rho_L; rho_L*u_L; E_L];
F_num = roe_flux(U_L, U_R, gamma);
% Apply flux
Minv = M1inv(:,:,e);
boundary_vector = Minv(:,end);
contribution = F_num * boundary_vector';
rhs(:,:,e) = rhs(:,:,e) + contribution;
else
% Internal right boundary
U_L = [U(1,end,e); U(2,end,e); U(3,end,e)];
U_R = [U(1,1,e+1); U(2,1,e+1); U(3,1,e+1)];
F_num = roe_flux(U_L, U_R, gamma);
% Apply to current element
Minv_current = M1inv(:,:,e);
boundary_vector_current = Minv_current(:,end);
contribution_current = F_num * boundary_vector_current';
rhs(:,:,e) = rhs(:,:,e) + contribution_current;
% Apply to right element
Minv_right = M1inv(:,:,e+1);
boundary_vector_right = Minv_right(:,1);
contribution_right = F_num * boundary_vector_right';
rhs(:,:,e+1) = rhs(:,:,e+1) - contribution_right;
end
end
end
%% Roe flux function
function F = roe_flux(U_L, U_R, gamma)
% Unpack variables
rho_L = U_L(1); u_L = U_L(2)/rho_L; E_L = U_L(3);
rho_R = U_R(1); u_R = U_R(2)/rho_R; E_R = U_R(3);
% Compute pressures
p_L = (gamma-1)*(E_L - 0.5*rho_L*u_L^2);
p_R = (gamma-1)*(E_R - 0.5*rho_R*u_R^2);
% Left and right fluxes
F_L = [rho_L*u_L; rho_L*u_L^2 + p_L; u_L*(E_L + p_L)];
F_R = [rho_R*u_R; rho_R*u_R^2 + p_R; u_R*(E_R + p_R)];
% Roe averages
sqrt_rho_L = sqrt(rho_L);
sqrt_rho_R = sqrt(rho_R);
sum_sqrt = sqrt_rho_L + sqrt_rho_R;
u_roe = (sqrt_rho_L*u_L + sqrt_rho_R*u_R)/sum_sqrt;
H_L = (E_L + p_L)/rho_L;
H_R = (E_R + p_R)/rho_R;
H_roe = (sqrt_rho_L*H_L + sqrt_rho_R*H_R)/sum_sqrt;
a_roe = sqrt((gamma-1)*(H_roe - 0.5*u_roe^2));
% Wave strengths
drho = rho_R - rho_L;
du = u_R - u_L;
dp = p_R - p_L;
dU = [drho; du; dp];
lambda = [u_roe-a_roe, u_roe, u_roe+a_roe];
% Entropy fix
epsilon = 0.1*a_roe;
alpha = abs(lambda);
for i = 1:3
if abs(lambda(i)) < epsilon
alpha(i) = (lambda(i)^2 + epsilon^2)/(2*epsilon);
end
end
% Numerical flux
F = 0.5*(F_L + F_R) - 0.5*max(alpha)*dU;
end
%% Gaussian quadrature
function [x, w] = get_gauss_quadrature(n)
beta = 0.5./sqrt(1-(2*(1:n-1)).^(-2));
T = diag(beta,1) + diag(beta,-1);
[V, D] = eig(T);
x = diag(D);
w = 2*V(1,:).^2;
[x, idx] = sort(x);
w = w(idx)';
end
%% Gauss-Lobatto quadrature
function x = get_gauss_lobatto_quadrature(n)
if n == 2
x = [-1; 1];
else
M = zeros(n-2);
for i = 1:n-3
M(i,i+1) = 0.5*sqrt(i*(i+2)/((i+0.5)*(i+1.5)));
M(i+1,i) = M(i,i+1);
end
[V, D] = eig(M);
x = [-1; diag(D); 1];
end
end
%% Lagrange basis evaluation
function [psi, dpsi] = evaluate_lagrange_basis(points_lagrange, points_evaluation)
n = length(points_lagrange);
m = length(points_evaluation);
psi = zeros(n, m);
dpsi = zeros(n, m);
% Compute weights
w = ones(1, n);
for i = 1:n
for j = 1:n
if i ~= j
w(i) = w(i)/(points_lagrange(i) - points_lagrange(j));
end
end
end
% Evaluate basis
for k = 1:m
x = points_evaluation(k);
for i = 1:n
% Compute product
product = 1;
for j = 1:n
if j ~= i
product = product * (x - points_lagrange(j));
end
end
psi(i,k) = w(i)*product;
% Compute derivative
sum_deriv = 0;
for j = 1:n
if j ~= i
term = 1;
for l = 1:n
if l ~= i && l ~= j
term = term * (x - points_lagrange(l));
end
end
sum_deriv = sum_deriv + term;
end
end
dpsi(i,k) = w(i)*sum_deriv;
end
end
end
分析上述代码并进行修正,给出修正后的代码