最近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スキームは非保存系のスキームでした。ということで、ある論文の方法を使って保存系に変形しました。線形移流方程式に対して、保存系、非保存系で解いた結果、同様の結果を得られたので、適用は出来そうです。
うまくいったら、結果を載せようと思います。
登録:
投稿 (Atom)