飛行ロボット大会も大きな大きな悔しさを残しながら終わったので(24日の日記参照)、shocktubeのプログラムを書くことにした。昨日の飲み会前に終わらせたいなと思っていたのに歯医者とかあって終わらず、結局今日になって書くのは終わった。結果を表示させてみたらやっぱりどう考えてもおかしい。これからバグ取り無間地獄。
1次元shocktubeをYee 2nd-order upwind TVD schemeで解いた。コードは以下の通り。計算は一瞬で終わった。ubuntuでtimeコマンドで時間測ってみたらだいたい0.3秒くらい。一回一回速度違うんだなと思った。コンパイラはg95。なんでgfortranじゃないかっていうとgfortranはVistaにインストールできなかったから。
- base.f90
program shocktube implicit none integer,parameter :: kiza_x=200 !kiza_xは空間刻み数で半分の位置が初期膜境界：偶数である必要がある real,parameter :: dt=0.000001 !dtは時間ステップ real,parameter :: cv=1007.0,r=287.0,gamma=1.4 !cvは定圧比熱、rは空気の気体定数,gammaは比熱比 real :: dx !dxは空間刻み距離 real :: q(3,kiza_x),e(3,kiza_x),u(kiza_x),h(kiza_x),ee(3,kiza_x-1),rho(kiza_x),p(kiza_x),t(kiza_x),ene(kiza_x),c(3,kiza_x) !real :: t(kiza_x) real :: right(3,3,kiza_x-1),left(3,3,kiza_x-1),b(2),alpha(3,kiza_x-1),sigma(3,kiza_x),g(3,kiza_x-1),gam(3,kiza_x),phi(3,kiza_x) real :: psi,sig,minmod integer :: i,j,k !q(物理量、場所)は保存量Q(密度、運動量、内部エネルギ) !e(物理量、場所)は流速ベクトルE（運動量、p+ρu^2,(e+p)u） !u(場所)は場所によるx方向速度 !h(場所)は全エンタルピ !ee(物理量、場所)は数値流束ベクトル。一つ右の値を表す !rho(場所)は密度ρ !p(場所)は静圧p !ene(場所)は全エネルギe !right()は右固有ベクトル !left()は左固有ベクトル !b()は固有ベクトル計算の定数 !alpha()は各特性線をまたいだ物理量の跳び !sigma()は計算に使うやつ !gはminmodが中に入っててαを比較 !gamはγの演算 !phiは流束ベクトル計算に使うベクトル !c()は固有値ベクトル !psiはΨの関数 !sigはσの関数 !minmodはminmod limiter real :: rho_av(kiza_x-1),u_av(kiza_x-1),h_av(kiza_x-1),c_av(kiza_x-1) !rho_av()はRoe平均密度 !u_av()はRoe平均速度 !h_av()はRoe平均全エンタルピ !c_av()はRoe平均特性速度 dx=0.2/kiza_x !dxは一個あたりの刻み幅 !print*, dx !計算開始 !do i=1,kiza_x p(1:kiza_x/2)=5.0*10**5 p(kiza_x/2+1:kiza_x)=1.0*10**5 !静圧の初期化 u(1:kiza_x)=0.0 !x方向速度の初期化 rho(1:kiza_x/2)=6.0 rho(kiza_x/2+1:kiza_x)=1.0 !密度の初期化 do i=1,kiza_x q(1,i)=rho(i) q(2,i)=rho(i)*u(i) t(i)=p(i)/(rho(i)*r) q(3,i)=rho(i)*(cv*t(i)+0.5*u(i)**2) ene(i)=q(3,i) ! print*, q(1,i),q(2,i),q(3,i) end do !保存量の初期化：3つめの全エネルギが合っているかあやしい。密度×定圧比熱×静温なんだけど。 !eneに全エネルギを代入してる。 do k=1,100 do i=1,kiza_x e(1,i)=rho(i)*u(i) e(2,i)=p(i)+rho(i)*u(i)**2 e(3,i)=(ene(i)+p(i))*u(i) ! print*, e(1,i),e(2,i),e(3,i) end do !流束ベクトルの計算 do i=1,kiza_x h(i)=ene(i)+p(i) end do !全エンタルピの計算 do i=1,kiza_x-1 rho_av(i)=sqrt(rho(i)*rho(i+1)) u_av(i)=(sqrt(rho(i))*u(i)+sqrt(rho(i+1))*u(i+1))/(sqrt(rho(i))+sqrt(rho(i+1))) h_av(i)=(sqrt(rho(i))*h(i)+sqrt(rho(i+1))*h(i+1))/(sqrt(rho(i))+sqrt(rho(i+1))) c_av(i)=sqrt((gamma-1)*(h_av(i)-0.5*u_av(i)**2)) !print*, rho_av(i),u_av(i),h_av(i),c_av(i) end do !Roe平均の密度、速度、エンタルピ,特性速度の計算 do i=1,kiza_x-1 c(1,i)=u_av(i)-c_av(i) c(2,i)=u_av(i) c(3,i)=u_av(i)+c_av(i) end do !固有値ベクトル do i=1,kiza_x-1 right(1,1:3,i)=1.0 right(2,1,i)=u_av(i)-c_av(i) right(2,2,i)=u_av(i) right(2,3,i)=u_av(i)+c_av(i) right(3,1,i)=h_av(i)-u_av(i)*c_av(i) right(3,2,i)=0.5*u_av(i)**2 right(3,3,i)=h_av(i)+u_av(i)*c_av(i) end do !右固有ベクトルの計算 do i=1,kiza_x-1 b(1)=(u_av(i)**2)/2.0*(gamma-1)/(c_av(i)**2) b(2)=(gamma-1)/(c_av(i)**2) left(1,1,i)=0.5*(b(1)+u_av(i)/c_av(i)) left(1,2,i)=-0.5*(1/c_av(i)+b(2)*u_av(i)) left(1,3,i)=0.5*b(2) left(2,1,i)=1-b(1) left(2,2,i)=b(2)*u_av(i) left(2,3,i)=-b(2) left(3,1,i)=0.5*(b(1)-u_av(i)/c_av(i)) left(3,2,i)=0.5*(1/c_av(i)-b(2)*u_av(i)) left(3,3,i)=0.5*b(2) end do !左固有ベクトルの計算 do i=1,kiza_x-1 !alpha(3,i)=q(1:3,i+1)-q(1:3,i) alpha(1:3,i)=matmul(left(1:3,1:3,i),(q(1:3,i+1)-q(1:3,i))) end do !物理量の跳びの計算。matmul()を使用 do i=1,kiza_x-1 sigma(1,i)=sig((u_av(i)-c_av(i)),psi(u_av(i)-c_av(i)),dt,dx) sigma(2,i)=sig(u_av(i),psi(u_av(i)),dt,dx) sigma(3,i)=sig((u_av(i)+c_av(i)),psi(u_av(i)+c_av(i)),dt,dx) end do !σの計算 do i=2,kiza_x-1 do j=1,3 g(j,i)=minmod(alpha(j,i),alpha(j,i-1)) end do end do !gの計算。最初の一つめの配列は未使用 do i=2,kiza_x-2 do j=1,3 if(alpha(j,i) /= 0) then gam(j,i)=sigma(j,i)*(g(j,i+1)-g(j,i))/alpha(j,i) else gam(j,i)=sigma(j,i)*0 end if end do end do !gam()の計算 do i=2,kiza_x-2 do j=1,3 phi(j,i)=sigma(j,i)*(g(j,i+1)-g(j,i))-psi(c(j,i)+gam(j,i))*alpha(j,i) end do end do !Φの計算 do i=2,kiza_x-2 ee(1:3,i)=0.5*((e(1:3,i+1)+e(1:3,i))+matmul(right(1:3,1:3,i),phi(1:3,i))) end do !数値流束ee()の計算 do i=3,kiza_x-2 q(1:3,i)=q(1:3,i)-dt/dx*(ee(1:3,i)-ee(1:3,i-1)) end do !次ステップのQの計算 open(10,file='result.txt') do i=1,kiza_x rho(i)=q(1,i) u(i)=q(2,i)/q(1,i) p(i)=(gamma-1)*(q(3,i)-0.5*q(1,i)*q(2,i)) t(i)=p(i)/(rho(i)*r) write(10,*) rho(i),' ',u(i),' ',p(i),' ',t(i) end do close(10) !q()から物理量を算出している end do end real function psi(z) real :: z,delta delta=10 if(abs(z) >= delta) then psi=abs(z) else psi=(z**2+delta**2)/(2.0*delta) end if end !Ψの関数 real function sig(z,psi,dt,dx) real :: z,psi,dt,dx sig=0.5*(psi-dt/dx*z**2) end !σの関数 real function minmod(x,y) real :: x,y if(abs(x)<abs(y)) then minmod=x else if(abs(x)>abs(y)) then minmod=y else minmod=0 end if end !minmodの関数
コメント