albada.m
代码语言：javascript复制figure(1), clf, hold on
x = 0.0:0.01:4.0; y = (x.*x x)./(1 x.*x); % Graph of van Albada limiter
plot(x,y)
x = 0.0:0.01:1.1; y = 0.5*(1 sqrt(2))*x; plot(x,y)
x = 0.0:0.01:4.0; y = 0.5*(1 sqrt(2))*ones(size(x)); plot(x,y)
y = ones(size(x)); plot(x,y)
x = 0.0:0.01:1.3; plot(x,x)
y = 0.0:0.01:1; x = ones(size(y)); plot(x,y,'--')
extrap.m
代码语言：javascript复制function[zL,zR] = extrap(z, limtype); % Slope-limited extrapolation
% Function called: limiter
J = length(z);
zL(1) = z(1);
for j = 2:J-1
b = (z(j 1) - z(j)); a = (z(j) - z(j-1));
zL(j) = z(j) 0.5*limiter(a, b, limtype)*(z(j 1)-z(j));
zR(j-1) = z(j) 0.5*limiter(b, a, limtype)*(z(j-1)-z(j));
end
zR(J-1) = z(J);
limiter.m
代码语言：javascript复制function y = limiter(a,b,limtype)
% Limiter function as defined in Sect. 4.8 and as used in Sect. 10.8
if b==0, y=0; else
if limtype == 2 % van Albada
y = (a^2 a*b)/(a^2 b^2);
elseif limtype == 1 % Minmod
y = max(0, min(a/b,1));
else % No limiting, no MUSCL
y = 0;
end
end
Liou_Steffen_MUSCL_scheme.m
代码语言：javascript复制clear all
global PRL CRL MACHLEFT gamma pleft pright rholeft rhoright uleft...
uright tend epsi lambda % lambda = dt/dx
% .....................Input............................
gamma = 1.4; % Ratio of specific heats
J = 48; % Number of grid cells
limtype = 2; % Enter 1 for minmod or 2 for van Albada or something else
% if MUSCL not wanted
gammab = 1/(gamma - 1); gam1 = gamma-1; gamgam = gamma^gamma;
problem_specification
h = 1/J; % Cell size
dt = lambda*h; % Time step
n = floor(tend/dt); % Number of time-steps
% Definition of grid numbering
% x=0 x=1
% grid |---o---|---o---|---o--- ... --|---o---|
% 1 1 2 2 3 J-1 J
xcenter = h*[1:J] - h/2; % Location of cell centers
press = zeros(size(xcenter)); % Preallocation of pressure,
rhoold = press; uold = press; % density, velocity,
rhonew = press; mnew = press; % momentum,
totenew = press; % total energy
for j = 1:length(xcenter) % Initial conditions
if xcenter(j) < 0.5, press(j) = pleft; rhoold(j) = rholeft; uold(j) = uleft;
else, press(j) = pright; rhoold(j) = rhoright; uold(j) = uright;
end
end
% toten = rho*(total energy)
totenold = rhoold.*(0.5*uold.*uold gammab*press./rhoold);
totenleft = totenold(1); totenright = totenold(J);
mold = rhoold.*uold; % m is momentum
c = sqrt(gamma*press./rhoold); % Sound speed
mach = uold./c;
enthalpy = 0.5*uold.*uold gammab*c.^2;
% Preallocation of auxiliary variables
machplus = mach; machminus = mach;
presplus = mach; presminus = mach;
% Preallocation of Liou-Steffen fluxes
flux1 = zeros(J-1,1); flux2 = flux1; flux3 = flux1;
% Preallocation of extrapolated states
Z1L = zeros(J-1,1); Z2L = Z1L; Z3L = Z1L;
Z1R = Z1L; Z2R = Z1L; Z3R = Z1L;
U1L = Z1L; U1R = Z1L; U2L = Z1L; U2R= Z1L;
U3R = Z1L; U3L = Z1L;
flopcount = flops; % Operations counter
t = 0;
for i = 1:n, t = t dt;
rhostar = rhoold; mstar = mold; totenstar = totenold;
rkalpha = 0.25; listeeulerstep
rkalpha = 1/3; listeeulerstep
rkalpha = 0.5; listeeulerstep
rkalpha = 1; listeeulerstep
rhonew = rhostar; mnew = mstar; totenew = totenstar;
uold = mnew./rhonew; press = gam1*(totenew - 0.5*mnew.*uold);
rhoold = rhonew; totenold = totenew; mold = mnew;
c = sqrt(gamma*press./rhoold); mach = uold./c;
end
flopcount = flops - flopcount;
entropy = log(press./rhoold.^gamma);
figure(1), clf
subplot(2,3,1),hold on,title('DENSITY','fontsize',14),plot(xcenter,rhonew,'o')
subplot(2,3,2),hold on,title('VELOCITY','fontsize',14),plot(xcenter,uold,'o')
subplot(2,3,3),hold on,title('PRESSURE','fontsize',14),plot(xcenter,press,'o')
subplot(2,3,4),hold on,title('MACHNUMBER','fontsize',14),plot(xcenter,mach,'o')
subplot(2,3,5),hold on,title('ENTROPY','fontsize',14),plot(xcenter,entropy,'o')
subplot(2,3,6), axis('off'), hold on, title('Liou-Steffen scheme','fontsize',14)
text(0,0.9,['lambda = ', num2str(lambda),' t = ',num2str(n*dt)])
text(0,0.75,' Primitive extrapolation'), text(0,0.6,'SHK Runge-Kutta')
if limtype == 0, s1 = 'No MUSCL';
elseif limtype == 1, s1 = 'MUSCL, minmod limiter';
else, s1 = 'MUSCL, van Albada limiter';
end
text(0,0.45,s1)
text(0,0.3,['flopcount = ',num2str(flopcount)])
Riemann % Plot exact solution
