发布一个基于PLIC-VOF方法的两相流求解器，详见 https://github.com/daidezhi/interPlicFoam

@izumi 关于数值计算和CFD中的松弛技术，写了一点，对你的算例有用没用凑合看看吧

@飞火流星jyj 不好意思，这两天太忙没顾得上，刚刚把代码发上来。链接地址：meshInfoFoam，它会将体单元的体积、中心和面单元的面积、中心、面积矢量、单位矢量全部写入初始文件夹0，具体结果可以参考链接meshInfoFoam output。代码里还测试了一些涉及单个网格元素操作的成员函数。

@飞火流星jyj

    //class: fvMesh
    Info << "\n-Class: fvMesh--------" << endl;

    //- Return the object registry - resolve conflict polyMesh/lduMesh.
        //  Type: virtual const objectRegistry &
    Info << mesh.thisDb() << endl;

    //- Return reference to name.
        //  Type: const word &
    Info << mesh.name() << endl;

    //- Return reference to boundary mesh.
        //  Type: const fvBoundaryMesh &
    mesh.boundary();

    //- Internal face owner.
        // Type: const labelUList &
    labelList owners(mesh.owner());

    //- Internal face neighbour.
        //  Type: const labelUList &
    labelList neighbours(mesh.neighbour());

    //- Return cell volumes.
        //  Type: const DimensionedField< scalar, volMesh > &
    Info << mesh.V() << endl;

    //- Return old-time cell volumes.
        //  Type: const DimensionedField< scalar, volMesh > &
    Info << mesh.V0() << endl;

    //- Return old-old-time cell volumes.
        //  Type: const DimensionedField< scalar, volMesh > &
    Info << mesh.V00() << endl;

    //- Return sub-cycle cell volumes.
        //  Type: tmp< DimensionedField< scalar, volMesh > >
    Info << mesh.Vsc() << endl;

    //- Return sub-cycle old-time cell volumes.
        //  Type: tmp< DimensionedField< scalar, volMesh > >
    Info << mesh.Vsc0() << endl;

    //- Return cell face area vectors.
        //  Type: const surfaceVectorField &
    Info << mesh.Sf() << endl;

    //- Return cell face area magnitudes.
        //  Type: const surfaceScalarField &
    Info << mesh.magSf() << endl;

    //- Return cell face motion fluxes.
        //  Type: const surfaceScalarField &
    Info << mesh.phi() << endl;

    //- Return cell centres as volVectorField.
        //  Type: const volVectorField &
    Info << mesh.C() << endl;

    //- Return face centres as surfaceVectorField.
        //  Type: const surfaceVectorField &
    Info << mesh.Cf() << endl;

    //- Return face deltas as surfaceVectorField.
        //  Type: tmp< surfaceVectorField >
    Info << mesh.delta() << endl;  

    Info << "----------------------\n" << endl;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //class: fvMesh

    pointField points(mesh.points());
    faceList faces(mesh.faces());
    cellList cells(mesh.cells());    
    


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //Class: point = vector
    Info << "\n-Class: point---------" << endl;
    point &pt(points[3]);
    Info << "pt = " << pt << endl;
    //- Return x component
    Info << "pt.x() = " << pt.x() << endl;
    //- Return y component
    Info << "pt.y() = " << pt.y() << endl;
    //- Return z component
    Info << "pt.z() = " << pt.z() << endl;
    Info << "----------------------\n" << endl;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //class: edge
    Info << "\n-Class: edges----------" << endl;
    edge eg(faces[1].faceEdge(1)), eg_(faces[1].faceEdge(3));
    label pt_1(eg.start()), &pt_2(eg.end());
    Info << "eg = " << eg << endl;
    //- Return start vertex label
    Info << "eg.start() = " << pt_1 << endl;
    //- Return end vertex label
    Info << "eg.end() = " << pt_2 << endl;
    //- Given one vertex, return the other
    Info << "eg.otherVertex(eg.end()) = " << eg.otherVertex(pt_2) << endl;
    //- Return common vertex
        //  - -1: no common vertex
    Info << "eg.commonVertex(eg_) = " << eg.commonVertex(eg_) << endl;
    //- Return reverse edge
    Info << "eg.reverseEdge() = " << eg.reverseEdge() << endl;
    //- Return centre (centroid)
    Info << "eg.centre(points) = " << eg.centre(points) << endl;
    //- Return the vector (end - start)
    Info << "eg.vec(points) = " << eg.vec(points) << endl;
    //- Return scalar magnitude
    Info << "eg.mag(points) = " << eg.mag(points) << endl;
    //- Return edge line
    Info << "eg.line(points) = " << eg.line(points) << endl;
    //- compare edges
        //  Returns:
        //  -  0: different
        //  - +1: identical
        //  - -1: same edge, but different orientation
    Info << "Foam::edge::compare(eg, eg.reverseEdge()) = " 
         << Foam::edge::compare(eg, eg.reverseEdge()) << endl;
    Info << "----------------------\n" << endl;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //class: face
    Info << "\n-Class: face----------" << endl;
    face &fe(faces[4]);
    Info << "fe = " << fe << endl;
    //- Return true if the face is empty
    Info << "fe.empty() = " << fe.empty() << endl;
    //- Return No. of points corresponding to this face
    Info << "fe.size() = " << fe.size() << endl;
    //- Return first point
    Info << "fe.first() = " << fe.first() << endl;
    //- Return last point
    Info << "fe.last() = " << fe.last() << endl;
    //- Return n-th point
    Info << "fe.operator[](0) = " << fe.operator[](0) << endl;
    //- Return the points corresponding to this face
    Info << "fe.points(points) = " << fe.points(points) << endl;
    //- Centre point of face
    Info << "fe.centre(points) = " << fe.centre(points) << endl;
    //- Calculate average value at centroid of face
    Info << "fe.average(points, points) = " 
         << fe.average(points, points) << endl;
    //- Magnitude of face area
    Info << "fe.mag(points) = " << fe.mag(points) << endl;
    //- Vector normal; magnitude is equal to area of face
    Info << "fe.normal(points) = " << fe.normal(points) << endl;
    //- Return face with reverse direction
        //  The starting points of the original and reverse face are identical.
    Info << "fe.reverseFace() = " << fe.reverseFace() << endl;
    //- Which vertex on face (face index given a global index)
        //  returns -1 if not found
    Info << "fe.which(1966) = " << fe.which(1966) << endl;
    //- Next vertex on face
    Info << "fe.nextLabel(1) = " << fe.nextLabel(1) << endl;
    //- Previous vertex on face
    Info << "fe.prevLabel(1) = " << fe.prevLabel(1) << endl;
    //- Return number of edges
    Info << "fe.nEdges() = " << fe.nEdges() << endl;
    //- Return edges in face point ordering,
        //  i.e. edges()[0] is edge between [0] and [1]
    Info << "fe.edges() = " << fe.edges() << endl;
    //- Return n-th face edge
    Info << "fe.faceEdge(1) = " << fe.faceEdge(1) << endl;
    //- Return the edge direction on the face
        //  Returns:
        //  -  0: edge not found on the face
        //  - +1: forward (counter-clockwise) on the face
        //  - -1: reverse (clockwise) on the face
    Info << "fe.edgeDirection(fe.faceEdge(1)) = " 
         << fe.edgeDirection(fe.faceEdge(1)) << endl;
    //- compare faces
        //   0: different
        //  +1: identical
        //  -1: same face, but different orientation
    Info << "Foam::face::compare(fe, fe) = " 
         << Foam::face::compare(fe, fe) << endl;
    Info << "----------------------\n" << endl;


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    //class: cell
    Info << "\n-Class: cell----------" << endl;
    cell &cl(cells[100]);
    Info << "cl = " << cl << endl;
    //- Return true if the cell is empty
    Info << "cl.empty() = " << cl.empty() << endl;
    //- Return No. of faces corresponding to this cell
    Info << "cl.size() = " << cl.size() << endl;
    //- Return first face
    Info << "cl.first() = " << cl.first() << endl;
    //- Return last face
    Info << "cl.last() = " << cl.last() << endl;
    //- Return n-th face
    Info << "cl.operator[](0) = " << cl.operator[](0) << endl;
    //- Return number of faces
    Info << "cl.nFaces() = " << cl.nFaces() << endl;
    //- Return labels of cell vertices
    Info << "cl.labels(faces) = " << cl.labels(faces) << endl;
    //- Return the cell vertices
    Info << "cl.points(faces, points) = " << cl.points(faces, points) << endl;
    //- Return cell edges
    Info << "cl.edges(faces) = " << cl.edges(faces) << endl;
    //- Returns cell centre
    Info << "cl.centre(points, faces) = " << cl.centre(points, faces) << endl;
    //- Returns cell volume
    Info << "cl.mag(points, faces) = " << cl.mag(points, faces) << endl;
    Info << "----------------------\n" << endl;

1. 对流项div(rhoPhi, U)和div(phi, nuTilda)使用高阶格式；
2. 网格，做一下网格尺寸无关性实验。

@yhdthu
不可压流动中驱动流体流动的是压力梯度，一般对流场中的一个特定位置设定参考位置和参考压力值。如果包含第一类压力边界条件，使用p_rgh进行设置更为方便。而且多相流中，使用p_rgh能使得计算更加稳定和准确，非得使用p的话，就得像fluent使用诸如body-force-weighted等特殊的压力插值格式计算压力梯度了。在不同的流体交界面附近，压力p是连续的，但是压力梯度是非连续的，因为两边的流体密度不一样。
参考高度的不同只会影响p的值，不会影响其梯度的值，因为它们之间相差一个梯度为零的等值压力场。
最后，体积力项包含在p和p_rgh的转换关系里，详细参考东岳流体文档 http://dyfluid.com/interFoam.html ，公式(22)-(25)。

@yhdthu
连续方程

动量插值

较新版本的Tecplot (例如2015版本)已经支持直接读取OpenFOAM计算结果。

if (Pstream::parRun())
{
    if (Pstream::master())
    {
        Info << "balabala..." << endl;
    }
}
else
{
     Info << "balabala..." << endl;
}

这样串行并行就都OK了，可能需要#include "Pstream.H"。希望能帮到你。

@金石为开 SIMPLE算法本身是为不可压流动中速度-压力的耦合问题设计出来的，这里温度场的求解没有这类耦合问题，用simple只是为了达到按时间推进和设置非正交修正的目的，在源代码里你把while (simple.loop())改成while (runTime.loop())或者while (runTime.run())效果都是一样的。如果不需要非正交修正，则只需把

while (simple.correctNonOrthogonal())
{
    solve
    (
        fvm::ddt(T) - fvm::laplacian(DT, T)     
    );
}

改成

solve
(
    fvm::ddt(T) - fvm::laplacian(DT, T)     
);

此时我们就没有引入simple的必要了，所以这里使用simple只是为了方便设置非正交修正次数罢了。
这里有一个算例，温度场开始随时间剧烈变化，最后趋于稳定
OpenFOAM求解器开发

http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.83.867&rep=rep1&type=pdf

https://www-users.cs.umn.edu/~saad/IterMethBook_2ndEd.pdf

@kisssinger 是用的虚拟机？还没遇到过你这种错误，蛮奇怪的

线程数不对吧，看你的例子不是应该是4么

@炎冰灬灬 再等两天，上周四单位突然加活了 [大哭]

@炎冰灬灬 肯定的啊

@炎冰灬灬 我这周做个demo教程，过了今天就不怎么忙了

@Hungryandfool

1. 代码搞完了，还在测试，最近很忙，每天下班只想睡觉。
2. 几何类VOF方法：
   a. IsoAdvector第一次提出是在2016年，我在OpenFOAM基础上开发PLIC-VOF始于14年下半年，不过开发进度很慢，那时刚开始学习OpenFOAM。有一个完成度很高的2D版本，大家有兴趣的话我会把这个代码整理整理，供大佬们批评指正，几何运算是基于CGAL库的；
   b. PLIC和IsoAdvector是两类不同的几何重构方法，具体区别在于如何表征混合单元内的interface。PLIC使用平面，IsoAdvector使用等值面。确定interface位置时PLIC求解平面方程$\vec{n} \vec{X} + D = 0$里的距离$D$，IsoAdvector求解等值面阈值。界面重构精度应该差不多，求解效率还没比较，现在很难说谁高谁低；
   c. 我们求解器的公开版本advection部分直接借用了isoAdvection的，做了一些改动适配PLIC方法。主要是因为懒，不想单独开发advection部分了，不过这也符合开源精神吧，只做自己需要做和感兴趣的部分；
   d. 目前PLIC方法的一个研究热点是多面体任意单元里interface的快速定位技术。这方面研究主要集中在卡塔赫纳理工大学的Dr. Lopez [López, Joaquín, et al. "A new volume conservation enforcement method for PLIC reconstruction in general convex grids." Journal of Computational Physics 316 (2016): 338-359] (他们去年开发了适配凹多面体的快速算法，不过我感觉这个学术意义更大一些，基本上网格生成都是避免凹单元的出现的)，LANL的Dr. Diot [Diot, Steven, and Marianne M. François. "An interface reconstruction method based on an analytical formula for 3D arbitrary convex cells." Journal of Computational Physics 305 (2016): 63-74] 和Dr. Shashkov [Dyadechko, Vadim, and Mikhail Shashkov. "Moment-of-fluid interface reconstruction." Los Alamos Report LA-UR-05-7571 (2005)] (MOF方法本质还是PLIC方法，不过在$\vec{n} \vec{X} + D = 0$的求解上都引入了全新概念)。
3. 会有算例教程，但是事情多，我尽量快一些，但是这些事只有我一个人在弄，所以会比较慢。

PS：我不大名鼎鼎哦，只是小卒，只想靠着CFD混口饭吃，学界和业界都太他妈卷了

@东岳

@nanxuan PLIC应该也不行。本来有计划把这个求解器适配所有OF版本的，但是精力实在有限，进度一直不咋滴。