clc;clear;
% Driver script for solving the 1D Euler equations
% Polynomial order used for approximation
N = 6;
Np = N + 1; % 每元素节点数
Nfaces = 2; % 1D元素面数
Nfp = 1; % 每个面的节点数
NODETOL = 1e-10;
% 定义网格范围
xmin = 0.0;
xmax = 1.0;
K = 250; % 元素数量
% Generate simple mesh
[Nv, VX, K, EToV] = MeshGen1D(xmin, xmax, K);
% 构建参考元素
r = JacobiGL(0, 0, N);
V = Vandermonde1D(N, r);
invV = inv(V);
Dr = Dmatrix1D(N, r, V);
% 创建表面积分项
LIFT = Lift1D(Np, Nfaces, Nfp, V);
% 构建所有节点的坐标
va = EToV(:,1)'; vb = EToV(:,2)';
x = ones(Np,1)*VX(va) + 0.5*(r+1)*(VX(vb)-VX(va));
% 计算几何因子
[rx, J] = GeometricFactors1D(x, Dr);
% 计算边缘节点掩码
fmask1 = find( abs(r+1) < NODETOL)';
fmask2 = find( abs(r-1) < NODETOL)';
Fmask = [fmask1; fmask2]';
Fx = x(Fmask(:), :);
% 构建表面法线和表面处的逆度量
[nx] = Normals1D(Nfp, Nfaces, K);
% 构建连接矩阵
[EToE, EToF] = Connect1D(EToV);
% 构建连接映射
[vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO] = BuildMaps1D(EToE, EToF, x, Np, Nfp, Nfaces, NODETOL, K, Fmask); % 修改返回值
% 计算 Fscale
Fscale = 1./(J(Fmask,:));
% 设置初始条件 -- Sod's problem
gamma = 1.4;
MassMatrix = inv(V')/V;
cx = ones(Np,1)*sum(MassMatrix*x,1)/2;
rho = ones(Np,K).*( (cx<0.5) + 0.125*(cx>=0.5));
rhou = zeros(Np,K);
Ener = ones(Np,K).*((cx<0.5) + 0.1*(cx>=0.5))/(gamma-1.0);
FinalTime = 0.2;
% Solve Problem
[rho,rhou,Ener] = Euler1D(rho, rhou, Ener, FinalTime, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, V, invV, Nfp, Nfaces); % 添加 x 参数
figure; plot(x, Ener); drawnow; pause(0.1);
figure;
subplot(3,1,1); plot(x, rho); title('Density');
subplot(3,1,2); plot(x, rhou./rho); title('Velocity');
subplot(3,1,3); plot(x, (gamma-1)*(Ener-0.5*rhou.^2./rho)); title('Pressure');
function [Nv, VX, K, EToV] = MeshGen1D(xmin, xmax, K)
Nv = K+1;
% Generate node coordinates
VX = (1:Nv);
for i = 1:Nv
VX(i) = (xmax-xmin)*(i-1)/(Nv-1) + xmin;
end
% read element to node connectivity
EToV = zeros(K, 2);
for k = 1:K
EToV(k,1) = k; EToV(k,2) = k+1;
end
end
function [x] = JacobiGL(alpha, beta, N)
x = zeros(N+1,1);
if (N==1)
x(1)=-1.0;
x(2)=1.0;
return;
end
[xint,w] = JacobiGQ(alpha+1, beta+1, N-2);
x = [-1, xint', 1]';
end
function [V1D] = Vandermonde1D(N, r)
V1D = zeros(length(r), N+1);
for j=1:N+1
V1D(:,j) = JacobiP(r(:), 0, 0, j-1);
end
end
function [Dr] = Dmatrix1D(N, r, V)
Vr = GradVandermonde1D(N, r);
Dr = Vr/V;
end
function [LIFT] = Lift1D(Np, Nfaces, Nfp, V)
Emat = zeros(Np, Nfaces*Nfp);
% Define Emat
Emat(1,1) = 1.0; Emat(Np,2) = 1.0;
% inv(mass matrix)*\s_n (L_i,L_j)_{edge_n}
LIFT = V*(V'*Emat);
end
function [rx, J] = GeometricFactors1D(x, Dr)
xr = Dr*x; J = xr; rx = 1./J;
end
function [nx] = Normals1D(Nfp, Nfaces, K)
nx = zeros(Nfp*Nfaces, K);
% Define outward normals
nx(1, :) = -1.0; nx(2, :) = 1.0;
end
function [EToE, EToF] = Connect1D(EToV)
Nfaces = 2;
% Find number of elements and vertices
K = size(EToV,1); TotalFaces = Nfaces*K; Nv = K+1;
% List of local face to local vertex connections
vn = [1,2];
% Build global face to node sparse array
SpFToV = spalloc(TotalFaces, Nv, 2*TotalFaces);
sk = 1;
for k=1:K
for face=1:Nfaces
SpFToV( sk, EToV(k, vn(face))) = 1;
sk = sk+1;
end
end
% Build global face to global face sparse array
SpFToF = SpFToV*SpFToV' - speye(TotalFaces);
% Find complete face to face connections
[faces1, faces2] = find(SpFToF==1);
% Convert face global number to element and face numbers
element1 = floor( (faces1-1)/Nfaces ) + 1;
face1 = mod( (faces1-1), Nfaces ) + 1;
element2 = floor( (faces2-1)/Nfaces ) + 1;
face2 = mod( (faces2-1), Nfaces ) + 1;
% Rearrange into Nelements x Nfaces sized arrays
ind = sub2ind([K, Nfaces], element1, face1);
EToE = (1:K)'*ones(1,Nfaces);
EToF = ones(K,1)*(1:Nfaces);
EToE(ind) = element2; EToF(ind) = face2;
end
function [vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO] = BuildMaps1D(EToE, EToF, x, Np, Nfp, Nfaces, NODETOL, K, Fmask) % 修改函数定义
% number volume nodes consecutively
nodeids = reshape(1:K*Np, Np, K);
vmapM = zeros(Nfp, Nfaces, K);
vmapP = zeros(Nfp, Nfaces, K);
for k1=1:K
for f1=1:Nfaces
% find index of face nodes with respect to volume node ordering
vmapM(:,f1,k1) = nodeids(Fmask(:,f1), k1);
end
end
for k1=1:K
for f1=1:Nfaces
% find neighbor
k2 = EToE(k1,f1); f2 = EToF(k1,f1);
% find volume node numbers of left and right nodes
vidM = vmapM(:,f1,k1); vidP = vmapM(:,f2,k2);
x1 = x(vidM); x2 = x(vidP);
% Compute distance matrix
D = (x1 -x2 ).^2;
if (D<NODETOL)
vmapP(:,f1,k1) = vidP;
end
end
end
vmapP = vmapP(:); vmapM = vmapM(:);
% Create list of boundary nodes
mapB = find(vmapP==vmapM); vmapB = vmapM(mapB);
% Create specific left (inflow) and right (outflow) maps
mapI = 1; mapO = K*Nfaces; vmapI = 1; vmapO = K*Np; % 定义 vmapI, mapI, vmapO, mapO
end
function [rho, rhou, Ener] = Euler1D(rho, rhou, Ener, FinalTime, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, V, invV, Nfp, Nfaces) % 添加 x 参数
% Parameters
CFL = 1.0; time = 0;
% Prepare for adaptive time stepping
mindx = min(x(2,:)-x(1,:));
% Limit initial solution
rho = SlopeLimitN(rho, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
rhou = SlopeLimitN(rhou, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
Ener = SlopeLimitN(Ener, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
% outer time step loop
while(time < FinalTime)
Temp = (Ener - 0.5*(rhou).^2./rho)./rho;
cvel = sqrt(gamma*(gamma-1)*Temp);
dt = CFL*min(min(mindx./(abs(rhou./rho)+cvel)));
if(time+dt > FinalTime)
dt = FinalTime-time;
end
% 3rd order SSP Runge-Kutta
% SSP RK Stage 1.
[rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho, rhou, Ener, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添加 vmapI, mapI, vmapO, mapO
rho1 = rho + dt*rhsrho;
rhou1 = rhou + dt*rhsrhou;
Ener1 = Ener + dt*rhsEner;
% Limit fields
rho1 = SlopeLimitN(rho1, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
rhou1 = SlopeLimitN(rhou1, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
Ener1 = SlopeLimitN(Ener1, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
% SSP RK Stage 2.
[rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho1, rhou1, Ener1, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添加 vmapI, mapI, vmapO, mapO
rho2 = (3*rho + rho1 + dt*rhsrho )/4;
rhou2 = (3*rhou + rhou1 + dt*rhsrhou)/4;
Ener2 = (3*Ener + Ener1 + dt*rhsEner)/4;
% Limit fields
rho2 = SlopeLimitN(rho2, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
rhou2 = SlopeLimitN(rhou2, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
Ener2 = SlopeLimitN(Ener2, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
% SSP RK Stage 3.
[rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho2, rhou2, Ener2, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces); % 添加 vmapI, mapI, vmapO, mapO
rho = (rho + 2*rho2 + 2*dt*rhsrho )/3;
rhou = (rhou + 2*rhou2 + 2*dt*rhsrhou)/3;
Ener = (Ener + 2*Ener2 + 2*dt*rhsEner)/3;
% Limit solution
rho = SlopeLimitN(rho, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
rhou = SlopeLimitN(rhou, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
Ener = SlopeLimitN(Ener, Np, V, invV, K, x, Dr); % 添加 x 和 Dr 参数
% Increment time and adapt timestep
time = time + dt;
end
end
function [rhsrho, rhsrhou, rhsEner] = EulerRHS1D(rho, rhou, Ener, N, Np, K, r, x, VX, Dr, LIFT, nx, Fx, Fscale, vmapM, vmapP, vmapB, mapB, vmapI, mapI, vmapO, mapO, Fmask, rx, J, gamma, Nfp, Nfaces) % 修改函数定义
% compute maximum velocity for LF flux
pres = (gamma-1.0)*(Ener - 0.5*(rhou).^2./rho);
cvel = sqrt(gamma*pres./rho); lm = abs(rhou./rho)+cvel;
% Compute fluxes
rhof = rhou; rhouf = rhou.^2./rho + pres; Enerf = (Ener + pres).*rhou./rho;
% Compute jumps at internal faces
drho = zeros(Nfp*Nfaces, K); drho(:) = rho(vmapM) - rho(vmapP);
drhou = zeros(Nfp*Nfaces, K); drhou(:) = rhou(vmapM) - rhou(vmapP);
dEner = zeros(Nfp*Nfaces, K); dEner(:) = Ener(vmapM) - Ener(vmapP);
drhof = zeros(Nfp*Nfaces, K); drhof(:) = rhof(vmapM) - rhof(vmapP);
drhouf = zeros(Nfp*Nfaces, K); drhouf(:) = rhouf(vmapM) - rhouf(vmapP);
dEnerf = zeros(Nfp*Nfaces, K); dEnerf(:) = Enerf(vmapM) - Enerf(vmapP);
LFc = zeros(Nfp*Nfaces, K); LFc(:) = max(lm(vmapP), lm(vmapM));
% Compute fluxes at interfaces
drhof(:) = nx(:).*drhof(:)/2.0 - LFc(:)/2.0.*drho(:);
drhouf(:) = nx(:).*drhouf(:)/2.0 - LFc(:)/2.0.*drhou(:);
dEnerf(:) = nx(:).*dEnerf(:)/2.0 - LFc(:)/2.0.*dEner(:);
% Boundary conditions for Sod's problem
rhoin = 1.000; rhouin = 0.0;
pin = 1.000; Enerin = pin/(gamma-1.0);
rhoout = 0.125; rhouout = 0.0;
pout = 0.100; Enerout = pout/(gamma-1.0);
% Set fluxes at inflow/outflow
rhofin = rhouin; rhoufin = rhouin.^2./rhoin + pin;
Enerfin = (pin/(gamma-1.0) + 0.5*rhouin^2/rhoin + pin).*rhouin./rhoin;
lmI = lm(vmapI)/2; nxI = nx(mapI); % 使用 vmapI 和 mapI
drhof (mapI) = nxI*(rhof (vmapI) - rhofin )/2.0 - lmI*(rho(vmapI) - rhoin);
drhouf(mapI) = nxI*(rhouf(vmapI) - rhoufin)/2.0 - lmI*(rhou(vmapI) - rhouin);
dEnerf(mapI) = nxI*(Enerf(vmapI) - Enerfin)/2.0 - lmI*(Ener(vmapI) - Enerin);
rhofout = rhouout; rhoufout = rhouout.^2./rhoout + pout;
Enerfout = (pout/(gamma-1.0) + 0.5*rhouout^2/rhoout + pout).*rhouout./rhoout;
lmO = lm(vmapO)/2; nxO = nx(mapO); % 使用 vmapO 和 mapO
drhof (mapO) = nxO*(rhof(vmapO) - rhofout)/2.0 - lmO*(rho (vmapO) - rhoout);
drhouf(mapO) = nxO*(rhouf(vmapO) - rhoufout)/2.0 - lmO*(rhou(vmapO) - rhouout);
dEnerf(mapO) = nxO*(Enerf(vmapO) - Enerfout)/2.0 - lmO*(Ener(vmapO) - Enerout);
% compute right hand sides of the PDE's
rhsrho = -rx.*(Dr*rhof) + LIFT*(Fscale.*drhof);
rhsrhou = -rx.*(Dr*rhouf) + LIFT*(Fscale.*drhouf);
rhsEner = -rx.*(Dr*Enerf) + LIFT*(Fscale.*dEnerf);
end
function ulimit = SlopeLimitN(u, Np, V, invV, K, x, Dr) % 添加 x 和 Dr 参数
% Compute cell averages
uh = invV*u; uh(2:Np,:)=0; uavg = V*uh; v = uavg(1,:);
% Apply slope limiter as needed.
ulimit = u; eps0 = 1.0e-8;
% find end values of each element
ue1 = u(1,:); ue2 = u(end,:);
% find cell averages
vk = v; vkm1 = [v(1), v(1:K-1)]; vkp1 = [v(2:K), v(K)];
% Apply reconstruction to find elements in need of limiting
ve1 = vk - minmod([(vk-ue1); vk-vkm1; vkp1-vk]);
ve2 = vk + minmod([(ue2-vk); vk-vkm1; vkp1-vk]);
ids = find(abs(ve1-ue1) > eps0 | abs(ve2-ue2) > eps0);
% Check to see if any elements require limiting
if(~isempty(ids))
% create piecewise linear solution for limiting on specified elements
uhl = invV*u(:,ids); uhl(3:Np,:)=0; ul = V*uhl;
% apply slope limiter to selected elements
ulimit(:,ids) = SlopeLimitLin(ul, x(:,ids), vkm1(ids), vk(ids), vkp1(ids), Np, Dr); % 添加 Dr 参数
end
end
function ulimit = SlopeLimitLin(ul, xl, vm1, v0, vp1, Np, Dr) % 添加 Dr 参数
% Compute various geometric measures
ulimit = ul; h = xl(Np,:) - xl(1,:);
x0 = ones(Np,1)*(xl(1,:) + h/2);
hN = ones(Np,1)*h;
% Limit function
ux = (2./hN).*(Dr*ul); % 使用传入的 Dr
ulimit = ones(Np,1)*v0 + (xl-x0).*(ones(Np,1)*minmod([ux(1,:); (vp1-v0)./h; (v0-vm1)./h]));
end
function mfunc = minmod(v)
m = size(v,1); mfunc = zeros(1, size(v,2));
s = sum(sign(v),1)/m;
ids = find(abs(s) == 1);
if(~isempty(ids))
mfunc(ids) = s(ids).*min(abs(v(:,ids)), [], 1);
end
end
function [x, w] = JacobiGQ(alpha, beta, N)
if (N == 0)
x(1) = -(alpha-beta)/(alpha+beta+2);
w(1) = 2;
return;
end
% Form symmetric matrix from recurrence.
J = zeros(N+1);
h1 = 2*(0:N) + alpha + beta;
J = diag(-1/2*(alpha^2-beta^2)./(h1+2)./h1) + ...
diag(2./(h1(1:N)+2).*sqrt((1:N).*((1:N)+alpha+beta).*...
((1:N)+alpha).*((1:N)+beta)./(h1(1:N)+1)./(h1(1:N)+3)),1);
if (alpha+beta < 10*eps)
J(1,1) = 0.0;
end
J = J + J';
% Compute quadrature by eigenvalue solve
[V, D] = eig(J);
x = diag(D);
w = (V(1,:)').^2 * 2^(alpha+beta+1)/(alpha+beta+1) * gamma(alpha+1) * ...
gamma(beta+1)/gamma(alpha+beta+1);
end
function [P] = JacobiP(x, alpha, beta, N)
xp = x;
dims = size(xp);
if (dims(2) == 1)
xp = xp';
end
PL = zeros(N+1, length(xp));
% Initial values P_0(x) and P_1(x)
gamma0 = 2^(alpha+beta+1)/(alpha+beta+1) * gamma(alpha+1) * ...
gamma(beta+1)/gamma(alpha+beta+1);
PL(1,:) = 1.0/sqrt(gamma0);
if (N == 0)
P = PL';
return;
end
gamma1 = (alpha+1)*(beta+1)/(alpha+beta+3)*gamma0;
PL(2,:) = ((alpha+beta+2)*xp/2 + (alpha-beta)/2)/sqrt(gamma1);
if (N == 1)
P = PL(N+1,:)';
return;
end
% Repeat value in recurrence.
aold = 2/(2+alpha+beta)*sqrt((alpha+1)*(beta+1)/(alpha+beta+3));
% Forward recurrence using the symmetry of the recurrence.
for i=1:N-1
h1 = 2*i + alpha + beta;
anew = 2/(h1+2)*sqrt( (i+1)*(i+1+alpha+beta)*(i+1+alpha)*...
(i+1+beta)/(h1+1)/(h1+3));
bnew = - (alpha^2-beta^2)/h1/(h1+2);
PL(i+2,:) = 1/anew*( -aold*PL(i,:) + (xp-bnew).*PL(i+1,:));
aold = anew;
end
P = PL(N+1,:)';
end
function [DVr] = GradVandermonde1D(N, r)
DVr = zeros(length(r), (N+1));
% Initialize matrix
for i=0:N
[DVr(:,i+1)] = GradJacobiP(r(:), 0, 0, i);
end
end
function [dP] = GradJacobiP(r, alpha, beta, N)
dP = zeros(length(r), 1);
if(N == 0)
dP(:,:) = 0.0;
else
dP = sqrt(N*(N+alpha+beta+1))*JacobiP(r(:), alpha+1, beta+1, N-1);
end
end
对上述代码分析,给出相对应的数学理论公式