clear all;
close all;
%Input Data
%--------------------------ONLY CHANGE THESE--------------------------%
nelem=21; %Number of Elements
nop=6; %Polynomial Order
time_final=1.5; %final time in revolutions
integration_points=1; %=1 for LGL and =2 for LG
integration_type=1; %=1 is inexact and =2 is exact
space_method='dg'; %cgc=CGc; cgd=CGd; dg=DG
%--------------------------ONLY CHANGE THESE--------------------------%
%For nop=2,4,6 works upto Time=10 (nop*nelem=120) dt=0.001.
kstages=4; %RK1, RK2, RK3, or RK4
dt=0.001;
nplots=10; %plotting variable - Number of Frames ploted
iplot_movie=0;
store_movie=0;
iplot_figures=1;
iplot_elements=1;
iplot_modes=1;
%alpha=2.0/3.0; %1=conservative, 0=non-conservative
alpha=1; %1=conservative, 0=non-conservative
limit=0; %=0 no limiting, =1 Krivodonova Limiter, =2 Bound-Preserving Limiter
icase=1; %=1 is a Gaussian with flat bottom;
xmu=0; %filtering strength: 1 is full strength and 0 is no filter
ifilter=0; %time-step frequency that the filter is applied.
filter_type=2; %=1 is Modal Hierarchical (no need for DSS)
%and =2 is regular Legendre (better for DG)
diss=1; %=1 dissipation (Rusanov Flux) and =0 no dissipation (Central Flux)
% visc=6.6e-4;
visc=0.0001;
tau=0e3;
flux_method='rusanov'; %rusanov,roe, or energy
LDG_alpha=0.5; LDG_beta=1.0-LDG_alpha;
%Viscous Operators
if (strcmp(space_method,'dg')>0)
elliptic_method='LDG';
else
elliptic_method='SIP';
end
%Store Constants
eps=1e-15;
ntime=time_final/dt;
iplot=round(ntime/nplots);
%Store Constants
ngl=nop + 1;
method_text = [space_method];
%Compute i,e -> I pointer
I=0;
Ipointer=zeros(ngl,nelem);
for e=1:nelem
for i=1:ngl
I=I+1;
Ipointer(i,e)=I;
end
end
%Compute Interpolation and Integration Points
[xgl,wgl]=legendre_gauss_lobatto(ngl);
if (integration_points == 1)
integration_text=['LGL'];
if (integration_type == 1)
noq=nop;
elseif (integration_type == 2)
noq=3*nop/2;
end
nq=noq + 1;
[xnq,wnq]=legendre_gauss_lobatto(nq);
elseif (integration_points == 2)
integration_text=['LG'];
noq=nop;
nq=noq + 1;
[xnq,wnq]=legendre_gauss(nq);
end
main_text=[space_method ':' integration_text ]; %连接两个字符串
%Compute Lagrange Polynomial and derivatives
[psi,dpsi] = lagrange_basis3(ngl,nq,xgl,xnq);
%Compute Filter Matrix
[f,vdm,vdm_inv] = filter_init(ngl,xgl,xmu,filter_type);
%Create Grid
[coord,intma_cg]=create_grid_dg(ngl,nelem,xgl,icase);
dx_min=coord(2,1)-coord(1,1);
%Form Global Matrix Pointer
npoin_cg=nop*nelem + 1;
npoin_dg=ngl*nelem;
intma=zeros(ngl,nelem);
if ( strcmp(space_method,'cgc') > 0 )
npoin=npoin_cg;
intma=intma_cg;
else %cgd or dg
npoin=npoin_dg;
ip=0;
for e=1:nelem
for i=1:ngl
ip=ip+1;
intma(i,e)=ip;
end
end
end
if ( strcmp(space_method,'dg') > 0 )
npoin_cg=npoin;
intma_cg=intma;
end
%Compute Exact Solution
time=0;
p_movie=zeros(npoin,nplots);
u_movie=zeros(npoin,nplots);
time_movie=zeros(nplots,1);
mass_movie=zeros(nplots,1);
energy_movie=zeros(nplots,1);
[qe] = exact_solution_dg(intma,coord,npoin,nelem,ngl,time,icase);
q_init=qe;
%Compute Initial Mass and Energy
[mass0,energy0] = compute_Mass_and_Energy(qe,intma,coord,wnq,psi,nelem,ngl,nq);
%Create Periodic BC Pointer
iperiodic=zeros(npoin,1);
for i=1:npoin
iperiodic(i)=i;
end
%Create Local/Element Mass
Mmatrix = create_Mmatrix(intma_cg,coord,npoin_cg,nelem,ngl,nq,wnq,psi);
%Initialize State Vector
qp=qe; q0=qe; q1=qe;
iframe=0;
%Initialize Arrays
rhs=zeros(npoin,2);
rhs_t=zeros(ngl,1);
rhs_rk4=zeros(npoin,2,4);
visc_elem=zeros(npoin,2);
%Compute RK Time-Integration Coefficients
[a0,a1,beta] = compute_ti_coefficients(kstages);
%Time Integration
for itime=1:ntime
time=time + dt;
u_max=-100000;
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
u_max=max(u_max, abs(qp(ip,1)) );
end
end
courant=u_max*dt/dx_min;
%disp(['itime = ',num2str(itime),' time = ', num2str(time),' courant = ', num2str(courant)]);
for ik=1:kstages
%Create RHS vector for Inviscid Operators
%-----------------------------------------------------------------------------%
%-------------------Student add your RHS functions here ----------------------%
%rhs = create_rhs_volume(qp,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,alpha);
%rhs = create_rhs_flux(rhs,qp,intma,nelem,ngl,diss,flux_method);
%-----------------------------------------------------------------------------%
%-----------------------------------------------------------------------------%
%Create RHS vector for Viscous Operators
if (strcmp(elliptic_method,'SIP') > 0 )
rhs = create_Lmatrix_SIPDG(rhs,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,qp,tau,visc);
elseif (strcmp(elliptic_method,'LDG') > 0 )
rhs = create_Lmatrix_LDG(rhs,Mmatrix,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,qp,LDG_alpha,LDG_beta,visc);
end
%Apply Communicator
%-----------------------------------------------------------------------------%
%-------------------Student add your RHS functions here ----------------------%
%rhs = create_rhs_communicator(rhs,space_method,intma,intma_cg,coord,Mmatrix,npoin,npoin_cg,nelem,ngl,nq,wnq,psi,dpsi,alpha,0);
%-----------------------------------------------------------------------------%
%-------------------Student add your RHS functions here ----------------------%
%Solve System
qp=a0(ik)*q0 + a1(ik)*q1 + dt*beta(ik)*rhs;
%Update
rhs_rk4(:,:,ik)=rhs(:,:);
q1=qp;
end %ik
%Do RK4 Solution
if (kstages == 4)
qp(:,:)=q0(:,:) + dt/6.0*( rhs_rk4(:,:,1) + 2*rhs_rk4(:,:,2)...
+ 2*rhs_rk4(:,:,3) + rhs_rk4(:,:,4) );
end
%Filter Solution
if (mod(itime,ifilter) == 0 && ifilter > 0)
qp = apply_filter_dg(qp,f,intma,nelem,ngl);
end
%Limit Solution
limit_element=zeros(nelem,1);
if limit == 0 || limit == 1
[limit_element,qp,qmodal] = limiter_nodal_v3(qp,vdm,vdm_inv,intma,nelem,ngl,limit);
elseif limit == 2
[limit_element,qp,qmodal] = limiter_nodal_bp(qp,q0,vdm,vdm_inv,intma,nelem,ngl);
end
%DSS Solution for CG after Limiting
if (limit > 0 || ( mod(itime,ifilter) == 0 && ifilter > 0 ) )
if ( strcmp(space_method,'cgd') > 0 || strcmp(space_method,'cgc') > 0) %CG
qp = apply_dss(qp,intma,intma_cg,coord,Mmatrix,npoin,npoin_cg,nelem,ngl,nq,wnq,psi,dpsi,alpha,0);
end
end
%Plot Modes
if (limit >= 0 && iplot_modes == 1 && mod(itime,iplot) == 0)
[test] = plot_modes(qmodal,Ipointer,limit_element,intma,coord,qp,nelem,ngl,nq,nop,noq,time,diss,main_text);
end
%Update Q
q0=qp;
qb=q0;
qb=0;
if (iplot_movie == 1 && mod(itime,iplot) == 0)
[p_movie,u_movie,time_movie,mass_movie,energy_movie,L2_norm,iframe] = store_movie_frames(qp,q0,intma,coord,wnq,psi,npoin,nelem,ngl,nq,p_movie,u_movie,time_movie,mass_movie,energy_movie,iframe,mass0,energy0,time,icase);
end
end %itime
if (iplot_movie == 1)
hmax=max(p_movie(:));
hmin=min(p_movie(:));
xmax=max(max(coord));
xmin=min(min(coord));
figure;
for i=1:iframe
subplot(3,1,1); %H Solution
% hmax=max(p_movie(:,i));
% hmin=min(p_movie(:,i));
% xmax=max(max(coord));
% xmin=min(min(coord));
for e=1:nelem
for j=1:ngl
jp=intma(j,e);
x(j)=coord(j,e);
y(j)=p_movie(jp,i);
end
plot_handle=plot(x,y,'r-','LineWidth',2);
hold on;
end
axis([xmin xmax hmin hmax]);
title_text=[main_text ', Ne = ' num2str(nelem) ', N = ' num2str(nop) ', Q = ' num2str(noq) ', Time = ' num2str(time_movie(i))];
title([title_text],'FontSize',18);
hold on;
%Plot Elements
for e=1:nelem
i1=intma(1,e);
i2=intma(ngl,e);
x1=coord(1,e);
x2=coord(ngl,e);
y1=p_movie(i1,i);
y2=p_movie(i2,i);
plot(x1,y1,'ko');
plot(x2,y2,'ko');
end
title([title_text],'FontSize',18);
xlabel('x','FontSize',18);
ylabel('h','FontSize',18);
set(gca, 'FontSize', 18);
hold off;
subplot(3,1,2); %L2 Norms
ymax=max(max(L2_norm));
ymin=min(min(L2_norm));
tmax=max(time_movie);
tmin=min(time_movie);
tt=time_movie(1:i);
h_norm=L2_norm(1,1:i);
plot_handle=plot(tt,h_norm,'b-','LineWidth',2);
u_norm=L2_norm(2,1:i);
hold on;
plot_handle=plot(tt,u_norm,'k-','LineWidth',2);
xlabel('Time','FontSize',18);
ylabel('L^2 Norm','FontSize',18);
set(gca, 'FontSize', 18);
legend('||h||_2','||u||_2');
axis([ tmin tmax ymin ymax]);
hold on;
subplot(3,1,3); %Mass Conservation
ymax=max(mass_movie);
ymin=min(mass_movie);
tmax=max(time_movie);
tmin=min(time_movie);
tt=time_movie(1:i);
yt=mass_movie(1:i);
plot_handle=plot(tt,yt,'r-','LineWidth',2);
xlabel('Time','FontSize',18);
ylabel('\Delta M','FontSize',18);
set(gca, 'FontSize', 18);
% axis([ tmin tmax ymin ymax]);
hold on;
tt=time_movie(1:i);
yt=energy_movie(1:i);
plot_handle=plot(tt,yt,'b-','LineWidth',2);
legend('\Delta M','\Delta E');
hold on;
M_i=getframe(gcf);
M(i)=M_i;
pause(0.01);
hold off;
end
if store_movie == 1
file_movie=[main_text '_e' num2str(nelem) '_p' num2str(nop) '_q' num2str(noq) '_diss' num2str(diss) '.avi'];
movie2avi(M,file_movie,'fps',5);
end
end
if (iplot_figures)
%Exact Solution
[qe] = exact_solution_dg(intma,coord,npoin,nelem,ngl,time,icase);
figure; %H Plot
hmax=max(qp(:,1));
hmin=min(qp(:,1));
xmax=max(coord(:));
xmin=min(coord(:));
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
x(i)=coord(i,e);
y(i)=qp(ip,1);
ye(i)=q_init(ip,1);
end
plot_handle=plot(x,y,'r-','LineWidth',2);
hold on;
plot_handle=plot(x,ye,'k:','LineWidth',2);
end
axis([xmin xmax hmin hmax]);
title_text=[main_text ' visc = ' num2str(visc) ', Ne = ' num2str(nelem) ', N = ' num2str(nop) ', Q = ' num2str(noq) ', Time = ' num2str(time)];
title([title_text],'FontSize',18);
hold on;
%Plot Elements
if (iplot_elements == 1)
for e=1:nelem
i1=intma(1,e);
i2=intma(ngl,e);
x1=coord(1,e);
x2=coord(ngl,e);
y1=qp(i1,1);
y2=qp(i2,1);
plot(x1,y1,'ko');
plot(x2,y2,'ko');
end
end
title([title_text],'FontSize',18);
xlabel('x','FontSize',18);
ylabel('u(x,t)','FontSize',18);
set(gca, 'FontSize', 18);
hold off;
end
%%
function [xgl,wgl] = legendre_gauss_lobatto(P)
p=P-1;
ph=floor( (p+1)/2 );
for i=1:ph
x=cos( (2*i-1)*pi/(2*p+1) );
for k=1:20
[L0,L0_1,L0_2]=legendre_poly(p,x);
%Get new Newton Iteration
dx=-(1-x^2)*L0_1/(-2*x*L0_1 + (1-x^2)*L0_2);
x=x+dx;
if (abs(dx) < 1.0e-20)
break
end
end
xgl(p+2-i)=x;
wgl(p+2-i)=2/(p*(p+1)*L0^2);
end
%Check for Zero Root
if (p+1 ~= 2*ph)
x=0;
[L0,L0_1,L0_2]=legendre_poly(p,x);
xgl(ph+1)=x;
wgl(ph+1)=2/(p*(p+1)*L0^2);
end
%Find remainder of roots via symmetry
for i=1:ph
xgl(i)=-xgl(p+2-i);
wgl(i)=+wgl(p+2-i);
end
end
%%
function [L0,L0_1,L0_2] = legendre_poly(p,x)
L1=0;L1_1=0;L1_2=0;
L0=1;L0_1=0;L0_2=0;
for i=1:p
L2=L1;L2_1=L1_1;L2_2=L1_2;
L1=L0;L1_1=L0_1;L1_2=L0_2;
a=(2*i-1)/i;
b=(i-1)/i;
L0=a*x*L1 - b*L2;
L0_1=a*(L1 + x*L1_1) - b*L2_1;
L0_2=a*(2*L1_1 + x*L1_2) - b*L2_2;
end
end
%%
function [psi,dpsi] = lagrange_basis3(P,Q,xgl,xnq)
%Perform Quadrature
for l=1:Q
xl=xnq(l);
%Construct Basis
for i=1:P
xi=xgl(i);
psi(i,l)=1;
dpsi(i,l)=0;
for j=1:P
xj=xgl(j);
%Basis
if (j ~= i)
psi(i,l)=psi(i,l)*(xl-xj)/(xi-xj);
end
ddpsi=1;
if (j ~= i)
for k=1:P
xk=xgl(k);
%Derivative of Basis
if (k ~=i && k ~=j)
ddpsi=ddpsi*(xl-xk)/(xi-xk);
end
end
dpsi(i,l)=dpsi(i,l) + ddpsi/(xi-xj);
end
end
end
end
end
%%
function [f,vdm,vdm_inv] = filter_init(P,xgl,xmu,filter_type)
%Constants
p=P-1;
ph=floor( (p+1)/2 );
alpha=17;
order=18;
%Initialize
leg=zeros(P,P);
leg2=zeros(P,P);
f=zeros(P,P);
%Compute Legendre Polynomial Matrix
for i=1:P
x=xgl(i);
for j=1:P
jj=j-1;
[L0,L0_1,L0_2]=legendre_poly(jj,x);
leg(i,j)=L0;
end
end
%Construct Hierarchical Modal Legendre Basis 分层模型
leg2=leg;
if filter_type == 1
for i=1:P
x=xgl(i);
leg2(i,1)=0.5*(1 - x);
leg2(i,2)=0.5*(1 + x);
for j=3:P
leg2(i,j)=leg(i,j) - leg(i,j-2);
end
end
leg=leg2; %Hierarchical Model Filter (no need for DSS)
end
%Compute Inverse
leg_inv=inv(leg);
%Store Vandermonde Matrices
vdm=leg;
vdm_inv=leg_inv;
%Compute Weight
ibeg=round(2*P/3);
for i=1:P
weight(i)=exp( - alpha*( (i-1)/p )^order );
end
weight;
%Construct 1D Filter Matrix
for i=1:P
for j=1:P
sum=0;
for k=1:P
sum=sum + leg(i,k)*weight(k)*leg_inv(k,j);
end %k
f(i,j)=xmu*sum;
end %j
f(i,i)=f(i,i) + (1-xmu);
end %i
end
%%
function [coord,intma] = create_grid_dg(ngl,nelem,xgl,icase)
%Set some constants
if (icase == 1)
xmin=0;
xmax=2;
elseif (icase == 6)
xmin=-3;
xmax=+3;
end
dx=(xmax-xmin)/nelem;
%Generate Grid Points
ip=1;
coord(1,1)=xmin;
for i=1:nelem
x0=xmin + (i-1)*dx;
coord(1,i)=x0;
intma(1,i)=ip;
for j=2:ngl
ip=ip + 1;
coord(j,i)=( xgl(j)+1 )*dx/2 + x0;
intma(j,i)=ip;
end %j
end %i
end
%%
function [qe] = exact_solution_dg(intma,coord,npoin,nelem,ngl,time,icase)
%Initialize
qe=zeros(npoin,2);
qb=zeros(npoin,1);
%Set some constants
xmin=min(min(coord)); %min(coord)返回每一列最小值
xmax=max(max(coord)); %max(max(coord))返回所有元素最大值
xm=0.5*(xmax + xmin);
xc=xm;
xl=xmax-xmin;
%Generate Grid Points
for e=1:nelem
for i=1:ngl
x=coord(i,e);
ip=intma(i,e);
if (icase == 1) %Gaussian IC with flat bottom
qe(ip,1)=sin(pi*(x-time)) + 0.01;
qe(ip,2)=0;
%if(x>=0)
% qe(ip,1)=1;
%else
% qe(ip,1)=-1;
%end
%qe(ip,2)=0;
end
end %i
end %e
end
%%
function [m1,m2] = compute_Mass_and_Energy(qp,intma,coord,wnq,psi,nelem,ngl,nq)
%Compute Mass and Energy
m1=0; m2=0;
for e=1:nelem
dx=coord(ngl,e)-coord(1,e);
jac=dx/2;
for l=1:nq
wq=wnq(l)*jac;
for j=1:ngl
jp=intma(j,e);
h=(qp(jp,1));
U=0.5*h^2;
m1=m1 + wq*(h)*psi(j,l); %质量是节点量 h 与形函数和权重的乘积，能量是动能 U 与形函数和权重的乘积
m2=m2 + wq*(U)*psi(j,l);
end
end
end
end
%%
function Mmatrix = create_Mmatrix(intma,coord,npoin,nelem,ngl,nq,wnq,psi)
%Initialize
Mmatrix=zeros(npoin,npoin);
for e=1:nelem
%Store Coordinates
for i=1:ngl
x(i)=coord(i,e);
end
dx=x(ngl)-x(1);
jac=dx/2;
%Do LGL Integration
for l=1:nq
wq=wnq(l)*jac;
for i=1:ngl
I=intma(i,e);
h_i=psi(i,l);
for j=1:ngl
J=intma(j,e);
h_j=psi(j,l);
Mmatrix(I,J)=Mmatrix(I,J) + wq*h_i*h_j;
end %j
end %i
end %l
end %e
end
%%
function [a0,a1,beta] = compute_ti_coefficients(kstages)
a0=zeros(kstages);
a1=zeros(kstages);
beta=zeros(kstages);
if kstages == 1 %RK1-SSP
a0(1)=1; a1(1)=0; beta(1)=1;
elseif kstages == 2 %RK2-SSP
a0(1)=1; a1(1)=0; beta(1)=1; %SSP
a0(2)=1.0/2.0; a1(2)=1.0/2.0; beta(2)=1.0/2.0;
% a0(1)=1; a1(1)=0; beta(1)=0.5; %Non-SSP
% a0(1)=1; a1(1)=0; beta(1)=1;
elseif kstages == 3 %RK3-SSP
a0(1)=1; a1(1)=0; beta(1)=1;
a0(2)=3.0/4.0; a1(2)=1.0/4.0; beta(2)=1.0/4.0;
a0(3)=1.0/3.0; a1(3)=2.0/3.0; beta(3)=2.0/3.0;
elseif kstages == 4 %RK4-Non-SSP
a0(1)=1; a1(1)=0; beta(1)=1.0/2.0;
a0(2)=1; a1(2)=0; beta(2)=1.0/2.0;
a0(3)=1; a1(3)=0; beta(3)=1.0;
a0(4)=1; a1(4)=0; beta(4)=0;
end
end
%%
function rhs = create_Lmatrix_SIPDG(rhs,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,q,tau,visc)
rhs_temp=zeros(npoin,1);
%VOLUME Term: IBP
rhs_temp=create_Lmatrix_IBP(rhs_temp,intma,coord,nelem,ngl,nq,wnq,dpsi,q);
%Flux
rhs_temp=create_Flux_SIPDG(rhs_temp,intma,coord,nelem,ngl,nq,psi,dpsi,q,tau);
%Update
rhs=rhs + visc*rhs_temp;
end
%%
function rhs = create_Lmatrix_IBP(rhs,intma,coord,nelem,ngl,nq,wnq,dpsi,q)
%Initialize
x=zeros(ngl,1);
for e=1:nelem
%Store Coordinates
for i=1:ngl
x(i)=coord(i,e);
end
dx=x(ngl)-x(1);
jac=dx/2;
ksi_x=2/dx;
for l=1:nq
wq=wnq(l)*jac;
%Form Derivative
q_x=0;
for j=1:ngl
J=intma(j,e);
q_x=q_x + dpsi(j,l)*ksi_x*q(J);
end
%Form RHS
for i=1:ngl
dhdx_i=dpsi(i,l)*ksi_x;
I=intma(i,e);
rhs(I)=rhs(I) - wq*dhdx_i*q_x;
end %i
end %l
end %e
end
%%
function rhs = create_Flux_SIPDG(rhs,intma,coord,nelem,ngl,nq,psi,dpsi,q,tau)
%FLUX
for e=1:nelem
eL=e;
eR=e+1;
if (eR > nelem)
eR=1;
end
%Left Values
dx_L=coord(ngl,eL)-coord(1,eL);
jac_L=dx_L/2;
ksi_x_L=2/dx_L;
IL=intma(ngl,eL);
dhdx_L=dpsi(ngl,nq)*ksi_x_L;
h_L=psi(ngl,nq);
q_L=q(IL);
dqdx_L=0;
for i=1:ngl
ip=intma(i,eL);
dqdx_L=dqdx_L + dpsi(i,nq)*q(ip)*ksi_x_L;
end
n_L=+1; %(pointing from Left -> Right)
%Right Values
dx_R=coord(ngl,eR)-coord(1,eR);
ksi_x_R=2/dx_R;
IR=intma(1,eR);
dhdx_R=dpsi(1,1)*ksi_x_R;
h_R=psi(1,1);
q_R=q(IR);
dqdx_R=0;
for i=1:ngl
ip=intma(i,eR);
dqdx_R=dqdx_R + dpsi(i,1)*q(ip)*ksi_x_R;
end
n_R=-1;
%Numerical Fluxes
q_star=0.5*(q_L + q_R);
mu_L=(ngl-1)*(ngl-1+1)/(dx_L)*tau; %Shabazi
mu_R=(ngl-1)*(ngl-1+1)/(dx_R)*tau;
mu=max(mu_L,mu_R);
dqdx_mean=0.5*( dqdx_L + dqdx_R);
dqdx_diss=mu*(q_R-q_L);
%Add to Left/Right Elements
rhs(IL)=rhs(IL) + n_L*h_L*dqdx_mean - 0*h_L*dqdx_diss + 0*n_L*dhdx_L*(q_L-q_star); %not valid for periodic
rhs(IR)=rhs(IR) + n_R*h_R*dqdx_mean - 0*h_R*dqdx_diss + 0*n_R*dhdx_R*(q_R-q_star);
end
end
%%
function rhs = create_Lmatrix_LDG(rhs,Mmatrix,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,q,alpha,beta,visc)
%Build Q: VOLUME
rhs_temp=create_Dmatrix(intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,q);
%Build Q: FLUX
rhs_temp=create_Flux_q(rhs_temp,intma,nelem,ngl,q,alpha,beta);
%Build Q using Lift Operator
Q=Mmatrix\rhs_temp;
%Build q: VOLUME
rhs_temp=create_Dmatrix(intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,Q);
%Build q: FLUX
rhs_temp=create_Flux_q(rhs_temp,intma,nelem,ngl,Q,beta,alpha);
%Add Viscous term
rhs = rhs + visc*rhs_temp;
end
%%
function rhs = create_Dmatrix(intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,q)
%Initialize
rhs=zeros(npoin,1);
x=zeros(ngl,1);
for e=1:nelem
%Store Coordinates
for i=1:ngl
x(i)=coord(i,e);
end
dx=x(ngl)-x(1);
jac=dx/2;
ksi_x=2/dx;
for l=1:nq
wq=wnq(l)*jac;
%interpolate q
q_e=0;
for j=1:ngl
J=intma(j,e);
q_e=q_e + psi(j,l)*q(J);
end
%Form RHS
for i=1:ngl
dhdx_i=dpsi(i,l)*ksi_x;
I=intma(i,e);
rhs(I)=rhs(I) - wq*dhdx_i*q_e;
end %i
end %l
end %e
end
%%
function rhs = create_Flux_q(rhs,intma,nelem,ngl,q,alpha,beta)
%Build Q: FLUX
for e=1:nelem
eL=e;
eR=e+1;
if (eR > nelem)
eR=1;
end
%LGL Integration
IL=intma(ngl,eL);
IR=intma(1,eR);
q_l=q(IL);
q_r=q(IR);
nvector=+1; %(pointing from Left -> Right)
qstar=alpha*q_l + beta*q_r;
%Add to Left
rhs(IL)=rhs(IL) + nvector*qstar;
%Add to Right
rhs(IR)=rhs(IR) - nvector*qstar;
end %ie
end
%%
function qp = apply_filter_dg(qp,f,intma,nelem,ngl)
%Initialize
h_e=zeros(ngl,1);
U_e=zeros(ngl,1);
hf=zeros(ngl,1);
Uf=zeros(ngl,1);
%Integrate Divergence of Flux
for e=1:nelem
%Store Local Solution
for i=1:ngl
ip=intma(i,e);
h_e(i)=qp(ip,1);
U_e(i)=qp(ip,2);
end
%Form Filtered Variable
hf=f*h_e;
Uf=f*U_e;
%Store It
for i=1:ngl
ip=intma(i,e);
qp(ip,1)=hf(i);
qp(ip,2)=Uf(i);
end
end %ie
end
%%
function [limit_element,qp,qmodal] = limiter_nodal_v3(qp,vdm,vdm_inv,intma,nelem,ngl,limit)
%Initialize
limit_element=zeros(nelem,1);
qmodal=zeros(ngl,nelem);
qtemp=zeros(ngl,nelem);
qnodal=zeros(ngl,1);
for k=2:-1:1
%Store Element-wise Solution
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
qtemp(i,e)=qp(ip,k);
end
end
%Compute Modal Solution
for e=1:nelem
qmodal(:,e)=vdm_inv(:,:)*qtemp(:,e);
end
%Limit Solution
if (limit > 0)
[limit_element,qmodal] = limiter_modal(qmodal,nelem,ngl);
end
%Compute Nodal Solution
for e=1:nelem
qnodal(:)=vdm(:,:)*qmodal(:,e);
for i=1:ngl
ip=intma(i,e);
qp(ip,k)=qnodal(i);
end
end
end %k
%Get Modal Solution for Plotting
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
qtemp(i,e)=qp(ip,1);
end
end
for e=1:nelem
qmodal(:,e)=vdm_inv(:,:)*qtemp(:,e);
end
end
%%
function [limit_element,qp] = limiter_modal(qp,nelem,ngl)
%Initialize
limit_element=zeros(nelem,1);
%Integrate Flux Terms
for ie=1:nelem
iel=ie-1;
ier=ie+1;
if (ie == 1)
iel=nelem;
end
if (ie == nelem)
ier=1;
end
iend=ngl;
iterations=0;
limit=1;
while limit == 1
iterations=iterations + 1;
%Get Coefficients
q_e=qp(iend,ie);
dq_p=qp(iend-1,ier) - qp(iend-1,ie);
dq_m=qp(iend-1,ie) - qp(iend-1,iel);
%Limit Solution
[qlimit,limit]=minmod(q_e,dq_p,dq_m);
if iterations == 1
limit_element(ie)=limit;
end
%Limit Solution
if limit == 1
qp(iend,ie)= qlimit;
iend=iend - 1;
%iend_min=round(ngl/2);
iend_min=1;
if iend == iend_min
limit = 0;
end
end %if
end %while
end %ie
end
%%
function [qlimit,limit] = minmod(q1,q2,q3)
tol=1e-16;
sign1=sign(q1);
sign2=sign(q2);
sign3=sign(q3);
q=zeros(3,1);
q(1)=abs(q1);
q(2)=abs(q2);
q(3)=abs(q3);
if (sign1 == sign2 && sign2 == sign3)
qlimit=sign1*min(q);
else
qlimit=0;
end
limit=1;
if ( abs(qlimit-q1) <= tol)
limit=0;
end
end
%%
function [limit_element,qp,qmodal] = limiter_nodal_bp(qp,q0,vdm,vdm_inv,intma,nelem,ngl)
%Initialize
limit_element=zeros(nelem,1);
qtemp=zeros(ngl,nelem);
qnodal=zeros(ngl,nelem);
qmodal=zeros(ngl,nelem);
qmodal_loc=zeros(ngl,1);
qlocal=zeros(ngl,1);
eps=1e-6;
tvector=zeros(3,1);
%Loop through variables
for k=1:2
qtemp=zeros(ngl,nelem);
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
qtemp(i,e)=q0(ip,k);
qnodal(i,e)=qp(ip,k);
end
end
%Global Extrema
global_max=max(max(qtemp));
global_min=min(min(qtemp));
%Loop through element
for e=1:nelem
%Store Local Space
qlocal(:)=qnodal(:,e);
%Map to Modal Space
qmodal_loc(:)=vdm_inv(:,:)*qlocal(:);
%Local Extrema
local_max=max(qlocal);
local_min=min(qlocal);
%Computer Limiter Weight
qmean=qmodal_loc(1);
t1_den=local_max - qmean;
t2_den=local_min - qmean;
tvector(1)=abs((global_max - qmean)/( (local_max - qmean) + eps ));
tvector(2)=abs((global_min - qmean)/( (local_min - qmean) + eps ));
tvector(3)=1;
theta=min(tvector);
if (theta < 1)
limit_element(e)=1;
end
for i=1:ngl
ip=intma(i,e);
qp(ip,k)=theta*qlocal(i) + (1-theta)*qmean;
end
end %e
end %k
%Map to Modal Space
qtemp=zeros(ngl,nelem);
for e=1:nelem
for i=1:ngl
ip=intma(i,e);
qtemp(i,e)=qp(ip,1);
end
end
for e=1:nelem
qmodal(:,e)=vdm_inv(:,:)*qtemp(:,e);
end
end
%%
function rhs = apply_dss(rhs,intma,intma_cg,coord,Mmatrix,npoin,npoin_cg,nelem,ngl,nq,wnq,psi,dpsi,alpha,ivolume_integrate)
%Initialize
rhs_continuous=zeros(npoin_cg,2);
%Volume Integrate
if (ivolume_integrate == 1)
%rhs_volume = compute_volume_integral(rhs,intma,coord,npoin,nelem,ngl,nq,wnq,psi);
rhs_volume = create_rhs_volume(rhs,intma,coord,npoin,nelem,ngl,nq,wnq,psi,dpsi,alpha);
else
rhs_volume=rhs;
end
%Gather Solution to make it Global
for e=1:nelem
for i=1:ngl
ip_cg=intma_cg(i,e);
ip=intma(i,e);
for k=1:2
rhs_continuous(ip_cg,k)=rhs_continuous(ip_cg,k) + rhs_volume(ip,k);
end
end %i
end %e
%Scatter Solution to make it Local
rhs=zeros(npoin,2);
rhs_continuous(:,1)=Mmatrix\rhs_continuous(:,1);
rhs_continuous(:,2)=Mmatrix\rhs_continuous(:,2);
for e=1:nelem
for i=1:ngl
ip_cg=intma_cg(i,e);
ip=intma(i,e);
for k=1:2
rhs(ip,k)=rhs_continuous(ip_cg,k);
end
end %i
end %e
end
%%
function [test] = plot_modes(qmodal,Ipointer,limit_element,intma,coord,qp,nelem,ngl,nq,nop,noq,time,diss,main_text)
test = 0;
II=Ipointer';
QQ(:,:)=qmodal';
icheck=[1:ngl];
%-------PLOT MODES------%
subplot(2,1,1);
%bar(II(1:nelem,icheck),QQ(1:nelem,icheck),'stacked');
bar(QQ(1:nelem,icheck));
%xlabel('Methods','FontSize',18);
%ylabel('Modes','FontSize',18);
%legend('0','1','2','3','4','5','6','7','8');
title_text=[main_text ': diss = ', num2str(diss) ': Ne = ' num2str(nelem) ', N = ' num2str(nop) ', Q = ' num2str(noq) ', Time = ' num2str(time)];
title([title_text],'FontSize',18);
qmin=min(min(qmodal));
qmax=max(max(qmodal));
axis([ 1 nelem qmin qmax]);
set(gca,'FontSize',18);
%hold on;
%-------PLOT Limiter------%
% subplot(3,1,2);
% bar(limit_element);
% qmin=min(limit_element);
% qmax=max(limit_element);
% qmin=0;
% qmax=1;
% axis([ 1 nelem qmin qmax]);
% %ylabel('Limiter','FontSize',18);
% set(gca,'FontSize',18);
%-------PLOT Solution------%
subplot(2,1,2)
qq=zeros(nq,nelem);
for ie=1:nelem
for i=1:ngl
ip=intma(i,ie);
x(i)=coord(i,ie);
y(i)=qp(ip,1);
end
plot_handle=plot(x,y,'r-');
set(plot_handle,'LineWidth',2);
hold on;
end
% axis([ -1 +1 -0.25 1.25]);
%ylabel('Solution','FontSize',18);
set(gca,'FontSize',18);
% M_i=getframe(gcf);
% M_modes(i)=M_i;
pause(0.1);
hold off;
test=1;
end
%%
function [p_movie,u_movie,time_movie,mass_movie,energy_movie,L2_norm,iframe] = store_movie_frames(qp,q0,intma,coord,wnq,psi,npoin,nelem,ngl,nq,p_movie,u_movie,time_movie,mass_movie,energy_movie,iframe,mass0,energy0,time,icase)
iframe=iframe + 1;
p_movie(:,iframe)=qp(:,1);
u_movie(:,iframe)=qp(:,2);
time_movie(iframe)=time;
[m1,m2] = compute_Mass_and_Energy(qp,intma,coord,wnq,psi,nelem,ngl,nq);
mass_movie(iframe)=m1/mass0;
energy_movie(iframe)=m2/energy0;
%Compute Norm
[qe] = exact_solution_dg(intma,coord,npoin,nelem,ngl,time,icase);
h_top=norm(q0(:,1)-qe(:,1),2);
h_bot=norm(qe(:,1),2);
u_top=norm(q0(:,2)-qe(:,2),2);
u_bot=norm(qe(:,2),2);
L2_norm(1,iframe)=h_top;
L2_norm(2,iframe)=u_top;
end
可以对上述代码所使用的函数进行删减，把代码中的burgers方程换位带有sod初值条件的Euler方程，且不含有精确解的表达式，给出修正后的代码
本文介绍了一个用于创建矩形域中规则网格的算法，并提供了MATLAB函数实现。该算法可以根据需要在X和Y方向上进行包裹，并通过调整参数实现网格中心及边缘的变形效果。
扫一扫
4032
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
