0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@rino_yume(rino yume)

Fortranの次を考える r001 ― Docker + Julia 1.13.1の数値計算環境を作る

0
Posted at

Fortranの次を考える r001 ― Docker + Julia 1.13.1の数値計算環境を作る

はじめに

数値計算ではやはりfortranに頼ってしまう今日この頃ですが、前回の記事は、河川・氾濫・流出などの過去の資産であるfortran数値計算を今後どうするか？どのような言語へ変換したり、実装していくかを考え、まずはJuliaを試してみることにしました。

今回はその最初の版として r001 を作ります。

まだ河川計算には入りません。最初に確認するのは、次の4点です。

  • Docker上でJulia 1.13.1を固定して動かせること
  • 配列をその場で更新できること
  • 差分法を使った小さな数値計算を書けること
  • テストで計算条件と保存性を確認できること

最初から1D河川や2D氾濫へ進むと、結果が合わないときに原因を切り分けにくくなります。

r001 では、今後使うことになる「配列更新」「時間刻み」「境界条件」「CFL条件」「保存量の確認」を、小さな計算で一度通しておきます。

1. 今回使うJulia

2026年10月5日時点で、Julia公式サイトのCurrent Stable Releaseは Julia 1.13.1 です。リリース日は2026年9月25日です。

Docker Official Imageにも 1.13.1-bookworm が用意されているため、今回はこのタグを固定して使います。

Julia 1.13.1
    ↓
Docker Official Image
    ↓
julia:1.13.1-bookworm

latest を使わず、バージョンとベースイメージを固定しておくことで、後から同じ環境を再現しやすくします。

参考:

2. r001で作るもの

今回の構成は次のようにしました。

20261006_j01.png

最初の版なので、外部Juliaパッケージは必須にしていません。

Julia本体と標準ライブラリだけで動かします。

これは機能を増やさないためというより、最初の検証で問題が起きたときに、原因をできるだけ少なくするためです。

3. ディレクトリ構成

今回のプロジェクトは次のようにしました。
参考のために、ZIPファイルを置いておきますので、ダウンロードして活用してください。

julia_hydraulic_r001/
├── Dockerfile
├── compose.yaml
├── Project.toml
├── README.md
├── src/
│   ├── JuliaHydraulicR001.jl
│   ├── array_ops.jl
│   ├── advection.jl
│   └── io.jl
├── examples/
│   ├── array_demo.jl
│   └── advection_demo.jl
├── test/
│   └── runtests.jl
├── benchmark/
│   └── benchmark.jl
├── scripts/
│   ├── run_all.jl
│   └── verify.sh
└── output/
    └── .gitkeep

src には計算本体を置き、examples には動作確認用の小さな実行例を置いています。

テストは test/runtests.jl に分けました。

今後 r002 で1D河道計算へ進むときも、この構成をそのまま使えます。

4. Dockerfile

Dockerfileはかなり小さくしています。

FROM julia:1.13.1-bookworm

WORKDIR /workspace

COPY Project.toml ./
COPY src ./src
RUN julia --project=. -e 'using JuliaHydraulicR001; println("JuliaHydraulicR001 precompile: OK")'

COPY . .

CMD ["julia", "--project=.", "scripts/run_all.jl"]

ビルド時に一度パッケージを読み込み、少なくともモジュールとして読み込める構成になっているかを確認する形にしています。

今回は外部依存パッケージがないので、パッケージ取得処理は入れていません。

5. Compose設定

compose.yaml は次のようにしています。

services:
  julia:
    build:
      context: .
      dockerfile: Dockerfile
    working_dir: /workspace
    user: "${LOCAL_UID:-1000}:${LOCAL_GID:-1000}"
    environment:
      HOME: /tmp/julia-home
      JULIA_DEPOT_PATH: /tmp/julia-depot
      JULIA_NUM_THREADS: "${JULIA_NUM_THREADS:-auto}"
    volumes:
      - .:/workspace
    command: ["julia", "--project=.", "scripts/run_all.jl"]

WSL2やLinuxでbind mountを使うと、コンテナがrootで作ったファイルがホスト側でもroot所有になることがあります。

今回は結果CSVをホスト側へ出力するため、UID/GIDをComposeへ渡せるようにしています。

実行前に次を設定します。

export LOCAL_UID=$(id -u)
export LOCAL_GID=$(id -g)

JULIA_DEPOT_PATH は /tmp 側へ分けています。

今回のように外部パッケージを使わない段階では大きな差はありませんが、プロジェクトディレクトリへJuliaのキャッシュを混ぜないようにしています。

6. Juliaプロジェクト

Project.toml は次の内容です。

name = "JuliaHydraulicR001"
uuid = "0d3f9f72-0f57-4b30-a775-5c26f42f7d2d"
authors = ["Qiita sample project"]
version = "0.1.0"

[compat]
julia = "1.13"

[extras]
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[targets]
test = ["Test"]

この段階では、計算ライブラリは追加していません。

まずJulia本体だけで、どこまで素直に数値計算を書けるかを確認します。

7. 配列更新を作る

最初に、将来の保存量更新を意識した小さな関数を作りました。

function add_source!(
    state::AbstractArray{T},
    source::AbstractArray{T},
    dt::Real,
) where {T<:AbstractFloat}
    axes(state) == axes(source) ||
        throw(DimensionMismatch("state and source must have identical axes"))
    dt >= 0 || throw(ArgumentError("dt must be non-negative"))

    dt_t = T(dt)
    for i in eachindex(state, source)
        state[i] += dt_t * source[i]
    end

    return state
end

やっていることは単純です。

stateⁿ⁺¹ = stateⁿ + Δt × source

将来、2D氾濫計算で水深や保存量を更新するときにも、同じような形の処理が多く出てきます。

関数名の末尾に ! を付けているのは、引数 state をその場で変更する関数であることを分かりやすくするためです。

8. 配列計算のサンプル

examples/array_demo.jl では、4×3の配列を使います。

using JuliaHydraulicR001

state = fill(1.0, 4, 3)
source = reshape(collect(0.0:0.1:1.1), 4, 3)
dt = 0.5

add_source!(state, source, dt)

display(state)
println(total_amount(state, 1.0))

ここでは、2次元配列をセル中心値に見立てています。

まだ水理計算ではありませんが、配列をその場で更新する基本形を確認できます。

9. 簡単な数値計算として1D線形移流を使う

次に、小さな差分計算を入れました。

対象は1次元線形移流方程式です。

∂u/∂t + c ∂u/∂x = 0

c は一定の移流速度です。

この式自体は河川のSaint-Venant方程式ではありません。

ただし、今後必要になる次の要素を確認できます。

  • 1次元格子
  • 隣接セル参照
  • 時間更新
  • 境界条件
  • CFL条件
  • 保存量

そのため、r001 の確認用として使うことにしました。

10. 一次風上差分

正の流速 c > 0 の場合、一次風上差分は次の形です。

uᵢⁿ⁺¹ = uᵢⁿ - σ (uᵢⁿ - uᵢ₋₁ⁿ)

ここで、

σ = c Δt / Δx

です。

今回の実装では、

|σ| <= 1

を確認します。

負の流速では、風上側が反対になるため i+1 側を参照します。

周期境界条件を使い、左端と右端を接続しています。

20261006_j02.png

この境界条件にしておくと、領域の外へ値が出ていかないため、総量保存を確認しやすくなります。

11. 1ステップ更新

計算本体は upwind_step! に分けました。

function upwind_step!(
    u_next::AbstractVector{T},
    u::AbstractVector{T},
    velocity::Real,
    dt::Real,
    dx::Real,
) where {T<:AbstractFloat}
    length(u_next) == length(u) ||
        throw(DimensionMismatch("u_next and u must have the same length"))

    sigma = Float64(velocity) * Float64(dt) / Float64(dx)

    abs(sigma) <= 1.0 + 16 * eps(Float64) ||
        throw(ArgumentError("CFL condition violated"))

    n = length(u)

    if sigma >= 0
        for i in eachindex(u)
            im1 = i == 1 ? n : i - 1
            u_next[i] = u[i] - sigma * (u[i] - u[im1])
        end
    else
        for i in eachindex(u)
            ip1 = i == n ? 1 : i + 1
            u_next[i] = u[i] - sigma * (u[ip1] - u[i])
        end
    end

    return u_next
end

u と u_next の2本を用意し、タイムステップごとに入れ替えます。

毎回新しい配列を作らず、用意した配列を再利用する形です。

12. 終了時刻に合わせて時間刻みを調整する

入力したCFL数から、まず最大時間刻みを求めます。

Δt_limit = CFL × Δx / |c|

ただし、この値をそのまま使うと、最後の1ステップだけ終了時刻を越える場合があります。

そこで、必要ステップ数を整数で決めてから、実際の時間刻みを再計算しています。

dt_limit = config.cfl * dx / abs(config.velocity)
nsteps = ceil(Int, config.final_time / dt_limit)
dt = config.final_time / nsteps

これで、指定した終了時刻にちょうど到達します。

また、再計算した dt は dt_limit 以下になるため、設定したCFL数を越えません。

13. 初期条件

初期条件にはガウス分布を使っています。

function gaussian_initial_condition(
    x::Real;
    center::Real = 0.30,
    width::Real = 0.06,
    amplitude::Real = 1.0,
)
    width > 0 || throw(ArgumentError("width must be positive"))

    return Float64(amplitude) *
        exp(-0.5 * ((Float64(x) - Float64(center)) / Float64(width))^2)
end

矩形波でも計算できますが、通常のサンプルでは滑らかな分布を使っています。

14. CSVへ保存する

結果は output/advection_result.csv に保存します。

列は次の3列です。

x,u_initial,u_final

外部CSVパッケージは使わず、標準のファイルI/Oだけで書いています。

open(path, "w") do io
    println(io, "x,u_initial,u_final")

    for i in eachindex(result.x, result.u0, result.u)
        println(io, "$(result.x[i]),$(result.u0[i]),$(result.u[i])")
    end
end

r001 ではこれで十分です。

データ量が増えたり、GeoPackageやZarrなどと連携する段階になったら、出力形式は別途整理します。

15. テスト

数値計算では、実行できるだけでは不十分です。

今回は次をテストしています。

テスト 確認内容
配列更新 state += dt * source が正しいか
配列サイズ 異なるサイズを渡した場合に停止するか
正流速 CFL=1で1セル右へ移動するか
負流速 CFL=1で逆方向へ移動するか
CFL条件 1を越えた場合に停止するか
総量 周期境界で初期・終了総量が一致するか
異常値 NaN / Infが出ていないか
入力条件 不正なセル数、流速、終了時刻を拒否するか

特に分かりやすいのが、CFL=1のテストです。

正の流速で、

[1, 2, 3, 4]

を1ステップ進めると、周期境界では、

[4, 1, 2, 3]

になります。

これは近似値ではなく、今回の離散式ではそのまま1セル移動するケースです。

u = [1.0, 2.0, 3.0, 4.0]
u_next = similar(u)

upwind_step!(u_next, u, 1.0, 1.0, 1.0)

@test u_next == [4.0, 1.0, 2.0, 3.0]

16. 総量保存を確認する

周期境界の線形移流では、今回の離散式で総量が保存されることを確認します。

離散的な総量は、

Σ uᵢ Δx

としています。

テストでは、初期値と終了値を、

@test isapprox(
    result.mass_final,
    result.mass0;
    rtol = 1e-12,
    atol = 1e-12,
)

で確認します。

ここで確認したいのは、「波形が完全に同じ形のまま移動すること」ではありません。

一次風上差分には数値拡散があるため、CFLが1でない場合には波形が少しなまります。

r001 では、その性質を隠さず、総量が保たれていることと、計算が安定していることを確認します。

17. 実行方法

まずUID/GIDを設定します。

export LOCAL_UID=$(id -u)
export LOCAL_GID=$(id -g)

イメージをビルドします。

docker compose build

Juliaのバージョンを確認します。

docker compose run --rm julia julia --version

テストは次です。

docker compose run --rm julia \
  julia --project=. -e 'using Pkg; Pkg.test()'

配列計算のサンプルは次です。

docker compose run --rm julia \
  julia --project=. examples/array_demo.jl

1D移流計算は次です。

docker compose run --rm julia \
  julia --project=. examples/advection_demo.jl

まとめてデモを実行する場合は、

docker compose run --rm julia

です。

18. 一括検証スクリプト

確認漏れを減らすために、scripts/verify.sh も用意しました。

./scripts/verify.sh

このスクリプトでは、順番に次を実行します。

Docker / Compose バージョン
        ↓
Docker build
        ↓
Julia version
        ↓
Pkg.test()
        ↓
配列サンプル
        ↓
1D移流サンプル
        ↓
簡易計測

途中でエラーが出た場合は、その時点で終了します。

19. 計算時間の確認

JuliaはJITコンパイルを使うため、最初の呼び出し時間をそのまま数値計算の実行時間として比較しないようにします。

benchmark/benchmark.jl では、最初に1回 upwind_step! を呼び出してから計測しています。

upwind_step!(u_next, u, velocity, dt, dx)

その後、同じカーネルを繰り返し実行します。

今回の簡易計測は、性能を正式に評価するためのものではありません。

r002 以降でFortran版と速度を比較する段階では、BenchmarkTools.jlなども含めて計測方法を整理します。

20. 今回の確認状況

ここは明確に分けておきます。

今回の配布物作成環境では、DockerとJuliaの実行環境が利用できなかったため、次の2点は実行済みとはしていません。

  • docker compose build
  • Julia上での Pkg.test()

一方で、次は確認しています。

  • Project.toml のTOML構文
  • compose.yaml のYAML構文
  • 正負流速に対する一次風上差分の更新式
  • CFL=1での1セル移動
  • サンプル条件でのCFL値
  • 周期境界での離散総量保存の数値ロジック

実機での最終確認は、同梱した scripts/verify.sh を使う形にしています。

未実行のものをPASSとして扱わないようにしました。

21. r001でやらないこと

今回の版では、次はまだ扱いません。

  • Saint-Venant方程式
  • Local Inertial法
  • 河道断面
  • Manning粗度
  • 上流流量境界
  • 下流水位境界
  • Node–Reach河道網
  • 合流・分岐
  • GPU

これらは次の段階です。

r001 で確認したかったのは、Juliaで数値計算コードを書くときの基本的な流れです。

22. 次はr002へ

次の r002 では、今回の線形移流モデルを卒業して、単一河道の1D計算へ進みます。

流れは次のように考えています。

20261006_j03.png

r002 では、既存のFortran版と同じ条件を与え、流量・水位・水収支を比較できる形にします。

ここから先が、本格的な水理計算になります。

23. まとめ

r001 では、Juliaを使った数値計算の土台を作りました。

今回入れたのは、

DockerでJulia 1.13.1を固定
        ↓
配列更新
        ↓
1D差分計算
        ↓
CFL確認
        ↓
保存量確認
        ↓
テスト
        ↓
CSV出力

までです。

内容としては小さいですが、次の河川計算へ進む前に確認しておきたい部分を分けて実装できました。

特に、今後も「計算できた」だけで終わらせず、境界条件、保存量、基準解との比較をテストとして残していきます。

次は r002 として、単一河道の1D河川計算をJuliaで実装します。

参考資料

0
1
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?