この文書では、WSL2のUbuntu24.04上でLBM（格子ボルツマン法）流体解析ツールXLBを活用する方法を紹介します。

「XLB（Accelerated Lattice Boltzmann）」は、Autodeskが公開している、GPUによる高速計算と機械学習との連携を重視したPython製LBM（格子ボルツマン法）ライブラリです。2次元・3次元の流体解析に対応し、標準的なBGKモデルや高レイノルズ数流れに適したKBCモデルなどを利用できます。
最大の特徴は、NVIDIA Warp、Google JAX、Autodesk Neonなどの計算バックエンドを利用できる点です。いずれも、PythonコードからGPUを活用した高速なLBM計算を比較的容易に実装できます。

以下の【01:JAXcpu導入方法と動作確認】では、バックエンドにJAXを CPUで利用する環境を構築し、例題lid_driven_cavity_2d.pyの動作確認まで進めました。

さらに【03:JAXgpu導入方法と動作確認】では、バックエンドにJAXを GPUで利用する環境を構築し、例題lid_driven_cavity_2d.pyの動作確認まで進めました。

この【04】では、【03】の手順で構築したJAXのGPU活用の環境をもとに、例題 flow_past_sphere_3d.pyを分析してゆきます。
https://github.com/Autodesk/XLB/blob/main/examples/cfd/flow_past_sphere_3d.py

XLBのGitHub公開例題の中で、最初に理解する例題としては lid_driven_cavity_2d.py（2次元 Lid-Driven Cavity）があり【02】で分析しました。次に３次元の例題で、本記事【04】では、この例題flow_past_sphere_3d.pyの分析を行います。これはXLBの基本やスクリプトの設定を学ぶための例題として、最適です。

例題「3次元球体後部の流れ」の概要

この例題は、3次元の流体中に球を置き、その周囲を流れる流れをXLBで解析する例題です。例題スクリプトでは、３次元格子、D3Q19モデル、JAXバックエンド、FP32計算精度、BGK衝突モデルを用い、球表面にはHalfway Bounce Back、入口には速度分布、出口には流出境界条件を設定しています。さらに、球に作用する抗力・揚力を計算し、速度分布を可視化します。

例題の概念図を、以下に示します。

以下の説明では、スクリプトのコードそのものは示さず、対象とする例題 flow_past_sphere_3d.py が「何を計算しているのか」を、物理モデル → LBM → 境界条件 → 球 → 抗力・揚力 → 出力の順に、文章と数式で詳しく説明します。

1. この例題で何を解析するのか

この例題では、3次元の流路の中に球形の固体を置き、その周囲を流体が通過する流れを解析します。

概念的には、上図のような状態です。流体が球に衝突すると、球の前面では流れがよどみ、球の側面では流れが加速し、球の後方では流れが剥離して、後流や渦が形成されます。

このような物体周りの流れは、流体力学やCFDにおける代表的な問題です。

このXLBの例題では、単に速度場を求めるだけではなく、球に作用する抗力（Drag）と揚力（Lift）まで計算するところが重要です。例題スクリプトでは、球表面に対する MomentumTransfer を使って力を求めています。

2. 計算領域は3次元

この例題 flow_past_sphere_3d.pyでは、2次元の lid_driven_cavity_2d.py とは異なり、3次元のLBM格子を使用します。このスクリプトでは、解析空間の格子の基準分割数 (512, 128, 128) を、基本例題のため縮小するために分割しています。
スクリプトでは、512//2のように「//」を用いていますが、これは単純な除算ではなく、整数除算として余りを切り捨てます。例えば5//2=2となります。

格子サイズは、以下の設定です。
$N_x=512/2=256$　$N_y=128/2=64$　$N_z=128/2=64$

したがって、$256\times64\times64$ の3次元格子を計算します。基準分割数 $512\times128\times128$ に比べて、格子数が 1/8 になっています。

3. LBMでは何を計算しているのか

この例題では、通常のNavier-Stokes方程式を直接離散化するのではなく、格子ボルツマン法（LBM） によって流れを計算しています。
LBMでは、各格子点において、流体の密度や速度を直接時間発展させるのではなく、$f_i(\mathbf{x},t)$ という 粒子分布関数（distribution function） を計算します。ここで (i) は粒子の移動方向を表します。

この例題では、格子の速度モデルとして D3Q19 を使用しています。

4. D3Q19とは

D3Q19は、今回の例題 flow_past_sphere_3d.py では、以下のように、指定されています。

velocity_set = xlb.velocity_set.D3Q19(
    precision_policy=precision_policy
)

名前を分解すると、以下の意味です。

D3 → 3次元（3 Dimensions）
Q19 → 19個の離散速度方向（19 discrete velocities）

具体的には、下図のような定義になっています。

厳密には、中心方向を含めて19個の離散速度があり、格子点では19個の分布関数 $f_0,f_1,\ldots,f_{18}$ を扱います。D3Q19は、3次元LBM（格子ボルツマン法）で使われる代表的な速度モデルです。

5. BGKモデルによる時間発展

この例題では、衝突モデルとしてBGKモデルを使用しています。
LBMの基本的な時間発展は、

$ f_i(\mathbf{x}+\mathbf{e}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)-\omega\left[f_i(\mathbf{x},t)-f_i^{eq}(\mathbf{x},t)\right]$

で表されます。

ここで、$f_i$：現在の分布関数、$f_i^{eq}$：平衡分布関数、$e_i$：D3Q19の離散速度、$\omega$：緩和パラメータです。この例題では、$\omega=1.6$ と設定されています。
つまり、この例題ではBGKモデルによって分布関数を平衡状態へ緩和させながら、流体の時間発展を計算しています。
基本式などの、概要は下図のとおりです。

6. 分布関数から密度と速度を求める

LBMで直接計算しているのは (f_i) ですが、私たちが知りたいのは、「密度 $(\rho)$・速度 $(\mathbf{u})$」です。
これらは分布関数から、$\rho=\sum_i f_i$ として密度を求め、$\rho\mathbf{u}=\sum_i f_i\mathbf{e}_i$ から速度を求めます。

例題スクリプトでは、この変換をXLBの Macroscopic という機能によって行っています。

7. 入口には放物線状の速度分布を与える

この例題の特徴の一つが、入口の速度条件です。入口では単純に、$u_x=u_{\max}$ とするのではなく、断面内で放物線状に変化する速度分布を与えています。
例題スクリプトでは、$u_x=u_{\max}\max(0,1-r^2)$ という形になっています。

ここで、$r^2=\left(\frac{2(y-y_c)}{H_y}\right)^2+\left(\frac{2(z-z_c)}{H_z}\right)^2$ です。
また、$u_y=0$、$u_z=0$ としています。これは円形または矩形流路におけるPoiseuille型の速度分布を模したものです。

8. 流路の壁は固定壁

計算領域の上面・下面・前面・背面は、流体が通過できない固定壁として扱われます。例題スクリプトでは、これらの壁にFullway Bounce Backを使用しています。
LBMにおけるBounce Backは、簡単に言えば、「壁に向かって進んできた分布を反対方向へ跳ね返す」 という考え方です。
これによって壁面におけるno-slip条件、$\mathbf{u}=0$ を表現します。

9. 球は固体として扱われる

この例題の中心となるのが球です。球の半径は、$r_{\mathrm{sphere}}=\frac{N_y}{12}$ として設定されています。
そして、格子上の各点について、$(x-x_c)^2+(y-y_c)^2+(z-z_c)^2<r^2$ を満たす点を球の内部として抽出しています。つまり、連続的なCAD形状を作っているのではなく、LBM格子上で球形の固体領域を構成しているわけです。

球表面には、Halfway Bounce Backを使用しています。したがって球も、$\mathbf{u}=0$ の固体壁として扱われます。

10. 出口では流体を外へ流す

流路の右側は出口です。入口とは異なり、出口では速度を直接指定していません。例題スクリプトでは、Extrapolation Outflowという境界条件を使用しています。これは、計算領域の内部の情報を利用して、流体が自然に計算領域の外へ流れ出るようにする境界条件です。

11. 球の後ろには後流が形成される

球の前方では、流体が球を避けて左右・上下に分かれます。球の後方では、流れが再び合流しますが、流れはすぐには一様な状態に戻りません。その結果、球の後ろには、 後流（wake） が形成されます。
流れの条件によっては、後流中に周期的な渦が形成されます。概念図のような現象になります。この分析では、条件を変えることで、単純な層流や周期的な渦を確認します。

この例題では速度場を計算しているので、このような球後方の流れ構造を可視化できます。

12. 球に作用する抗力と揚力を計算する

この例題の大きな特徴がここです。単に、「流体の速度を計算して終わり」ではありません。
球表面における運動量の変化から、球に作用する力を求めています。XLBでは MomentumTransfer を使用しています。

球に作用する力を、$\mathbf{F}=(F_x,F_y,F_z)$ とします。例題スクリプトでは、流れに平行な$F_x$ を抗力方向の力、流れに直交する$F_z$ を揚力方向の力として取り出しています。

13. 抗力係数

さらに、抗力係数$C_D$ を、単なる力 $F_x$ ではなく、無次元化した抗力係数を計算しています。
例題スクリプトでは、$C_D=\frac{2F_D}{u_{\max}^2 A}$ という形で計算しています。
ここで、$F_D$：抗力、$u_{\max}$：代表速度、$A$：球の断面積、です。

球の投影断面積は、$A=\pi r^2$ として計算されています。したがって、$C_D=\frac{2F_D}{u_{\max}^2\pi r^2}$ となります。これは流体力学で一般的に使われる抗力係数です。

14. 揚力係数

同様に、揚力についても、揚力係数 $C_L=\frac{2F_L}{u_{\max}^2 A}$ として無次元化しています。つまり、この例題では最終的に、$C_D$と、$C_L$という、球周りの流れを評価する重要な指標を得ることができます。

15. 時間発展

計算は1回だけ行うのではありません。この例題では、$N_{\mathrm{step}}=10000$ステップまでLBMの時間発展を行います。各ステップで、$f_i^n \rightarrow f_i^{n+1}$ という更新を行います。

概念的には、以下を10,000回繰り返します。

初期状態 → Collision → Streaming → Boundary
→ 次の時刻 → Collision → Streaming → Boundary  → ...

16. 1000ステップごとに後処理する

解析結果のなかで、すべてのステップで画像を保存するのではなく、1000ステップごとに後処理を行います。
後処理では、以下の処理を行います。

1. 球に作用する力を計算 → 2. CD を計算 → 3. CL を計算 → 4. 密度を計算
 → 5. 速度を計算 → 6. 速度の大きさを計算 → 7. 画像を保存

17. 速度の大きさを計算する

速度は、$\mathbf{u}=(u_x,u_y,u_z)$ として得られます。その大きさは、$|\mathbf{u}|=\sqrt{u_x^2+u_y^2+u_z^2}$ です。例題スクリプトでも、この式によって u_magnitude を計算しています。

18. この例題で最終的に得られるもの

このプログラムの出力を整理すると、大きく3種類あります。
① 流れ場：$u_x,\quad u_y,\quad u_z$ および $|\mathbf{u}|$ です。
② 密度場：$\rho$ も計算されています。
③ 球に作用する力：$F_D,\quad F_L$から、$C_D,\quad C_L$を計算します。
したがって、以下のような構成になっています。

              flow_past_sphere_3d
        ┌──────────────┼──────────────┐
     流れ場           密度         流体力
   ux, uy, uz                     Drag / Lift
       |u|                         CD / CL
      画像

19. この例題の本質

この flow_past_sphere_3d.py は、XLBを理解するうえでかなり良い発展例題です。
lid_driven_cavity_2d.py では、「LBMで流れを計算する」 ことが中心でした。
一方、この例題では、「LBMで3次元の物体周り流れを計算し、物体に作用する力まで評価する」 ところまで進んでいます。

全体を数式でまとめると、$f_i \rightarrow \text{BGK} \rightarrow \text{Streaming} \rightarrow \text{Boundary} \rightarrow f_i^{n+1}$ を繰り返し、$f_i \rightarrow (\rho,\mathbf{u})$ を求め、さらに球表面での運動量交換から、$\mathbf{F} \rightarrow F_D,F_L \rightarrow C_D,C_L$ を求める、という流れです。

つまり、この例題を理解すると、「XLBはLBMを計算するだけのライブラリではなく、境界条件、物体周りの流れ、物体に作用する流体力まで扱えるCFD基盤である」ということが見えてきます。
なお、例題スクリプトでは、速度場については save_image() による中央 (y) 断面の速度大きさの画像保存が実装されています。

例題スクリプトの具体的な解説

例題スクリプトの処理内容・構成は基本的に変更せず、学習しやすいように ①〜⑫の12ブロックに分け、それぞれに「何をしているか」「なぜ必要か」「物理的な意味」を詳しくコメントとして挿入した1本のスクリプトにまとめます。

① 計算条件を決める
② XLBを初期化
③ 256 × 64 × 64 のGridを作る
④ 入口・出口・壁・球を定義
⑤ 入口速度を設定
⑥ 境界条件を設定
⑦ BGK + D3Q19 のStepperを作る
⑧ 密度・速度を計算する仕組みを作る
⑨ 球に働く力を計算する仕組みを作る
⑩ 時間発展：Collision（BGK）＋ Boundary Condition ＋ Streaming
⑪ CD・CL・速度分布を計算
⑫ 画像を保存

import xlb
from xlb.compute_backend import ComputeBackend
from xlb.precision_policy import PrecisionPolicy
from xlb.grid import grid_factory
from xlb.operator.stepper import IncompressibleNavierStokesStepper
from xlb.operator.boundary_condition import (
    FullwayBounceBackBC,
    HalfwayBounceBackBC,
    RegularizedBC,
    ExtrapolationOutflowBC,
)
from xlb.operator.macroscopic import Macroscopic
from xlb.utils import save_image
import warp as wp
import numpy as np
import jax.numpy as jnp
import time


# ============================================================================
# ① ライブラリの読み込み
# ============================================================================
#
# このブロックでは、XLBによるLBM計算に必要な機能を読み込みます。
# XLB：
#   格子ボルツマン法（LBM）の計算を行うためのライブラリ
# ComputeBackend：
#   XLBでどの計算バックエンドを使用するか指定する
#   今回は JAX を使用
# PrecisionPolicy：
#   浮動小数点数の精度を指定する
# grid_factory：
#   LBM計算を行う3次元格子（Grid）を作成する
# IncompressibleNavierStokesStepper：
#   非圧縮Navier-Stokes方程式に対応したLBMの時間発展を行う
# Boundary Condition：
#   入口、壁、球表面、出口などの境界条件を設定する
# Macroscopic：
#   LBMの分布関数から密度 rho と速度 u を計算する
# save_image：
#   計算結果を画像として保存する
# Warp：
#   XLB内部で使用される可能性があるため読み込んでいる
# NumPy / JAX：
#   格子インデックスや速度分布などの配列処理に使用する
# time：
#   計算時間を測定する
# ============================================================================


# ============================================================================
# ② シミュレーション条件の設定
# ============================================================================
#
# omega：
#   BGK衝突モデルの緩和パラメータです。
#   BGKモデルでは、分布関数を平衡分布へ緩和させます。
#   omega = 1.6 が今回の設定です。
# --------------------------------------------------------------------------
#
# grid_shape：
#   計算領域の格子数を指定します。
#       (512 // 2, 128 // 2, 128 // 2)
#   Pythonの // は整数除算です。
#       512 // 2 = 256　128 // 2 = 64
#   したがって実際の格子サイズは
#       (256, 64, 64)
#   となります。
#   x方向：256　y方向： 64　z方向： 64
#   したがって総格子点数は
#       256 × 64 × 64 = 1,048,576  です。
# --------------------------------------------------------------------------
#
# compute_backend：
#   XLBの計算バックエンドとしてJAXを指定しています。
#   したがって、この例ではXLBの主要な計算を
#       XLB → JAX
#   という構成で実行します。
# --------------------------------------------------------------------------
#
# precision_policy：
#   FP32FP32を指定しています。
#   今回は単精度（32bit浮動小数点）を中心とした計算設定です。
# --------------------------------------------------------------------------
#
# velocity_set：
#   D3Q19を指定しています。
#   D3Q19とは、D3 = 3次元　Q19 = 19個の離散速度方向 という意味です。
#   つまり、各格子点に19方向の粒子分布関数
#       f_0, f_1, ..., f_18
#   を持たせてLBM計算を行います。
# --------------------------------------------------------------------------
#
# u_max：
#   入口での最大流速です。
#       u_max = 0.04
#   ここではLBMの格子単位で指定されています。
# --------------------------------------------------------------------------
#
# num_steps：
#   時間発展を10,000ステップ行います。
# --------------------------------------------------------------------------
#
# post_process_interval：
#   1,000ステップごとに計算結果を取り出して、
#   抗力・揚力・速度分布などを計算します。
# ============================================================================

omega = 1.6

grid_shape = (
    512 // 2,
    128 // 2,
    128 // 2,
)

compute_backend = ComputeBackend.JAX

precision_policy = PrecisionPolicy.FP32FP32

velocity_set = xlb.velocity_set.D3Q19(
    precision_policy=precision_policy,
    compute_backend=compute_backend,
)

u_max = 0.04

num_steps = 10000

post_process_interval = 1000


# ============================================================================
# ③ XLBの初期化
# ============================================================================
#
# XLB全体に対して、
#   ・使用する速度モデル
#   ・標準バックエンド
#   ・標準精度
# を登録します。
# velocity_set：
#   D3Q19
# default_backend：
#   JAX
# default_precision_policy：
#   FP32FP32
# これ以降のXLBの各オブジェクトが、
# この設定を基準として動作します。
# ============================================================================

xlb.init(
    velocity_set=velocity_set,
    default_backend=compute_backend,
    default_precision_policy=precision_policy,
)


# ============================================================================
# ④ 計算格子（Grid）の作成
# ============================================================================
#
# grid_factory()を使って、実際のLBM計算用のGridを作ります。
# 今回は、grid_shape = (256, 64, 64)
# なので、
#       x方向 256
#       y方向  64
#       z方向  64
# また、compute_backend=ComputeBackend.JAX
# なので、JAXバックエンド用のGridとして作成されます。
# このgridオブジェクトを利用して、
#   ・入口
#   ・出口
#   ・壁
#   ・球
# などの格子位置を後から指定します。
# ============================================================================

grid = grid_factory(
    grid_shape,
    compute_backend=compute_backend,
)


# ============================================================================
# ⑤ 計算領域と球の形状を定義
# ============================================================================
#
# ここでは、計算領域の境界と、流路中央に置く球の位置を定義します。
#
# --------------------------------------------------------------------------
# 5-1. 計算領域の境界
# --------------------------------------------------------------------------
#
# bounding_box_indices()から計算領域の境界格子を取得します。
# boxには、
#   left
#   right
#   bottom
#   top
#   front
#   back
# などの境界情報が入ります。
# ============================================================================

box = grid.bounding_box_indices()

# エッジ部分を除いた境界を取得
box_no_edge = grid.bounding_box_indices(remove_edges=True)

# 左側を入口にする
inlet = box_no_edge["left"]

# 右側を出口にする
outlet = box_no_edge["right"]


# --------------------------------------------------------------------------
# 5-2. 流路の壁を定義
# --------------------------------------------------------------------------
#
# 流路の上下・前後の4面を壁として設定します。
#   bottom：下壁
#   top   ：上壁
#   front ：手前側
#   back  ：奥側
# x方向のleft/rightは入口・出口なので、ここでは壁にしません。
# ============================================================================

walls = [
    box["bottom"][i]
    + box["top"][i]
    + box["front"][i]
    + box["back"][i]
    for i in range(velocity_set.d)
]

# 重複する格子インデックスを削除
walls = np.unique(
    np.array(walls),
    axis=-1
).tolist()


# --------------------------------------------------------------------------
# 5-3. 球のサイズと位置を定義
# --------------------------------------------------------------------------
#
# 球の半径を流路のy方向サイズから決めています。
#       sphere_radius = grid_shape[1] // 12
# 今回、grid_shape[1] = 64 なので、
#       sphere_radius = 64 // 12 = 5
# となります。
# ※Pythonの整数除算なので、厳密には5です。
# --------------------------------------------------------------------------
#
# 球の中心位置：
#       x = grid_shape[0] // 6
#       y = grid_shape[1] // 2
#       z = grid_shape[2] // 2
# つまり、
#   x方向：入口から少し入った位置
#   y方向：流路中央
#   z方向：流路中央
# に球を配置します。
# ============================================================================

sphere_radius = grid_shape[1] // 12

x = np.arange(grid_shape[0])
y = np.arange(grid_shape[1])
z = np.arange(grid_shape[2])

X, Y, Z = np.meshgrid(
    x,
    y,
    z,
    indexing="ij"
)

# 球の内部に存在する格子点を検索
indices = np.where(
    (X - grid_shape[0] // 6) ** 2
    + (Y - grid_shape[1] // 2) ** 2
    + (Z - grid_shape[2] // 2) ** 2
    < sphere_radius**2
)

# D3Q19の各方向について球表面の境界位置を作成
sphere = [
    tuple(indices[i].tolist())
    for i in range(velocity_set.d)
]


# ============================================================================
# ⑥ 入口速度分布（Poiseuille流）を定義
# ============================================================================
#
# ここでは入口から流入する速度を設定します。
# 入口速度は一様速度ではなく、断面方向に放物線状に変化する
# Poiseuille流の速度分布を使用します。
# 中央で速度が最大となり、壁に近づくほど速度が小さくなります。
# 概念的には、
#
#       壁
#       ↓
#       0
#       │      ／＼
#       │    ／    ＼
#       │  ／        ＼
#       │／            ＼
#       0────────────────→
#             入口速度
# のような分布です。
# 今回は3次元流れなので、y-z断面上で速度分布を作ります。
# ============================================================================

def bc_profile():

    # y方向とz方向の流路高さ
    H_y = float(grid_shape[1] - 1)
    H_z = float(grid_shape[2] - 1)


    # ----------------------------------------------------------------------
    # 6-1. JAXバックエンドの場合
    # ----------------------------------------------------------------------
    #
    # JAXを使用する場合は、JAX配列として入口速度分布を作ります。
    # ----------------------------------------------------------------------

    if compute_backend == ComputeBackend.JAX:

        def bc_profile_jax():

            y = jnp.arange(grid_shape[1])
            z = jnp.arange(grid_shape[2])

            Y, Z = jnp.meshgrid(
                y,
                z,
                indexing="ij"
            )

            # 断面中心からの距離
            y_center = Y - (H_y / 2.0)
            z_center = Z - (H_z / 2.0)

            # 中心からの距離を正規化
            r_squared = (
                (2.0 * y_center / H_y) ** 2.0
                + (2.0 * z_center / H_z) ** 2.0
            )

            # x方向速度
            #
            #       u_x = u_max (1 - r²)
            #
            # r²が1を超えた場合は速度を0にします。
            u_x = u_max * jnp.maximum(
                0.0,
                1.0 - r_squared
            )

            # y方向速度は0
            u_y = jnp.zeros_like(u_x)

            # z方向速度も0
            u_z = jnp.zeros_like(u_x)

            # 3方向の速度をまとめる
            return jnp.stack([
                u_x,
                u_y,
                u_z
            ])

        return bc_profile_jax


    # ----------------------------------------------------------------------
    # 6-2. Warpバックエンドの場合
    # ----------------------------------------------------------------------
    #
    # 今回のcompute_backendはJAXなので、この部分は実際には使用されません。
    # ただし、XLBの例題としてJAX以外のバックエンドにも対応できるよう、
    # Warp版の入口速度分布も定義されています。
    # ----------------------------------------------------------------------

    else:

        wp_dtype = precision_policy.compute_precision.wp_dtype

        H_y = wp_dtype(grid_shape[1] - 1)
        H_z = wp_dtype(grid_shape[2] - 1)

        two = wp_dtype(2.0)


        @wp.func
        def bc_profile_warp(index: wp.vec3i):

            # y,z方向の格子インデックス
            y = wp_dtype(index[1])
            z = wp_dtype(index[2])

            # 断面中心からの距離
            y_center = y - (H_y / two)
            z_center = z - (H_z / two)

            # 正規化された距離
            r_squared = (
                (two * y_center / H_y) ** two
                + (two * z_center / H_z) ** two
            )

            # 放物線状速度分布
            return wp.vec(
                wp_dtype(u_max)
                * wp.max(
                    wp_dtype(0.0),
                    wp_dtype(1.0) - r_squared
                ),
                length=1
            )

        return bc_profile_warp


# ============================================================================
# ⑦ 境界条件を設定
# ============================================================================
#
# ここでは流路の各境界に物理的な条件を与えます。
# 今回は4種類の境界条件を使用しています。
#   ① 入口       ：RegularizedBC
#   ② 流路壁     ：FullwayBounceBackBC
#   ③ 出口       ：ExtrapolationOutflowBC
#   ④ 球表面     ：HalfwayBounceBackBC
# --------------------------------------------------------------------------
# 入口
# --------------------------------------------------------------------------
#
# velocity境界条件として、先ほど作ったPoiseuille速度分布を指定します。
# ============================================================================

bc_left = RegularizedBC(
    "velocity",
    profile=bc_profile(),
    indices=inlet
)


# --------------------------------------------------------------------------
# 入口速度を単純に一定値にすることも可能
# 以下は今回使用していません。
# bc_left = RegularizedBC(
#     "velocity",
#     prescribed_value=(u_max, 0.0, 0.0),
#     indices=inlet
# )
# --------------------------------------------------------------------------


# --------------------------------------------------------------------------
# 流路壁
# --------------------------------------------------------------------------
#
# Fullway Bounce Back境界条件です。
# 壁に向かって進んだ分布関数を反射させることで、
# 壁面で流体が通過しない条件を表現します。
# ============================================================================

bc_walls = FullwayBounceBackBC(
    indices=walls
)


# --------------------------------------------------------------------------
# 出口
# --------------------------------------------------------------------------
#
# ExtrapolationOutflowBCを使います。
# 流体が計算領域の右側から自然に流出できるようにします。
# ============================================================================

bc_outlet = ExtrapolationOutflowBC(
    indices=outlet
)


# --------------------------------------------------------------------------
# 球表面
# --------------------------------------------------------------------------
#
# 球の表面にはHalfway Bounce Backを使用します。
# 流体粒子が球に衝突した際に反射させることで、
# 固体球の存在をLBM格子上で表現します。
# これによって、
#       流体 → ● 球 → 後流
# という球まわりの流れを計算できます。
# ============================================================================

bc_sphere = HalfwayBounceBackBC(
    indices=sphere
)


# すべての境界条件をリストにまとめる
boundary_conditions = [
    bc_walls,
    bc_left,
    bc_outlet,
    bc_sphere
]


# ============================================================================
# ⑧ LBMの時間発展計算を準備
# ============================================================================
#
# IncompressibleNavierStokesStepperは、
# LBMによる非圧縮Navier-Stokes流れの時間発展を担当します。
# collision_type="BGK"
# により、衝突モデルとしてBGKモデルを使用します。
# BGKモデルでは、分布関数を局所平衡分布へ緩和させます。
# その後、Streaming（並進）によって各方向へ分布関数を移動させます。
# prepare_fields()では、時間発展計算に必要な値を準備します。
# f_0, f_1：
#   時間発展で交互に使用する2つの分布関数バッファ
# bc_mask：
#   境界条件を適用するための情報
# missing_mask：
#   境界などによって不足する分布関数の情報
# ============================================================================

stepper = IncompressibleNavierStokesStepper(
    grid=grid,
    boundary_conditions=boundary_conditions,
    collision_type="BGK",
)

f_0, f_1, bc_mask, missing_mask = stepper.prepare_fields()


# ============================================================================
# ⑨ 巨視的物理量（密度・速度）の計算を準備
# ============================================================================
#
# LBMでは、基本的に分布関数 f_i を計算します。
# しかし、可視化や流体力学的な解析では、
#       密度 rho
#       速度 u
# が必要です。
# そこでMacroscopicを使用して、
#       rho = Σ f_i
#       rho u = Σ f_i e_i
# から密度と速度を計算します。
# 今回はD3Q19とJAXを使って巨視的物理量を計算します。
# ============================================================================

macro = Macroscopic(
    compute_backend=ComputeBackend.JAX,
    precision_policy=precision_policy,
    velocity_set=xlb.velocity_set.D3Q19(
        precision_policy=precision_policy,
        compute_backend=ComputeBackend.JAX
    ),
)


# --------------------------------------------------------------------------
# XLBの配列をJAX配列へ変換するための変換器
# --------------------------------------------------------------------------
#
# 後処理ではJAXの計算機能を利用するため、
# 必要に応じて分布関数をJAX配列へ変換します。
# ============================================================================

to_jax = xlb.utils.ToJAX(
    "populations",
    velocity_set.q,
    grid_shape
)


# ============================================================================
# ⑩ 球に作用する力（抗力・揚力）の計算を準備
# ============================================================================
#
# 球に流体からどれだけの力が作用しているかを計算します。
# MomentumTransfer：
#   境界との運動量交換から力を求める方法です。
# 今回は球表面の境界条件bc_sphereに対して
# MomentumTransferを適用します。
# これによって球に作用する力を計算できます。
# boundary_force[0]
#   → x方向の力
#   → 流れ方向の力なので「抗力 Drag」
# boundary_force[2]
#   → z方向の力
#   → 今回はこれを「揚力 Lift」として扱います。
# ============================================================================

from xlb.operator.force.momentum_transfer import MomentumTransfer

momentum_transfer = MomentumTransfer(
    bc_sphere,
    compute_backend=compute_backend
)


# --------------------------------------------------------------------------
# 球の正面投影面積
# --------------------------------------------------------------------------
#
# 抗力係数・揚力係数を計算するために、
# 球の流れ方向から見た断面積
#       A = π r²
# を計算します。
# ============================================================================

sphere_cross_section = np.pi * sphere_radius**2


# ============================================================================
# ⑪ 後処理（Post Processing）
# ============================================================================
#
# ここでは、一定時間ごとに計算結果を取り出して、
#   ・抗力
#   ・揚力
#   ・抗力係数 CD
#   ・揚力係数 CL
#   ・密度
#   ・速度
#   ・速度の大きさ
# を計算します。
#
# また、速度分布を画像として保存します。
# ============================================================================

def post_process(step, f_0, f_1):

    # ----------------------------------------------------------------------
    # Warp側の処理が完了するまで同期
    # ----------------------------------------------------------------------

    wp.synchronize()


    # ----------------------------------------------------------------------
    # 11-1. 球に作用する力を計算
    # ----------------------------------------------------------------------
    #
    # Momentum Transfer法によって球に作用する力を計算します。
    # ======================================================================

    boundary_force = momentum_transfer(
        f_0,
        f_1,
        bc_mask,
        missing_mask
    )


    # x方向の力を抗力とする
    drag = boundary_force[0]

    # z方向の力を揚力とする
    lift = boundary_force[2]


    # ----------------------------------------------------------------------
    # 11-2. 抗力係数 CD と揚力係数 CLを計算
    # ----------------------------------------------------------------------
    #
    # 抗力係数：
    #       CD = 2 FD / (rho U² A)
    # という形が一般的です。
    # このプログラムでは、正規化されたLBM量として
    #       CD = 2 * drag / (u_max² A)
    # として計算しています。
    # 同様に、
    #       CL = 2 * lift / (u_max² A)
    # としています。
    # ======================================================================

    cd = (
        2.0
        * drag
        / (u_max**2 * sphere_cross_section)
    )

    cl = (
        2.0
        * lift
        / (u_max**2 * sphere_cross_section)
    )

    print(
        f"CD={cd}, CL={cl}"
    )


    # ----------------------------------------------------------------------
    # 11-3. JAX配列へ変換
    # ----------------------------------------------------------------------
    #
    # f_0がJAX配列でない場合は、JAX配列へ変換します。
    # その後のMacroscopic計算をJAXで行うためです。
    # ======================================================================

    if not isinstance(f_0, jnp.ndarray):

        f_0 = to_jax(f_0)

        wp.synchronize()


    # ----------------------------------------------------------------------
    # 11-4. 密度 rho と速度 u を計算
    # ----------------------------------------------------------------------
    #
    # 分布関数f_iから、
    #       rho
    #       u = (u_x, u_y, u_z)
    # を計算します。
    # ======================================================================

    rho, u = macro(f_0)


    # ----------------------------------------------------------------------
    # 11-5. 境界セルを除去
    # ----------------------------------------------------------------------
    #
    # 可視化・後処理では、外側の境界セルを除きます。
    # ======================================================================

    u = u[
        :,
        1:-1,
        1:-1,
        1:-1
    ]

    rho = rho[
        :,
        1:-1,
        1:-1,
        1:-1
    ]


    # ----------------------------------------------------------------------
    # 11-6. 速度の大きさを計算
    # ----------------------------------------------------------------------
    #
    # 速度ベクトル
    #       u = (u_x, u_y, u_z)
    # の大きさ
    #       |u|
    #       = sqrt(u_x² + u_y² + u_z²)
    # を計算します。
    # 速度分布をカラー画像として表示するときに使用します。
    # ======================================================================

    u_magnitude = jnp.sqrt(
        u[0] ** 2
        + u[1] ** 2
        + u[2] ** 2
    )


    # ----------------------------------------------------------------------
    # 11-7. 後処理用の物理量をまとめる
    # ----------------------------------------------------------------------

    fields = {
        "u_magnitude": u_magnitude,
        "u_x": u[0],
        "u_y": u[1],
        "u_z": u[2],
        "rho": rho[0],
    }


    # ----------------------------------------------------------------------
    # 11-8. y方向中央断面の速度分布を画像保存
    # ----------------------------------------------------------------------
    #
    # y方向の中央
    #       y = grid_shape[1] // 2
    # の断面を取り出します。
    # つまり、
    #       y = 64 // 2 = 32
    # の断面です。
    # ここで保存されるのは、
    # 球を通過する流れの中央断面における速度の大きさです。
    #
    # 球の後ろには速度の低い領域や後流が形成されるため、
    # 球まわりの流れの様子を画像で確認できます。
    # ======================================================================

    save_image(
        fields["u_magnitude"][
            :,
            grid_shape[1] // 2,
            :
        ],
        timestep=step
    )

    print(
        f"Post-processed step {step}: "
        f"Saved u_magnitude slice at "
        f"y={grid_shape[1] // 2}"
    )


# ============================================================================
# ⑫ LBM時間発展ループ
# ============================================================================
#
# ここが実際の流体計算を行う中心部分です。
# LBMでは、時間ステップごとに
#       ① 衝突（Collision）
#       ② 境界条件
#       ③ 並進（Streaming）
# などを含むLBMの時間発展を行います。
# XLBでは、それらをstepper()がまとめて処理します。
# ============================================================================

start_time = time.time()


for step in range(num_steps):


    # ----------------------------------------------------------------------
    # 12-1. LBMを1ステップ進める
    # ----------------------------------------------------------------------
    #
    # stepper()によって、
    #       BGK衝突
    #       ↓
    #       境界条件
    #       ↓
    #       Streaming
    # などを含むLBMの時間発展が実行されます。
    # omega：
    #   BGKモデルの緩和パラメータ
    # step：
    #   現在の時間ステップ
    # ----------------------------------------------------------------------

    f_0, f_1 = stepper(
        f_0,
        f_1,
        bc_mask,
        missing_mask,
        omega,
        step
    )


    # ----------------------------------------------------------------------
    # 12-2. バッファを入れ替える
    # ----------------------------------------------------------------------
    #
    # f_0とf_1を交換します。
    #
    # LBMの時間発展では、現在の分布関数から次の時間の分布関数を
    # 計算するため、2つのバッファを交互に使用します。
    # イメージ：
    #       f_0 → 計算 → f_1
    #       次のステップ
    #       f_1 → 計算 → f_0
    # という「ダブルバッファ」の構造です。
    # ----------------------------------------------------------------------

    f_0, f_1 = f_1, f_0


    # ----------------------------------------------------------------------
    # 12-3. 一定ステップごとに後処理
    # ----------------------------------------------------------------------
    #
    # 1,000ステップごと、または最後のステップで
    # post_process()を実行します。
    # ここで、
    #       CD
    #       CL
    #       速度分布
    #       密度
    # などを計算・保存します。
    # ----------------------------------------------------------------------

    if (
        step % post_process_interval == 0
        or step == num_steps - 1
    ):

        post_process(
            step,
            f_0,
            f_1
        )


        # ------------------------------------------------------------------
        # 12-4. 計算時間を表示
        # ------------------------------------------------------------------

        end_time = time.time()

        elapsed = end_time - start_time

        print(
            f"Completed step {step}. "
            f"Time elapsed for "
            f"{post_process_interval} steps: "
            f"{elapsed:.6f} seconds."
        )

        # 次の計測区間の開始時刻を更新
        start_time = time.time()


# ============================================================================
# プログラム終了
# ============================================================================
#
# このプログラムによって、
#       3次元流路
#          ↓
#       入口から流入
#          ↓
#       球に衝突
#          ↓
#       球表面でBounce Back
#          ↓
#       球の後方に流れが形成
#          ↓
#       出口から流出
# という3次元球まわり流れのLBMシミュレーションを行います。
# また、一定時間ごとに
#       CD（抗力係数）
#       CL（揚力係数）
#       速度分布画像
# を出力します。
# この例ではcollision_type="BGK"としているため、
# XLBの時間発展にBGK衝突モデルが使われています。
# ============================================================================

例題スクリプトの実行の方法と出力

以上の解説を理解して、例題flow_past_sphere_3d.pyを実行してみます。
作業用ディレクトリXLBgpuWorkのなかに、例題ディレクトリFLOW3Dを作り、ここにリポジトリの例題ファイルflow_past_sphere_3d.pyをコピーします。

(XLBjaxGPU) dalab@dalab2312nb01:~$ cd ~/XLBjaxgpuWork/
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork$ ls
LDC2D  XLB-git
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork$ mkdir FLOW3D
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork$ ls
FLOW3D  LDC2D  XLB-git
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork$ cd FLOW3D/
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ cp ~/XLBjaxgpuWork/XLB-git/examples/cfd/flow_past_sphere_3d.py .
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ ls
flow_past_sphere_3d.py

実行時間を計測するために、コマンドtimeを用いて、pythonにより実行します。

(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ time python flow_past_sphere_3d.py
Single-GPU support is available: 1 GPU detected.
registered bc RegularizedBC_8601947340207 with id 1
registered bc FullwayBounceBackBC_8601947340222 with id 2
registered bc ExtrapolationOutflowBC_8601931827006 with id 3
registered bc HalfwayBounceBackBC_8601931819604 with id 4
WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC!
Warp 1.17.0 initialized:
   CUDA Toolkit 12.9, Driver 13.2
   Devices:
     "cpu"      : "x86_64"
     "cuda:0"   : "NVIDIA GeForce RTX 3070 Ti Laptop GPU" (8 GiB, sm_86, mempool enabled)
   Kernel cache:
     /home/dalab/.cache/warp/1.17.0
CD=0.0, CL=0.0
Post-processed step 0: Saved u_magnitude slice at y=32
Completed step 0. Time elapsed for 1000 steps: 1.243107 seconds.
CD=3.488799571990967, CL=0.0
：
Post-processed step 9999: Saved u_magnitude slice at y=32
Completed step 9999. Time elapsed for 1000 steps: 10.956343 seconds.

real    1m58.211s
user    0m22.224s
sys     0m6.927s

解析ログより、以下が確認できます。

「Single-GPU support is available」より１つのGPUを利用しています
4種類の境界条件のid設定と、no-slipの扱いが指摘されています
GPUは「"cuda:0":"NVIDIA GeForce RTX 3070 Ti Laptop GPU" (8 GiB, sm_86, mempool enabled)」です
最後のtimeコマンドの出力を見ると、処理時間は全体で118秒になりました
最後の段階で「抗力係数：CD=4.17169713973999, 揚力係数：CL=0.006799831986427307」です

下図のように、GPUが動作していることが分かります。

解析結果として、以下のような１１ステップのpng画像ファイルが作られます。これはY方向中央断面の速度分布です。
スクリプトの定義に従い、0000から9999までの1000ステップ間隔の出力になります。

(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ ls
flow_past_sphere_3d.py        flow_past_sphere_3d_3000.png  flow_past_sphere_3d_7000.png
flow_past_sphere_3d_0000.png  flow_past_sphere_3d_4000.png  flow_past_sphere_3d_8000.png
flow_past_sphere_3d_1000.png  flow_past_sphere_3d_5000.png  flow_past_sphere_3d_9000.png
flow_past_sphere_3d_2000.png  flow_past_sphere_3d_6000.png  flow_past_sphere_3d_9999.png

画像ファイルは、以下の内容です。このスクリプトの設定を分析すると、レイノルズ数は約１０なので、層流として球周りの整然としたながれになっているので、どのステップでも同じような状態です。

計算環境はGPU利用として構築しており、環境変数でも「JAX_PLATFORMS=cuda」として設定しています。
この設定を一時的に変更することで、以下のように、CPUによる解析も行うことができます。

(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ JAX_PLATFORMS=cpu time python flow_past_sphe
re_3d.py
No GPU support is available; CPU fallback will be used.
：
5886.01user 1481.45system 15:34.05elapsed 788%CPU (0avgtext+0avgdata 2269432maxresident)k
873663inputs+168outputs (384major+1068557230minor)pagefaults 0swaps

この場合には、「No GPU support is available」としてGPUは利用していません。代わりに、以下のように、「Intel Corei9-12900H」により計算が行われています。
解析時間は、「5886.0秒」なので、GPU利用の場合50倍高速であることが分かります。

例題スクリプトの修正と改良

カルマン渦の分析

このスクリプトの設定からレイノルズ数を分析すると、約10であり層流となっていました。最初の概念図のように、カルマン渦を再現することを考えます。

ここは、単に u_max を大きくするだけではなく、球周り流れのレイノルズ数 $Re_D$ を、意図的に設定するように設定します。今回の目的が「球の後方に周期的なカルマン渦（渦放出）を発生させる」ことなら、まず $Re_D \approx 300$ を目標にすることをおすすめします。

球周り流れの研究では、$Re_D$ が約210で軸対称性が崩れ、約270～275を超えると非定常流れ・周期的な渦放出が現れるとされています。特に $Re=300$ では周期的な渦放出が報告されています。

【1. 現在の設定】
スクリプトでは、$\omega=1.6, u_{\max}=0.04$ で、球の半径は、$r=\frac{64}{12}=5.333$、
したがって、$D=2r=10.667$ です。

D3Q19では、$c_s^2=\frac13$ なので、LBMの動粘性係数は、$\nu=c_s^2\left(\frac{1}{\omega}-\frac12\right)$ です。

現在の $\omega=1.6$ では、$nu=0.04167$ となり、$Re_D=\frac{u_{\max}D}{\nu}$ から、$Re_D=\frac{0.04\times10.667}{0.04167}\approx10.2$ でした。

つまり、現在の条件ではカルマン渦が周期的に放出される領域からかなり離れています。

【2. 目標を (Re_D=300) にする】
今回、$Re_D=300$ を目標にします。すると、$Re_D=\frac{UD}{\nu}$ なので、$\nu=\frac{UD}{Re_D}$ です。現在の球径 $D=10.667$ をそのまま使い、$U=0.08$ とすると、$\nu=\frac{0.08\times10.667}{300}$ で、$nu\approx0.002844$ となります。

この $\nu$ に対応する $\omega$ は、$\nu=\frac13\left(\frac1\omega-\frac12\right)$ なので、$\omega=\frac{1}{3\nu+0.5}$ より、$\omega\approx1.9664$ となります。
ただし、この条件は $\omega$ が2にかなり近いため、数値的な安定性の面ではあまり余裕がありません。

【3. そこで、球を少し大きくする】
例えば、$r=8$ として、$D=16$ にします。そして、$U=0.08$ とします。このとき $Re_D=300$ に必要な粘性係数は、$\nu=\frac{0.08\times16}{300} = 0.004267$ です。

これに対応する緩和パラメータは、$\omega=\frac{1}{3(0.004267)+0.5}$ となり、$\omega\approx1.9501$ です。こちらの方が、現在の球径のまま $Re=300$ にするより少し余裕があります。

【4. 修正する主要パラメータ】
したがって、まず次のように変更するのがよいです。

omega = 1.9501
u_max = 0.08
num_steps = 20000
post_process_interval = 100

そして球の半径を、以下とします。

sphere_radius = grid_shape[1] // 8

すると、$r=64/8=8$ なので、$D=16$ になります。そして、$Re_D=\frac{0.08\times16}{0.004267}\approx300$ です。つまり、$Re_D\approx300$ を狙えます。

【5. なぜ num_steps を増やすのか】
ここも非常に重要です。カルマン渦は定常状態ではありません。

定常流 → 不安定化 → 渦が発生 → 渦が下流へ流れる → 次の渦 → 次の渦 → ...

という時間的な変化を観察する必要があります。したがって、num_steps = 10000よりも、num_steps = 20000程度にして、時間変化を見る方がよいでしょう。

【6. post_process_intervalも重要】
現在は、post_process_interval = 1000です。カルマン渦を画像で追跡するなら、1000ステップは少し粗い可能性があります。最初は、post_process_interval = 100程度にして、「step 0 → step 100 →...」と細かく保存することをおすすめします。

【7. スクリプトの修正の手順】
以下の手順で、オリジナルのスクリプトをとして保存して、層流の解析画像を削除してから、カルマン渦用のスクリプトを作り、環境変数をGPU設定に戻して実行してみます。

(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ cp flow_past_sphere_3d.py flow_past_sphere_3d-org.py
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ rm *.png
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ ls
flow_past_sphere_3d-org.py  flow_past_sphere_3d.py
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ vi flow_past_sphere_3d.py
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ JAX_PLATFORMS=cuda
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ python flow_past_sphere_3d.py
Single-GPU support is available: 1 GPU detected.
registered bc RegularizedBC_8241249715775 with id 1
registered bc FullwayBounceBackBC_8241249715781 with id 2
registered bc ExtrapolationOutflowBC_8241249734627 with id 3
registered bc HalfwayBounceBackBC_8241249715232 with id 4
WARNING! Assuming no-slip condition for BC type = HalfwayBounceBackBC!
Warp 1.17.0 initialized:
   CUDA Toolkit 12.9, Driver 13.2
   Devices:
     "cpu"      : "x86_64"
     "cuda:0"   : "NVIDIA GeForce RTX 3070 Ti Laptop GPU" (8 GiB, sm_86, mempool enabled)
   Kernel cache:
     /home/dalab/.cache/warp/1.17.0
CD=0.0, CL=0.0
Post-processed step 0: Saved u_magnitude slice at y=32
Completed step 0. Time elapsed for 100 steps: 1.278503 seconds.
CD=2.2415826320648193, CL=-0.09965435415506363
:

試してみたのですが、これではカルマン渦が再現できませんでした。また解析結果が多すぎるので、さらに以下の設定に調整しました。

u_max = 0.10
num_steps = 10000
post_process_interval = 100
sphere_radius = grid_shape[1] // 6

改めて実行してみます。

(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ vi flow_past_sphere_3d.py
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ rm *.png
(XLBjaxGPU) dalab@dalab2312nb01:~/XLBjaxgpuWork/FLOW3D$ python flow_past_sphere_3d.py

以下のように、カルマン渦が確認できました。