「CuthillMckee法によるOpenFOAMのDIC前処理に関するスレッド並列化」を発表してきました

はじめまして、藤原と申します。 私はフィックスターズのインターンに2週間参加し、2018年4月21日に弊社で開催された第72回オープンCAE勉強会＠関東（流体など）【大崎】に成果発表をしてきました。 インターン後の今もアルバイトとして、プロジェクトに継続して参加しています。

この記事では特にプレゼンでは紹介しにくいソースコードに重点をおいて、発表内容の紹介をしたいと思います。

動機

弊社では2017年12月に、OpenCAEシンポジウム2017@名古屋 講演会で、『OpenFOAMのスレッド並列化に向けて基礎検討』(富岡稔、伊藤優希、丸石崇史、吉藤尚生)についての講演を行いました。 この講演の内容を要約すると、

従来のOpenFOAMの疎行列の格納形式「lduMatrix」では、並列化が不可能である。 そのため、OpenFOAMのスレッド並列を実現するためには、少なくとも格納形式を並列化できるものに変更しなければならない。 そこで、弊社ではlduMatrixの代わりに、「CSR形式」を使用することで実際にスレッド並列化するための見通しを立てることに成功した。 しかし、CSR形式を用いてもDIC前処理 など、一部の処理ではデータ依存性により並列化することが難しく、Cuthill-Mckee法やMulticoloring法などを用いて工夫する必要がある。

という内容でした。 このブログではCuthill-Mckee法を実装し、DIC前処理を実際にスレッド並列で動作することを確認したので、そのアルゴリズムと実装方法についてまとめさせていただきます。

DIC前処理を並列化するために

先述の発表にあるとおり、DIC前処理ではデータの依存性があるので、行毎に並列すると、ある一つのデータに関して計算中に書き込みと読み出しが同時に実行され、結果が不安定になる可能性があります。 計算結果が不安定にならないようにするために、以下の二種類の方法があります。

• 依存し合うデータが同じ色にならないように塗り分けて、色ごとに並列で計算する「Multicoloring法」
• データの依存関係をグラフで表現し、並列化する前と収束性は変わらない範囲で並列計算をする「Cuthill-Mckee法」

Cuthill-Mckee法(後述)とMulticoloring法の特徴を以下の表にまとめます。

今回はバンド幅が狭い点(キャッシュヒット率が高い)点と収束性が変わらない点に注目し、Cuthill-Mckee法を実装してみました。

Cuthill-Mckee法

アルゴリズム

もともと、Cuthill-Mckee法は並列化するためのアルゴリズムではなく、行列のfill-inやバンド幅を減らすためのアルゴリズムです。 Cuthill-Mckee法のアルゴリズムは、おおまかに以下の4ステップで構成されます。

1. 行列を隣接行列(グラフ)として捉える。
2. 最も次数が少ないノードを見つけ、そこを幅優先探索のスタート地点とする。
3. 幅優先探索をし、スタート地点からの距離を「レベル」と名付ける。
4. レベルをベースに行列を並び替える

このアルゴリズムだけでは、同じレベル内に隣接ノード(依存しあっている関係)が存在し、まだ並列化できません。 今回は、「レベルに新しくノードを追加する時に、すでに同じレベル内で隣接しているようなノードがある場合は追加しない」という処理を追加し、隣接ノードを排除しました。 詳細は中島研吾先生のスライドを参照してください。

検証用コード

いきなりOpenFOAM上に実装するのは苦なので、検証用コードにCuthill-Mckee法の関数を実装してみました。

• 行列として見る場合: i は行 j は列。
• グラフとして見る場合: i 行目 j 列目に値があれば、それは、i -> j に辺が伸びている

と捉えると読みやすいと思います。

結果

以下に行列を並び替えた結果を示しています。Multicoloring法に比べ、Cuthill-Mckee法はバンド幅が狭くなっている事がわかります。

original

Multicoloring法

Cuthill-Mckee法

収束性

Cuthill-Mckee法は並列化前と収束性は変わっておらず、Multicoloring法に比べ良くなっていることがわかります。

OpenFOAMへの実装

DICPreconditionerのコンストラクタでCuthill-Mckee法のレベルを計算し、DICPreconditioner::preconditionで求められたレベルを元に並列計算を行う。

という実装にしました。

もともとOpenFOAMに実装されていたpreconditionの前進代入に関しては、以下のようになっており、

for (label face=0; face<nFaces; face++)
{
    wAPtr[uPtr[face]] -= rDPtr[uPtr[face]]*upperPtr[face]*wAPtr[lPtr[face]];
}

弊社でCSR形式に書き換えることで以下のようなコードになりました。

for (label i=0; i<n; i++)
{
    for (label index = csr.offset[i]; index < csr.offset[i+1]; index++)
    {
        const auto j = csr.column[index];
        if(j < i)
        {
            wAPtr[i] -= rDPtr[i] * csr.data[index] * wAPtr[j];
        }
    }
}

ここでwAPtr[i]とwAPtr[j]が相互依存し、並列化が出来なかったため、Cuthill-Mckee法を導入した結果が、以下のソースコードになります。

for (auto level = decltype(maxLevel)(0); level < (maxLevel); level++)
{
    const auto rowStart = levelOffset[level];
    const auto rowEnd = levelOffset[level + 1];
    #pragma omp for
    for (auto rowIndex = rowStart; rowIndex < rowEnd; rowIndex++)
    {
        const auto i = row[rowIndex];
        const auto colStart = csr.offset[i];
        const auto colEnd = csr.offset[i + 1];
        for (auto colIndex = colStart; colIndex < colEnd; colIndex++)
        {
            const auto j = csr.column[colIndex];
            if (j < i)
            {
                wAPtr[i] -= rDPtr[i] * csr.data[colIndex] * wAPtr[j];
            }
        }
    }
}

row[]には、レベル順にソートされた行番号が入っており、levelOffset[]はレベルが変わるrowのindexを表しています。

たとえば、

であれば、 row は [1(レベル1),0(レベル2),2(レベル2),4(レベル2),3(レベル3)] になり、 levelOffset は [0,1,4,5] が格納されます。 ループを回す都合上、配列の終わりを表す 5 は入れておいたほうが楽なので入っています。

rowStartとrowEndは同じレベル内のrowのindexを表し、同じレベルであればどのような順番で計算してもiとjの間は依存関係がないため、

wAPtr[i] -= rDPtr[i] * csr.data[colIndex] * wAPtr[j];

のコードも、レベル内で並列計算しても安全に動作する。 という流れになっています。

結果、Cuthill-Mckee法を用いることでデータの依存関係があるDIC前処理を、安全に並列計算ができるようになりました。

今後の予定

まだOpenFOAMにはDIC前処理しかスレッド並列化を実装しておらず、またプロファイラも実装していないため、従来使用されてきたフラットMPIとの実行時間の比較が難しい状況になっています。 今後はDIC前処理以外の部分にもスレッド並列を適用し、ソルバ全体のスレッド並列化を試みつつ、プロファイラを実装し、個々の処理の実行時間の比較をしていきたいと考えています。

また、効果測定後にスレッド並列化したOpenFOAMのコードも公開する予定です。

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください。