等間隔WENO

WENOは格子間をラグランジュ補間して、スムーズインジケーターで修正する手法らしいです。ただ、あくまで私の理解であり、正しい認識を提供するものではありません。

河川の流れを対象とし、RANSを用いた計算だと、WENOのような高級な手法を用いることはありませんが、勉強がてらやってみました。

汚いですが、ソースも晒しておきます。たぶん間違いがあると思いますので、あくまで参考ということで。

あとは不等間隔の場合、係数が異なるように思いますので、座標変換して解くような手法（格子間隔を一定として変換）でない限り、修正が必要かと思います。QUICKと同様ですね。

#advection #################################################
# WENO  ########################################
############################################################
#import#####################################################
from pylab import *
def calflux1(u,f,flm,flp):
alpha = 0.0
for i in range(1,imx-1):
alpha_tmp = u
if abs(alpha_tmp) > abs(alpha):
alpha = alpha_tmp
for i in range(1,imx-1):
flm[i] = 0.5*(f[i]*u - alpha*f[i])
flp[i] = 0.5*(f[i]*u + alpha*f[i])
flm[0] = 0.0
flm[imx-1] = 0.0
flp[0] = 0.0
flp[imx-1] = 0.0
def calflux2(i,omegam,omegap,flm,fip,fl2m,fl2p):
print i
fltmpp = [0]
fltmpm = [0]
for k in range(0,3):
fltmpp = fltmpp + [0]
fltmpm = fltmpm + [0]

fltmpp[0] =  1./3.*flp[i-2] - 7./6.*flp[i-1] + 11./6.*flp[i  ]
fltmpp[1] = -1./6.*flp[i-1] + 5./6.*flp[i  ] +  1./3.*flp[i+1]
fltmpp[2] =  1./3.*flp[i  ] + 5./6.*flp[i+1] -  1./6.*flp[i+2]
fl2p[i] = 0.0
for k in range(0,3):
fl2p[i] = fl2p[i] + omegap[k]*fltmpp[k]

fltmpm[1] = -1./6. * flm[i-1] + 5./6. * flm[i  ] +  2./6. * flm[i+1]
fltmpm[2] =  2./6. * flm[i  ] + 5./6. * flm[i+1] -  1./6. * flm[i+2]
fltmpm[3] = 11./6. * flm[i+1] - 7./6. * flm[i+2] +  2./6. * flm[i+3]
fl2m[i] = 0.0
for k in range(0,3):
fl2m[i] = fl2m[i] + omegam[k]*fltmpm[k]
in_ff = open("test.txt")
dx=1.0
dt=0.2
u=0.2
c1 = 1./10.
c2 = 6./10.
c3 = 3./10.
eps = 1.0e-10
print "-"*40
print "u  =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
f= [0]
fn=  [0]
fini=[0]
fl=  [0]
fr=  [0]
fex = [0]
x=   [0]
abm= [0]
xl=  [0]
flm = [0]
flp = [0]
flux = [0]
flux1 =[0]
flux2 =[0]
flux3 =[0]
fl2m = [0]
fl2p = [0]
i=1
for x in in_ff:
print i
a = x.split(" ")
for k in range(0,1):
print a[k]
xl=xl+[float(a[0])]
f=f+[float(a[1])]
fini=fini+[float(a[1])]
fn=fn+[float(a[1])]
fl=fl+[float(a[1])]
fr=fr+[float(a[1])]
i=i+1
fex = fex+[float(a[1])]
flux = flux + [0]
flux1 = flux1 + [0]
flux2 = flux2 + [0]
flux3 = flux3 + [0]
flm = flm + [0]
flp = flp + [0]
fl2m = fl2m + [0]
fl2p = fl2p + [0]

imx=i
adm=[0.0]
alpham = [0]
alphap = [0]
omegam = [0]
omegap = [0]


for k in range(0,3):
alpham = alpham + [0]
alphap = alphap + [0]
omegam = omegam + [0]
omegap = omegap + [0]
print "*"*40
print "xl,f"
for i in range(0,imx-1):
print float(xl[i]),"   ",float(f[i])
# set parameter ####################################
# 1st order upwind ####################################
for n in range(0,1000):
calflux1(u,f,flm,flp)

# WENO
for i in range(3,imx-4):

isp1 = 13./12.*(flp[i-2]-2.0*flp[i-1]+flp[i])*(flp[i-2]-2.0*flp[i-1]+flp[i])+1.0/4.0*(flp[i-2]-4.0*flp[i-1]+3.0*flp[i])*(flp[i-2]-4.0*flp[i-1]+3.0*flp[i])
isp2 = 13./12.*(flp[i-1]-2.0*flp[i]+flp[i+1])*(flp[i-1]-2.0*flp[i]+flp[i+1])+1.0/4.0*(flp[i-1]-flp[i+1])*(flp[i-1]-flp[i+1])
isp3 = 13./12.*(flp[i]-2.*flp[i+1]+flp[i+2])*(flp[i]-2.*flp[i+1]+flp[i+2])+1./4.0*(3.*flp[i]-4.*flp[i+1]+flp[i+2])*(3.*flp[i]-4.*flp[i+1]+flp[i+2])

ism1 = 13./12.*(flm[i+1]-2.0*flm[i+2]+flm[i+3])*(flm[i+1]-2.0*flm[i+2]+flm[i+3])+1.0/4.0*(3.0*flm[i+1]-4.0*flm[i+2]+flm[i+3])*(3.0*flm[i+1]-4.0*flm[i+2]+flm[i+3])
ism2 = 13./12.*(flm[i]-2.0*flm[i+1]+flm[i+2])*(flm[i]-2.0*flm[i+1]+flm[i+2])+1.0/4.0*(flm[i]-flm[i+2])*(flm[i]-flm[i+2])
ism3 = 13./12.*(flm[i-1]-2.0*flm[i  ]+flm[i+1])*(flm[i-1]-2.0*flm[i]+flm[i+1])+1./4.0*(flm[i-1]-4.*flm[i]+3.*flm[i+1])*(flm[i-1]-4.*flm[i]+3.*flm[i+1])
#
alpham[0] = c1/(eps + ism1)/(eps+ism1)
alpham[1] = c2/(eps + ism2)/(eps+ism2)
alpham[2] = c3/(eps + ism3)/(eps+ism3)
#
alphap[0] = c1/(eps + isp1)/(eps+isp1)
alphap[1] = c2/(eps + isp2)/(eps+isp2)
alphap[2] = c3/(eps + isp3)/(eps+isp3)
#
sigmaalpham = alpham[0] +alpham[1] +alpham[2]
sigmaalphap = alphap[0] +alphap[1] +alphap[2]
#
omegam[0] = alpham[0] / sigmaalpham
omegam[1] = alpham[1] / sigmaalpham
omegam[2] = alpham[2] / sigmaalpham

omegap[0] = alphap[0] / sigmaalphap
omegap[1] = alphap[1] / sigmaalphap
omegap[2] = alphap[2] / sigmaalphap

calflux2(i,omegam,omegap,flm,flp,fl2m,fl2p)
fl2m[2] = 0.0
fl2m[1] = 0.0
fl2m[0] = 0.0
fl2p[2] = 0.0
fl2p[1] = 0.0
fl2p[0] = 0.0

fl2m[imx-1] = 0.0
fl2m[imx-2] = 0.0
fl2m[imx-3] = 0.0
fl2p[imx-1] = 0.0
fl2p[imx-2] = 0.0
fl2p[imx-3] = 0.0
# EULER ##############################
for i in range(3,imx-4):#0~imx-1
fn[i]=f[i]-dt*(fl2p[i]-fl2p[i-1]+fl2m[i]-fl2m[i-1])/(dx)
#shift#########################################
for i in range(2,imx-2):
f[i]=fn[i]
print fn
###################################################
# Exact ###########################################
for i in range(0,imx-1):
if i>=50 and i<=60:
fex[i] = 1.0
else:
fex[i] = 0.0
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
ff.write(str(xl[i]))
ff.write("    ")
ff.write(str(f[i]))
ff.write("\n")
ff.close()
#output
#for i in range(0,imx-1):
# print f[i]
print "-"*100
print xl
print f
#####################graph#########################
y1=f
y2=fini
y3 = fex
plot(xl,y1,label = 'WENO')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('WENO')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('WENO')
show()
###################################################