HartenYee

汚いですし、直したいところはたくさんありますが、
晒すことにします。

HartenYeeの2nd Order Upwindです。
おそらく数値計算ハンドブックにも載ってるような記憶がありますが、
基本的には藤井先生の本をトレースしたものです。
「test.txt」は０〜１００の格子番号iに対して10<i<20でf＝1.0、その他で０としています。

あとは[pylab]を入れていればそのまま描画されるはずです。
ない場合は、import以下を消して最後のグラフ描画部分も消して下さい。

数年前のプログラムに、修正を施しただけですので、
もっと、きれいに書けます（もっとOOPな感じで）が、まあ、残しておきます（笑）

OOPならC++、Rubyで書いてますが、Rubyが一番オブジェクト指向にしっくりきます。

ただ、科学技術計算ならPythonが一番使いやすいと思います。描画は特に。

ちなみにHartenYeeのスキームはあまり利用されていないようですが、
Non-MUSCL系だと好きだったりします。

「中央差分に修正流束を加算する」というスキームの形になってるのが理解しやすいです。



#advection #################################################
# Harten Yee Upwind ########################################
############################################################
#import#####################################################
from pylab import *


# Harten Yee ##############################################
def cal_fai(z):
eps = 0.0001
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 minmod(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"
#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 ################################
f= [0]
fn=  [0]
fini=[0]
fl=  [0]
fr=  [0]
x=   [0]
abm= [0]
xl=  [0]
f_tilde = [0.0]
limiter = [0.0]
fex=[0]
df = [0]
gi = [0]
gamma = [0]
fai = [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])]
f_tilde = f_tilde+[0.0]
fex = fex+[float(a[1])]
limiter = limiter+[0.0]
df = df + [0.0]
gi = gi +[0.0]
gamma = gamma + [0.0]
fai = fai + [0.0]
i=i+1
imx=i
print "imx=",imx
print "*"*40
print "xl,f"
for i in range(0,imx):
print i,"  ",float(xl[i]),"   ",float(f[i])
# Time integration ################################
for n in range(0,1000):
#LAX-WEDNROFF##############################
for i in range(0,imx-1):
df[i] = f[i+1] - f[i]
for i in range(2,imx-1):
gi[i] = minmod(df[i],df[i-1])
for i in range(3,imx-2):
if df[i]!=0.0:
#print i,df[i],gi[i],gi[i+1]
gamma[i] = cal_sigma(u,dt,dx)*(gi[i+1]-gi[i])/df[i]
else:
gamma[i] = 0.0
for i in range(2,imx-1):
fai[i] = cal_sigma(u,dt,dx)*(gi[i]+gi[i+1])-cal_fai(u+gamma[i])*df[i]
for i in range(1,imx-3): # Numerical Flux
CS = 0.5*(u* f[i+1] + u*f[i] ) # CS:Centered Space
CF = 0.5*fai[i]       # CF:Correct Flux
f_tilde[i] = CS + CF       # Numerical FLUX
#print i,limit
for i in range(2,imx-2):#0~imx-1
fn[i] = f[i] - dt/dx*(f_tilde[i] - f_tilde[i-1])
#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 = '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()
###################################################