【Navier-Stokes方程式に対する数値計算手法】サイトマップ
はじめに
Navier-Stokes方程式を数値的に解く手法には、非常に多くのバリエーションがあり、どれから勉強してよいのか?何が違うのか?と悩む方も多いことでしょう。私もその一人でした。そこで、Navier-Stokes方程式に対する数値計算手法の本を山ほど買い込みだんだん勉強して参りました。その成果は当ブログの記事となっています。それらの記事をまとめたのがこのページです。
Navier-Stokes方程式に対する数値計算手法には、大きく分けるとMAC法系列の解法とSIMPLE法系列の解法があります。なのでMAC法とSIMPLE法を理解すれば、他の方法は簡単に理解できると思います。まだ書いている途中ですが、だんだんアップデートしていくのでよろしくお願いします。
MAC法
全ての始まりです。
フラクショナルステップ法
SIMPLE法
完全陰解法
SIMPLEC法
SIMPLE法を改良した方法です。
SIMPLER法
SIMPLE法を改良した方法です。
SIMPLEST法
SIMPLE法を改良した方法です。
PISO法
流れ関数ー渦度法
勉強中です。二次元の流れにしか適用できません。計算機の性能が今ほど優れていなかった時代によく用いられていました。当時の教科書には大体載っていましたが、現在では論文や本で見ることはほとんどありません。
文献まとめ
文献まとめを作成中です。なるべく網羅的なものにしたいです。
- 作者: 河村哲也
- 出版社/メーカー: 朝倉書店
- 発売日: 2014/03/26
- メディア: 単行本
- この商品を含むブログ (1件) を見る
- 作者: 桑原邦郎,河村哲也
- 出版社/メーカー: 朝倉書店
- 発売日: 2005/03/01
- メディア: 単行本
- クリック: 6回
- この商品を含むブログを見る
- 作者: 河村哲也,平野広和,池川昌弘,川原睦人,登坂宣好,数値流体力学編集委員会
- 出版社/メーカー: 東京大学出版会
- 発売日: 1995/02
- メディア: 単行本
- この商品を含むブログを見る
【英語】"Sapiens: A Brief History of Humankind"『サピエンス全史』【9冊目】
どんな本か?
Yuval Noah Harari 著 "Sapiens: A Brief History of Humankind" を読み終わりました。日本語の題は『サピエンス全史』です。ホリエモンも絶賛していました。この本は、人類はどうして発展することができたのか、そして今後人類はどうなっていくのか、を歴史上の例をふんだんに用いながら解説した名著です。歴史にとどまらず、あらゆる分野の知識を縦横無尽に使いこなしており、Jared Diamond の "Guns, Germs, and Steel: The Fates of Human Societies" を読んだ時のような知的興奮を覚えました。
Sapiens: A Brief History of Humankind
- 作者: Yuval Noah Harari
- 出版社/メーカー: Vintage
- 発売日: 2015/04/30
- メディア: ペーパーバック
- この商品を含むブログ (4件) を見る
- 作者: ユヴァル・ノア・ハラリ,柴田裕之
- 出版社/メーカー: 河出書房新社
- 発売日: 2016/09/08
- メディア: 単行本
- この商品を含むブログ (50件) を見る
- 作者: ユヴァル・ノア・ハラリ,柴田裕之
- 出版社/メーカー: 河出書房新社
- 発売日: 2016/09/08
- メディア: 単行本
- この商品を含むブログ (24件) を見る
- 作者: ユヴァル・ノア・ハラリ
- 出版社/メーカー: 河出書房新社
- 発売日: 2016/09/16
- メディア: Kindle版
- この商品を含むブログ (15件) を見る
Guns, Germs, and Steel: The Fates of Human Societies
- 作者: Jared Diamond
- 出版社/メーカー: W W Norton & Co Inc
- 発売日: 2017/03/07
- メディア: ペーパーバック
- この商品を含むブログを見る
要点
3つの革命によりホモ・サピエンスが発展してきたことが解説されます。3つの革命とは、「認知革命」、「農業革命」、「科学革命」です。認知革命により、ホモ・サピエンスは虚構を信じることができるようになり、互いに協力することができるようになりました。その例が、国家、法人、通貨、宗教などです。どれも実際には存在していないものです。農業革命により、生産に余剰が生じ、その余剰を充てることにより、専門的な職業が誕生しました。また、人口が増加しました。科学革命により、ホモ・サピエンスは自らの力を増大させることができることに気付きました。進歩と成長により資本主義が生まれます。科学、資本主義、帝国の相互作用によりホモ・サピエンスは発展し、その結果として現在の我々の世界があります。そして、今後ホモ・サピエンスは、この力(遺伝子工学)を使ってホモ・サピエンスをホモ・サピエンスでない新たな人、superhuman(超人)を生み出していくことが示唆され、この本は終わります。
英文の難易度
原著はヘブライ語とのことですが、英語版も著者が書いているそうです。著者はヘブライ大学の教授で歴史学者です。もちろん校閲は入っているでしょうが、ノン・ネイティブなのにこんなに流暢な英文を書けるとは驚きです。英文の構文は素直で読みやすいです。ただ、単語のレベルは非常に高いです。現在私の語彙力は20,000前後あるはずですが、それでも1ページあたり3語ほど未知語が出てきます。しかし、逆に言えば単語の勉強にもなるということです。内容が面白いことは言うまでもありません。どんな単語が出てくるかについては、今後別の記事でまとめる予定です。
おわりに
Yuval Noah Harari 著 "Sapiens: A Brief History of Humankind" を読みました。今は Jared Diamond "Upheaval: How Nations Cope with Crisis and Change" を読んでいます。日本語訳はまだないのですが『激変』とすればよいでしょうか。この本は、フィンランド、チリ、日本、ドイツ、オーストラリア、インドネシア、アメリカの7つの国を取り上げ、それぞれの国がどのように危機的な状況を乗り越えたかが分析されます。今後の世界を生きていくうえで非常に有益な内容です。どうも、「選択的に変化すること」がキーワードですね。こちらも読み終えたらまとめていきます。
Upheaval: How Nations Cope with Crisis and Change
- 作者: Jared Diamond
- 出版社/メーカー: Allen Lane
- 発売日: 2019/05/07
- メディア: ペーパーバック
- この商品を含むブログを見る
【英語】『ハイレベル実戦英文法』Section 1
はじめに
最近英語を受験生に教える機会があったのですが、あまりに英文法問題がわからなくて恥ずかしくなりました。なので、英文法を復習していこうと思います!今回から手元にあった、『ハイレベル実戦英文法』という本を1 Sectionずつ読んでいこうと思います。
- 作者: 猿谷宣弘
- 出版社/メーカー: ベレ出版
- 発売日: 2017/08/19
- メディア: 単行本
- この商品を含むブログを見る
この本は、「英語の新聞を読むことができるようになるための英語の文法事項をまとめたもの」とのことなのでちょうどよいです。私たちが忘れている文法事項も出てくると思います。
Section 1 -ing の用法 進行形、名詞化、形容詞化
ingの用法が例文とともにまとめてあります。気になった点を例文とともにまとめておきます。
- 代名詞が動名詞の意味上の主語になる場合には所有格でも目的格でもよい
なんとなく所有格しか認識していませんでした。目的格でもよいのですね。つまり
What are the advantages and disadvantages of his not keeping a driver's license?
What are the advantages and disadvantages of him not keeping a driver's license?
のようにhisでもhimでもよいのです。
forの後は普通の名詞がくるものだと思い込んでいました。つまり
People often blame parents for their children's being rude in a public place.
のように所有格にできるんですね。
おわりに
自分は英文法を理解していると思い込んでいたことがわかりました…頑張って復習します!!!
【数値流体力学】勾配を制限するか、それとも流束を制限するか
MUSCL法など勾配制限関数(slope limiter)を用いる方法は、再構築した単調な値を用いて数値流束を計算するので単調な流束(1次風上差分、HLLなど)を使う必要があります。一方、高次の流束を使う場合は、解の単調性を保つために、流束自体を流束制限関数で制限する必要があります。
【有限体積法】1次元移流方程式をMUSCL法で解く C++コード付き
今回は1次元移流方程式を1次の風上差分+MUSCL法で解いていきます。MUSCL法とは、セル内で一定だった未知数のプロファイルを1次以上の関数で補外しセル界面での値を計算し、その値を用いて数値流速を計算することによって高次精度を達成する方法です。その際に、数値振動が発生しないように勾配制限関数(slope limiter)が用いられます。今回はminmod関数を使いました。
1次元移流方程式とは
のような一階の偏微分方程式のことをいいます。ここで、 は溶質の濃度、 は風などによる移流速度です。
今回の計算条件は、移流速度が正 で、初期条件は領域の左端の矩形波(長方形の波)です。比較のために普通の風上差分で解いた結果も載せます。
これが計算結果です。青がMUSCLで、赤が1次の風上差分です。
MUSCLのほうが初期プロファイルをよく保っているのがわかります。また、1次の風上差分の数値粘性の大きさもわかります。こんなに拡散されているのですね。
では、計算コードです。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <string> #include <sstream> using namespace std; //flux for advection equation// inline double Fluxadv(double V, double uL, double uR) { return max(V,0.0)*uL-max(-V,0.0)*uR; } //flux limiter minmod// inline double minmod(double d1, double d2) { if(d1*d2<=0) { return 0.0; } else if(abs(d1)<abs(d2)) { return d1; } else { return d2; } } //reconstruction// inline void reconst(int Node, double dt, double dx, double u[], double uiL[], double uiR[]) { //inner cells// for(int i=1;i<Node-1;i++) { //jump// double dul=u[i]-u[i-1]; double dur=u[i+1]-u[i]; //slope with slope limiter// double k=0.0; double b=(3.0-k)/(1.0-k); double dulLim=minmod(dur,b*dul); double durLim=minmod(dul,b*dur); double delu=(1.0+k)/2.0*dulLim+(1.0-k)/2.0*durLim; //boundary extrapolation// uiL[i]=u[i]-0.5*delu; uiR[i]=u[i]+0.5*delu; } //BC// //i=0// uiL[0]=uiL[1]; uiR[0]=uiR[1]; //i=Node-1// uiL[Node-1]=uiL[Node-2]; uiR[Node-1]=uiR[Node-2]; } inline void output(double u[], int Node, double dx) { stringstream ss; string name; ofstream fo; static int count=0; ss<<count; name=ss.str(); name="u_" + name + ".txt"; fo.open(name.c_str ()); for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<u[i]<<endl; } fo<<endl; fo<<endl; count+=1; } int main() { const int Cell=400; const int Node=Cell+1; int div=100;//for output double LL=1.0; double dx=LL/(Node-1); double dt=0.001; double TotalTime=10; double V=0.1; double u[Node]; double ub[Node]; double un[Node]; double FL; double FR; double uiL[Node]; double uiR[Node]; //IC// for(int i=0;i<Node;i++) { u[i]=0.0; if((i>=80)&&(i<=120)) u[i]=1.0; } //IC output// output(u,Node,dx); //time step// for(int k=1;double(k)*dt<=TotalTime;k++) { //reconstruction// reconst(Node,dt,dx,u,uiL,uiR); //inner cells// for(int i=1;i<Node-1;i++) { //FR=Fluxadv(V,u[i],u[i+1]); //FL=Fluxadv(V,u[i-1],u[i]); FR=Fluxadv(V,uiR[i],uiL[i+1]); FL=Fluxadv(V,uiR[i-1],uiL[i]); un[i]=u[i]-dt/dx*(FR-FL); } //BC// //left boundary// un[0]=un[1]; //right boundary// un[Node-1]=un[Node-2]; //update// for(int i=0;i<Node;i++) { u[i]=un[i]; } //output// if(k%div==0) { output(u,Node,dx); cout<<"t="<<double(k)*dt<<endl; } } return 0; }
こちらではToro流のslope limiterを入れてあります。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <string> #include <sstream> using namespace std; //flux for advection equation// inline double Fluxadv(double V, double uL, double uR) { return max(V,0.0)*uL-max(-V,0.0)*uR; } //flux limiter superbee// inline double superbee(double d1, double d2, double delta, double omega) { double zeta; if(d1*d2<=0) { zeta=0.0; } else { if(d1/d2<=0.5) { zeta=2.0*d1/d2; } else if(d1/d2<=1.0) { zeta=1.0; } else { double zetaR=2.0/(1.0-omega+(1.0+omega)*d1/d2); zeta=min(zetaR,d1/d2); zeta=min(zeta,2.0); } } return zeta*delta; } //reconstruction// inline void reconst(int Node, double dt, double dx, double u[], double uiL[], double uiR[]) { //inner cells// for(int i=1;i<Node-1;i++) { //jump// double dul=u[i]-u[i-1]; double dur=u[i+1]-u[i]; //slope// double omega=0.0; double delu=(1.0+omega)/2.0*dul+(1.0-omega)/2.0*dur; //slope limiter// delu=superbee(dul,dur,delu,omega); //boundary extrapolation// uiL[i]=u[i]-0.5*delu; uiR[i]=u[i]+0.5*delu; } //BC// //i=0// uiL[0]=uiL[1]; uiR[0]=uiR[1]; //i=Node-1// uiL[Node-1]=uiL[Node-2]; uiR[Node-1]=uiR[Node-2]; } inline void output(double u[], int Node, double dx) { stringstream ss; string name; ofstream fo; static int count=0; ss<<count; name=ss.str(); name="u_" + name + ".txt"; fo.open(name.c_str ()); for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<u[i]<<endl; } fo<<endl; fo<<endl; count+=1; } int main() { const int Cell=800; const int Node=Cell+1; int div=100;//for output double LL=1.0; double dx=LL/(Node-1); double dt=0.001; double TotalTime=10; double V=0.1; double u[Node]; double ub[Node]; double un[Node]; double FL; double FR; double uiL[Node]; double uiR[Node]; //IC// for(int i=0;i<Node;i++) { u[i]=0.0; if((i>=80)&&(i<=160)) u[i]=1.0; } //IC output// output(u,Node,dx); //time step// for(int k=1;double(k)*dt<=TotalTime;k++) { //reconstruction// reconst(Node,dt,dx,u,uiL,uiR); //inner cells// for(int i=1;i<Node-1;i++) { //FR=Fluxadv(V,u[i],u[i+1]); //FL=Fluxadv(V,u[i-1],u[i]); FR=Fluxadv(V,uiR[i],uiL[i+1]); FL=Fluxadv(V,uiR[i-1],uiL[i]); un[i]=u[i]-dt/dx*(FR-FL); } //BC// //left boundary// un[0]=un[1]; //right boundary// un[Node-1]=un[Node-2]; //update// for(int i=0;i<Node;i++) { u[i]=un[i]; } //output// if(k%div==0) { output(u,Node,dx); cout<<"t="<<double(k)*dt<<endl; } } return 0; }
【浅水流方程式】MUSCL-Hancock法とHLL法でダム崩壊問題を解いてみました 動画とC++コード付き
【浅水流方程式】サイトマップ(ここから関連記事が探せます)
http://mathlang.hatenablog.com/entry/2018/12/18/213650
今回は、HLL法にMUSCL-Hancock法を組み合わせて高次精度化したスキームでダム崩壊問題を解きました。
普通のHLL法が時間・空間方向に対して1次精度なのにたいして、MUSCL-Hancock法を組み合わせたHLL法は、時間・空間方向に対して2次精度になります!普通のGodunov法(有限体積法)では各セルでの未知数は定数としますが、MUSCL-Hancock法では、各セルにおいて未知数は線型に変化するとします。こうすることにより、より高精度な数値流束を求めることが可能となります。また、衝撃波近傍での数値振動を回避するために、slope limiter(勾配制限関数)の一種であるsuperbee limiterを採用しています。
解いている方程式は、摩擦、勾配を無視した1次元浅水流方程式(Shallow Water Equations)です。
参考文献としてはToroの"Shock-Capturing Methods for Free-Surface Shallow Flows"のpp.205-210です。
Shock-Capturing Methods for Free-Surface Shallow Flows
- 作者: Eleuterio F. Toro
- 出版社/メーカー: Wiley
- 発売日: 2001/03/23
- メディア: ハードカバー
- この商品を含むブログを見る
こちらが計算結果です。青い線が水深です。
より少ない格子数で、普通のHLL法よりも衝撃波がシャープに捉えられています。
本当にそうだろうか?と疑問を持つ方のために比較動画を用意しました。セル数は100です。青がMUSCL-Hancock法を入れたHLL法で、赤がただのHLL法です。
やはりMUSCLのほうが衝撃波の先端がシャープです。しかも、膨張波も角まで捉えられているし、途中の形状もより正確です。
では計算コードです。Toro流のslope limiterが入っています。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <sstream> using namespace std; const double g=9.8; //celerity// inline double a(double h) { return sqrt(g*h); } //function for TRRS// inline double gk(double h0, double hlr) { return sqrt(g/2.0*(h0+hlr)/(h0*hlr)); } //Riemann Solver Based on Exact Depth Positivity// inline void EDPRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double aL=a(hL); double aR=a(hR); double uL=qL/hL; double uR=qR/hR; hst=(hL+hR)/2.0-(uR-uL)*(hL+hR)/4.0/(aL+aR); ust=(uL+uR)/2.0-(hR-hL)*(aL+aR)/(hL+hR); } //Two-Rarefaction Riemann Solver// inline void TRRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double aL=a(hL); double aR=a(hR); double uL=qL/hL; double uR=qR/hR; hst=1.0/g*((aL+aR)/2.0+(uL-uR)/4.0)*((aL+aR)/2.0+(uL-uR)/4.0); ust=(uL+uR)/2.0+aL-aR; } //Two-Shock Riemann Solver// inline void TSRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double uL=qL/hL; double uR=qR/hR; TRRS(hst,ust,hL,hR,qL,qR); //EDPRS(hst,ust,hL,hR,qL,qR); double hcon=hst; hst=(gk(hst,hL)*hL+gk(hst,hR)*hR+uL-uR)/(gk(hst,hL)+gk(hst,hR)); ust=(uL+uR)/2.0+((hcon-hR)*gk(hcon,hR)-(hcon-hL)*gk(hcon,hL))/2.0; } //function for HLL// inline double qk(double hst, double hlr) { double qk=1.0; if(hst>hlr) qk=sqrt((hst+hlr)*hst/(hlr*hlr)/2.0); return qk; } //left wave estimation// inline double SLeft(double hst, double hL, double qL) { double uL=qL/hL; return uL-a(hL)*qk(hst,hL); } //right wave estimation// inline double SRight(double hst, double hR, double qR) { double uR=qR/hR; return uR+a(hR)*qk(hst,hR); } //flux for continuity equation// inline double Fluxcon(double h, double q) { return q; } //flux for momentum equation// inline double Fluxmom(double h, double q) { return q*q/h+g*h*h/2.0; } //hll flux// inline void HLLFlux(double FLR[], double hL, double hR, double qL, double qR) { double hst,ust; //estimastion of hst and ust// TRRS(hst,ust,hL,hR,qL,qR); //TSRS(hst,ust,hL,hR,qL,qR); //EDPRS(hst,ust,hL,hR,qL,qR); //wave speed estimation// double SL=SLeft(hst,hL,qL); double SR=SRight(hst,hR,qR); if(SL>=0) { FLR[0]=Fluxcon(hL,qL); FLR[1]=Fluxmom(hL,qL); } else if(SR<=0) { FLR[0]=Fluxcon(hR,qR); FLR[1]=Fluxmom(hR,qR); } else { FLR[0]=(SR*Fluxcon(hL,qL)-SL*Fluxcon(hR,qR)+SR*SL*(hR-hL))/(SR-SL); FLR[1]=(SR*Fluxmom(hL,qL)-SL*Fluxmom(hR,qR)+SR*SL*(qR-qL))/(SR-SL); } } inline void output(double h[], double z[], int Node, double dx) { stringstream ss; string name; ofstream fo; static int count=0; ss<<count; name=ss.str(); name="h_" + name + ".txt"; fo.open(name.c_str ()); for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<h[i]+z[i]<<endl; } fo<<endl; fo<<endl; for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<z[i]<<endl; } count+=1; } //flux limiter superbee// inline double superbee(double d1, double d2, double delta, double omega) { double zeta; if(d1*d2<=0) { zeta=0.0; } else { if(d1/d2<=0.5) { zeta=2.0*d1/d2; } else if(d1/d2<=1.0) { zeta=1.0; } else { double zetaR=2.0/(1.0-omega+(1.0+omega)*d1/d2); zeta=min(zetaR,d1/d2); zeta=min(zeta,2.0); } } return zeta*delta; } //reconstruction// inline void reconst(int Node, double dt, double dx, double h[], double q[], double hiLn[], double hiRn[], double qiLn[], double qiRn[]) { double hiL[Node]; double hiR[Node]; double qiL[Node]; double qiR[Node]; //inner cells// for(int i=1;i<Node-1;i++) { //jump// double dhl=h[i]-h[i-1]; double dhr=h[i+1]-h[i]; double dql=q[i]-q[i-1]; double dqr=q[i+1]-q[i]; //slope// double omega=0.0; double delh=(1.0+omega)/2.0*dhl+(1.0-omega)/2.0*dhr; double delq=(1.0+omega)/2.0*dql+(1.0-omega)/2.0*dqr; //slope limiter// delh=superbee(dhl,dhr,delh,omega); delq=superbee(dql,dqr,delq,omega); //boundary extrapolation// hiL[i]=h[i]-0.5*delh; hiR[i]=h[i]+0.5*delh; qiL[i]=q[i]-0.5*delq; qiR[i]=q[i]+0.5*delq; //time evolution of boundary extrapolated values// hiLn[i]=hiL[i]-0.5*dt/dx*(Fluxcon(hiR[i],qiR[i])-Fluxcon(hiL[i],qiL[i])); hiRn[i]=hiR[i]-0.5*dt/dx*(Fluxcon(hiR[i],qiR[i])-Fluxcon(hiL[i],qiL[i])); qiLn[i]=qiL[i]-0.5*dt/dx*(Fluxmom(hiR[i],qiR[i])-Fluxmom(hiL[i],qiL[i])); qiRn[i]=qiR[i]-0.5*dt/dx*(Fluxmom(hiR[i],qiR[i])-Fluxmom(hiL[i],qiL[i])); } //BC// //i=0// hiLn[0]=hiLn[1]; hiRn[0]=hiRn[1]; qiLn[0]=qiLn[1]; qiRn[0]=qiRn[1]; //i=Node-1// hiLn[Node-1]=hiLn[Node-2]; hiRn[Node-1]=hiRn[Node-2]; qiLn[Node-1]=qiLn[Node-2]; qiRn[Node-1]=qiRn[Node-2]; } int main() { const int n=200;//number of cells const int Node=n+1; double L=1.0; double dx=L/double(n); double dt=0.0001; double TotalTime=0.15; int div=20;//for output double q[Node]; double h[Node]; double qn[Node]; double hn[Node]; double FL[2],FR[2]; double InithL=1.0; double InithR=0.1; double z[Node]; double I[Node]; double mn=0.01; //for reconstruction// double hiLn[Node]; double hiRn[Node]; double qiLn[Node]; double qiRn[Node]; //IC// for(int i=0;i<Node;i++) { q[i]=0.0; h[i]=InithL; } for(int i=Node/2;i<Node;i++) { q[i]=0.0; h[i]=InithR; } //bottom// for(int i=0;i<Node;i++) { z[i]=0.0; double a=0.1-3.0*(float(i)*dx-0.5)*(float(i)*dx-0.5); if(a>0.0) z[i]=a; z[i]=0.0; //cout<<z[i]<<endl; } //slope// for(int i=1;i<Node-1;i++) { I[i]=-(z[i]-z[i-1])/dx; } //at boundary// I[0]=-(z[1]-z[0])/dx; I[Node-1]=-(z[Node-1]-z[Node-2])/dx; //IC output// output(h,z,Node,dx); //time step// for(int k=1;double(k)*dt<=TotalTime;k++) { //reconstruction// reconst(Node,dt,dx,h,q,hiLn,hiRn,qiLn,qiRn); //HLL// for(int i=1;i<Node-1;i++) { //left flux, between i-1 and i//i-1R,iL HLLFlux(FL,hiRn[i-1],hiLn[i],qiRn[i-1],qiLn[i]); //right flux, between i and i+1//iR,i+1L HLLFlux(FR,hiRn[i],hiLn[i+1],qiRn[i],qiLn[i+1]); hn[i]=h[i]-dt/dx*(FR[0]-FL[0]); qn[i]=q[i]-dt/dx*(FR[1]-FL[1]);//+dt*g*h[i]*(I[i]-mn*mn*q[i]*abs(q[i])/pow(h[i],10.0/3.0)); } //BC// //left boundary// hn[0]=hn[1]; qn[0]=qn[1]; //right boundary// hn[Node-1]=hn[Node-2]; qn[Node-1]=qn[Node-2]; //update// for(int i=0;i<Node;i++) { q[i]=qn[i]; h[i]=hn[i]; } //output// if(k%div==0) { output(h,z,Node,dx); cout<<"t="<<double(k)*dt<<endl; } } return 0; }
こちらは一般的なMUSCLに書き直したほうです。両者はほぼ一致します。
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <sstream> using namespace std; const double g=9.8; //celerity// inline double a(double h) { return sqrt(g*h); } //function for TRRS// inline double gk(double h0, double hlr) { return sqrt(g/2.0*(h0+hlr)/(h0*hlr)); } //Riemann Solver Based on Exact Depth Positivity// inline void EDPRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double aL=a(hL); double aR=a(hR); double uL=qL/hL; double uR=qR/hR; hst=(hL+hR)/2.0-(uR-uL)*(hL+hR)/4.0/(aL+aR); ust=(uL+uR)/2.0-(hR-hL)*(aL+aR)/(hL+hR); } //Two-Rarefaction Riemann Solver// inline void TRRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double aL=a(hL); double aR=a(hR); double uL=qL/hL; double uR=qR/hR; hst=1.0/g*((aL+aR)/2.0+(uL-uR)/4.0)*((aL+aR)/2.0+(uL-uR)/4.0); ust=(uL+uR)/2.0+aL-aR; } //Two-Shock Riemann Solver// inline void TSRS(double &hst, double &ust, double hL, double hR, double qL, double qR) { double uL=qL/hL; double uR=qR/hR; TRRS(hst,ust,hL,hR,qL,qR); //EDPRS(hst,ust,hL,hR,qL,qR); double hcon=hst; hst=(gk(hst,hL)*hL+gk(hst,hR)*hR+uL-uR)/(gk(hst,hL)+gk(hst,hR)); ust=(uL+uR)/2.0+((hcon-hR)*gk(hcon,hR)-(hcon-hL)*gk(hcon,hL))/2.0; } //function for HLL// inline double qk(double hst, double hlr) { double qk=1.0; if(hst>hlr) qk=sqrt((hst+hlr)*hst/(hlr*hlr)/2.0); return qk; } //left wave estimation// inline double SLeft(double hst, double hL, double qL) { double uL=qL/hL; return uL-a(hL)*qk(hst,hL); } //right wave estimation// inline double SRight(double hst, double hR, double qR) { double uR=qR/hR; return uR+a(hR)*qk(hst,hR); } //flux for continuity equation// inline double Fluxcon(double h, double q) { return q; } //flux for momentum equation// inline double Fluxmom(double h, double q) { return q*q/h+g*h*h/2.0; } //hll flux// inline void HLLFlux(double FLR[], double hL, double hR, double qL, double qR) { double hst,ust; //estimastion of hst and ust// TRRS(hst,ust,hL,hR,qL,qR); //TSRS(hst,ust,hL,hR,qL,qR); //EDPRS(hst,ust,hL,hR,qL,qR); //wave speed estimation// double SL=SLeft(hst,hL,qL); double SR=SRight(hst,hR,qR); if(SL>=0) { FLR[0]=Fluxcon(hL,qL); FLR[1]=Fluxmom(hL,qL); } else if(SR<=0) { FLR[0]=Fluxcon(hR,qR); FLR[1]=Fluxmom(hR,qR); } else { FLR[0]=(SR*Fluxcon(hL,qL)-SL*Fluxcon(hR,qR)+SR*SL*(hR-hL))/(SR-SL); FLR[1]=(SR*Fluxmom(hL,qL)-SL*Fluxmom(hR,qR)+SR*SL*(qR-qL))/(SR-SL); } } inline void output(double h[], double z[], int Node, double dx) { stringstream ss; string name; ofstream fo; static int count=0; ss<<count; name=ss.str(); name="h_" + name + ".txt"; fo.open(name.c_str ()); for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<h[i]+z[i]<<endl; } fo<<endl; fo<<endl; for(int i=0;i<Node;i++) { fo<<dx*float(i)<<" "<<z[i]<<endl; } count+=1; } //flux limiter minmod// inline double minmod(double d1, double d2) { if(d1*d2<=0) { return 0.0; } else if(abs(d1)<abs(d2)) { return d1; } else { return d2; } } //reconstruction// inline void reconst(int Node, double dt, double dx, double h[], double q[], double hiLn[], double hiRn[], double qiLn[], double qiRn[]) { double hiL[Node]; double hiR[Node]; double qiL[Node]; double qiR[Node]; //inner cells// for(int i=1;i<Node-1;i++) { //jump// double dhl=h[i]-h[i-1]; double dhr=h[i+1]-h[i]; double dql=q[i]-q[i-1]; double dqr=q[i+1]-q[i]; //slope limiter// double k=0.0; double b=(3.0-k)/(1.0-k); double dhlLim=minmod(dhr,b*dhl); double dhrLim=minmod(dhl,b*dhr); double dqlLim=minmod(dqr,b*dql); double dqrLim=minmod(dql,b*dqr); double omega=0.0; double delh=(1.0+omega)/2.0*dhlLim+(1.0-omega)/2.0*dhrLim; double delq=(1.0+omega)/2.0*dqlLim+(1.0-omega)/2.0*dqrLim; //boundary extrapolation// hiL[i]=h[i]-0.5*delh; hiR[i]=h[i]+0.5*delh; qiL[i]=q[i]-0.5*delq; qiR[i]=q[i]+0.5*delq; //time evolution of boundary extrapolated values// hiLn[i]=hiL[i]-0.5*dt/dx*(Fluxcon(hiR[i],qiR[i])-Fluxcon(hiL[i],qiL[i])); hiRn[i]=hiR[i]-0.5*dt/dx*(Fluxcon(hiR[i],qiR[i])-Fluxcon(hiL[i],qiL[i])); qiLn[i]=qiL[i]-0.5*dt/dx*(Fluxmom(hiR[i],qiR[i])-Fluxmom(hiL[i],qiL[i])); qiRn[i]=qiR[i]-0.5*dt/dx*(Fluxmom(hiR[i],qiR[i])-Fluxmom(hiL[i],qiL[i])); } //BC// //i=0// hiLn[0]=hiLn[1]; hiRn[0]=hiRn[1]; qiLn[0]=qiLn[1]; qiRn[0]=qiRn[1]; //i=Node-1// hiLn[Node-1]=hiLn[Node-2]; hiRn[Node-1]=hiRn[Node-2]; qiLn[Node-1]=qiLn[Node-2]; qiRn[Node-1]=qiRn[Node-2]; } int main() { const int n=200;//number of cells const int Node=n+1; double L=1.0; double dx=L/double(n); double dt=0.0001; double TotalTime=0.15; int div=20;//for output double q[Node]; double h[Node]; double qn[Node]; double hn[Node]; double FL[2],FR[2]; double InithL=1.0; double InithR=0.1; double z[Node]; double I[Node]; double mn=0.01; //for reconstruction// double hiLn[Node]; double hiRn[Node]; double qiLn[Node]; double qiRn[Node]; //IC// for(int i=0;i<Node;i++) { q[i]=0.0; h[i]=InithL; } for(int i=Node/2;i<Node;i++) { q[i]=0.0; h[i]=InithR; } //bottom// for(int i=0;i<Node;i++) { z[i]=0.0; double a=0.1-3.0*(float(i)*dx-0.5)*(float(i)*dx-0.5); if(a>0.0) z[i]=a; z[i]=0.0; //cout<<z[i]<<endl; } //slope// for(int i=1;i<Node-1;i++) { I[i]=-(z[i]-z[i-1])/dx; } //at boundary// I[0]=-(z[1]-z[0])/dx; I[Node-1]=-(z[Node-1]-z[Node-2])/dx; //IC output// output(h,z,Node,dx); //time step// for(int k=1;double(k)*dt<=TotalTime;k++) { //reconstruction// reconst(Node,dt,dx,h,q,hiLn,hiRn,qiLn,qiRn); //HLL// for(int i=1;i<Node-1;i++) { //left flux, between i-1 and i//i-1R,iL HLLFlux(FL,hiRn[i-1],hiLn[i],qiRn[i-1],qiLn[i]); //right flux, between i and i+1//iR,i+1L HLLFlux(FR,hiRn[i],hiLn[i+1],qiRn[i],qiLn[i+1]); hn[i]=h[i]-dt/dx*(FR[0]-FL[0]); qn[i]=q[i]-dt/dx*(FR[1]-FL[1]);//+dt*g*h[i]*(I[i]-mn*mn*q[i]*abs(q[i])/pow(h[i],10.0/3.0)); } //BC// //left boundary// hn[0]=hn[1]; qn[0]=qn[1]; //right boundary// hn[Node-1]=hn[Node-2]; qn[Node-1]=qn[Node-2]; //update// for(int i=0;i<Node;i++) { q[i]=qn[i]; h[i]=hn[i]; } //output// if(k%div==0) { output(h,z,Node,dx); cout<<"t="<<double(k)*dt<<endl; } } return 0; }