数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【浅水流方程式】ダム崩壊問題(dam break problem)の解析解を計算するC++コード

【浅水流方程式】サイトマップ(ここから関連記事が探せます)

浅水流方程式にはいくつか解析解が出る問題があり、テスト問題としてスキームの精度を調べるのに活用されています。その中で一番有名なのがダム崩壊問題(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

Shock-Capturing Methods for Free-Surface Shallow Flows

#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;
}

こちらもどうぞ