OpenFOAMでセルごとに割り当てるプロセッサの番号を指定する

動機

decomposeParでの分割方法はsimple、scotch、hierarchical、manualの4通りがある。
しかし前者3者で分割を行った時に、cyclicAMI面でうまく計算できないことがあった。
このためなるべく綺麗に分割したいと考え、manualで分割するためのutilityを作成した。

もしかしたらプロセッサを跨ぐ面が減れば大規模計算になった時にノード間通信料がへるかも。
ゼロ割エラーなどの場合はどの場所で止まったかのチェックに使えたりもする。

バージョンは4.1と5.xで動作確認済み。

やったこと

manualで指定するための設定ファイルを作成するようなutilityを作ってみた。
https://github.com/inabower/OpenFOAM4.1-utilities/tree/master/functionProcessor
下図はtestCase1について4つの方法で分割したそれぞれの結果を示している。色は割り当てられるプロセッサの番号を示している。
左3つは標準の分割方法。右が今回の方法。

拡大

外側と内側で分割してみた。
またinletとoutletはそれぞれ同一のプロセッサ内で計算させるようにしてみた。

作ったものの解説

functionProcessor

　実態としては実は何もしない。
空の"0/cellDist"と"constant/cellDecomposition"を作成するutility。
functionObject機能を通してその中身を編集する。

cellDistとcellDecompositionはどちらもどのセルがどのプロセッサに割り当てられるかを示している。
cellDistは上図のようにparaViewなどで表示するためのもの。
cellDecompositionはdecomposeParを実行するときに読み込むファイル。

例えば上githubのcase1では"system/constrolDict"の中は以下のようになっている。

system/controlDict

...

functions
{
    processor
    {
        libs ("libutilityFunctionObjects.so");
        type  coded;
        name  processor;

        codeEnd
        #{
            Info << "coded Function Start" << endl;

            volScalarField& dist = const_cast<volScalarField&>
            (mesh().lookupObject<volScalarField>("cellDist"));

            forAll(dist, i)
            {
                const vector& C = mesh().C()[i];
                const scalar& x = C[0];
                const scalar& y = C[1];
                const scalar& z = C[2];

                scalar& d = dist[i];
                const scalar r = 0.0254;
                const scalar RR = 0.1;

                d = 4;
                if(x<0.3)
                {
                  d=0;
                }
                else if(x<0.5)
                {
                  if(Foam::sqr(y)+Foam::sqr(z)<Foam::sqr(0.02))
                  {
                    d=1;
                  }
                  else
                  {
                    d=2;
                  }
                }
                else if(z<0.1)
                {
                  if
                  ( // equation of torus
                    Foam::sqr
                    (
                      Foam::sqrt(
                        Foam::sqr(x-0.5)+Foam::sqr(z-0.1)
                      ) - RR
                    )
                    + Foam::sqr(y) < Foam::sqr(0.02)
                  ) {
                    d = 1;
                  }
                  else
                  {
                    d = 2;
                  }
                }
                else if(z < 0.2)
                {
                  if(Foam::sqr(x-0.6)+Foam::sqr(y)<Foam::sqr(0.02))
                  {
                    d=1;
                  }
                  else
                  {
                    d=2;
                  }
                }
                else if(x<0.7)
                {
                  if
                  ( // equation of torus
                    Foam::sqr
                    (
                      Foam::sqrt(
                        Foam::sqr(x-0.7)+Foam::sqr(z-0.2)
                      ) - RR
                    )
                    + Foam::sqr(y) < Foam::sqr(0.02)
                  ) {
                    d = 1;
                  }
                  else
                  {
                    d = 2;
                  }
                }
                else
                {
                  d = 3;
                }
            }
        #};
    }
}

まずcellDistを読み込む。
編集したいのでconst_castをする。

volScalarField& dist = const_cast<volScalarField&>
(mesh().lookupObject<volScalarField>("cellDist"));

その後各セル毎のループの中で
中心座標をx,y,zに格納し、プロセッサ番号をdとしている。

const vector& C = mesh().C()[i];
const scalar& x = C[0];
const scalar& y = C[1];
const scalar& z = C[2];

scalar& d = dist[i];

そしてむかし習った数学を思い出しながらゴリゴリと場合分けしていく。

if(x<0.3)
{
  d=0;
}
else if(x<0.5)
{
  if(Foam::sqr(y)+Foam::sqr(z)<Foam::sqr(0.02))
  {
    d=1;
  }
  else
  {
    d=2;
  }
}

このように、controlDict内に分割方法を記載することで任意にどのセルをどのプロセッサに計算させるかを指定できる。
ただ何千並列になった場合にはこんなことやってられないと思うので、ワークステーション規模までの人におすすめ。
32並列くらいであればできなくもなかった。