% This Matlab script solves the one-dimensional convection % equation using a finite volume algorithm. The % discretization forward Euler in time. clear all; close all; clc; clf; red=[224,32,40]/255; green=[0,209,42]/255; blue=[26,46,250]/255; black=[0,0,0]/255; orange=[255,152,51]/255; magenta=[255,19,255]/255; global a; a=1; Nx = 128; L=2; x = linspace(0,L,Nx); CFL = 0.5; s.dx = L/(Nx-1); s.dt = CFL*s.dx/abs(a); cr=1:(Nx-1); cl=circshift(cr,1); cll=circshift(cr,2); crr=circshift(cr,-1); s.faces={}; for i=1:Nx s.faces{i}.Ul=0; s.faces{i}.Ur=0; s.faces{i}.l=1; s.faces{i}.n=[1;0;0]; s.faces{i}.nl=i-1; s.faces{i}.nr=i; if i==Nx s.faces{i}.cll=cll(1); s.faces{i}.cl=cl(1); s.faces{i}.cr=cr(1); s.faces{i}.crr=crr(1); else s.faces{i}.cll=cll(i); s.faces{i}.cl=cl(i); s.faces{i}.cr=cr(i); s.faces{i}.crr=crr(i); end if i==Nx s.faces{i}.nr=-1; end if i==1 s.faces{i}.nl=-1; end end s.cells={}; for i=1:(Nx-1) s.cells{i}.Q=0; s.cells{i}.V=s.dx; s.cells{i}.x=x(i)+s.dx/2; xc(i)=s.cells{i}.x; end Q=zeros(Nx-1,1); Q(x>=0&x<=0.25)=1; Qu=Q; Qglf=Q; Qllf=Q; Qlw=Q; Qroe=Q; Qhoroe=Q; figure(1); hold on; p1=plot(xc,Qglf,'-','Color',red,'DisplayName','Upwind','LineWidth',2); p2=plot(xc,Qglf,'--','Color',green,'DisplayName','Global Lax-Friedrich (GLF)','LineWidth',2); p3=plot(xc,Qllf,'-.','Color',blue,'DisplayName','Local Lax-Friedrich (LLF)','LineWidth',2); p4=plot(xc,Qlw,'-','Color',black,'DisplayName','Lax-Wendroff (LW)','LineWidth',2); p5=plot(xc,Qroe,'--','Color',orange,'DisplayName','Roe','LineWidth',2); p6=plot(xc,Qhoroe,'-','Color',magenta,'DisplayName','High-order Roe','LineWidth',2); hold off; ylim([-0.25,1.25]); xlim([0,L]); grid on; legend show; xlabel('x (m)'); set(gca,'fontweight','bold'); t = 0; tend=0.7; while(t0 F=Fl; else F=Fr; end end %Global Lax-Friedrichs if method==1 w=dx/dt*delta; F=0.5*(Fl+Fr)-0.5*w; end %Local Lax-Friedrichs if method==2 %w=max(|F'(ul)|,|F'(ur)|) wl=abs(a*nx); wr=abs(a*nx); w=max(wl,wr)*delta; F=0.5*(Fl+Fr)-0.5*w; end %Lax-Wendroff if method==3 A=a; %w=dt/dx*A*(Fr-Fl)=dt/dx*A^2*delta w=dt/dx*A^2*delta; F=0.5*(Fl+Fr)-0.5*w; end %Roe if method==4 A=a; %Jacobian A is a since f'(u)=a. Remember that f(u)=au w=abs(A)*delta; F=0.5*(Fl+Fr)-0.5*w; end end % Limiter function %%%%%%%%%%%%%%%%%%%%%% function phi=limiter(r,type) switch type case 'minmod' phi=max([0, min([r,1])]); case 'monotonizedcentral' phi=max([0, min([2*r,0.5*(r+1),2])]); case 'osher' phi=max([0,min(r,1.5)]); case 'sweby' phi=max([0, min([1.5*r,1]),min([r,1.5])]); case 'superbee' phi=max([0,min(2*r,1),min(r,2)]); case 'vanalbada1' phi=(r^2+r)/(r^2+1); case 'vanalbada2' phi=2*r/(r^2+1); case 'vanleer' phi=(r+abs(r))/(1+abs(r)); case 'zero' phi=0; otherwise phi=1; end end