最近PDLを使おうかと思ってます。
いや、仕事ではないですが、せっかくperlの勉強をしたいと思っているので。
自由度が高いけど、書き手によってコードのエレガント差が変わる部分が大きい言語は結構好きかもしれません。
学習レベルによってコードが変わるのは、自分の成長度がわかっていいですね。
業務では嫌われる要因でしょうが・・・。
モダンPerlの新版はいつ出るんでしょうか。
結構楽しみにしてます。
あと、Perl入学式に参加して、
「書き方は自由です」
というのは、勇気づけられました。
もちろん、技巧をこらしても良い、ということではなく、
「汚く書いても動くし、技巧をこらしても動くけど、ちゃんと考えて書いてね」
と言われてる気がします。
それはいわば、自然言語に近いのかなと。
文章の流れを考えて、くどく「私は」「私は」と書くのは、
読みづらいし、大事な所で主語が抜けるのは、本当に意味が分からない。
話し言葉と文章では異なりますし。
そういった意味では、使い手のレベルに依存し、
高度なコードも、低度なコードも可能なので、自然言語に似ていると感じます。
さすが、ラリーさんです。
今後も発展して欲しいと思いますし、無くなることはないでしょう。
Linuxにはほとんど入ってますし。
死ぬまでにどれだけ理解出来るか、楽しみです。
2014年6月18日水曜日
カルマンフィルタ
足立先生の「カルマンフィルタの基礎」を勉強させていただきます。
初心者にとっては非常に良い本だと思ってます。
その中で、例題プログラムがあり、MATLABなのですが、
ビンボーなため、持っておりません。
そこでpythonで書いてみました。
最初の例題は、スカラー変数なので、若干弄っております。
コメントはしません。
色々おこられそうなので。。。
本当に勉強になります。
pl.show()
初心者にとっては非常に良い本だと思ってます。
その中で、例題プログラムがあり、MATLABなのですが、
ビンボーなため、持っておりません。
そこでpythonで書いてみました。
最初の例題は、スカラー変数なので、若干弄っております。
コメントはしません。
色々おこられそうなので。。。
本当に勉強になります。
import numpy as np
import pylab as pl
import random
def kf(A,B,Bu,C,Q,R,u,y,xhat,P):
xhatm = A*xhat+Bu*u
Am = A;Bm=B;Cm=C;I =1
Pm = A*P*Am + B*Q*Bm
G = Pm*C/(C*Pm*Cm+R)
xhat_new = xhatm + G*(y-C*xhatm)
P_new = (I-G*C)*Pm
return xhat_new,P_new,G
A = 1.
b = 1.
c = 1.
Q = 1.
R = 5.
N = 300
x = np.zeros(N)
y = np.zeros(N)
v = np.zeros(N)
w = np.zeros(N)
xhat = np.zeros(N)
t = np.zeros(N)
for i in range(N):
v[i] = random.normalvariate(0,Q)
w[i] = random.normalvariate(0,R)
t[i] = i
#print v[i],w[i]
y[0] = c*x[0] + w[0]
for k in range(1,N):
x[k] = A*x[k-1]+b*v[k-1]
y[k] = c*x[k] + w[k]
P = 0
xhat[0] = 0.0
for k in range(1,N):
xhat[k],P,G = kf(A,b,0,c,Q,R,0,y[k],xhat[k-1],P)
pl.plot(t,y,label='measured')
pl.plot(t,x,label='ture')
pl.plot(t,xhat,label='estimate')
pl.grid(True)
pl.legend()
pl.savefig('kf')
pl.show()
とりあえずうまくいってるようです。
2014年6月15日日曜日
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()
###################################################
もうちょっと美しく出来れば良いんですけど・・・。
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()
###################################################
2014年5月11日日曜日
2014年3月22日土曜日
Fortran
新しい会社では全く使わなくなったFortran＋本の話です。
「数値計算のためのFrotran90/95プログラム入門」は以前の会社に入って良く読みました。京都大学の牛島先生の本です。
特定の分野では77は使われています。
ゼロから書く場合は90/95へのシフトがなされているようです。
大学の時は77が多くて、
・Windowsかつ小さなプログラムならfcpad＋salfordで解く
・Macならg77で解く
だった気がします。
Fortranは科学技術計算に特化していますが、文法は分っても、「自分が書きたいプログラム」の書き方の参考になるものは少なかった気がします。77だとそんなこと気にすることもないですが・・・。common文は使わない、とかimplicit realは使わないくらいでしょうか。
当時、90/95の書籍はあったのですが、流体の計算をする場合の書き方みたいなものを提示してある書籍かつインストール方法やgnuplotを使った可視化まで書いてある本はあまりなかった記憶があります。この本は自習用として良いのではないかとすぐ購入した覚えがあります。
結果、モジュールやサブルーチンの使い方、プログラム規模に応じた使い分け、浅水流方程式の特性曲線法を使ったプログラムまで掲載されており、手垢塗れになりました。付箋もいっぱい貼ってあります。
77はGOTO文を使わないとLOOPを抜け出せない、文番号をつかったFORMAT文、1行72文字まで等の制限がありますが、90/95は自由形式になり、それなりに使いやすくなってます。配列も使いやすくなってます。
Fortranを使う人間は、あまりプログラムがメインでない技術者が多く、文法とか書き方にこだわりすぎても仕方ありません。Cだとメモリに気を使わないとダメですし。
その点、Fortranはちょうど良いかと思います。フリーのライブラリも多くありますし。
フリーのコンパイラならgfortran、商用ならインテルコンパイラやPGIコンパイラがあります。Linuxかつ非商用目的なら、インテルコンパイラはフリーで使えます。
あとFortran一番の特徴は計算が速いことです。これが利用される目的として一番大きいですかね。 時間の数値積分のための繰り返し回数が多くなると、計算終わるまでに時間がかかります。計算が速いのは重要です。
「数値計算のためのFrotran90/95プログラム入門」は以前の会社に入って良く読みました。京都大学の牛島先生の本です。
特定の分野では77は使われています。
ゼロから書く場合は90/95へのシフトがなされているようです。
大学の時は77が多くて、
・Windowsかつ小さなプログラムならfcpad＋salfordで解く
・Macならg77で解く
だった気がします。
Fortranは科学技術計算に特化していますが、文法は分っても、「自分が書きたいプログラム」の書き方の参考になるものは少なかった気がします。77だとそんなこと気にすることもないですが・・・。common文は使わない、とかimplicit realは使わないくらいでしょうか。
当時、90/95の書籍はあったのですが、流体の計算をする場合の書き方みたいなものを提示してある書籍かつインストール方法やgnuplotを使った可視化まで書いてある本はあまりなかった記憶があります。この本は自習用として良いのではないかとすぐ購入した覚えがあります。
結果、モジュールやサブルーチンの使い方、プログラム規模に応じた使い分け、浅水流方程式の特性曲線法を使ったプログラムまで掲載されており、手垢塗れになりました。付箋もいっぱい貼ってあります。
77はGOTO文を使わないとLOOPを抜け出せない、文番号をつかったFORMAT文、1行72文字まで等の制限がありますが、90/95は自由形式になり、それなりに使いやすくなってます。配列も使いやすくなってます。
Fortranを使う人間は、あまりプログラムがメインでない技術者が多く、文法とか書き方にこだわりすぎても仕方ありません。Cだとメモリに気を使わないとダメですし。
その点、Fortranはちょうど良いかと思います。フリーのライブラリも多くありますし。
フリーのコンパイラならgfortran、商用ならインテルコンパイラやPGIコンパイラがあります。Linuxかつ非商用目的なら、インテルコンパイラはフリーで使えます。
あとFortran一番の特徴は計算が速いことです。これが利用される目的として一番大きいですかね。 時間の数値積分のための繰り返し回数が多くなると、計算終わるまでに時間がかかります。計算が速いのは重要です。
2014年2月10日月曜日
関数プログラミング
プログラミングの書き方はいくつかあると思いますが、
私のこれまでの感覚だと、手続き型とオブジェクト指向型がメインだったように思います。
最近の書籍、雑誌、ネットでは「関数プログラミング」というキーワードが多くあるとおもいます。
私が関数プログラミングに触れ始めたのはEmacsからです。Emacsを使いこなそうとすると、Elispを使わなければなりません。いや、使わなくてもEmacsは使えるんですが、これがないなら他のエディタを使った方が良いです。
Elispを見て、それまでの命令型と全く違って面食らいました・・・。
かっこまみれは良いとして、関数型の思想を勉強してないので、頭が大混乱・・・。
そこから関数型言語としてCommonLispとか、Schemeとかあるんだなと。
ハッカーと画家も買って読みました。読み物として面白いです。
が、その威力が未だに見えてきません。変数書き換え不可！なぜ？
ある程度大規模な科学技術計算（私は流体解析メイン）はFortran、C/C++で書かれることが多く、空間離散点が数百万、時間積分も長時間やる場合は、Intel、PGIのコンパイラが多いようです。
中規模であれば、javaやC#なんかも使われているようです。
小規模であったり、統計解析や定常計算等であれば、Excel、Python、Ruby、R等が使われることが多いようです。
さて、自分のフィールドである科学技術計算に関数プログラミングを使えないかと考えたのですが、なかなか難しいです・・・。
その理由の第一は関数プログラミングを十分理解出来ていないことです。これは本読んで、コーディングすれば少しずつ解決しそうです。
第二に関数プログラミングが適用されている他のフィールドで求められている機能について理解出来ていないことです。これはなかなか難しそうですが、第一の理由を解決する過程で何か答えの断片を得られそうです。
その上で、第一、第二を踏まえて科学技術計算にどのように適用するのかを考えようと思います。
一つの機能としては、並列化への適用は可能性がありそうです。科学技術計算をやる上で、「解析時間の短縮」は大きなウェートです。
可読性向上については、厳しいような気がします。計算屋は命令型しか勉強しないですし、オブジェクト指向すら知らない人が多いです。実際、オブジェクト指向で書くと、遅くなることが多いです。
あとは、保守管理、過去の遺産の利用等々、いくつかの側面についても検討する必要があると思います。
ちなみに、関数プログラミングについては、Twitterで、おそらくその道のプロ（NY先生）であろう方からの返信があって、考え始めました。資料もご提示いただき、少し頭が晴れてきました。このページを見ているとは思いませんが、お礼申し上げます。ありがとうございます。
また、考えがまとまったら書こうと思います。
私のこれまでの感覚だと、手続き型とオブジェクト指向型がメインだったように思います。
最近の書籍、雑誌、ネットでは「関数プログラミング」というキーワードが多くあるとおもいます。
私が関数プログラミングに触れ始めたのはEmacsからです。Emacsを使いこなそうとすると、Elispを使わなければなりません。いや、使わなくてもEmacsは使えるんですが、これがないなら他のエディタを使った方が良いです。
Elispを見て、それまでの命令型と全く違って面食らいました・・・。
かっこまみれは良いとして、関数型の思想を勉強してないので、頭が大混乱・・・。
そこから関数型言語としてCommonLispとか、Schemeとかあるんだなと。
ハッカーと画家も買って読みました。読み物として面白いです。
が、その威力が未だに見えてきません。変数書き換え不可！なぜ？
ある程度大規模な科学技術計算（私は流体解析メイン）はFortran、C/C++で書かれることが多く、空間離散点が数百万、時間積分も長時間やる場合は、Intel、PGIのコンパイラが多いようです。
中規模であれば、javaやC#なんかも使われているようです。
小規模であったり、統計解析や定常計算等であれば、Excel、Python、Ruby、R等が使われることが多いようです。
さて、自分のフィールドである科学技術計算に関数プログラミングを使えないかと考えたのですが、なかなか難しいです・・・。
その理由の第一は関数プログラミングを十分理解出来ていないことです。これは本読んで、コーディングすれば少しずつ解決しそうです。
第二に関数プログラミングが適用されている他のフィールドで求められている機能について理解出来ていないことです。これはなかなか難しそうですが、第一の理由を解決する過程で何か答えの断片を得られそうです。
その上で、第一、第二を踏まえて科学技術計算にどのように適用するのかを考えようと思います。
一つの機能としては、並列化への適用は可能性がありそうです。科学技術計算をやる上で、「解析時間の短縮」は大きなウェートです。
可読性向上については、厳しいような気がします。計算屋は命令型しか勉強しないですし、オブジェクト指向すら知らない人が多いです。実際、オブジェクト指向で書くと、遅くなることが多いです。
あとは、保守管理、過去の遺産の利用等々、いくつかの側面についても検討する必要があると思います。
ちなみに、関数プログラミングについては、Twitterで、おそらくその道のプロ（NY先生）であろう方からの返信があって、考え始めました。資料もご提示いただき、少し頭が晴れてきました。このページを見ているとは思いませんが、お礼申し上げます。ありがとうございます。
また、考えがまとまったら書こうと思います。
2014年2月9日日曜日
KKスキーム
河村・桑原スキームについてです。
1984年のLESに関する論文だったかと思います。
数値解析のうち、移流項の解法は以下の流れで勉強しました（私の場合）。
・中央差分（必ず発散）
・1次精度風上差分（解けるけど解が鈍る）
→中央差分の方が精度が良さそうだが解けない。なぜなら移流項（対流項）は移流速度に乗って情報が上流から下流に伝搬するからである。1次精度風上差分の場合、意図しない所で、拡散的な項（数値粘性）が入っているため解が鈍る。
・LaxWendroff法（2次精度だが振動が入る）
・MacCormack法（2次精度だが振動が入る）
→1次精度風上差分では数値粘性が入っていたため、安定的に解けた。同様の考え方で人工粘性を入れることで、振動を小さくする。
・KKスキーム
→3次精度を実現するために4次精度の中央差分＋4次の拡散項を入れる。振動は残るが、河村先生の論文以降、実現象解明へ適用されてきた実績がある。
以降、QUICK、TVD（MUSCL系、Non-MUSCL系）、CIP（オリジナル、有理関数補間、CSL系）の勉強をしました。個人的にはCIP-CSLR1が好きですが、既存のコードの書き換えは結構大変です。マルチモーメントなので・・・。
その後は、小松先生の6Point Scheme、牛島先生のQSIスキーム、ENO、WENO等もかじってます。
今回、気が向いて、KKスキームについてちょっと勉強し直しました。
ベンチマークで使われる2D-Cavity Flowに適用しようと考えてます。
非圧縮性流体の解析なので、MAC系の解法を使って、移流項にKKスキームをと・・。
KKスキームは非保存系のスキームでした。ということで、ある論文の方法を使って保存系に変形しました。線形移流方程式に対して、保存系、非保存系で解いた結果、同様の結果を得られたので、適用は出来そうです。
うまくいったら、結果を載せようと思います。
1984年のLESに関する論文だったかと思います。
数値解析のうち、移流項の解法は以下の流れで勉強しました（私の場合）。
・中央差分（必ず発散）
・1次精度風上差分（解けるけど解が鈍る）
→中央差分の方が精度が良さそうだが解けない。なぜなら移流項（対流項）は移流速度に乗って情報が上流から下流に伝搬するからである。1次精度風上差分の場合、意図しない所で、拡散的な項（数値粘性）が入っているため解が鈍る。
・LaxWendroff法（2次精度だが振動が入る）
・MacCormack法（2次精度だが振動が入る）
→1次精度風上差分では数値粘性が入っていたため、安定的に解けた。同様の考え方で人工粘性を入れることで、振動を小さくする。
・KKスキーム
→3次精度を実現するために4次精度の中央差分＋4次の拡散項を入れる。振動は残るが、河村先生の論文以降、実現象解明へ適用されてきた実績がある。
以降、QUICK、TVD（MUSCL系、Non-MUSCL系）、CIP（オリジナル、有理関数補間、CSL系）の勉強をしました。個人的にはCIP-CSLR1が好きですが、既存のコードの書き換えは結構大変です。マルチモーメントなので・・・。
その後は、小松先生の6Point Scheme、牛島先生のQSIスキーム、ENO、WENO等もかじってます。
今回、気が向いて、KKスキームについてちょっと勉強し直しました。
ベンチマークで使われる2D-Cavity Flowに適用しようと考えてます。
非圧縮性流体の解析なので、MAC系の解法を使って、移流項にKKスキームをと・・。
KKスキームは非保存系のスキームでした。ということで、ある論文の方法を使って保存系に変形しました。線形移流方程式に対して、保存系、非保存系で解いた結果、同様の結果を得られたので、適用は出来そうです。
うまくいったら、結果を載せようと思います。
2014年2月4日火曜日
2014年1月24日金曜日
河川のシミュレーション!
河村先生の書籍です。
現在は絶版になっております。
大学の時に購入しました。
解析の基本的な説明の後、1次元、2次元、3次元の解析について記述されています。
偏微分方程式については移流拡散方程式、ポアソン方程式の解法を、
また一般座標系の説明があります。
河川に焦点を当てた（浅水流方程式に焦点を当てた）書籍は他にないような気がします。
サンプルプログラムもありますが、実用性は低いように思います。
また、流砂量公式については記述がありません。
ただし、勉強し始めた人にとっては有用な本の一つであることは間違いないと感じました。最近は大学レベルでも1からプログラムを書くことは少なくなっているようですし、
あくまでソフロウェアのユーザーであっても、知っておいた方が良いように思います。
個人的にはKKスキームについて記述されていないのが、ちょっと残念です。
日本人が作った数少ない移流項の解法ですからね。
（河村先生のその他の本には記述されています。河川だから省いたのかなと邪推しています）
日本の技術書は薄い本が多いように感じます。
ゼロから勉強しようにも、技術のバックグラウンドや基礎方程式について記述されているだけの本が多いです。基本的には基礎方程式は複雑なことが多く、複雑な方程式を解くには数値計算に頼らざるを得ないことが多いです。そのためには数値解析の基本的事項およびそれを実用的に使う方法まで書くには厚くすべきだと思ってます。
その点、海外の書籍は厚く、読むのが大変ですが、独学に向いています。
日本でもそのような流れにしていただきたいです。
ちなみに河村先生は多くの書籍を出されています。数値計算の基本的な所を記述した教科書になるものが多いです。私はファンなので買ってしまいます・・・。
エクセルVBAを使った解析について書かれているものもありますので、興味のある人は是非！
現在は絶版になっております。
大学の時に購入しました。
解析の基本的な説明の後、1次元、2次元、3次元の解析について記述されています。
偏微分方程式については移流拡散方程式、ポアソン方程式の解法を、
また一般座標系の説明があります。
河川に焦点を当てた（浅水流方程式に焦点を当てた）書籍は他にないような気がします。
サンプルプログラムもありますが、実用性は低いように思います。
また、流砂量公式については記述がありません。
ただし、勉強し始めた人にとっては有用な本の一つであることは間違いないと感じました。最近は大学レベルでも1からプログラムを書くことは少なくなっているようですし、
あくまでソフロウェアのユーザーであっても、知っておいた方が良いように思います。
個人的にはKKスキームについて記述されていないのが、ちょっと残念です。
日本人が作った数少ない移流項の解法ですからね。
（河村先生のその他の本には記述されています。河川だから省いたのかなと邪推しています）
日本の技術書は薄い本が多いように感じます。
ゼロから勉強しようにも、技術のバックグラウンドや基礎方程式について記述されているだけの本が多いです。基本的には基礎方程式は複雑なことが多く、複雑な方程式を解くには数値計算に頼らざるを得ないことが多いです。そのためには数値解析の基本的事項およびそれを実用的に使う方法まで書くには厚くすべきだと思ってます。
その点、海外の書籍は厚く、読むのが大変ですが、独学に向いています。
日本でもそのような流れにしていただきたいです。
ちなみに河村先生は多くの書籍を出されています。数値計算の基本的な所を記述した教科書になるものが多いです。私はファンなので買ってしまいます・・・。
エクセルVBAを使った解析について書かれているものもありますので、興味のある人は是非！
履歴
年明け一発目です。
せっかくなので、これまでのプログラミング履歴をまとめておこうと思います。
順番は適当です。
せっかくなので、これまでのプログラミング履歴をまとめておこうと思います。
順番は適当です。
- C/C++
- Fortran
- python
- ruby
- perl
- xojo
- VBA
- VB.NET
- java
- ELisp
- AWK
2013年12月18日水曜日
C++
C++はマルチパラダイム言語ですが、
最近、再勉強をしています。
楽しいです。とりあえず、微分方程式を解くためのアルゴリズムを中心に書いていますが、本当に色々な書き方があるんだなと。
本も数冊購入しました。一冊は本当に初学者の本です。なぜなら、C++の全体のイメージを思い出すことが必要だと思ったからです。新しい発見でしたが、比較的初学者の本を読んで、ちょっと説明が違うな、と認識できるところがありました。少しは成長していたんだなと。
別の一冊はSTL中心です。STLは本当にコードが短くなるだけでなく、デバッグすべき項目を減らせるなという印象です。世の中の、どの分野で、どの程度使われているかわかりませんが・・・。加えてBoostの知識も入れておこうと思います。
最後の一冊はWindowsに特化したMFC中心の本です。
今さら感が強いですが・・・。
VC++はC++とは全然違うものというのが再認識されました・・・。
ライブラリを含めると、C++という言語はあまりにも広すぎます。
プラットフォームに依存したVC++だと別物です。
なかなか学習の仕方が難しいですが、今のところ、やりたいことは微分方程式を解くことなので、いろいろな書き方を試したいと思います。
最近、再勉強をしています。
楽しいです。とりあえず、微分方程式を解くためのアルゴリズムを中心に書いていますが、本当に色々な書き方があるんだなと。
本も数冊購入しました。一冊は本当に初学者の本です。なぜなら、C++の全体のイメージを思い出すことが必要だと思ったからです。新しい発見でしたが、比較的初学者の本を読んで、ちょっと説明が違うな、と認識できるところがありました。少しは成長していたんだなと。
別の一冊はSTL中心です。STLは本当にコードが短くなるだけでなく、デバッグすべき項目を減らせるなという印象です。世の中の、どの分野で、どの程度使われているかわかりませんが・・・。加えてBoostの知識も入れておこうと思います。
最後の一冊はWindowsに特化したMFC中心の本です。
今さら感が強いですが・・・。
VC++はC++とは全然違うものというのが再認識されました・・・。
ライブラリを含めると、C++という言語はあまりにも広すぎます。
プラットフォームに依存したVC++だと別物です。
なかなか学習の仕方が難しいですが、今のところ、やりたいことは微分方程式を解くことなので、いろいろな書き方を試したいと思います。
2013年12月12日木曜日
全然書けてない・・・。
ようやく落ち着いたので、書き始めようと思います。
次の会社に向けて、C++の復習と統計の勉強を少しずつしています。
・C++について
プログラミング言語の学習は文法をマスターすることではありません。
その言語がどのような思想で設計され、与えられる課題に対して、どのようにアプローチすべきか？というような「思考方法」を学ぶことだと捉えています。
C++のようないろんな書き方がある場合、つまりは思考方法が複数あることを意味しています。そのような混沌としたような言語は非常に難しく感じます。
ただし、それが非常に面白いです。
結局、C言語の思考方法に振り返ることになってます（笑）
でもコーディングする人間が忘れてはならないことは、
「このコードが何をするために使われるか？」
だと思います。
・統計
「統計学は最・・・」とかいう本が出てますが、読んでません。
個人的には、統計学で説明できる事象とは何か？に興味があります。
特に時系列分析やスペクトル解析でしょうか。
統計学はあくまで考え方の一つであり、統計学に准ずれば正解、ということではありません。如何せん、現象論は考えていないからです。ただ、現象としては捉えられるけど、メカニズムは不明、という状況において、何らかしら定量的なアプローチをする方法としては有効であるように思います。すべての現象をニュートン力学や微分方程式で表現できる訳ではない（いや、ほとんどは出来るんですが、その場合、時空間解像度を相当確保しなければなりません）ですし、表現できたとしても、実現象を再現できるほど精度がある訳ではありません。
そういった現象と結果を結びつけたり、そのためのツールとして非常に有効な学問だと思います。平均値、相関、標準偏差とかです。
でも、中学や高校で確率の勉強が、工場の品質管理に使われるなんて想像もつきませんでした（笑）。鵜呑みにして勉強した結果、あ、そういうことか、と納得したのが大学です。
さらに言えば、水文学なんかもそうです。母数推定や関数へのフィッティング等。。。
怖いのは、ツールを使えばなんか答えがでる、ということです。
勉強していない人間がテキトーなことをして、定量的に評価することほど怖いことはありません。
統計学にはそういう恐ろしさも孕んでます。
次の会社に向けて、C++の復習と統計の勉強を少しずつしています。
・C++について
プログラミング言語の学習は文法をマスターすることではありません。
その言語がどのような思想で設計され、与えられる課題に対して、どのようにアプローチすべきか？というような「思考方法」を学ぶことだと捉えています。
C++のようないろんな書き方がある場合、つまりは思考方法が複数あることを意味しています。そのような混沌としたような言語は非常に難しく感じます。
ただし、それが非常に面白いです。
結局、C言語の思考方法に振り返ることになってます（笑）
でもコーディングする人間が忘れてはならないことは、
「このコードが何をするために使われるか？」
だと思います。
・統計
「統計学は最・・・」とかいう本が出てますが、読んでません。
個人的には、統計学で説明できる事象とは何か？に興味があります。
特に時系列分析やスペクトル解析でしょうか。
統計学はあくまで考え方の一つであり、統計学に准ずれば正解、ということではありません。如何せん、現象論は考えていないからです。ただ、現象としては捉えられるけど、メカニズムは不明、という状況において、何らかしら定量的なアプローチをする方法としては有効であるように思います。すべての現象をニュートン力学や微分方程式で表現できる訳ではない（いや、ほとんどは出来るんですが、その場合、時空間解像度を相当確保しなければなりません）ですし、表現できたとしても、実現象を再現できるほど精度がある訳ではありません。
そういった現象と結果を結びつけたり、そのためのツールとして非常に有効な学問だと思います。平均値、相関、標準偏差とかです。
でも、中学や高校で確率の勉強が、工場の品質管理に使われるなんて想像もつきませんでした（笑）。鵜呑みにして勉強した結果、あ、そういうことか、と納得したのが大学です。
さらに言えば、水文学なんかもそうです。母数推定や関数へのフィッティング等。。。
怖いのは、ツールを使えばなんか答えがでる、ということです。
勉強していない人間がテキトーなことをして、定量的に評価することほど怖いことはありません。
統計学にはそういう恐ろしさも孕んでます。
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()
###################################################
河川の流れを対象とし、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」は０〜１００の格子番号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()
###################################################
晒すことにします。
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()
###################################################
2013年10月10日木曜日
MUSCL
MUSCLをVBAで書いてみました。
Minmod関数でリミッターを掛けてます。
緑色はセル境界の内挿した値です。
河川ならこの程度の精度があれば十分だと解釈してます。
（個人的な意見です）
ただし、移流項に対してのみTVD条件を満たす（システム系で解かずに、運動方程式の移流項に対してのみTVD条件を満たす）場合、他の項による不連続部に対して、リミッターがかかってしまうのではないかと気になっています。
そうすると、熱移動と流れの数値解析や数値流体工学に議論されているハイブリット法が良いのか？と思いますが、アレは拡散項ありきなんですよね。河川だと拡散項が省略されることが多く、生成項が効いてきます。ハイブリッド法だと1次風上のみになってしまいます。
ここら辺の議論については非常に興味深いです。誰も議論しないのだろうか？いや、精度的には意味がないのか？
おそらく河川の流れの空間刻み幅は小さくしてもせいぜい数mなので、高次精度のスキームは必ず必要だと思ってます。
勉強します。
2013年9月29日日曜日
空間２次精度
空間２次精度についての話です。
河川の流れであれば、
1次精度風上差分では不十分でも、せいぜい２次精度があれば十分なようです。
ダムブレークのような問題を解く必要はあまりないですから。
２次精度だとベースをオリジナルのLaxWendroff法（2STEPではない）を用い、FluxLimiterを施すことで、大局２次精度を維持できます。
CIP法を用いた非保存系のスキームも使われているようですが、既存のプログラムがある場合、書き換えるのが煩雑です。
また、CIP法だとわずかに数値振動が残りますがTVDだと1次精度となって振動を抑えます。
TVDだとMUSCLも使われているようですが、TVD-LaxWendroffは使われることは少ないようです。
精度の高い乱流モデルを後々組み込むことを考えると、MUSCLという選択もわるくないと思います。
以上の話は、ナビエストークス方程式（or レイノルズ方程式）のうちの移流項に対する解法であり、
運動方程式、連続式をシステム系としてとらえた解法ではありません。なぜなら、生成項となる摩擦項や地形勾配項があるためであり、
C-Propertyを満たすことが難しいためです。
アカデミックの世界でも研究されている方はいるようですが、マニアックなまでに研究されている方はいないようです。
すでに解決された問題だと捉えられているのでしょう。
まあ、分っている人は少ないと思いますが・・・。
流体の数値計算法について（２）
藤井先生の本についての続編です。
一つ面白い記述に
「最近の多くの実用的な計算結果を見る限りでは、一般的には2次精度を維持することが必要といえる」
というのがあります。
乱流の解析では2次精度では不十分な場合がありますが、多くは2次精度程度である程度事足りるとの指摘です。
河川の流れの解析のように、平面2次元流れの解析であれば、この指摘は当てはまると思われます。
平面2次元流れは、どの程度の流れの解像度を得ることが必要かで乱流モデルを選択する必要がありますが、多くの場合は、ブシネスク近似を使って、ゼロ方程式を用いたモデルが多いようです。
①格子形状
格子形状は境界適合座標系を使い、スカラー量とベクトル量を互い違いに配置するスタッガード格子を採用します。
②解析手法
有限体積法を用いることが多いですが、モデルによっては有限差分法が多いです。有限要素法は滅多にありませんが、アカデミックな世界では実施されていることもあります。
また、移流項に対しては、1次風上差分が非常に多いはずです。藤井先生の本を読まれる方は「は？」と思われるかもしれませんが、実用的な解析では、計算が飛ばないことが非常に重要です。
はっきり言って、解析についてあまり知らなくても、計算結果を確実に得ることが求められます。２次精度以上を確保しようとすると、不安定になる可能性が多くなりますので、1次風上が多いのです。
一方、移流項以外は、２次精度の中央差分が多いです。各項の大きさについての議論は別にして、移流項で3次精度以上を確保しようとすると、他の項も同様に3次精度以上を確保する必要があります。
したがって、移流項に対しては大局的空間２次精度程度を確保することが適切のようです。
さて、藤井先生の本に戻りますが、前述の指摘がある程度的を得ています。
分野が違えど、面白い指摘だなと思った次第です。
なお、ブシネスク近似のうちk−εやLESを用いる場合はこの指摘は当てはまりません。
乱流モデルは、小生の間隔では3次精度以上が必要のようです。
一つ面白い記述に
「最近の多くの実用的な計算結果を見る限りでは、一般的には2次精度を維持することが必要といえる」
というのがあります。
乱流の解析では2次精度では不十分な場合がありますが、多くは2次精度程度である程度事足りるとの指摘です。
河川の流れの解析のように、平面2次元流れの解析であれば、この指摘は当てはまると思われます。
平面2次元流れは、どの程度の流れの解像度を得ることが必要かで乱流モデルを選択する必要がありますが、多くの場合は、ブシネスク近似を使って、ゼロ方程式を用いたモデルが多いようです。
①格子形状
格子形状は境界適合座標系を使い、スカラー量とベクトル量を互い違いに配置するスタッガード格子を採用します。
②解析手法
有限体積法を用いることが多いですが、モデルによっては有限差分法が多いです。有限要素法は滅多にありませんが、アカデミックな世界では実施されていることもあります。
また、移流項に対しては、1次風上差分が非常に多いはずです。藤井先生の本を読まれる方は「は？」と思われるかもしれませんが、実用的な解析では、計算が飛ばないことが非常に重要です。
はっきり言って、解析についてあまり知らなくても、計算結果を確実に得ることが求められます。２次精度以上を確保しようとすると、不安定になる可能性が多くなりますので、1次風上が多いのです。
一方、移流項以外は、２次精度の中央差分が多いです。各項の大きさについての議論は別にして、移流項で3次精度以上を確保しようとすると、他の項も同様に3次精度以上を確保する必要があります。
したがって、移流項に対しては大局的空間２次精度程度を確保することが適切のようです。
さて、藤井先生の本に戻りますが、前述の指摘がある程度的を得ています。
分野が違えど、面白い指摘だなと思った次第です。
なお、ブシネスク近似のうちk−εやLESを用いる場合はこの指摘は当てはまりません。
乱流モデルは、小生の間隔では3次精度以上が必要のようです。
2013年9月16日月曜日
書籍
「プログラミング言語AWK」
AWKの本です。最近、USP Labから出版されたものです。
もちろん著者はAWKの開発者であり、AWK版のK&Rと考えてもらえれば良いとおもいます。
ただ、出版といっても、絶版しては再販されている書籍であり、
需要がある書籍と言えましょう。
読んだ感想ですが、古くさいです。
ただ、K&Rよりも実用的であると感じます。
スクリプト言語であり、ワンライナーによるフィルタをかけるだけであれば、序盤を読めばかなり書けるようになると思います。
まあ、敢えてAWKを使う人はあまりいないと思います。今であれば、Perl、Python、Rubyを使えば比較的簡単にできます。
それがテキスト処理だとしてもです。
ただし、プログラミングの歴史を学ぶ点においては、知っておいて損はないかと。
加えて、「プログラミングが、なぜ、今の形になったのか？」ということを理解する一つの助けとなる気がします。
やはり、周辺知識は独立して使うことが出来るものではないですが、メインの知識と組み合わさって威力を発揮するものだと思います。
等間隔QUICKEST
# 久々にpython
for i in range(3,imx-3):#0~imx-1
c = u*dt/dx #
if u >0:
fn[i] = f[i] - c/6*(2*f[i+1]+3*f[i]-6*f[i-1]+f[i-2]) \
+ c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
- c*c*c/6 *(f[i+1]-3*f[i]+3*f[i-1]-f[i-2])
else:
fn[i] = f[i] - c/6*(2*f[i-1]+3*f[i]-6*f[i+1]+f[i+2]) \
+ c*c/2*(f[i+1]-2*f[i]+f[i-1]) \
- c*c*c/6 *(f[i-1]-3*f[i]+3*f[i+1]-f[i+2])
等間隔にしか使えないですが・・・。
一応、保存則を使っています。
明日中には不等間隔＋TVD化されたQUICKを晒したいです。
登録:
投稿 (Atom)