Harten-Yeeをちょっと修正

以前のものを多少修正しました。
もうちょっと美しく出来れば良いんですけど・・・。
1次風上も残してますが、２０行くらい少なくなってます。


#advection #################################################
# Harten Yee Upwind ########################################
############################################################
#import#####################################################
from pylab import *
def cal_fai(z):
    eps = 0.0000000001
    if abs(z) > eps:
        return abs(z)
    else:
        return (z*z+eps*eps)/(2*eps)
def cal_sigma(z,dt,dx):
    return 0.5*(cal_fai(z)-dt/dx*z*z)
def minmod2(x,y):
    if abs(x)<=abs(y) and x*y >=0.0:
        return x
    elif abs(x)>abs(y) and x*y >= 0.0:
        return y
    elif x*y<0:
        return 0.0
    else:
        print "Error minmode2"
# Harten Yee ##############################################
class Cell:
    def __init__(self,dxl,u,f,fn,fl,fr):
        self.dxl   = dxl
        self.u     = u
        self.f     = f
        self.fn    = fn
        self.fl    = fl
        self.fr    = fr
        self.df    = 0.0
        self.gi    = 0.0
        self.gamma = 0.0
        self.fai   = 0.0
    def upwind(self,dt,cellup,celldn):
        ultmp = 0.5*(u + cellup.u)
        urtmp = 0.5*(celldn.u + u)
        if self.u > 0:
            self.fl = ultmp*cellup.f
            self.fr = urtmp*self.f
        else:
            self.fl = ultmp*self.f
            self.fr = urtmp*celldn.f
    def caldf(self,celldn):
        self.df = celldn.f - self.f
    def calgi(self,cellup):
        self.gi = minmod2(self.df,cellup.df)
    def calgamma(self,dt,celldn):
        if self.df != 0:
            self.gamma = cal_sigma(self.u,dt,self.dxl)*(celldn.gi - self.gi)/self.df
        else:
            self.gamma = 0.0
    def calfai(self,celldn):
        self.fai = cal_sigma(self.u,dt,self.dxl)*(self.gi+celldn.gi)-cal_fai(self.u+self.gamma)*self.df
    def yee(self,dt,cellup,celldn):
        csdn = 0.5*(celldn.u*celldn.f+self.u*self.f)
        csup = 0.5*(self.u*self.f + cellup.u*cellup.f)
        cfdn = 0.5*self.fai
        cfup = 0.5*cellup.fai
        self.fr = csdn+cfdn
        self.fl = csup+cfup
    def euler(self,dt):
        self.fn = self.f - dt/self.dxl*(self.fr - self.fl)
    def shift(self):
        self.f = self.fn
#set parameters ############################################
in_ff = open('test.txt')
imx=100
dx =  1.0
dt =  0.2
u  =  0.2
print "-"*40
print "u  =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
cells = []
fini = []
i = 0
for x in in_ff:
    print i
    a = x.split("\t")
    for k in range(0,1):
        print a[k]
    cells = cells + [Cell(float(dx),float(u),float(a[1]),float(a[1]),0.0,0.0)]
    fini = fini + [float(a[1])]
    i=i+1
imx=i
print "imx=",imx
print "*"*40
print "xl,f"
for i in range(0,imx):
    print i,"  ",cells[i].f,cells[i].u
# Time integration ################################
for n in range(0,1000):
    #Harten Yee##############################
    for i in range(0,imx-1):
        cells[i].caldf(cells[i+1])
    for i in range(2,imx-1):
        cells[i].calgi(cells[i-1])
    for i in range(3,imx-2):
        cells[i].calgamma(dt,cells[i+1])
    for i in range(2,imx-1):
        cells[i].calfai(cells[i+1])
    for i in range(1,imx-3):
        cells[i].yee(dt,cells[i-1],cells[i+1])
    for i in range(2,imx-2):
        cells[i].euler(dt)
    #shift#########################################
    for i in range(2,imx-2):
        cells[i].shift()
# Exact ###########################################
fout = [0]
fex = [0]
xl  = [0]
for i in range(0,imx-1):
    if i>=49 and i<=59:
        fex = fex +[1.0]
    else:
        fex = fex +[0.0]
    fout = fout + [float(cells[i].f)]
    xl = xl + [float(i)]
#output############################################
ff=open("output.dat","w")
for i in range(0,imx-1):
    ff.write(str(i))
    ff.write("    ")
    ff.write(str(fout[i]))
    ff.write("\n")
ff.close()
#####################graph#########################
y1=fout
y2=fini
y3=fex
print len(fout),len(fini),len(fex)
plot(xl,y1,label = 'Harten-Yee')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('Hatren-Yee')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('Harten-Yee')
show()
###################################################