吸収境界条件の実装について/About implementation of absorbing boundary conditions

[msomeya1]

お世話になっております。地震研M2の染矢です。

線形分散波をシミュレーションするために、吸収境界条件をONにした計算を試みています（nxa=nya=60, apara=0.002）。 すると、30分を過ぎたあたりから計算領域の右下で波が爆発し始め、やがて振幅が1e6を超えて、計算が停止します。

おそらく、スポンジ（ダンパー）の定義に問題があると思います。 以下は、mod_rwg.f90のサブルーチンmake_abcからの抜粋です。

do j = 1, nlat
    do i = 1, nlon
        if(i < nxa) then
            tmp = exp(-((apara*(nxa-i       ))**2))
        else if(i > (nlon-nxa+1)) then
            tmp = exp(-((apara*(i-nlon+nxa-1))**2))
        else if(j < nya) then
            tmp = exp(-((apara*(nya-j       ))**2))
        else if(j > (nlat-nya+1)) then
            tmp = exp(-((apara*(j-nlat+nya-1))**2))
        else
            tmp = 1.0d0
        end if
    abc(i,j) = tmp
    end do
end do

nlon=100, nlat=150, nxa=nya=30, apara=0.002としてabcを描いてみると、四隅のスポンジが滑らかに接続されていないことがわかります。

先ほどの波が爆発した箇所は、ちょうど不連続にabcが飛ぶ箇所に対応していると思います。

そこで、make_abcを次のように改良してみます。

do j = 1, nlat
    jj = min(j, nlat+1-j)
    do i = 1, nlon
        ii = min(i, nlon+1-i)
        if(ii <= nxa .and. jj > nya) then
            tmp = exp(-((apara*(nxa-ii))**2))
        else if(ii > nxa .and. jj <= nya) then
            tmp = exp(-((apara*(nya-jj))**2))
        else if(ii <= nxa .and. jj <= nya) then
            tmp = exp(-((apara*(nxa-ii))**2+(apara*(nya-jj))**2))
        else
            tmp = 1.0d0
        end if
        abc(i,j) = tmp
    end do
end do

このコードでは、四辺（四隅以外）での処理はこれまでと同じですが、四隅ではユークリッド距離に応じて減衰の仕方を決めることにより、スポンジを滑らかに接続します。

この場合、最後まで爆発することなく、計算は正常に終了します。

以上、ご確認いただければ幸いです。

[jagurs-admin]

染谷さま

Bjどうもありがとうございます。 次の改†の•rに取り入れようと思います。ありがとうございました。大‰浃辘郡い扦埂 引きAきよろしくおŠいします。

Rˆ

---

差出人: 'msomeya1' via jagurs_admin @.> 送信日•r: Sunday, June 9, 2024 4:03:23 PM 宛先: jagurs-admin/jagurs @.> CC: Subscribed @.***> 件名: [jagurs-admin/jagurs] 吸…Ь辰缣跫Œg装について/About implementation of absorbing boundary conditions (Issue #20)

お世’になっております。地震研M2の染矢です。

€形分散波をシミュレ`ションするために、吸…Ь辰缣跫ONにした‹算を‡みています（nxa=nya=60, apara=0.002）。 すると、30分を^ぎたあたりから‹算I域の右下で波が爆kし始め、やがて振幅が1e6を超えて、‹算が停止します。

unstable.gif (view on web)https://github.com/jagurs-admin/jagurs/assets/118101078/2f77ddd3-4e45-48d4-85e9-84009e53af21

おそらく、スポンジ（ダンパ）の定xに†–}があると思います。 以下は、mod_rwg.f90のサブルチンmake_abcからの’i‚です。

do j = 1, nlat do i = 1, nlon if(i < nxa) then tmp = exp(-((apara*(nxa-i ))2)) else if(i > (nlon-nxa+1)) then tmp = exp(-((apara*(i-nlon+nxa-1))*2)) else if(j < nya) then tmp = exp(-((apara(nya-j ))2)) else if(j > (nlat-nya+1)) then tmp = exp(-((apara*(j-nlat+nya-1))**2)) else tmp = 1.0d0 end if abc(i,j) = tmp end do end do

nlon=100, nlat=150, nxa=nya=30, apara=0.002としてabcを描いてみると、四隅のスポンジが滑らかに接Aされていないことがわかります。

abc1.png (view on web)https://github.com/jagurs-admin/jagurs/assets/118101078/25aec620-b3dd-446e-8144-eea68660883e

先ほどの波が爆kしたw所は、ちょうど不BAにabcがwぶw所にŒ辘筏皮い毪人激い蓼埂

そこで、make_abcを次のように改良してみます。

do j = 1, nlat jj = min(j, nlat+1-j) do i = 1, nlon ii = min(i, nlon+1-i) if(ii <= nxa .and. jj > nya) then tmp = exp(-((apara*(nxa-ii))2)) else if(ii > nxa .and. jj <= nya) then tmp = exp(-((apara*(nya-jj))*2)) else if(ii <= nxa .and. jj <= nya) then tmp = exp(-((apara(nxa-ii))2+(apara*(nya-jj))**2)) else tmp = 1.0d0 end if abc(i,j) = tmp end do end do

このコドでは、四x（四隅以外）での„I理はこれまでと同じですが、四隅ではユクリッド距xに辘袱œp衰の仕方を›Qめることにより、スポンジを滑らかに接Aします。

abc2.png (view on web)https://github.com/jagurs-admin/jagurs/assets/118101078/1341b704-3c29-436f-be80-ed5166dca434

このˆ龊稀⒆钺幛蓼潜kすることなく、‹算は正常にK了します。

stable.gif (view on web)https://github.com/jagurs-admin/jagurs/assets/118101078/4c5c1bf3-ecb2-4f24-bc99-9dad455499c3

以上、ご_Jいただければ幸いです。

― Reply to this email directly, view it on GitHubhttps://github.com/jagurs-admin/jagurs/issues/20, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AEKCJZA36PVGJHFCRLLYVYDZGP43XAVCNFSM6AAAAABJAUNTXWVHI2DSMVQWIX3LMV43ASLTON2WKOZSGM2DEMBWGI2TSMQ. You are receiving this because you are subscribed to this thread.Message ID: @.***>

To unsubscribe from this group and stop receiving emails from it, send an email to @.**@.>.

[msomeya1]

馬場先生

ありがとうございます。 （一部文字化けしてしまっているのですが、次回改訂時に取り入れて頂けるということは読むことができました）

もう一つ、計算の初期条件について質問があります。 JAGURSでは、波高の初期値η(x,y,t=0)が与えられたとき（初期流束M(x,y,t=0), N(x,y,t=0)は0とします）、その後の時刻の波動場を計算することはできると思います。 （初期地殻変動ファイルとしてη(x,y,t=0)を与えることで実現可能） では、η(x,y,t=0)だけでなく、non-zeroの初期流束M(x,y,t=0), N(x,y,t=0)が与えられたとき、その後の時刻の波動場を計算するにはどうしたらよいでしょうか。 （leap-frogやstaggered-gridであるため、格子点や時刻を半分ずらして用意する必要はありますが、この点は事前にクリアできているものとします） ※最も簡単な線形長波での計算を想定しています

私が思いつく最も簡単な方法は、計算のリスタート機能を利用（乱用？）することです。 リスタートのためには計算終了時点での波高や流束などの情報が必要ですので、 逆に与えたい初期条件ファイルをリスタート情報ファイルに変換することができれば、リスタート機能によって任意のη, M, Nから計算を開始できるのではないか？ということです。

このようなことは可能でしょうか？ また、可能である場合、注意点などはありますでしょうか？

重ねての質問となり大変恐縮ですが、ご確認いただけますと幸いです。 よろしくお願いいたします。

染矢

[jagurs-admin]

Hi Someya-san,

Yes, I think you can use restart function for that purpose. Also, you can directly e-mail me to avoid garbled characters. @.**@.>

Regards,

-- Baba, 088-656-9721 (ext. 4231)

[msomeya1]

Thank you for the immediate reply. I will try the restart function.

Also, you can directly e-mail me to avoid garbled characters.

I see. Next time I ask a question, I will use email.

Someya