OpenFOAM | CFD WITH A MISSION

OpenFOAM のベクトル演算

OpenFOAM には，ベクトルに対する様々な演算子が用意されています．次のプログラムを実行して，それぞれの演算子の意味を確認しましょう．

Test Program

#include "vector.H"
#include "IOstreams.H"

using namespace Foam;

int main()
{
    vector a(0.2, -0.51, 0.5);
    Info<< "cmptMax(a) = " << cmptMax(a) << endl;
    Info<< "cmptMin(a) = " << cmptMin(a) << endl;
    Info<< "cmptSum(a) = " << cmptSum(a) << endl;
    Info<< "cmptAv(a) = " << cmptAv(a) << endl;
    Info<< "cmptProduct(a) = " << cmptProduct(a) << endl;
    Info<< "cmptMag(a) = " << cmptMag(a) << endl;

    vector b(1.0, 2.0, 3.0);
    vector c(-3.0, 4.0, 1.0);
    Info<< "cmptMultiply(b, c) = " << cmptMultiply(b, c) << endl;
    Info<< "max(b, c) = " << max(b, c) << endl;
    Info<< "min(b, c) = " << min(b, c) << endl;

    return 0;
}

#include "vector.H"

#include "IOstreams.H"

using namespace Foam;

int main()

{

vector a(0.2, -0.51, 0.5);

Info<< "cmptMax(a) = " << cmptMax(a) << endl;

Info<< "cmptMin(a) = " << cmptMin(a) << endl;

Info<< "cmptSum(a) = " << cmptSum(a) << endl;

Info<< "cmptAv(a) = " << cmptAv(a) << endl;

Info<< "cmptProduct(a) = " << cmptProduct(a) << endl;

Info<< "cmptMag(a) = " << cmptMag(a) << endl;

vector b(1.0, 2.0, 3.0);

vector c(-3.0, 4.0, 1.0);

Info<< "cmptMultiply(b, c) = " << cmptMultiply(b, c) << endl;

Info<< "max(b, c) = " << max(b, c) << endl;

Info<< "min(b, c) = " << min(b, c) << endl;

return 0;

}

Results

cmptMax(a) = 0.5
cmptMin(a) = -0.51
cmptSum(a) = 0.19
cmptAv(a) = 0.0633333
cmptProduct(a) = -0.051
cmptMag(a) = (0.2 0.51 0.5)
cmptMultiply(b, c) = (-3 8 3)
max(b, c) = (1 4 3)
min(b, c) = (-3 2 1)

cmptMax(a) = 0.5

cmptMin(a) = -0.51

cmptSum(a) = 0.19

cmptAv(a) = 0.0633333

cmptProduct(a) = -0.051

cmptMag(a) = (0.2 0.51 0.5)

cmptMultiply(b, c) = (-3 8 3)

max(b, c) = (1 4 3)

min(b, c) = (-3 2 1)

Meaning of the Operations

cmptMax

$$cmptMax({\bf a}) = max (a_{1},\;a_{2},\;a_{3})$$

cmptMin

$$cmptMin({\bf a}) = min (a_{1},\;a_{2},\;a_{3})$$

cmptSum

$$cmptSum({\bf a}) = a_{1} + a_{2} + a_{3}$$

cmptAv

$$cmptAv({\bf a}) = \frac{a_{1} + a_{2} + a_{3}}{3}$$

cmptProduct

$$cmptProduct({\bf a}) = a_{1} \times a_{2} \times a_{3}$$

cmptMag

$$cmptMag({\bf a}) = (|a_{1}|\;\;|a_{2}|\;\;|a_{3}|)$$

cmptMultiply

$$cmptMultiply({\bf b},\;{\bf c}) = (b_{1}\times c_{1}\;\;b_{2}\times c_{2}\;\;b_{3}\times c_{3})$$

max

$$max({\bf b},\;{\bf c}) = (max(b_{1},\;c_{1})\;\;max(b_{2},\;c_{2})\;\;max(b_{3},\;c_{3}))$$

min

$$min({\bf b},\;{\bf c}) = (min(b_{1},\;c_{1})\;\;min(b_{2},\;c_{2})\;\;min(b_{3},\;c_{3}))$$

これらのベクトル演算子は，VectorSpaceI.H ファイルで実装されています．

境界に隣接するセルからなる cellZone の生成

ある境界に隣接するセルのみから構成される cellZone は，setSet ユーティリティを使用して，以下のような手順で生成することができます．

readline>faceSet ductwall new patchToFace duct
readline>cellSet layer1Cells new faceToCell ductwall owner
readline>cellZoneSet layer1Zone new setToCellZone layer1Cells

readline>faceSet ductwall new patchToFace duct

readline>cellSet layer1Cells new faceToCell ductwall owner

readline>cellZoneSet layer1Zone new setToCellZone layer1Cells

１行目のコマンドで，境界名が duct の境界フェイスから，ductwall という名前の faceSet を生成しています．２行目のコマンドでは，この faceSet に属するフェイスの owner セルから構成される cellSet を生成し，layer1Cells と名付けています．最後に３行目のコマンドで，この cellSet から，layer1Zone という名前の cellZone を生成しています．

ここで，２行目のコマンドで使用している faceToCell オプションは，指定した faceSet に属する face の neighbour セル，owner セル，またはその両方のセル (any) から構成される cellSet を生成することができます．また，all オプションを指定した場合には，全てのフェイスが指定した faceSet に属するセルから構成される cellSet を生成します．境界フェイスの場合，owner セルしか存在しないので，owner と any のどちらを指定しても同じ cellSet が生成されます．

Usage: faceToCell <faceSet> neighbour|owner|any|all

Select cells that are the owner|neighbour|any of the faces in the faceSet 
or where all faces are in the faceSet

Usage: faceToCell <faceSet> neighbour|owner|any|all

Select cells that are the owner|neighbour|any of the faces in the faceSet

or where all faces are in the faceSet

実装は，faceToCell.C ファイルで確認できます．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "faceToCell.H"
#include "polyMesh.H"
#include "faceSet.H"

#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(faceToCell, 0);
    addToRunTimeSelectionTable(topoSetSource, faceToCell, word);
    addToRunTimeSelectionTable(topoSetSource, faceToCell, istream);

    template<>
    const char* Foam::NamedEnum
    <
        Foam::faceToCell::faceAction,
        4
    >::names[] =
    {
        "neighbour",
        "owner",
        "any",
        "all"
    };
}


Foam::topoSetSource::addToUsageTable Foam::faceToCell::usage_
(
    faceToCell::typeName,
    "\n    Usage: faceToCell <faceSet> neighbour|owner|any|all\n\n"
    "    Select cells that are the owner|neighbour|any"
    " of the faces in the faceSet or where all faces are in the faceSet\n\n"
);

const Foam::NamedEnum<Foam::faceToCell::faceAction, 4>
    Foam::faceToCell::faceActionNames_;


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

void Foam::faceToCell::combine(topoSet& set, const bool add) const
{
    // Load the set
    faceSet loadedSet(mesh_, setName_);


    // Handle owner/neighbour/any selection
    forAllConstIter(faceSet, loadedSet, iter)
    {
        const label faceI = iter.key();

        if ((option_ == OWNER) || (option_ == ANY))
        {
            const label cellI = mesh_.faceOwner()[faceI];

            addOrDelete(set, cellI, add);
        }

        if (mesh_.isInternalFace(faceI))
        {
            if ((option_ == NEIGHBOUR) || (option_ == ANY))
            {
                const label cellI = mesh_.faceNeighbour()[faceI];

                addOrDelete(set, cellI, add);
            }
        }
    }

    // Handle all selection.
    if (option_ == ALL)
    {
        // Count number of selected faces per cell.

        Map<label> facesPerCell(loadedSet.size());

        forAllConstIter(faceSet, loadedSet, iter)
        {
            const label faceI = iter.key();
            const label own = mesh_.faceOwner()[faceI];

            Map<label>::iterator fndOwn = facesPerCell.find(own);

            if (fndOwn == facesPerCell.end())
            {
                facesPerCell.insert(own, 1);
            }
            else
            {
                fndOwn()++;
            }

            if (mesh_.isInternalFace(faceI))
            {
                label nei = mesh_.faceNeighbour()[faceI];

                Map<label>::iterator fndNei = facesPerCell.find(nei);

                if (fndNei == facesPerCell.end())
                {
                    facesPerCell.insert(nei, 1);
                }
                else
                {
                    fndNei()++;
                }
            }
        }

        // Include cells that are referenced as many times as they have faces
        // -> all faces in set.
        forAllConstIter(Map<label>, facesPerCell, iter)
        {
            const label cellI = iter.key();

            if (iter() == mesh_.cells()[cellI].size())
            {
                addOrDelete(set, cellI, add);
            }
        }
    }
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //

void Foam::faceToCell::combine(topoSet& set, const bool add) const

{

// Load the set

faceSet loadedSet(mesh_, setName_);

// Handle owner/neighbour/any selection

forAllConstIter(faceSet, loadedSet, iter)

{

const label faceI = iter.key();

if ((option_ == OWNER) || (option_ == ANY))

{

const label cellI = mesh_.faceOwner()[faceI];

addOrDelete(set, cellI, add);

}

if (mesh_.isInternalFace(faceI))

{

if ((option_ == NEIGHBOUR) || (option_ == ANY))

{

const label cellI = mesh_.faceNeighbour()[faceI];

addOrDelete(set, cellI, add);

}

// Handle all selection.

if (option_ == ALL)

{

// Count number of selected faces per cell.

Map<label> facesPerCell(loadedSet.size());

forAllConstIter(faceSet, loadedSet, iter)

{

const label faceI = iter.key();

const label own = mesh_.faceOwner()[faceI];

Map<label>::iterator fndOwn = facesPerCell.find(own);

if (fndOwn == facesPerCell.end())

{

facesPerCell.insert(own, 1);

}

else

{

fndOwn()++;

}

if (mesh_.isInternalFace(faceI))

{

label nei = mesh_.faceNeighbour()[faceI];

Map<label>::iterator fndNei = facesPerCell.find(nei);

if (fndNei == facesPerCell.end())

{

facesPerCell.insert(nei, 1);

}

else

{

fndNei()++;

}

// Include cells that are referenced as many times as they have faces

// -> all faces in set.

forAllConstIter(Map<label>, facesPerCell, iter)

{

const label cellI = iter.key();

if (iter() == mesh_.cells()[cellI].size())

{

addOrDelete(set, cellI, add);

}

OpenFOAM では，それぞれのセルのフェイスの両側に位置する２つのセルを，owner セルと neighbour セルと呼んでいます（境界フェイスの場合は，neighbour セルのみ存在します）．ここで，セル番号が小さい方のセルが owner セルです．ここら辺の話題については，こちらのスライドで触れています．

OpenFOAM －空間の離散化と係数行列の取り扱い（Spatial Discretization and Coefficient Matrix）－ from Fumiya Nozaki

neighbour のスペルちょっと違和感ありますか？OpenFOAM Code Style Guide の最後の Orthography の項でも言及されているように，OpenFOAM はもともと，英国産のソフトウェアなので，今でも英国式のスペリングが好まれています．他には，centre なども頻出しますね．

cellSource を使用した cellZone におけるセル値の演算処理

OpenFOAM の Function Objects 機能を活用することで，結果処理の効率化が可能です．その１つである cellSource を使用すると，ユーザーが指定した cellZone (セルの集合) に属するセルの変数値に対して，セル体積を重みとした加重平均や体積積分などの処理を計算実行中に行って，値をモニタリングすることが可能です．例えば，porous という名前の cellZone において体積を重みとした温度の加重平均値を計算したい場合には，ケースの system/controlDict ファイルに下記のような設定を行います．

functions
{
    averageT
    {
        type            cellSource;
        functionObjectLibs ("libfieldFunctionObjects.so");
        enabled         true;
        log             true;
        valueOutput     true;
        source          cellZone;
        sourceName      porous;
        operation       volAverage;
        weightField     none;
        fields
        (
            T
        );
    }
}

functions

{

averageT

{

type cellSource;

functionObjectLibs ("libfieldFunctionObjects.so");

enabled true;

log true;

valueOutput true;

source cellZone;

sourceName porous;

operation volAverage;

weightField none;

fields

(

);

}

実行する演算処理を operation に指定します．選択可能な 11 個のオプションが，cellSource.C ファイルで確認できます．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "cellSource.H"
#include "fvMesh.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    template<>
    const char* NamedEnum<fieldValues::cellSource::sourceType, 2>::names[] =
    {
        "cellZone",
        "all"
    };


    template<>
    const char* NamedEnum<fieldValues::cellSource::operationType, 11>::names[] =
    {
        "none",
        "sum",
        "sumMag",
        "average",
        "weightedAverage",
        "volAverage",
        "weightedVolAverage",
        "volIntegrate",
        "min",
        "max",
        "CoV"
    };

    namespace fieldValues
    {
        defineTypeNameAndDebug(cellSource, 0);
        addToRunTimeSelectionTable(fieldValue, cellSource, dictionary);
    }
}

template<>

const char* NamedEnum<fieldValues::cellSource::operationType, 11>::names[] =

{

"none",

"sum",

"sumMag",

"average",

"weightedAverage",

"volAverage",

"weightedVolAverage",

"volIntegrate",

"min",

"max",

"CoV"

};

それぞれの演算オプションの実装は，cellSourceTemplates.C ファイルで確認できます．

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "cellSource.H"
#include "volFields.H"

// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //

template<class Type>
bool Foam::fieldValues::cellSource::validField(const word& fieldName) const
{
    typedef GeometricField<Type, fvPatchField, volMesh> vf;

    if (obr_.foundObject<vf>(fieldName))
    {
        return true;
    }

    return false;
}


template<class Type>
Foam::tmp<Foam::Field<Type>> Foam::fieldValues::cellSource::setFieldValues
(
    const word& fieldName,
    const bool mustGet
) const
{
    typedef GeometricField<Type, fvPatchField, volMesh> vf;

    if (obr_.foundObject<vf>(fieldName))
    {
        return filterField(obr_.lookupObject<vf>(fieldName));
    }

    if (mustGet)
    {
        FatalErrorInFunction
            << "Field " << fieldName << " not found in database"
            << abort(FatalError);
    }

    return tmp<Field<Type>>(new Field<Type>(0.0));
}


template<class Type>
Type Foam::fieldValues::cellSource::processValues
(
    const Field<Type>& values,
    const scalarField& V,
    const scalarField& weightField
) const
{
    Type result = Zero;
    switch (operation_)
    {
        case opSum:
        {
            result = sum(values);
            break;
        }
        case opSumMag:
        {
            result = sum(cmptMag(values));
            break;
        }
        case opAverage:
        {
            result = sum(values)/values.size();
            break;
        }
        case opWeightedAverage:
        {
            result = sum(weightField*values)/sum(weightField);
            break;
        }
        case opVolAverage:
        {
            result = sum(V*values)/sum(V);
            break;
        }
        case opWeightedVolAverage:
        {
            result = sum(weightField*V*values)/sum(weightField*V);
            break;
        }
        case opVolIntegrate:
        {
            result = sum(V*values);
            break;
        }
        case opMin:
        {
            result = min(values);
            break;
        }
        case opMax:
        {
            result = max(values);
            break;
        }
        case opCoV:
        {
            Type meanValue = sum(values*V)/sum(V);

            const label nComp = pTraits<Type>::nComponents;

            for (direction d=0; d<nComp; ++d)
            {
                scalarField vals(values.component(d));
                scalar mean = component(meanValue, d);
                scalar& res = setComponent(result, d);

                res = sqrt(sum(V*sqr(vals - mean))/sum(V))/mean;
            }

            break;
        }
        default:
        {
            // Do nothing
        }
    }

    return result;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

template<class Type>

Type Foam::fieldValues::cellSource::processValues

(

const Field<Type>& values,

const scalarField& V,

const scalarField& weightField

) const

{

Type result = Zero;

switch (operation_)

{

case opSum:

{

result = sum(values);

break;

}

case opSumMag:

{

result = sum(cmptMag(values));

break;

}

case opAverage:

{

result = sum(values)/values.size();

break;

}

case opWeightedAverage:

{

result = sum(weightField*values)/sum(weightField);

break;

}

case opVolAverage:

{

result = sum(V*values)/sum(V);

break;

}

case opWeightedVolAverage:

{

result = sum(weightField*V*values)/sum(weightField*V);

break;

}

case opVolIntegrate:

{

result = sum(V*values);

break;

}

case opMin:

{

result = min(values);

break;

}

case opMax:

{

result = max(values);

break;

}

case opCoV:

{

Type meanValue = sum(values*V)/sum(V);

const label nComp = pTraits<Type>::nComponents;

for (direction d=0; d<nComp; ++d)

{

scalarField vals(values.component(d));

scalar mean = component(meanValue, d);

scalar& res = setComponent(result, d);

res = sqrt(sum(V*sqr(vals - mean))/sum(V))/mean;

}

break;

}

default:

{

// Do nothing

}

return result;

}

ここで，V はセル体積を表しているので，volIntegrate オプションは，体積積分を計算します．volAverage は，この体積積分値を cellZone の合計体積で割ることで，体積を重みとした加重平均値を計算します．これに対して，average は，単なる算術平均を計算します．また，名前に weighted が付いているオプションでは，加重平均の重みに使用する変数を weightField に指定します．sumMag オプションで使用されている cmptMag 演算子についてはこちらの記事で取り上げています．

選択したオプションが，意図した演算ではなかったなんてことが無いように，１度はしっかり確認しておきたいですね．

References

体積平均値の計算

Structure of src Directory in OpenFOAM part1

OpenFOAM の各種ライブラリのソースコードは，src ディレクトリに設置されています．例えば，2015年のオープン CAE シンポジウムの講習会で取り上げた境界条件のライブラリ libfiniteVolume のソースコードは，src/finiteVolume/ 以下に設置されています．第1回目の今回は，以下の4つのディレクトリを取り上げます．

fvOptions

既存のソルバーのソースコードを直接編集することなく支配方程式にソース項（ポーラスメディアの抵抗ソースや熱交換器の体積熱ソースなど）を追加するなどのカスタマイズを行うことができる fvOptions 機能のライブラリ libfvOptions のソースコードが設置されています．各ソルバーが fvOptions 機能に対応しているかどうかは，ソルバーのソースコードで確認できます．下記は，simpleFoam の例です．

#include "fvCFD.H"

#include "singlePhaseTransportModel.H"

#include "turbulentTransportModel.H"

#include "simpleControl.H"

#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])

{

#include "setRootCase.H"

#include "createTime.H"

#include "createMesh.H"

simpleControl simple(mesh);

#include "createFields.H"

#include "createMRF.H"

#include "createFvOptions.H"

#include "initContinuityErrs.H"

fvOptions で使用可能なタイプのリストを，こちらのページにまとめています．

renumber

セル番号の付け替え（Ordering）を行うための手法（Cuthill-McKee，Sloan アルゴリズムなど）のソースコードが設置されています．renumberMesh ユーティリティは，このライブラリ librenumberMethods を使用して，セル番号の付け替えにより係数行列のバンド幅を小さくすることで，連立方程式の反復解法の計算時間の短縮化を図ります．デフォルトの手法である Cuthill-McKee 以外の手法を使用するためには，設定ファイルを用意し，実行時に -dict オプションでそのパスを指定します．

CuthillMcKee: デフォルトの手法
Sloan: boost C++ ライブラリの Sloan アルゴリズムを使用
manual
random
spring
structured

最初の２つの手法は，特別なパラメータ等のコントロールなしに使用することができます．それぞれ，バンド幅を小さくするのか，profile を小さくするのか等の特徴があります．

dynamicFvMesh

OpenFOAM では，メッシュの移動変形機能を Dynamic Mesh という名前で呼んでいます．並進・回転などの剛体運動，１次元方向の伸び縮み運動など，トポロジー変化を伴わないメッシュの移動変形のライブラリのソースコードが設置されています．次のスライドで詳しく説明しています．

[slideshare id=31811977&doc=openfoamdynamicmesh20140302-140302022414-phpapp01&w=500]

topoChangerFvMesh

トポロジー変化を伴うメッシュの移動変形のライブラリのソースコードが設置されています．

上記のリンクのように，OpenFOAM のリリースノートには，新機能や機能変更に関して比較的詳しい説明が記載されています．自分が使用しているバージョンに対応する情報を収集するためにも，良くチェックするようにしたいですね．

その2に続きます．