03. 1 次元移流と周期境界

解く方程式

\[\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0, \quad a > 0 \text{ は定数}\]

  • 時間積分：陽的オイラー法

  • 空間離散化：1 次精度の風上差分。a > 0 なので、左隣との差を使います。

\[u_i^{n+1} = u_i^n - \frac{a \Delta t}{\Delta x}\left(u_i^n - u_{i-1}^n\right)\]

時間刻みは、クーラン数 \(C = a \Delta t / \Delta x\) が 0 より大きく 1 以下になるように選びます。 ここでは \(C = 0.37\) を使います。

格子とゴーストセル

物理セル NX 個の左右に、1 点ずつゴーストセルを置きます。 風上差分は左隣 u[i-1] しか読まないので左の 1 点で足りますが、MISO と同じく左右対称に置きます。 周期境界では、毎ステップ、ゴーストセルを反対側の物理セルの値で埋めてから更新します。

index   0  |  1  ...  NX  |  NX+1
    ghost   physical cells  ghost

実装する内容

  1. 長さ NX + 2 の配列 u と u_new を、それぞれ Array1D<double, Backend> として確保します。

  2. 物理セルの初期値は周期的なガウシアンとします。 \(L = 1\)、\(x_0 = 0.5\)、\(\sigma = 0.05\) とし、物理セル i の座標を x = (i - 1) * L / NX とします。 \(d = x - x_0\) を周期の幅 \(L\) で折り返して \(-L/2 \le d \le L/2\) にしてから、\(\exp(-d^2 / 2\sigma^2)\) を計算します。 折り返しには util::fmod(d, L)（d を L で割った余り。結果は -L と L の間）を使い、余りが L/2 より大きければ L を引き、-L/2 より小さければ L を足します。

  3. この初期値の式は、厳密解との比較でも使うので関数にします。 ラムダ式の中から呼ぶので、__host__ __device__ inline double gaussian(double x) のように __host__ __device__ を付けておきます（CPU ビルドでは何もしない指定です）。

  4. 各ステップの前に、左のゴーストセル u[0] へ最後の物理セル u[NX] を、右のゴーストセル u[NX + 1] へ最初の物理セル u[1] を写します。 この処理も、Range1D{0, 1}（1 回だけ実行する範囲）の for_each で書きます。

  5. 物理セルの範囲 Range1D{1, NX + 1} の for_each で、更新前の u を読み、更新後の値を u_new へ書きます。

  6. ステップの最後に std::swap(u, u_new) で 2 つの配列の役割を入れ替えます。 std::swap には <utility> を読み込みます。

Array1D はコピーできませんが、ムーブはできます（第 4 章）。 std::swap はムーブを使って中身のポインタを入れ替えるので、要素はコピーしません。 入れ替えの後、u という名前の変数は前の u_new のメモリを持っています。 View は取り出したときのポインタを持ち続けるので、入れ替えの後は u.view() で View を取り直します。 ゴーストセルを埋める処理と更新の処理を、それぞれ配列を引数に取る関数にして、関数の中で View を取り出すと、取り直しを忘れません。

更新前の配列へ途中で書き込まないので、各セルは同じ時刻の隣の値を読めます（第 7 章の「実行順に頼らない」）。

結果の確認

まず NX = 256 で 1 ステップだけ計算し、同じ入力とゴーストセルからホスト側の for 文で求めた値と比べます。 初期の振幅 1 を基準に、最大絶対誤差が 1e-14 以下であることを確かめます。

次に時間を進め、初期分布を距離 \(a t\) だけ平行移動した厳密解と比べます。 厳密解は、初期値の式の \(x\) を \(x - a t\) に置き換え、同じ周期の折り返しを使って求めます。 相対 L2 誤差は、物理セルだけで次のように求めます。

\[E = \sqrt{\frac{\sum_i (u_i-u_{\mathrm{exact},i})^2}{\sum_i u_{\mathrm{exact},i}^2}}\]

分子と分母の和は、前の課題の reduce で求められます（reduce を 2 回呼びます）。 NX を 256、512、1024 と増やし、\(C\) と終了時刻を揃えて誤差の変化を調べます。 たとえばステップ数を 200、400、800 とすれば、同じ時刻まで進みます。 この方法は 1 次精度なので、十分細かい格子では、格子幅を半分にすると誤差もおよそ半分になります。

CPU で 1 ステップの比較と時間発展を確かめたら、第 8 章の MPI へ進めます。 GPU で同じ計算を動かす準備は、第 9 章で学びます。

GPU で実行する段階

CUDA でも同じ入力で、1 ステップの比較と時間発展を行います。 ラムダ式の中から呼ぶ gaussian には __host__ __device__ が必要です。 FMA などにより、CPU と GPU では丸めが異なる場合があります。 完全一致を一律に求めず、上の許容誤差と、格子を細かくしたときの誤差の変化で確かめます。