時間積分 - Kimmy's Lab.

テキスト「空力音の直接数値計算を支える技術」（暗黙の型宣言著）の補足課題3

7章まで読んだらやってみよう！

課題

(1) 時間積分にEuler法を用いて移流方程式を計算する．

初期値はテキスト式(7.11)とし，波が境界を出ない時間まで計算する．
境界では勾配ゼロ(df/dx=0)の境界条件を与える．
- 式7.11は境界でほぼ値がゼロ(f=0)だが，境界条件としてf=0を無理やり与えると境界で値が発散する．
- テキストでは周期境界条件を用いて，波が2周する時間分計算しているが，実際の空力音の数値計算では周期境界条件を用いることができるケースというのは限られるため，ここでは周期境界条件を用いず勾配ゼロの境界条件を与える．
誤差も計算して，実効的な精度を確認する．
テキストの以下プログラム参照．
- List7.2(移流方程式メインルーチン)
- List7.3(Euler法による時間積分)
空間微分は，補足課題2を参考に，以下を用いる．
- グローバルスキーム・・・6次精度
- 境界1点内側・・・5次精度
- 境界・・・4次精度片側差分

(2) 同様に，TVD Runge-Kutta法（以下，TVDRK）を用いて移流方程式を計算する．

List7.14(TVDRK法による時間積分)．

解答例

まずはEuler法．テキストを写経して，以下のファイルを作成する．

main.f90(List7.3参照)
module_TimeIntegration.f90(List7.4参照)

main.f90を以下のように修正する．（修正箇所は「!---」, ただのコメントは「!」）

    real(8), parameter :: EndT = 10d0 !--- 計算時間を，波が境界に到達しない時間に設定
    :
    !--- 境界条件用パラメータを追加
    real(8), parameter :: df1 = 0d0
    real(8), parameter :: dfN = 0d0

   !--- 3重対角行列を指定
    type(TriDiagonalMatrix) :: G

    :

    Initialize : block
        
        :

        !--- 両端で勾配ゼロになるようオフセットする(3.8章参照)
        f%var(1) = f%var(2)
        f%var(Nx)= f%var(Nx-1)
    end block Initialize

    :

    do n = 0, Nt-1
        Euler : block
            call compute1stDerivative(f%dx, f%var, dx, G, df1, dfN)  !--- 境界条件用パラメータdf1, dfNを追加
            call integrateEuler(f%var, f%dx, dt, conv)
        end block Euler
    end do

    !--- 誤差計算・出力用ブロックを追加
    ResultOut : block
        integer, parameter :: Nx_err = nint(wavelength/dx) ! 波長を誤差計算範囲とする
        integer :: i, i_wavecenter_f, i_wavecenter_f0
        real :: err

        err = 0d0

        i_wavecenter_f0= nint(Nx/2d0) + 1                                           ! 初期の波f0の中心位置
        i_wavecenter_f = i_wavecenter_f0 + nint(conv*dt*(Nt-1)/dx)    ! 伝播した波fの中心位置

        open(1, file='FileOut.txt', status='replace')
    
        do i=1,Nx
            write(1,*)  dble(i-1)*dx - Lx/2d0, f0%var(i), f%var(i)             ! x, f0, f を出力
        enddo

        do i=-Nx_err, Nx_err
            err = err + abs(f%var(i_wavecenter_f+i)-f0%var(i_wavecenter_f0+i))  ! 初期の波と伝播した波の振幅差を計算
        enddo
        print *, "MinAveError:", err/(Nx_err*2+1) ! 累積誤差を格子点数で割って，最小平均誤差を出力

        close(1)
    end block ResultOut

    :

    contains

    subroutine compute1stDerivative(d_f_dx, f, dx, G, GraditntB1, GraditntBN) !--- 境界条件用パラメータを追加
        use TypeTriDiagonalMatrix !--- 3重対角行列を指定
        implicit none

        :

        ! 境界条件用パラメータを追加
        real(8), intent(in   ), optional :: GraditntB1
        real(8), intent(in   ), optional :: GraditntBN

        call construct1stDerivativeVector(d_f_dx, f, dx, G, GraditntB1, GraditntBN)    !--- 境界条件用パラメータを追加
        call solve(G, d_f_dx)
    end subroutine compute1stDerivative

end program main

空間微分には，補足課題2で作成した以下のファイルを用いる．

module_CompactScheme.f90(List3.4参照)
TypeTriDiagonalMatrix.f90(List4.2参照)

以上でEuler法のプログラムは完成．

パラメータ「CFL」をさまざまに変化させて，誤差の変化を確認する．確認結果は後述．

次に，TVDRK法．main.f90のEuler法用ブロックを，TVDRK法用ブロックに変更．

main.f90

    :
    do n = 0, Nt-1
        TVDRK : block
            call cacheTVDRK(f%var)            ! fnの値を保持
            do step=1,3
                call compute1stDerivative(f%dx, f%var, dx, G, df1, dfN)
                call integrateTVDRK(f%var, f%dx, dt, step, conv)   ! TVDRKのサブルーチンをcall
            enddo
        end block TVDRK
    end do
    :

さらに，module_TimeIntegration.f90にテキストList7.14のサブルーチンcacheTVDRKとintegrateTVDRKを追加する．

テキストにも書いてあるように，Euler法からTVDRK法への変更は簡単♪

TVDRK法は，Euler法から得られた結果と現時間ステップの結果に重みづけして，次の時間ステップの値を決めている，と解釈できる．

パラメータ「CFL」をさまざまに変化させて，誤差の変化を確認した結果は図1のとおり．

f:id:kimmy123:20200531221009p:plain — 図1 移流方程式誤差

Euler法はslope 2（CFLを1/10にすると誤差が1/100）よりも傾きが浅い．近似直線を引いてみると，約1.1次精度であった．

TVDRK法(Nx = 400)は，0.1<CFL<0.6でほぼslope3に沿っているので約3次精度． 0.6<CFLでslope3より誤差が大きいのは，格子が荒すぎる影響．

CFL<0.1で誤差の減少がサチっているのは？

Nx=400では0.00001付近でサチっているのに対し，格子を2倍荒くしたNx=200では誤差が0.001付近でサチっている． CFL=0.002では，Nx=400の誤差は，Nx=200の誤差の約 $2^{-6}$ 倍であった．空間微分のグローバルスキームが6次精度なので，格子を2倍にすることで誤差が $2^{-6}$ 倍に減少している．つまり，誤差の減少がサチるのは，空間微分の精度の影響．