OpenFOAM のソースコードを眺めていると,fvc::reconstruct という演算子を目にします.例えば,buoyantSimpleFoam の UEqn.H ではこの演算子が下記のように使用されています.
| 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |     if (simple.momentumPredictor())     {         solve         (             UEqn          ==             fvc::reconstruct             (                 (                   - ghf*fvc::snGrad(rho)                   - fvc::snGrad(p_rgh)                 )*mesh.magSf()             )         );         fvOptions.correct(U);     } | 
これは,具体的にどのような演算を行っているのでしょうか?こんな時は,ソースコードを確認しましょう.ソースコードを参照して実装方法をすぐに確認できるのが,オープンソースの1つのメリットですね.この演算子のソースファイルは,fvcReconstruct.C です.このファイルの 76 行目に演算の定義が記述されています.
| 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 | template<class Type> tmp <     GeometricField     <         typename outerProduct<vector,Type>::type, fvPatchField, volMesh     > > reconstruct (     const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf ) {     typedef typename outerProduct<vector, Type>::type GradType;     const fvMesh& mesh = ssf.mesh();     surfaceVectorField SfHat(mesh.Sf()/mesh.magSf());     tmp<GeometricField<GradType, fvPatchField, volMesh>> treconField     (         new GeometricField<GradType, fvPatchField, volMesh>         (             IOobject             (                 "volIntegrate("+ssf.name()+')',                 ssf.instance(),                 mesh,                 IOobject::NO_READ,                 IOobject::NO_WRITE             ),             inv(surfaceSum(SfHat*mesh.Sf()))&surfaceSum(SfHat*ssf),             extrapolatedCalculatedFvPatchField<GradType>::typeName         )     );     treconField.ref().correctBoundaryConditions();     return treconField; } | 
私は一見しただけでは何を行っているのか分かりませんでした.ヘッダーファイル fvcReconstruct.H も確認してみます.
| 24 25 26 27 28 29 30 31 | InNamespace     Foam::fvc Description     Reconstruct volField from a face flux field. SourceFiles     fvcReconstruct.C | 
冒頭の Description の項目に,Reconstruct volField from a face flux field という記述があります.最初に挙げた buoyantSimpleFoam での使用例を確認すると,fvc::reconstruct(引数) の引数は,フェイス(セル界面)で値をもっています(face flux field).具体的には,
- snGrad(rho) は,フェイスの法線方向の密度の勾配を,
- mesh.magSf() は,フェイスの面積を
それぞれ計算するので,両方とも値はフェイスで保持しています.この例では,OpenFOAM の用語を使用すると,surfaceScalarField を引数をとして受け取り,セル中心で値をもつ volVectorField を返しています.
| 次回予告 | 
次回の投稿では,簡単な例でこの演算子を書き下して,数式を使用した説明にチャレンジしたいと思います.最終的には,重力のような体積力を考慮した計算において,圧力の変数として \(p\) ではなく \(p_{rgh}\) を使用することの数値計算上の利点に関する説明につなげられたらと考えています.では!