【浅水流方程式】サイトマップ（ここから関連記事が探せます）
浅水流方程式にはいくつか解析解が出る問題があり、テスト問題としてスキームの精度を調べるのに活用されています。その中で一番有名なのがダム崩壊問題（dam break problem）です。真ん中にダムがあり、ダムの左右には水深の違う水面があります。水路は勾配が無いものとします。時刻0の瞬間にダムを取り去った時に水がどのように動くか、という問題です。
今回は、この問題の解析解を計算するC++コードを載せます。なぜ解析解なのに計算コードを載せるのか？実はこの場合の解析解は陽的には書けず、ある非線型の方程式の解になっています。コードではこれをニュートン法で解いています。その後水深と流速をテキストファイルに出力しています。
コードはToroの"Shock-Capturing Methods for Free-Surface Shallow Flows"の5章と6章をもとに書いています。wet bed、dry bedともに計算できます。初期水深が正で、後にdry bedがあらわれる場合も大丈夫です。
Shock-Capturing Methods for Free-Surface Shallow Flows
- 作者: Eleuterio F. Toro
- 出版社/メーカー: Wiley
- 発売日: 2001/03/23
- メディア: ハードカバー
- この商品を含むブログを見る
#include <iostream> #include <cmath> #include <fstream> #include <iomanip> #include <string> #include <sstream> using namespace std; const double g=9.8; inline double a(double h) { return sqrt(g*h); } enum output { excel=0, gnuplot=1, }; enum filename { time=0, number=1, }; enum category { wet_Lshock_Rshock=0, wet_Lrarefaction_Rshock=1, wet_Lshock_Rrarefaction=2, wet_Lrarefaction_Rrarefaction=3, dry_R=4, dry_L=5, dry_Generation=6, }; inline void output_excel(double h, double u, ofstream &fh, ofstream &fu) { fh<<h<<" "; fu<<u<<" "; } inline void output_gnuplot(double t, double x, double h, double u, ofstream &fh, ofstream &fu) { fh<<t<<" "<<x<<" "<<h<<endl; fu<<t<<" "<<x<<" "<<u<<endl; } inline double F(double h, double hstar) { if(hstar<=h)//rarefaction { return 2.0*(sqrt(g*hstar)-sqrt(g*h)); } else//shock { return (hstar-h)*sqrt(g/2.0*(hstar+h)/hstar/h); } } inline double dF(double h, double hstar) { if(hstar<=h)//rarefaction { return g/a(h); } else//shock { double gk=sqrt(g/2.0*(hstar+h)/hstar/h); return gk-g*(hstar-h)/(4.0*hstar*hstar*gk); } } inline void newton_method(double hl, double hr, double &hstar, double du) { for(int i=0;i<10000;i++) { double nhstar=hstar-(F(hl,hstar)+F(hr,hstar)+du)/(dF(hl,hstar)+dF(hr,hstar)); hstar=nhstar; } } inline void wave(double ul, double ur, double ustar, double hl, double hr, double hstar, ofstream &fk, double &shl, double &stl, double &sl, double &shr, double &str, double &sr, double u_crit, double du, category &wave_category) { int left,right; if(hr==0)//dry bed { shl=ul-a(hl); stl=ul+2.0*a(hl); fk<<"left rarefaction dry bed with head speed "<<shl<<" and tail speed "<<stl<<endl; wave_category=dry_R; } else if(hl==0) { shr=ur+a(hr); str=ur-2.0*a(hr); fk<<"right rarefaction dry bed with head speed "<<shr<<" and tail speed "<<str<<endl; wave_category=dry_L; } else if(u_crit<=du) { shl=ul-a(hl); stl=ul+2.0*a(hl); shr=ur+a(hr); str=ur-2.0*a(hr); fk<<"generation of dry bed"<<endl; fk<<"left rarefaction dry bed with head speed "<<shl<<" and tail speed "<<stl<<endl; fk<<"right rarefaction dry bed with head speed "<<shr<<" and tail speed "<<str<<endl; wave_category=dry_Generation; } else//wet bed { //left wave// if(hstar<=hl)//rarefaction { shl=ul-a(hl); stl=ustar-a(hstar); fk<<"left rarefaction with head speed"<<" "<<shl<<" "<<"and tail speed"<<" "<<stl<<endl; left=1; } else//shock { double ql=sqrt((hstar+hl)*hstar/hl/hl/2.0); sl=ul-a(hl)*ql; fk<<"left shock with shock speed"<<" "<<sl<<endl; left=0; } //right wave// if(hstar<=hr)//rarefaction { shr=ur+a(hr); str=ustar+a(hstar); fk<<"right rarefaction with head speed"<<" "<<shr<<" "<<"and tail speed"<<" "<<str<<endl; right=1; } else//shock { double qr=sqrt((hstar+hr)*hstar/hr/hr/2.0); sr=ur+a(hr)*qr; fk<<"right shock with shock speed"<<" "<<sr<<endl; right=0; } if((left==0)&&(right==0)) wave_category=wet_Lshock_Rshock; if((left==1)&&(right==0)) wave_category=wet_Lrarefaction_Rshock; if((left==0)&&(right==1)) wave_category=wet_Lshock_Rrarefaction; if((left==1)&&(right==1)) wave_category=wet_Lrarefaction_Rrarefaction; } } inline void exact_solution(double hl, double hr, double ul, double ur, double hstar, double ustar, double sl, double sr, double shl, double stl, double shr, double str, output output_case, category wave_category, filename filename_case) { double dt=0.01; double dx=0.01; double xl1,xl2,xr1,xr2; //initial condition// switch(output_case) { case excel: { ofstream fh,fu; fh.open("excel_h.txt"); fu.open("excel_u.txt"); for(double x=-5.0;x<=5.0;x=x+dx) { fh<<x<<" "; fu<<x<<" "; } fh<<endl; fu<<endl; for(double x=-5.0;x<=5.0;x=x+dx) { if(x<=0.0) { output_excel(hl,ul,fh,fu); } if(x>0.0) { output_excel(hr,ur,fh,fu); } } fh<<endl; fu<<endl; for(double t=dt;t<=0.8;t=t+dt)//time { for(double x=-5.0;x<=5.0;x=x+dx)//space { switch(wave_category) { case(wet_Lshock_Rshock): { xl1=sl*t; xr1=sr*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xr1)) output_excel(hstar,ustar,fh,fu); if(x>xr1) output_excel(hr,ur,fh,fu); break; } case(wet_Lrarefaction_Rshock): { xl1=shl*t; xl2=stl*t; xr1=sr*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if((xl2<x)&&(x<=xr1)) output_excel(hstar,ustar,fh,fu); if(xr1<x) output_excel(hr,ur,fh,fu); break; } case(wet_Lshock_Rrarefaction): { xl1=sl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xr1)) output_excel(hstar,ustar,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if(x>xr2) output_excel(hr,ur,fh,fu); break; } case(wet_Lrarefaction_Rrarefaction): { xl1=shl*t; xl2=stl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if((xl2<x)&&(x<=xr1)) output_excel(hstar,ustar,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if(x>xr2) output_excel(hr,ur,fh,fu); break; } case(dry_R): { xl1=shl*t; xl2=stl*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if(xl2<x) output_excel(hr,ur,fh,fu); break; } case(dry_L): { xr1=str*t; xr2=shr*t; if(x<=xr1) output_excel(hl,ul,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if(xr2<x) output_excel(hr,ur,fh,fu); break; } case(dry_Generation): { xl1=shl*t; xl2=stl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_excel(hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if((xl2<x)&&(x<=xr1)) output_excel(0,0,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_excel(hh,uu,fh,fu); } if(xr2<x) output_excel(hr,ur,fh,fu); break; } } } fh<<endl; fu<<endl; } fh.close(); fu.close(); break; } case gnuplot: { ofstream fh,fu; fh.open("gnuplot_h_0.txt"); fu.open("gnuplot_u_0.txt"); for(double x=-5.0;x<=5.0;x=x+dx) { if(x<=0.0) { fh<<0<<" "<<x<<" "<<hl<<endl; fu<<0<<" "<<x<<" "<<ul<<endl; } if(x>0.0) { fh<<0<<" "<<x<<" "<<hr<<endl; fu<<0<<" "<<x<<" "<<ur<<endl; } } fh.close(); fu.close(); int count=0; for(double t=dt;t<=0.8;t=t+dt)//time { ofstream fh,fu; switch(filename_case) { case time: { stringstream ss; string nameh,nameu; ss<<t; nameh=ss.str(); nameh="gnuplot_h_" + nameh + ".txt"; fh.open(nameh.c_str ()); nameu=ss.str(); nameu="gnuplot_u_" + nameu + ".txt"; fu.open(nameu.c_str ()); break; } case number: { count+=1; stringstream ss; string nameh,nameu; ss<<count; nameh=ss.str(); nameh="gnuplot_h_" + nameh + ".txt"; fh.open(nameh.c_str ()); nameu=ss.str(); nameu="gnuplot_u_" + nameu + ".txt"; fu.open(nameu.c_str ()); break; } } for(double x=-5.0;x<=5.0;x=x+dx)//space { switch(wave_category) { case(wet_Lshock_Rshock): { xl1=sl*t; xr1=sr*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xr1)) output_gnuplot(t,x,hstar,ustar,fh,fu); if(x>xr1) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(wet_Lrarefaction_Rshock): { xl1=shl*t; xl2=stl*t; xr1=sr*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if((xl2<x)&&(x<=xr1)) output_gnuplot(t,x,hstar,ustar,fh,fu); if(xr1<x) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(wet_Lshock_Rrarefaction): { xl1=sl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xr1)) output_gnuplot(t,x,hstar,ustar,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if(x>xr2) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(wet_Lrarefaction_Rrarefaction): { xl1=shl*t; xl2=stl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if((xl2<x)&&(x<=xr1)) output_gnuplot(t,x,hstar,ustar,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if(x>xr2) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(dry_R): { xl1=shl*t; xl2=stl*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if(xl2<x) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(dry_L): { xr1=str*t; xr2=shr*t; if(x<=xr1) output_gnuplot(t,x,hl,ul,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if(xr2<x) output_gnuplot(t,x,hr,ur,fh,fu); break; } case(dry_Generation): { xl1=shl*t; xl2=stl*t; xr1=str*t; xr2=shr*t; if(x<=xl1) output_gnuplot(t,x,hl,ul,fh,fu); if((xl1<x)&&(x<=xl2)) { double aa=(ul+2.0*a(hl)-x/t)/3.0; double hh=aa*aa/g; double uu=(ul+2.0*a(hl)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu);; } if((xl2<x)&&(x<=xr1)) output_gnuplot(t,x,0,0,fh,fu); if((xr1<x)&&(x<=xr2)) { double aa=(-ur+2.0*a(hr)+x/t)/3.0; double hh=aa*aa/g; double uu=(ur-2.0*a(hr)+2.0*x/t)/3.0; output_gnuplot(t,x,hh,uu,fh,fu); } if(xr2<x) output_gnuplot(t,x,hr,ur,fh,fu); break; } } } fh.close(); fu.close(); } break; } } } int main() { double hl=1.0; double hr=0.0; double ul=0.0; double ur=0.0; double du=ur-ul; double shl,stl,sl,shr,str,sr; double ustar,hstar; double u_crit=2.0*(a(hl)+a(hr)); //output with excel, gnuplot// output output_case=gnuplot; //7 patterns// category wave_category; //filename for gnuplot_each, time or number// filename filename_case=number; ofstream fk; fk.open("answer.txt"); if((hl>0)&&(hr>0))//wet bed { if(u_crit>du)//complete wet bed { //initial depth by two-rarefaction approximation// hstar=1.0/g*((a(hl)+a(hr))/2.0-du/4.0)*((a(hl)+a(hr))/2.0-du/4.0);//initial value //newton method// newton_method(hl,hr,hstar,du); fk<<"hstar"<<" "<<hstar<<endl; ustar=(ul+ur)/2.0+(F(hr,hstar)-F(hl,hstar))/2.0; fk<<"ustar"<<" "<<ustar<<endl; } } wave(ul,ur,ustar,hl,hr,hstar,fk,shl,stl,sl,shr,str,sr,u_crit,du,wave_category); exact_solution(hl,hr,ul,ur,hstar,ustar,sl,sr,shl,stl,shr,str,output_case,wave_category,filename_case); return 0; }
こちらもどうぞ
