2013年10月20日日曜日

RCIP


CIPは非常に優秀な移流方程式に対する解法ですが、不連続部に若干のオーバーシュート、アンダーシュートが残ります。

これは補間関数に3次関数を選択したためであり、これを避けるために、別の関数を使えば良いじゃん!というのが、RCIP(Rational CIP)基本的な考え方です。有理関数を使えばオーバーシュートアンダーシュートを抑えられます。

ソースを晒しておきます。まあ、汚いです。

#advection #################################################
# RCIP ########### ########################################
############################################################
#import#####################################################
from pylab import * 
#set parameters ############################################
in_ff = open('test.txt')
imx=100
dx=1.0
dt=0.2
u=0.2
alpha = 1.0
eps = 1.0e-15
print "-"*40
print "u  =",u
print "dx =",dx
print "dt =",dt
print "-"*40
############################################################
#initial condition ################################
f= [0]
fn=  [0]
df = [0]
dfn = [0]
fini=[0]
fl=  [0]
fr=  [0]
x=   [0]
abm= [0]
xl=  [0]
f_tilde = [0]
fex = [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]
fex = fex + [0]
df = df + [0]
dfn = dfn + [0]
i=i+1
imx=i
adm=[0.0]
for i in range(1,imx):
adm=adm+[0.0]
print "*"*40
print "xl,f"
for i in range(0,imx-1):
print float(xl[i]),"   ",float(f[i])
# set parameter ####################################
for n in range(0,1000):
###############################
for i in range(1,imx-1):
if u>=0.0:
ip = i - 1
dd = -dx
else:
ip = i + 1
dd = dx
xi = -u*dt
S  = (f[ip]-f[i])/dd
if S - f[ip] >0.0:
isgn = 1.0
else:
isgn = -1.0
TS = f[ip]-S
BB = ((abs(S-df[i])+eps)/(abs(S-df[ip])+eps)-1.0)/dd

a3 = (df[i]-S+(df[ip]-S)*(1.0+alpha*BB*dd))/(dd*dd)
a2 = S*alpha*BB+(S-df[i])/dd-a3*dd
a1 = df[i] +f[i]*alpha * BB
a0 = f[i]
fn[i]  = (a3*xi*xi*xi+a2*xi*xi+a1*xi+a0)/(1.0+alpha*BB*xi)
dfn[i] = (3.0*a3*xi*xi+2.0*a2*xi+a1)/(1.0+alpha*BB*xi)-(alpha*BB)*(fn[i])/(1.0+alpha*BB*xi)
#shift#########################################
for i in range(1,imx-1):
f[i]=fn[i]
df[i]=dfn[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 = 'RCIP')
plot(xl,y2,label = 'INITIAL')
plot(xl,y3,label = 'EXACT')
xlabel('x')
ylabel('f')
title('RCIP')
axis([-0,100.0,-0.2,1.2])
legend()
grid(True)
savefig('RCIP')
show()
###################################################


等間隔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()
###################################################

2013年10月16日水曜日

HartenYee

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

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

2013年10月10日木曜日

MUSCL


MUSCLをVBAで書いてみました。
Minmod関数でリミッターを掛けてます。

緑色はセル境界の内挿した値です。

河川ならこの程度の精度があれば十分だと解釈してます。
(個人的な意見です)

ただし、移流項に対してのみTVD条件を満たす(システム系で解かずに、運動方程式の移流項に対してのみTVD条件を満たす)場合、他の項による不連続部に対して、リミッターがかかってしまうのではないかと気になっています。

そうすると、熱移動と流れの数値解析や数値流体工学に議論されているハイブリット法が良いのか?と思いますが、アレは拡散項ありきなんですよね。河川だと拡散項が省略されることが多く、生成項が効いてきます。ハイブリッド法だと1次風上のみになってしまいます。

ここら辺の議論については非常に興味深いです。誰も議論しないのだろうか?いや、精度的には意味がないのか?

おそらく河川の流れの空間刻み幅は小さくしてもせいぜい数mなので、高次精度のスキームは必ず必要だと思ってます。

勉強します。