0
0

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の次を考える r002 ― Juliaで単一河道1Dを実装し、REAL*4 Fortranと比較する

0
Posted at

Fortranの次を考える r002 ― Juliaで単一河道1Dを実装し、REAL*4 Fortranと比較する

前回の r001 では、Julia 1.13.1をDocker上で使うための最小環境を作り、配列更新と1次元線形移流まで確認できる構成にしました。

今回はその次の段階として、単一河道の1D水理計算へ進みます。

ただし、Juliaで新しく計算できたというだけでは確認として不十分です。

手元にはFortranで積み上げてきた数値計算があります。今後Juliaへ移していくのであれば、同じ条件を与えて、少なくとも

  • 流量
  • 水位
  • 水収支

を比較できる状態を先に作っておいた方が安心です。

そこで r002 では、単一河道に範囲を絞り、JuliaとFortranが同じ入力ファイルを読み、同じ離散式で計算し、全保存時刻の結果を比較できるようにしました。

r002 の目的は、まず1本の河道で比較の仕組みを固めることです。

1. r002でやること

今回の流れは次のようにしました。

20261013_f01.png

比較条件をコード内に二重に書かず、case/single_reach_case.csv を正本にします。

FortranもJuliaも、この同じCSVを読みます。

さらに、元になった既存の single_reach_backwater.json も case/single_reach_backwater_source.json として同梱し、SHA-256を固定しています。preflight.py では、共有CSVの主要条件がこの元JSONと一致していることも検査します。

これで、

実は時間刻みが違っていた

Manning粗度係数が片方だけ違っていた

といった比較以前の問題を減らせます。

2. 比較に使う単一河道

今回の条件は、これまで背水・逆流の最小確認に使っていた単一河道ケースを基にしています。

UP
│
│  L = 1000 m
│  5 segments
│
▼
OUT

主要条件は次のとおりです。

項目 値
計算時間 600 s
時間刻み 2 s
保存間隔 20 s
河道延長 1000 m
分割数 5
上流河床高 1.0 m
下流河床高 0.5 m
河床幅 10.0 m
法勾配 H:V 1.0
Manning n 0.03
初期水位 3.0 m
上流流量 0.0 m³/s
下流水位 3.0 → 4.5 m

下流水位は600秒かけて線形に上昇させます。

上流から流量を与えないため、このケースでは下流側の水位上昇によって背水が伝わり、負流量が発生するかどうかを確認できます。

20261013_f02.png

3. 今回の1D計算

r002 では、Local Inertial型の式を使います。

連続式は、

∂A/∂t + ∂Q/∂x = 0

です。

  • A: 流水断面積 [m²]
  • Q: 流量 [m³/s]

運動量式は、移流加速度を省略して、

∂Q/∂t + g A ∂H/∂x + g n² Q|Q| / (A R^(4/3)) = 0

としています。

  • H: 水位 [m]
  • R: 径深 [m]
  • n: Manning粗度係数

今回の目的は新しい水理式を提案することではなく、FortranとJuliaの比較経路を作ることです。

そのため、計算式は意図的に小さくしています。

4. 水位と流量を交互に配置する

計算では、水位をセル中心、流量をセル境界に置きます。

5セルの場合、概念的には次の配置です。

Q0    H1    Q1    H2    Q2    H3    Q3    H4    Q4    H5    Q5
|-----●-----|-----●-----|-----●-----|-----●-----|-----●-----|

上流端は Q0 を流量境界として与えます。

下流端は、外側に指定水位を置き、最後のセルとの水面勾配から Q5 を計算します。

20261013_f03.png

5. 台形断面

断面は台形として扱います。

水深を h、河床幅を b、法勾配を z とすると、流水断面積は

A = b h + z h²

です。

潤辺は、

P = b + 2 h sqrt(1 + z²)

なので、径深は

R = A / P

です。

Julia側では、次のように独立した関数にしています。

function trapezoid_area(depth::Float32, bottom_width::Float32, side_slope::Float32)
    return bottom_width * depth + side_slope * depth * depth
end

function hydraulic_radius(depth::Float32, bottom_width::Float32, side_slope::Float32)
    area = trapezoid_area(depth, bottom_width, side_slope)
    perimeter = bottom_width + 2f0 * depth * sqrt(1f0 + side_slope * side_slope)
    return area / perimeter
end

面積から水深へ戻す処理も解析式で実装しています。

水収支を断面積で更新したあと、水位へ戻すときに使います。

6. REAL*4とFloat32を合わせる

今回の比較でかなり大事にした点です。

Fortran側では、

use iso_fortran_env, only: real32
real(real32) :: stage, q, dt

としています。

Julia側も、

Float32

に統一しました。

つまり、

Fortran REAL*4相当
        ↓
real(real32)

Julia
        ↓
Float32

です。

Juliaだけ Float64 にしてしまうと、単純な移植比較なのか、浮動小数点精度を変えた影響なのかを切り分けにくくなります。

r002 ではまず精度をそろえます。

7. ディレクトリ構成

プロジェクトは次の構成にしました。

参考のために、ZIPファイルを置いておきますので、ダウンロードして活用してください。

julia_hydraulic_r002/
├── Dockerfile
├── compose.yaml
├── Project.toml
├── README.md
├── VALIDATION.md
├── qiita_r002.md
│
├── case/
│   ├── single_reach_case.csv
│   ├── single_reach_backwater_source.json
│   ├── comparison_tolerance.csv
│   └── README.md
│
├── src/
│   ├── JuliaHydraulicR002.jl
│   ├── config.jl
│   ├── geometry.jl
│   ├── solver.jl
│   ├── io.jl
│   └── comparison.jl
│
├── fortran/
│   └── reference_solver.f90
│
├── reference/
│   ├── fortran_timeseries.csv
│   ├── fortran_summary.csv
│   ├── fortran_compiler.txt
│   └── README.md
│
├── examples/
│   ├── run_julia_case.jl
│   └── compare_fortran_julia.jl
│
├── test/
│   └── runtests.jl
│
├── tools/
│   └── independent_float32_check.py
│
├── scripts/
│   ├── preflight.py
│   ├── build_fortran.sh
│   ├── run_fortran.sh
│   ├── run_julia.sh
│   ├── compare.sh
│   ├── validate_reference.sh
│   ├── run_all.sh
│   └── verify.sh
│
├── build/
└── output/

reference/fortran_timeseries.csv は、配布時に実際にgfortranで計算して作った基準値です。

8. 共通入力CSV

条件は次のような単純な key,value CSVにしました。

key,value
duration_s,600
dt_s,2
save_interval_s,20
gravity_mps2,9.80665
min_depth_m,0.02
length_m,1000
n_segments,5
bed_upstream_m,1.0
bed_downstream_m,0.5
bottom_width_m,10.0
side_slope_hv,1.0
manning_n,0.03
initial_stage_m,3.0
upstream_flow_m3s,0.0
downstream_stage_start_m,3.0
downstream_stage_end_m,4.5

FortranとJuliaの両方がこのファイルを読みます。

将来、比較ケースを増やすときも、計算条件をプログラムに埋め込まなくて済みます。

9. 流量更新

Julia側の流量更新部分は次の形です。

function _update_face_q(
    q_old::Float32,
    stage_left::Float32,
    stage_right::Float32,
    depth_left::Float32,
    depth_right::Float32,
    distance::Float32,
    cfg::SingleReachConfig,
)
    hface = 0.5f0 * (depth_left + depth_right)
    area = trapezoid_area(hface, cfg.bottom_width_m, cfg.side_slope_hv)
    radius = hydraulic_radius(hface, cfg.bottom_width_m, cfg.side_slope_hv)

    slope_wse = (stage_right - stage_left) / distance

    numerator = q_old - cfg.dt_s * cfg.gravity_mps2 * area * slope_wse

    friction =
        cfg.dt_s * cfg.gravity_mps2 * cfg.manning_n * cfg.manning_n * abs(q_old) /
        (max(area, 1f-9) * max(radius, 1f-9)^(4f0 / 3f0))

    return numerator / (1f0 + friction)
end

摩擦項に Q|Q| を使う形なので、流量が負になっても同じ式で扱えます。

このケースでは、下流水位を上昇させるため、途中から負流量が現れます。

10. 連続式で水位を更新する

各セルでは、流入量と流出量から断面積を更新します。

離散化すると、

A(i) new = A(i) old + dt / dx × (Qin - Qout)

です。

Juliaでは次の部分です。

area_new = area_old + cfg.dt_s / dx * (q_new[i] - q_new[i + 1])

その後、台形断面積から水深を戻し、河床高を加えて水位にします。

断面積
  ↓
水深
  ↓
河床高 + 水深
  ↓
水位

11. 水収支

今回の比較では、水位と流量だけでなく水収支も保存します。

初期貯留量を V0、時刻 t の貯留量を V(t)、境界からの累積正味流入量を B(t) とすると、

Residual(t) = V(t) - V0 - B(t)

としました。

境界正味流入量は、流量の正方向を上流から下流として、

B += dt × (Qup - Qdown)

です。

下流から逆流してくる場合には Qdown < 0 になるため、河道内への流入量として正しく加算されます。

12. Fortran側を実際に動かす

Fortran側は REAL(real32) で実装しています。

ビルドは、

gfortran \
  -std=f2008 \
  -O2 \
  -Wall -Wextra \
  -fcheck=all -fbacktrace \
  fortran/reference_solver.f90 \
  -o build/fortran_reference

です。

今回の成果物作成環境では、

GNU Fortran (Debian 14.2.0-19) 14.2.0

で実際にビルドしています。

計算も600秒まで実行しました。

主要結果は次のとおりでした。

項目 Fortran REAL*4結果
最終上流側セル水位 4.50278378 m
最終下流側セル水位 4.49196148 m
最小流量 -61.1819839 m³/s
最終下流流量 -51.5203323 m³/s
初期貯留量 27582.5000 m³
最終貯留量 51456.2891 m³
累積境界正味流入量 23873.8066 m³
最終水収支残差 -0.017578125 m³
最大絶対水収支残差 0.0458984375 m³

REAL*4 なので、水収支が機械精度レベルのゼロになるわけではありません。

一方で、約2.39万m³の境界正味流入に対して、最大残差は約0.046m³でした。

この値はJulia側との比較でも確認します。

13. 全保存時刻を比較する

比較は最終時刻だけではありません。

保存間隔が20秒なので、600秒計算では31時刻あります。

比較対象は、

31時刻
  ×
  ├─ 水位 5セル
  ├─ 流量 6面
  ├─ 河道内貯留量
  ├─ 境界正味流入量
  └─ 水収支残差

です。

Julia側では compare_timeseries() が2つのCSVを読み、列ごとの最大絶対差を求めます。

passed, rows = compare_timeseries(
    "reference/fortran_timeseries.csv",
    "output/julia_timeseries.csv",
)

結果は、

output/comparison.csv

へ保存します。

14. 比較許容値

最初の比較許容値は次のようにしました。

分類 最大絶対差
水位 0.0001 m
流量 0.001 m³/s
貯留量・水収支 0.1 m³
時刻 0.0001 s

これは水理モデルそのものの最終的な品質基準ではありません。

r002 での目的は、同じFloat32/REAL*4計算を別言語へ移したときの回帰比較です。

実モデルの許容誤差については、対象モデルや観測精度を含めて別途決める必要があります。

15. 独立したFloat32再計算も入れる

今回の成果物を作っている環境にはJulia本体が入っていません。

そのため、Julia実行について「PASS」とはしていません。

ただし、Fortran側の式とJulia側へ転記した式に大きな食い違いがないかを見るため、NumPy Float32 だけで同じ離散式を独立実装した検査を行いました。

この検査はJuliaの代替ではありません。

結果は、Fortran出力に対して、

水位       最大差 約 5e-9 m
流量       最大差 約 5e-8 m3/s
貯留量     最大差 約 5e-5 m3
境界累積量 最大差 約 4.7e-5 m3

でした。

これで少なくとも、共有入力とREAL*4 Fortran基準値、離散式の転記については整合していることを確認しています。

Juliaについては、次のDocker検証を実行して初めて確定します。

16. Dockerで検証する

DockerfileはJulia 1.13.1を基準にし、その中へgfortranを追加します。

つまり、FortranとJuliaを同じコンテナ環境で実行できます。

20261013_f04.png

WSL2上では、プロジェクトディレクトリで、

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

./scripts/verify.sh

を実行します。

verify.sh は、

Docker build
    ↓
Julia / gfortran version確認
    ↓
Fortran基準値の再生成・byte比較
    ↓
Julia Pkg.test()
    ↓
Julia本計算
    ↓
Fortran / Julia 全時刻比較

まで実行します。

17. Fortran基準値も固定する

reference/fortran_timeseries.csv は、今回実際にgfortranで計算した結果です。

validate_reference.sh では、Fortranを再ビルド・再計算したあと、

cmp reference/fortran_timeseries.csv output/fortran_timeseries.csv
cmp reference/fortran_summary.csv output/fortran_summary.csv

で基準値とbyte比較します。

これにより、Julia側の比較以前にFortran基準値が変わっていないかを確認できます。

18. Juliaのテスト

test/runtests.jl では、次を確認します。

  • 台形断面の A -> h 逆変換
  • 共通入力CSVの読込
  • 下流水位境界の補間
  • 有限値で計算できること
  • 背水により負流量が発生すること
  • 上流側水位が初期値より上昇すること
  • 最大水収支残差が0.1m³未満であること
  • Fortran REAL*4基準CSVとの全時刻比較

最後の項目が r002 の中心です。

19. 今回確認できたことと、まだ確認していないこと

ここは分けておきます。

今回の作成環境で確認済み

  • Fortranソースのgfortranビルド
  • REAL(real32) による600秒計算
  • 31保存時刻の結果CSV生成
  • 水収支計算
  • Fortran基準CSVの再実行byte一致
  • 独立Float32実装との数値照合

この作成環境では未確認

  • Julia 1.13.1による実計算
  • Pkg.test()
  • FortranとJuliaの実比較PASS
  • Docker Compose build/run

未実行項目はPASS扱いにしません。

WSL2側で ./scripts/verify.sh を通した結果を、Julia版の最終確認とします。

20. 既存Fortran全体との関係

今回同梱したFortranコードは、r002 の単一河道比較だけに範囲を絞った参照実装です。

既存のNILIM2.0 Fortranには、河川以外にも、

  • 2D氾濫
  • 1D–2D交換
  • 下水道
  • 構造物
  • 横流入
  • 複数河道・Node–Reach

などがあります。

それらを一度にJuliaへ移すと、差が出たときの原因を追いにくくなります。

そのため、今回は、

r001
配列・差分計算
   ↓
r002
単一河道1D + Fortran比較
   ↓
r003
Node–Reach
   ↓
r004
横流入
   ↓
r005
降雨流出

という順番を維持します。

20261013_f05.png

21. まとめ

r002 では、ようやく水理計算に入りました。

今回作ったのは大きな河川モデルではありません。

1本の河道を5分割しただけです。

ただし、今後の移植で必要になる、

共通入力
   ↓
Fortran REAL*4
   ↓
基準結果

共通入力
   ↓
Julia Float32
   ↓
計算結果

2つを全時刻で比較

という流れを先に作ることができました。

新しい言語へ移すときは、コードが動くことよりも、これまでの計算と何が同じで、何が変わったのかを追えることの方が重要です。

次の r003 では、この単一河道を Node–Reach へ広げます。

合流・分岐へ進む前に、まずNodeとReachの保存量・流量の受け渡しを整理し、今回と同じようにFortran側と比較できる形を維持していきます。

参考

0
0
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
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?