;demo_advection
;Written by Thomas Wiegelmann
;11.09.2008
;Email: wiegelmann@mps.mpg.de
;This program contains several explicit
;methods to solve equations in conservative form:
;du/dt=-F(u)
;Advection equation
;Invicid Burgers equation
;Diffusion equation
;Diffusive Burgers equation
;
;
;Euler method (It is unstable for hyperbolic equations)
;Leap-Frog 2. and 4th order (Unstable for parabolic Eq.)
;Upwind method 1. order
;Lax method (very diffusive)
;Lax Wendroff method
;
;Remark: In practice Leap-Frog and Lax-Wendroff
;are the most important of these methods
;Upwind methods play some role in shock-capturing method
;See PDE-lecture for details
;
FUNCTION rho_0, x
;Gauss profile
x0=3.0
rho_bg=0.1
rho0=rho_bg+1.0*exp(-(x-x0)^2/2)
return,rho0
END
;
FUNCTION rho_1, x
;Rectangular profile
;Square wave test for advection equation
x0=5
rho0=0.0*x+0.1
u=where((x ge x0-1.) and (x le x0+1.))
rho0[u]=1.1
return,rho0
END
;
FUNCTION compute_F, utest
common eq_type, eq_type
;Flux for all methods
;It is of course not efficients to compute F always
;for all methods, where only one is used,
;but this program is aimed only for demonstrations
v0=0.1
;Advection equation
F1=v0*utest
;Invicsid Burgers-Equation (a nonlinear Equation)
F2=0.5*v0*utest*utest
; diffusion equation or diffusion term for Burgers
dx=10.0/n_elements(utest)
D=0.01
F3=-D*(shift(utest,-1)-shift(utest,1))/(2*dx)
CASE eq_type OF
   1: F0=F1
   2: F0=F2
   3: F0=F3
   4: F0=F2+F3
ENDCASE
return,F0
END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Main program
common eq_type, eq_type
xpixel=450 & ypixel=400
;window,xsize=xpixel,ysize=ypixel
mydevice = !D.NAME

;device,set_resolution=[xpixel,ypixel]
eq_type=1;
print,'Which Equation to solve?'
print,'1=Advection, 2=inviscid Burgers equation (shocks)'
print,'3=Diffussion Eq., 4=Diffusive Burgers Eq.'
read,eq_type
CASE eq_type OF
   1: title1='Advection Equation'
   2: title1='Inviscid Burgers Eq.'
   3: title1='Diffusion Equation'
   4: title1='Diffusive Burgers Eq.'
ENDCASE
prof=1
print,'Initial profile: 1=Gauss,  2=Square box'
read,prof
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Initial conditions
nx=100
xmax=10.0
xmin=0.0
x=xmax*findgen(nx)/(nx)-xmin
x0=5
rho0=rho_0(x) ; Gauss-shape
if prof eq 2 then $
rho0=rho_1(x) ; Square box
dt=0.5; iteration time step
;dt=1.0
dx=x[1]-x[0]; equidistant grid
rho=rho0 & rho1=rho0 & rho2=rho0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Here to choose methods
;
;Methods: 0=Euler,
        ; 1=Leap-Frog 2th order,
        ; 2=Leap-Frog 4th order;
        ; 3=upwind 1th order
        ; 4=Lax,
        ; 5=Lax-Wendroff
method=5;
print,'Which method?'
print,'0=Euler,1=Leap-Frog 2th order, 2=Leap-Frog 4th order'
print,'3=upwind 1th order, 4=Lax, 5=Lax-Wendroff'
read,method
xi=0.00; viscosity (only for Leap-Frog)
if (method eq 1) or (method eq 2) then begin
print,'Input viscosity 0.0<=xi<1.0'
read,xi
endif
;
;
for it=0, 200 do begin
t=it*dt
F=compute_F(rho)
;
;Euler-scheme
if method eq 0 then $
rho1=rho-(dt/(2*dx))*(shift(F,-1)-shift(F,1))
;
;Leap-Frog scheme
;2th order
if method eq 1 then $
rho1=rho2-(dt/dx)*(shift(F,-1)-shift(F,1))
;4th order
if method eq 2 then begin
rho1=rho2-(dt/(6*dx))*$
(-shift(F,-2)+8*shift(F,-1)-8*shift(F,1)+shift(F,2))
endif
;viscosity
;get rid of numerical wiggles
;
ave_rho=(shift(rho1,1)+shift(rho1,-1))/2
if (method eq 1) or (method eq 2) then $
rho1=(1.0-xi)*rho1+xi*ave_rho


;
;Upwind scheme 1. order
if method eq 3 then $
rho1=rho-(dt/(dx))*(F-shift(F,1))
;
;Lax Method
if method eq 4 then $
rho1=0.5*(shift(rho,1)+shift(rho,-1))-(dt/(2*dx))*(shift(F,-1)-shift(F,1))

;
;Lax-Wendroff
if method eq 5 then begin
;half timestep rho_h, Lax-step
rho_h=0.5*(shift(rho,-1)+rho)-(dt/(2.0*dx))*(shift(F,-1)-F)
F_h=compute_F(rho_h)
rho1=rho-(dt/dx)*(shift(F_h,0)-shift(F_h,1))
endif
;
;Update density
rho2=rho
rho=rho1
;set_plot,'z'
plot,x,rho,xstyle=1,ystyle=1,xrange=[xmin,xmax],$
yrange=[-0.2,1.2],xtitle='x',ytitle='rho(x)',thick=2,$
charsize=1.5,title=title1

oplot,x,rho0,linestyle=1
xyouts,3,-0.05,'t='+strcompress(string(t)),size=2
;

wait,0.05
endfor
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;




end