Please refresh the page if equations are not rendered correctly.
---------------------------------------------------------------

场数据操作基础

cpp中，函数必须先定义后使用。但是为了提高代码的可读性，可以只在程序的开头对函数进行声明，如以下代码块的第7行，函数的具体定义则放在文件的最后。

#include "fvCFD.H"

// This is a function declaration; this method will calculate some scalar value
// given the current time, location in space x, and a reference point x0. The
// function also accepts a scaling factor, scale.
// The actual implementation, or definition, is given at the end of this code.
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale);

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

}

1. 读取字典文件（transportProperties为例）

    Info << "Reading transportProperties\n" << endl;

    IOdictionary transportProperties
    (
        IOobject
        (
            "transportProperties", // name of the dictionary
            runTime.constant(), // location in the case - this one is in constant
            mesh, // needs the mesh object reference to do some voodoo - unimportant now
            IOobject::MUST_READ_IF_MODIFIED, // the file will be re-read if it gets modified during time stepping
            IOobject::NO_WRITE // read-only
        )
    );

    // Create a scalar constant for kinematic viscosity by reading the value from the dictionary.
    dimensionedScalar nu  // 定义有量纲标量
    (
        "nu", // name of the variable
        dimViscosity,     // 定义量纲
        transportProperties.lookup("nu") // takes the value from the dictionary and returns it, passing it to the object constructor as an argument
    );

以上代码块的第20行（dimViscosity）定义的是动力粘度的单位。可以在终端中输入以下命令查看它是如何定义的：

grep -r dimViscosity $FOAM_SRC/OpenFOAM/

你应该能得到下面的返回值：

/opt/openfoam8/src/OpenFOAM/dimensionSet/dimensionSets.C: const dimensionSet dimViscosity(dimArea/dimTime);
/opt/openfoam8/src/OpenFOAM/dimensionSet/dimensionSets.C: const dimensionSet dimDynamicViscosity(dimDensity*dimViscosity);
/opt/openfoam8/src/OpenFOAM/dimensionSet/dimensionSets.H: extern const dimensionSet dimViscosity;

返回值中给出dimViscosity动力粘度的量纲是面积/时间（m^2/s）。在dimensionSets.C文件中，我们可以进一步查看长度、时间、面积及其组合量纲是如何定义的：

        const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0);
        const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0);
        const dimensionSet dimArea(sqr(dimLength));
        const dimensionSet dimViscosity(dimArea/dimTime);*/

当然，与其使用dimViscosity，你也可以使用以下代码直接定义量纲：

        // This is what gets used here. But, an alternative would be to type in the units directly:
        dimensionSet(0,2,-1,0,0,0,0),

2. 读取场数据（以时间步文件夹中的压力p和速度U文件为例）

这一节的代码以时间步文件夹中的压力p和速度U文件为例，展示了数据的场读取过程，具体读取的文件所在时间步文件夹由system/controlDict (i.e. latestTime, startTime, etc.)指定。在OpenFOAM中，速度、压力等场数据是存储在单元中心的，因此变量初始化是采用的是vol<Type>Field。

These read the fields p and U from the time folders, as specified in system/controlDict (i.e. latestTime, startTime, etc.)

    Info<< "Reading field p\n" << endl;
    volScalarField p // note that pressure is a scalar field
    (
        IOobject
        (
            "p", // name of the field
            runTime.timeName(), // current time name, i.e. the time folder to read from
            mesh,
            IOobject::MUST_READ, // always gets imported, will throw an error if the field is missing
            IOobject::AUTO_WRITE // will get saved automatically when the controlDict parameters will request it
        ),
        mesh // initialises the field to match the size of the mesh with default (0) values
    );

读取速度的方法与读取压力场数据一致，不同的是速度为矢量，因此声明时应选用volVectorField而不是volScalarField。

    Info<< "Reading field U\n" << endl;
    volVectorField U // note that velocity is a vector field
    (
        IOobject
        (
            "U",
            runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh
    );

3

让我们首先定义一个向量作为计算的原点（originVector），其值在程序执行过程中不会改变。然后我们计算距离原点最远处的那个网格单元中心点到原点的距离。具体实现时，首先计算所有网格的中心点（采用mesh.C() 获取）到原点的距离，然后取最大值。

Let us define a vector whose values will not change over the course of the program execution.

下面代码块的第6行使用max()方法取一个数组的最大值。使用dimensionedVector("x0",dimLength,originVector)定义一个量纲为长度的矢量。使用mag()计算一组矢量中的每个矢量的长度。这里有一个疑问，为什么要使用.value()方法将距离转换为无量纲数呢？

    const vector originVector(0.05,0.05,0.005);

    // Calculate the distance from the origin to the cell centre furthest away.
    // In Python, this is equivalent to: np.sqrt(np.sum((x0-x)**2))
    // The .value() method is called to convert from a dimensionedScalar to a regular scalar.
    const scalar rFarCell = max(
        mag(dimensionedVector("x0",dimLength,originVector)-mesh.C())  //
        ).value(); // convert to dim-less scalar

从以下代码开始，将根据模拟的需要进行时间步进循环，直到触发终止条件跳出循环。

This part of the code performs time stepping loops as long as required by the simulation.

    Info<< "\nStarting time loop\n" << endl;

    // This will increment the current time automatically
    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        // Loop over all cells in the mesh and calculate the pressure value.
        for (label cellI=0; cellI<mesh.C().size(); cellI++)
        {
            // cellI describes a series of integers, each corresponding to an index of an individual cell in the grid.

            // Call the method and compute p.
            // Note how mesh.C() and p elements are accessed by using the [] operators, as in a regular C array.
            // .value() is also called to convert the time to a dim-less scalar
            p[cellI] = calculatePressure(runTime.time().value(), mesh.C()[cellI], originVector, rFarCell);
        }

实际编程中，对所有单元逐个遍历并计算会严重降低代码的执行效率。更好的办法是直接对数组进行操作，正如前面计算距离原点最远的单元中心到原点的距离rFarCell时所作的那样。为了完整性和说明某些基本特征，这里也演示了这种迭代方法，但一般来说，除非绝对必要，否则不鼓励采用这种方法。

另外，我们也可以操作边界面的场数据，会在后面的教程中讲到。 一旦我们有了压力场数据p，就可以使用fvc::grad(p)计算压力梯度，进一步得到速度场U了。此处需要注意的是，压力梯度grad(p)的单位不是m/s。所以我们乘以一个单位时间变量dimensionedScalar("tmp",dimTime,1.)，使之成为如此。当然，这只是为了说明这个问题。

NOTE: a more appropriate way to calculate things in OpenFOAM is through performing operations on field objects and not iterating cell-by-cell, where possible. How to do this has been shown above, where rFarCell is being computed. The iterative approach has been presented for completeness and to illustrate certain basic features, but is, generally, discouraged, unless absolutely necessary.

NOTE: it is also possible to interact with boundary face values, but this will be addressed in a separate tutorial. Next we calculate the gradient of p and substitute for U. Note that the units of grad(p) are not m/s, so we multiply by a unit time variable to make them so. This is just for illustration purposes, of course.

        // The result will point either towards or away from the origin, depending on the sign of the time varying "pressure".
        U = fvc::grad(p)*dimensionedScalar("tmp",dimTime,1.);

        // If requested by controlDict, save the fields to files.
        runTime.write();

3. 压力函数

下面函数中，首先将单元中心点\vec{x}到原点\vec{x_0}的距离使用前文计算得到的最远端距离rFarCell进行归一化，得到取值在[0,1]区间内的半径值r。然后获取半径的倒数。最终，压力采用以下函数求得：
p=\frac{sin \left( 2 \pi t \right)}{r}
其中，t取值为runTime.time().value()。

// definition of the custom function
scalar calculatePressure(scalar t, vector x, vector x0, scalar scale)
{
    // Calculates the distance between the base point and x, which is given by the magnitude (hence the name mag)
    // of the vector between x0 and x. The value is scaled by the passed factor, with the intention of making
    // it vary between 0 and 1.
    scalar r (mag(x-x0)/scale);

    // Calculate the inverse of r and apply a limiter to avoid dividing by zero.
    scalar rR (1./(r+1e-12));

    // definition of a frequency
    scalar f (1.);

    // Return a sinusoidally varying pressure with maximum at x0.
    // Note how we call the OpenFOAM sin method by referring to the Foam namespace.
    // This is done to differentiate between the native cpp implementation of a method with the same name
    // and thus avoid an ambiguous expression.
    return Foam::sin(2.*Foam::constant::mathematical::pi*f*t)*rR;
}

到此为止，场数据的基本操作指令讲解完了。你可以使用ParaView对压力场p和速度场U进行可视化。

Best to visualize the results by plotting p iso-contours with range (-10,10) and applying a glyph filter to the U field in Paraview.

4. 案例文件下载