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?

@htymjun

CUDA Fortranで始めるGPU加速CFD — Sod衝撃波管をKEEPスキームで解く3パターン比較

0
Posted at

はじめに

計算流体力学 (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 との組み合わせが定番

参考資料

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?