はじめに

計算流体力学 (CFD) の格子計算は，各格子点が独立して同じ演算を繰り返す構造を持つ．これはGPUが得意とする「大量のデータに同じ命令を並列適用するSIMT (Single Instruction Multiple Threads)」の典型例であり，うまく実装すれば CPU のみの実装に比べて 10 倍以上の高速化が期待できる．

CUDA Fortran は NVIDIA HPC SDK に含まれる nvfortran コンパイラが提供する Fortran 拡張で，既存の Fortran コードに GPU 計算を導入できるのが大きな特徴である．CUDA C/C++ と同様の概念（スレッド階層・device 属性・カーネル起動）を Fortran 構文で書けるため，科学技術計算の資産を活かしやすい．

この記事では，圧縮性流れの標準的なベンチマーク問題であるSod衝撃波管問題を題材に，3種類の GPU 実装バリアントを通じて CUDA Fortran の基本を解説する．

バリアント	ディレクトリ	GPU の使い方	メモリレイアウト
① AoS ディレクティブ	`keep_directive_AoS/`	`!$cuf kernel do`	AoS: `q(3, nx)`
② SoA ディレクティブ	`keep_directive_SoA/`	`!$cuf kernel do`	SoA: `q(nx, 3)`
③ 明示的カーネル	`keep_kernel/`	`attributes(global)`	SoA: `q(nx, 3)`

3つのコードはすべて同じ物理問題を解く．異なるのは GPU プログラミングのスタイルとメモリレイアウトだけで，それがコードのシンプルさとパフォーマンスにどう影響するかが本記事のメインテーマである．

サンプルコードは以下のリポジトリに公開している:
https://github.com/htymjun/CFCF

問題設定: Sod衝撃波管

物理的な設定

Sod衝撃波管 (Sod shock tube) は，長さ $L_x$ の 1 次元管の中央に隔膜を置き，左右で圧力・密度を変えた初期状態から始める問題だ．解析解が存在するため，数値スキームの精度検証に広く使われる．

変数	左側 ($x \leq L_x/2$)	右側 ($x > L_x/2$)
密度 $\rho$	$\rho_l = 1.293\ \mathrm{kg/m^3}$	$0.125,\rho_l$
速度 $u$	$0$	$0$
圧力 $p$	$\rho_l R T_{lr}$	$0.1,\rho_l R T_{lr}$

隔膜を取り除くと，衝撃波・接触不連続・膨張波の 3 つの波が伝播する．

支配方程式: 1次元圧縮性 Navier-Stokes 方程式

保存形式の 1 次元圧縮性 Navier-Stokes 方程式：

$$\frac{\partial \mathbf{Q}}{\partial t} + \frac{\partial (\mathbf{E}_c - \mathbf{E}_v)}{\partial x} = 0$$

保存変数ベクトル $\mathbf{Q}$，対流フラックス $\mathbf{E}$，粘性フラックス $\mathbf{E}_v$：

$$
\mathbf{Q} = \begin{pmatrix} \rho \ \rho u \ E \end{pmatrix}, \quad
\mathbf{E}_c = \begin{pmatrix} \rho u \ \rho u^2 + p \ (E + p) u \end{pmatrix}, \quad
\mathbf{E}_v = \begin{pmatrix} 0 \ \tau \ \tau u - q_T \end{pmatrix}
$$

ここで $E = \rho e + \frac{1}{2}\rho u^2$ は全エネルギー，$p = (\gamma-1)(E - \frac{1}{2}\rho u^2)$ は状態方程式から計算する．

粘性応力と熱流束はSutherlandの粘性則に従う温度依存粘性 $\mu(T)$ を使って：

$$\tau = \frac{4}{3}\mu \frac{\partial u}{\partial x}, \qquad q_T = -\frac{\mu C_p}{\mathrm{Pr}} \frac{\partial T}{\partial x}$$

$$\mu(T) = \mu_0 \left(\frac{T_0 + S}{T + S}\right)\left(\frac{T}{T_0}\right)^{3/2}$$

KEEPスキーム (対流フラックス)

通常の中心差分で対流項を離散化すると，運動エネルギーが非物理的に増大して計算が発散することがある．KEEP (Kinetic Energy and Entropy Preserving) スキームは，離散レベルで運動エネルギーとエントロピーの保存性を厳密に満たすように対流フラックスを設計した手法で，数値安定性が高い．

界面 $j+\frac{1}{2}$ での対流フラックス定義（$C$ は質量フラックスの補助量）：

$$C_{j+1/2} = \frac{1}{4}(\rho_j + \rho_{j+1})(u_j + u_{j+1})$$

$$E_{\rho u, j+1/2} = \frac{1}{2} C_{j+1/2}(u_j + u_{j+1}) + \frac{1}{2}(p_j + p_{j+1})$$

$$E_{E, j+1/2} = {C_{j+1/2} \cdot \frac{R}{2(\gamma-1)}(T_j + T_{j+1})} + {\frac{1}{2} C_{j+1/2}, u_j u_{j+1}} + {\frac{1}{2}(u_j p_{j+1} + u_{j+1} p_j)}$$

時間積分は 1 次精度のオイラー法，CFL = 0.1 で安定性を確保している．

CUDA Fortranの基礎知識 (最小限)

コードを読む前に，CUDA Fortranで必要な概念を最低限押さえておく．

GPUスレッド階層

GPU 上の計算はスレッドの集合として実行される：

Grid（グリッド）
└─ Block（ブロック） × 複数
   └─ Thread（スレッド） × 複数

CFD では典型的に「1スレッドが1格子点（または1界面）を担当する」設計にする．スレッド階層の制御がパフォーマンスを左右する鍵になる．

`cudafor` モジュールと `device` 属性

use cudafor

real(dp), device :: x_gpu(nx)              ! GPU グローバルメモリに確保
real(dp), allocatable, device :: y_gpu(:)  ! 動的確保も可

device 属性を付けた変数は GPU のグローバルメモリに置かれる．CPU/GPU 間のメモリ転送は代入演算子 = で書ける（内部では cudaMemcpy が呼ばれる）：

q_gpu = q    ! CPU → GPU  (H2D)
q = q_gpu    ! GPU → CPU  (D2H)

`!$cuf kernel do` ディレクティブ

既存の Fortran ループをほぼそのまま GPU 並列化できる最も簡単な方法：

!$cuf kernel do(1)<<<grid_size, block_size>>>
do j = 1, nx
  a(j) = b(j) + c(j)   ! GPU 上で並列実行される
enddo

<<<grid_size, block_size>>> はグリッドサイズとブロックサイズの指定．合計スレッド数 = grid_size × block_size．

`attributes(global)` 明示的カーネル

より細かい制御が必要なときは，GPU カーネル関数として明示的に書く：

attributes(global) subroutine my_kernel(n, a, b, c)
  integer, intent(in), value :: n        ! スカラーは value 属性が必須
  real(dp), intent(in),  device :: a(n), b(n)
  real(dp), intent(out), device :: c(n)
  integer :: j
  j = (blockIdx%x - 1) * blockDim%x + threadIdx%x   ! グローバルスレッドインデックス
  if (j > n) return                                   ! 余分なスレッドの境界チェック
  c(j) = a(j) + b(j)
end subroutine my_kernel

! 呼び出し
call my_kernel<<<nx/128, 128>>>(n, a, b, c)

スカラー引数への value 属性と境界チェック if (j > n) return を忘れないこと．グリッドサイズを切り上げると余分なスレッドが生まれるため，境界チェックがないとメモリ違反になる．

コード構成と全体像

ディレクトリ構成

shock_tube_keep/
├── keep_directive_AoS/    # バリアント ①
│   ├── parameters.f90
│   ├── solver.f90
│   ├── main.f90
│   └── Makefile
├── keep_directive_SoA/    # バリアント ②
│   └── ...（同構成）
└── keep_kernel/           # バリアント ③
    └── ...（同構成）

各バリアントは parameters.f90 → solver.f90 → main.f90 の 3 ファイルで構成される（Fortran はモジュール依存順にコンパイルする必要がある）．

parameters.f90（3バリアント共通）

物理定数・格子パラメータをすべてここで定義する：

module parameters
  implicit none
  integer,  parameter :: dp   = kind(1.0d0)  ! 倍精度 kind 値
  integer,  parameter :: nx   = 4096         ! 格子点数
  integer,  parameter :: nt   = 2500         ! 時間ステップ数
  real(dp), parameter :: gamma = 1.4d0
  real(dp), parameter :: Rgas  = 287.03d0    ! 気体定数 [J/(kg·K)]
  real(dp), parameter :: Pr    = 0.72d0      ! プラントル数
  real(dp), parameter :: mu0   = 1.716d-5   ! 参照粘性 [Pa·s]
  real(dp), parameter :: T0    = 273.2d0    ! 参照温度 [K]
  real(dp), parameter :: S     = 111d0      ! サザーランド定数 [K]
  real(dp), parameter :: Re    = 25000d0    ! レイノルズ数
  real(dp), parameter :: Tlr   = 300.d0     ! 初期温度 [K]
  real(dp), parameter :: CFL   = 0.1d0
  real(dp), parameter :: rhol  = 1.293d0    ! 左側初期密度 [kg/m^3]
  ! 派生量
  real(dp), parameter :: Cp   = gamma * Rgas / (gamma - 1.d0)
  real(dp), parameter :: a    = sqrt(Rgas * 300d0)          ! 音速 [m/s]
  real(dp), parameter :: Lx   = Re * mu0 / (1.293d0 * a)   ! 領域長さ [m]
  real(dp), parameter :: dx   = Lx / dble(nx-1)
  real(dp), parameter :: dtdx = CFL * dx / (a * dx)         ! = CFL/a
end module parameters

dtdx = CFL/a は時間更新 q = q - dtdx * (flux_right - flux_left) で使う $\Delta t / \Delta x$ の値だ（dx が約分されている点に注意）．

main.f90 の構造（バリアント①を例に）

program main
  use cudafor
  use parameters, only : dp, nx, Rgas, gamma, Lx, rhol, a, Tlr
  use solver
  implicit none
  real(dp) :: q(3,nx)                         ! CPU メモリ上の保存変数
  real(dp), allocatable, device :: q_gpu(:,:) ! GPU メモリ上の保存変数

  call init(nx, q)       ! CPU 上で初期条件設定

  q_gpu = q              ! (1) CPU → GPU 転送
  call time_step(nx, q_gpu)   ! (2) GPU 上で 2500 ステップ計算
  q = q_gpu              ! (3) GPU → CPU 転送

  ! 結果を無次元化してファイル出力
  open(10, file='keep.dat')
  do n = 1, nx
    ...
    write(10,"(4es16.8e3)") x/Lx, rho/rhol, u/a, p/(rhol*Rgas*Tlr)
  enddo
end program main

CPU↔GPU 転送はそれぞれ1回だけ．時間積分の 2500 ステップはすべて GPU 上で完結する設計になっている．

メモリレイアウト: AoS vs SoA

3変数（$\rho$, $\rho u$, $E$）×格子点数の 2D 配列をどちらの軸で並べるかで，2つのレイアウトが生まれる．

AoS (Array of Structures): `q(3, nx)`

変数番号が先，格子点番号が後：

q(1,1)  q(2,1)  q(3,1) | q(1,2)  q(2,2)  q(3,2) | q(1,3) ...
←── 格子点1の3変数 ───→ ← 格子点2の3変数 ──→

隣接格子点 j と j+1 の同じ変数（例: q(1,j) と q(1,j+1)）はメモリ上で 3要素離れている（ストライドアクセス）．

SoA (Structure of Arrays): `q(nx, 3)`

格子点番号が先，変数番号が後：

q(1,1)  q(2,1)  q(3,1)  ...  q(nx,1) | q(1,2)  q(2,2) ...
←──────────── 変数1 (ρ) の全格子点 ──────────→

隣接格子点 j と j+1 の同じ変数（例: q(j,1) と q(j+1,1)）はメモリ上で隣接している．

GPUでの有利・不利

GPU では同一 Warp（32スレッドのグループ）が連続するメモリ番地にアクセスすると Coalesced アクセス となり，メモリ帯域を最大限に利用できる．逆に飛び飛びのアドレスに分散すると Non-coalesced となり，性能が大きく低下する．（H100等のL1, L2 cacheが強いGPUにおいてはこの限りではないかも…）

レイアウト	隣接格子点アクセス	Coalesced？
AoS: `q(3, nx)`	`q(1,j)` と `q(1,j+1)` が 3要素おき	Non-coalesced（ストライド）
SoA: `q(nx, 3)`	`q(j,1)` と `q(j+1,1)` が隣接	Coalesced

本コードのカーネルは，Warp 内で $j, j+1, j+2, \ldots$ を並列処理するため，SoA の方がメモリ効率の面で有利だ．

バリアント①: AoS ディレクティブ

solver.f90 の eev サブルーチン

フラックス計算を行う eev サブルーチンが GPU 処理の中心だ：

subroutine eev(nx, q, e, ev)
  integer, intent(in)             :: nx
  real(dp), intent(inout), device :: q(3,nx)        ! AoS レイアウト
  real(dp), intent(inout), device :: e(3,nx-1), ev(3,nx-1)
  real(dp), device :: qq(3,nx), T(nx), mu(nx)
  real(dp) c, m, g, ie, ke, Lp, tw, twu, kT
  integer :: j

  ! ステップ1: 保存変数 → 原始変数（密度, 速度, 圧力, 温度, 粘性係数）
  !$cuf kernel do(1)<<<32,4>>>
  do j = 1, nx
    qq(1,j) = q(1,j)                                           ! ρ
    qq(2,j) = q(2,j) / qq(1,j)                                ! u = ρu/ρ
    qq(3,j) = (gamma-1)*(q(3,j) - 0.5*qq(1,j)*qq(2,j)**2)    ! p = (γ-1)(E - ½ρu²)
    T(j)    = qq(3,j) / (qq(1,j) * Rgas)                      ! T = p/(ρR)
    mu(j)   = mu0 * ((T0+S)/(T(j)+S)) * (T(j)/T0)**1.5d0     ! サザーランド則
  enddo

  ! ステップ2: KEEPスキームと粘性フラックス（界面 j+1/2 を計算）
  !$cuf kernel do(1)<<<32,4>>>
  do j = 1, nx-1
    c  = 0.25d0 * (qq(1,j)+qq(1,j+1)) * (qq(2,j)+qq(2,j+1))  ! 質量フラックス補助量
    m  = 0.5d0  *  c * (qq(2,j)+qq(2,j+1))                    ! 運動量（対流部）
    g  = 0.5d0  * (qq(3,j)+qq(3,j+1))                         ! 圧力項
    ie = c * 0.5d0 * Rgas*(T(j)+T(j+1)) / (gamma-1.d0)        ! 内部エネルギー輸送
    ke = 0.5d0  *  c * qq(2,j) * qq(2,j+1)                    ! 運動エネルギー輸送
    Lp = 0.5d0  * (qq(2,j)*qq(3,j+1) + qq(2,j+1)*qq(3,j))    ! 圧力仕事

    tw  = 4d0*(mu(j)+mu(j+1))*(-qq(2,j)+qq(2,j+1)) / (6d0*dx) ! 粘性応力
    twu = 0.5d0 * tw * (qq(2,j)+qq(2,j+1))                    ! 粘性散逸
    kT  = Cp*(mu(j)+mu(j+1))*(-T(j)+T(j+1)) / (Pr*2d0*dx)    ! 熱流束

    e(1,j) = c;       e(2,j) = m + g;   e(3,j) = ie + ke + Lp  ! 対流フラックス
    ev(1,j) = 0.d0;  ev(2,j) = tw;     ev(3,j) = twu + kT      ! 粘性フラックス
  enddo
end subroutine eev

<<<32, 4>>> で合計 128 スレッドを起動している．nx = 4096 に対して 128 スレッドでは各スレッドが複数の格子点を担当する（コンパイラがループを自動分割する）．

AoSのデメリット: ステップ2で qq(2,j) と qq(2,j+1) を参照するとき，Warp の各スレッド（$j$, $j+1$, ..., $j+31$）がそれぞれ qq(2,*) にアクセスするが，AoS では変数番号が先のため隣接格子点が 3要素おきに並んでいる．これが Non-coalesced アクセスを引き起こす．

バリアント②: SoA ディレクティブ

AoS との違いはインデックスの順序だけで，物理計算は完全に同一だ．

変更点の比較

! ===== AoS (keep_directive_AoS) =====
real(dp), device :: qq(3,nx), T(nx), mu(nx)

!$cuf kernel do(1)<<<32,4>>>
do j = 1, nx
  qq(1,j) = q(1,j)
  qq(2,j) = q(2,j) / qq(1,j)
  qq(3,j) = (gamma-1)*(q(3,j) - 0.5*qq(1,j)*qq(2,j)**2)
  T(j) = qq(3,j) / (qq(1,j) * Rgas)
  mu(j) = mu0 * ((T0+S)/(T(j)+S)) * (T(j)/T0)**1.5d0
enddo

! ===== SoA (keep_directive_SoA) =====
real(dp), device :: rho(nx), u(nx), p(nx), T(nx), mu(nx)

!$cuf kernel do(1)<<<32,4>>>
do j = 1, nx
  rho(j) = q(j,1)
  u(j)   = q(j,2) / rho(j)
  p(j)   = (gamma-1)*(q(j,3) - 0.5*rho(j)*u(j)**2)
  T(j)   = p(j) / (rho(j) * Rgas)
  mu(j)  = mu0 * ((T0+S)/(T(j)+S)) * (T(j)/T0)**1.5d0
enddo

SoA では原始変数を別々の 1D 配列 rho(nx), u(nx), ... に格納する．こうすると Warp 内で同じ変数の連続格子点にアクセスするとき，アドレスが隣接するため Coalesced アクセスになる．

フラックス計算部分も同様にインデックスが変わる：

! AoS: e(1,j), e(2,j), e(3,j) → e(変数番号, 界面番号)
e(1,j) = c;  e(2,j) = m + g;  e(3,j) = ie + ke + Lp

! SoA: e(j,1), e(j,2), e(j,3) → e(界面番号, 変数番号)
e(j,1) = c;  e(j,2) = m + g;  e(j,3) = ie + ke + Lp

ディレクティブの書き方 !$cuf kernel do(1)<<<32,4>>> はまったく変わらない．

バリアント③: 明示的CUDAカーネル

SoA ディレクティブとの違いは，フラックス計算部分を attributes(global) カーネル関数 として明示的に書くことだ．これにより「1スレッドが1界面を担当する」きめ細かな設計が可能になる．

kernel サブルーチン

attributes(global) subroutine kernel(nx, rho, u, p, T, mu, e, ev)
  integer, intent(in), value    :: nx
  real(dp), intent(in),  device :: rho(nx), u(nx), p(nx), T(nx), mu(nx)
  real(dp), intent(out), device :: e(nx-1,3), ev(nx-1,3)
  real(dp) c, m, g, ie, ke, Lp, tw, twu, kT
  integer :: j, jt

  jt = threadIdx%x
  j  = (blockIdx%x - 1) * blockDim%x + jt    ! グローバルスレッドインデックス
  if (j > nx-1) return                         ! 境界チェック

  ! 界面 j+1/2 のフラックスを計算（SoA バリアントと同一の数式）
  c  = 0.25d0 * (rho(j)+rho(j+1)) * (u(j)+u(j+1))
  m  = 0.5d0  *  c * (u(j)+u(j+1))
  g  = 0.5d0  * (p(j)+p(j+1))
  ie = c * 0.5d0 * Rgas*(T(j)+T(j+1)) / (gamma-1.d0)
  ke = 0.5d0  *  c * u(j) * u(j+1)
  Lp = 0.5d0  * (u(j)*p(j+1) + u(j+1)*p(j))

  tw  = 4d0*(mu(j)+mu(j+1))*(-u(j)+u(j+1)) / (6d0*dx)
  twu = 0.5d0 * tw * (u(j)+u(j+1))
  kT  = Cp*(mu(j)+mu(j+1))*(-T(j)+T(j+1)) / (Pr*2d0*dx)

  e(j,1) = c;       e(j,2) = m + g;   e(j,3) = ie + ke + Lp
  ev(j,1) = 0.d0;  ev(j,2) = tw;     ev(j,3) = twu + kT
end subroutine kernel

eev サブルーチンからの呼び出し

subroutine eev(nx, q, e, ev)
  ...
  ! ステップ1: 原始変数抽出はディレクティブのまま
  !$cuf kernel do(1)<<<32,4>>>
  do j = 1, nx
    rho(j) = q(j,1);  u(j) = q(j,2) / rho(j)  ! ... 省略 ...
  enddo

  ! ステップ2: フラックス計算は明示的カーネルで
  call kernel<<<nx/128, 128>>>(nx, rho, u, p, T, mu, e, ev)
end subroutine eev

<<<nx/128, 128>>> = <<<32, 128>>> で 4096 スレッドを起動．nx - 1 = 4095 個の界面に対して 1 スレッド 1 界面で対応する（1スレッドだけ仕事なしになる）．

スレッドインデックスの計算イメージ

blockIdx%x=1, blockDim%x=128 のとき:
  threadIdx%x=1   → j = (1-1)*128 + 1 = 1    （界面 1+1/2）
  threadIdx%x=128 → j = (1-1)*128 + 128 = 128

blockIdx%x=2 のとき:
  threadIdx%x=1   → j = (2-1)*128 + 1 = 129
  ...
blockIdx%x=32 のとき:
  threadIdx%x=128 → j = (32-1)*128 + 128 = 4096  → if (j > nx-1) return で処理をスキップ

3バリアントの比較まとめ

項目	① AoS ディレクティブ	② SoA ディレクティブ	③ 明示的カーネル
メモリレイアウト	`q(3, nx)`	`q(nx, 3)`	`q(nx, 3)`
GPU 記述スタイル	`!$cuf kernel do`	`!$cuf kernel do`	`attributes(global)`
Coalesced アクセス	非 Coalesced	Coalesced	Coalesced
スレッド数 (フラックス)	128 (1スレッド≒32格子点)	128 (1スレッド≒32格子点)	4096 (1スレッド=1界面)
書きやすさ	既存 Fortran に近い	インデックス変更のみ	スレッドindex等が増える
最適化の自由度	低い（コンパイラ依存）	低い（コンパイラ依存）	高い（Shared Memory等も可）

どれを選ぶべきか

とにかく動かしたい・既存コードを流用したい → ① AoS ディレクティブ
変更量が最小で，まず動かして確認するのに適している．
パフォーマンスを上げたい・メモリレイアウトを最適化したい → ② SoA ディレクティブ
AoS から SoA への変更は機械的なインデックス順の入れ替えだけ．Coalesced アクセスになるため，メモリ帯域律速のケースでは速度向上が見込める．
完全なスレッド制御が必要・Shared Memory を使いたい → ③ 明示的カーネル
1スレッド1界面の設計で GPU の並列性を最大限活用できる．将来的に Shared Memory キャッシュや原子演算を使う場合もこちらを選ぶ．

ビルドと実行方法

環境要件

NVIDIA HPC SDK >= 24（nvfortran コンパイラ）
CUDA 対応 GPU
環境モジュール nvhpc-openmpi3 のロード

ビルド・実行

module load nvhpc-openmpi3   # 環境モジュールをロード

cd shock_tube_keep/keep_directive_AoS   # または他の2つ
make
./sod

実行すると以下が出力される：

keep.dat — 無次元化された解（列: $x/L_x$, $\rho/\rho_l$, $u/a$, $p/(\rho_l R T_{lr})$）
CPU→GPU転送，計算，GPU→CPU転送それぞれの計算時間

Makefileのコンパイルフラグ

FFLAGS = -cuda -acc -fast -mfma -Mchkptr -Minfo -gpu=ptxinfo,rdc,lto
SRC    = parameters.f90 solver.f90 main.f90

フラグ	意味
`-cuda`	CUDA Fortran を有効化
`-acc`	OpenACC を有効化（`!$cuf kernel do` に必要）
`-fast`	高レベル最適化（`-O2` + ベクトル化など）
`-mfma`	Fused Multiply-Add (FMA) 命令を使う
`-Minfo`	最適化情報を標準エラーに出力
`-gpu=cc**`	必要に応じてCompute Capabilityを指定
`-gpu=ptxinfo`	Stack frameやRegister spill等の情報を確認
`-gpu=rdc,lto`	Relocatable Device Code + リンク時最適化

おわりに

本記事では，Sod衝撃波管問題を例に CUDA Fortran の 3 つの実装パターンを比較した．

ディレクティブ (!$cuf kernel do) vs 明示的カーネル (attributes(global)): 記述の容易さと最適化の自由度のトレードオフ
AoS vs SoA: インデックスの順序がメモリアクセスパターンに直結し，Coalesced アクセスの有無がパフォーマンスに影響する

さらに学ぶために

このコードはシンプルな構成になっているが，さらに高性能な実装に向けた発展方向として：

Shared Memory の活用: 同一 Block 内で隣接格子点の値をキャッシュすれば，グローバルメモリアクセスを削減できる
高次精度スキーム: 現在は 1 次精度のオイラー時間積分．4 次 Runge-Kutta などへの拡張が考えられる
多次元化: 2D/3D への拡張では do(2), do(3) ディレクティブや dim3 型グリッドが必要
MPI + CUDA Fortran: 複数 GPU への分散並列化では MPI との組み合わせが定番

参考資料

CUDA Fortran Programming Guide (NVIDIA)
CUDA Fortran Interfaces Reference (NVIDIA)
Ruetsch, G. & Fatica, M. (2024). CUDA Fortran for Scientists and Engineers, Elsevier.
Kuya, Y. Totani, K. & Kawai, S. (2018). Kinetic energy and entropy preserving schemes for compressible flows by split convective forms, Journal of Computational Physics.