2015-04-12

fieldAverage 使用说明

fieldAverage是 OpenFOAM 中的一种 functionObject，用来计算时均值。其基本用法是作为一个 function object 放在 controlDict 文件中，运行 solver 的同时计算指定场的时均值，以下是一个示例：

functions
{
    fieldAverage1
    {
        type            fieldAverage;
        functionObjectLibs ( "libfieldFunctionObjects.so" );
        outputControl   outputTime;
        fields
        (
            Ua
            {
                 mean        on;
                 prime2Mean  off;
                 base        time;
            }
            Ub
            {
                 mean        on;
                 prime2Mean  off;
                 base        time;
            }
        );
    }
}

但是，在有效地使用fieldAverage之前，有一些问题需要澄清，其中最重要的一个是：fieldAverage到底在对哪个时间段求的时间平均？本文通过对fieldAverage的源码进行分析，试图厘清这些细节，并给出一个可靠的fieldAverage使用说明。涉及到的源文件包括：fieldAverage.H, fieldAverage.C, fieldAverageTemplates.C, fieldAverageItem.C，位于目录$WM_PROJECT_DIR/src/postProcessing/functionObjects/field/fieldAverage下。

首要问题：`fieldAverage`对哪个时间区间进行时均计算？

上面这段代码是从我运行的一个算例的 controlDict 中摘出来的。算例运行过程中发现，每一个时间步都要调用计算时均值相关的代码，并屏幕上的输出也会打印出Calculating averages，而每一个输出数据文件里，都有*Mean数据，其中*表示在controlDict中定义了需要求解时均值的某个场。
经过摸索发现，跟时均值求解密切相关的一段代码位于fieldAverageTemplates.C，下面是从中摘出来的一段核心代码：

182 template<class Type>
183 void Foam::fieldAverage::calculateMeanFieldType(const label fieldI) const
184 {
185  const word& fieldName = faItems_[fieldI].fieldName();
187  if (obr_.foundObject<Type>(fieldName))
188  {
189    const Type& baseField = obr_.lookupObject<Type>(fieldName);
190 
191    Type& meanField = const_cast<Type&>
192    (
193    obr_.lookupObject<Type>(faItems_[fieldI].meanFieldName())
194    );
195 
196    scalar dt = obr_.time().deltaTValue();
197    scalar Dt = totalTime_[fieldI];
198 
199    if (faItems_[fieldI].iterBase())
200    {
201      dt = 1.0;
202      Dt = scalar(totalIter_[fieldI]);
203    }
204 
205    scalar alpha = (Dt - dt)/Dt;
206    scalar beta = dt/Dt;
207 
208    if (faItems_[fieldI].window() > 0)
209    {
210      const scalar w = faItems_[fieldI].window();
211 
212      if (Dt - dt >= w)
213      {
214        alpha = (w - dt)/w;
215        beta = dt/w;
216      }
217    }
219    meanField = alpha*meanField + beta*baseField;
220  }
221 }

以下是我根据理解整理的代码解析：

变量baseField定义为当前时间步计算得到的场的值(Ua)，meanField定义为上一个时间步对应的场的时均值(UaMean)。
dt定义为时间步长，Dt定义为当前所在的时间。注意199行，iterBase()这个函数定义为当base (见controlDict) 为ITER时，返回true。此时，dt定义为1，意义很显然，当以时间步数为基准的时候，下一步和上一步当然是只相差1个时间步；Dt则被定义为当前所在的时间步(第xx步)。所以，Dt和dt的真实含义取决于base。fieldAverageItem.C中给base赋了默认值为ITER。base可以指定为time，具体见下文。
alpha 定义为(Dt-dt)/Dt, beta 定义为 dt/Dt。然后当前时间步的meanField定义为alpha*meanField + beta*baseField。举例说明，假设dt = 1, Dt = 20, 假设上一时间步的时均值为meanField_pre，则新时刻时均值为
1
19/20*meanField_pre + 1/20*baseField) = (19*meanField_pre + baseField)/20
19*meanField_pre的含义很明显，前19步的时均值乘以总时间步数，也就是前 19 步 Field 值的加和，再加上当前时间步的值baseField，就是前 20 步的 Field 值的和。用这个值除以20，得到新时刻的时均值meanField。由此可见，在这种情况下，每个时间步输出的时均值是从开始时刻到当前时间步的Field值对总时间的平均。
很多时候，我们并不想得到从开始到结束整个模拟时段内的平均值，而希望得到指定时段内的平均值，这时，可以通过指定window变量的值来达到目的。代码208行，注意到这里出现了一个新函数window。可以查到window函数返回值默认是-1。但可以在controlDict中base下面定义变量windows的值，当window函数返回值大于0时，208-217将被执行。接着上一点的假设，dt = 1, 定义并初始化为window函数的返回值，假设为20。则：
- 当 Dt 小于或等于20时，不满足Dt - dt >= w，此时时均值将按照上一点叙述的进行计算。
- 当 Dt大于20以后，假设Dt = 21，此时 alpha = (w-dt)/w = (20-1)/20，bata = dt/w = 1/19。令meanField_pre为Dt = 20 时的时均值，则Dt = 21时的时均值
  1
  meanField = 19/20 * meanField_pre + 1/20 * baseField = (19 * meanField_pre + baseField)/20
对比发现，如果window = -1(默认值)，则Dt = 21的时均值应该是20/21 * meanField_pre + 1/21 * baseField。而定义了windows = 20以后，此处的19 * meanField_pre 可以理解为Dt = 21这一步之前的19步的Field值的和，加上baseField，则为Dt = 21 以及其之前19步，共20步的Field值的和。再除以20，则为 Dt = 21以及其之前19步这20步的时均值。所以，当设定window为一正整数w时，输出的时均值是当前时间步以及其之前w-1步，这w步内Field的时均值。

附上一个更一般化的示例：

fieldAverage1
{
    type fieldAverage;
    functionObjectLibs ("libfieldFunctionObjects.so");
    resetOnRestart true;
    resetOnOutput false;
    startTime     0.1; // 开始计算时均值的时间
    endTime       1.0; // 停止计算时均值的时间
    outputControl   outputTime;
    fields
    (
        U
        {
            mean            on;
            prime2Mean      on;
            base            time; //以物理时间为基础来计算平均，而不是时间步数。
            window          10.0;
            windowName      w1; //optional
        }
        p
        {
            mean            on;
            prime2Mean      on;
            base            time;
        }
    );
}

resetOnRestart的值决定当solver继续运行时，是否要读取最近一个时间步的meanField的值来计算接下来时刻的时均值；resetOnOutput，顾名思义，是否要在每一次输出到文件以后重置meanField的值。这两个开关的默认值都是false。
mean这个开关的含义无需多言，计算公式如下：
$$ \overline {x} = \frac{1}{N} \sum \limits _{i=1}^{N} x_i $$
prime2Mean的计算公式如下：
$$ \overline{x’ ^{\,2}} = \frac{1}{N}\displaystyle\sum\limits_{i=1}^N (x_i - \overline{x})^2 $$
所以，如果prime2Mean为on，mean必须为on。

此外，如果计算已经结束，controlDict 中定义的 function 仍可以用execFlowFunctionObjects来执行。只是，这样的运行只能利用输出到文件的数据来进行计算了。举例说，假如时间步长是 0.001s，每 0.1s 输出一次，那么同样是1-2s的时均值，solver运行过程中求解的公式是:

1	( field.at(1.001) +field.at(1.002) + ... + field.at(2.0) )/1000

而solver运行完以后利用execFlowFunctionObjects计算的时均值应该是：

1	( field.at(1.1) + field.at(1.2) + ... + field.at(2) )/10

注意，这个结果没有经过直接的验证，是我根据原理推演的结果。
有时候，运行execFlowFunctionObjects会报错说找不到phi，这时加上-noFlow选项，就不会报错了。

小结

base用来指定作时间平均的基础，是基于时间步数(ITER)还是物理时间(time);
window用来作平均的时间段的长度，如果不设定，则求的是从开始到当前时间这个时间段的平均值。window的数值的实际含义依base而定，如果base是ITER，则window=20表示当前步及其前 19 个时间步从 20 个时间步内的平均，而如果base是time,则表示的是 20s 内的平均。

错误更正：
2015.09.05: 在 controlDict 里， base 只能是 time 或 iteration，而不能是上文中说的 ITER， ITER 只是代码中使用的一个代称。感谢“启之” 指出！

补充:
2015.11.26：上文没有提 prime2Mean 的计算方法，今天仔细看了一下，记录如下。
prime2Mean 的公式上面提到了，计算代码需要结合三个个函数来看，分别是 fieldAverage.C 文件里的 calcAverages 函数（事实上，上文提到的计算时均值的函数 calculateMeanFieldType 也是在这个函数中调用的）以及 fieldAverageTemplates.C 中定义的 calculatePrime2MeanFields 和 addMeanSqrToPrime2Mean 两个函数。我们看 calcAverages 看起

void Foam::fieldAverage::calcAverages()
{
   ......
   ......
    Info<< "    Calculating averages" << nl;

    addMeanSqrToPrime2Mean<scalar, scalar>();
    addMeanSqrToPrime2Mean<vector, symmTensor>();

    calculateMeanFields<scalar>();
    calculateMeanFields<vector>();
    calculateMeanFields<sphericalTensor>();
    calculateMeanFields<symmTensor>();
    calculateMeanFields<tensor>();

    calculatePrime2MeanFields<scalar, scalar>();
    calculatePrime2MeanFields<vector, symmTensor>();
  ......
  ......
}

可以看到，计算时均值 mean ，其实只要调用 calculateMeanFields 函数就可以了，但是计算 prime2Mean 则分成两部分：分别是在计算 mean 之前调用 addMeanSqrToPrime2Mean ，以及在计算 mean 之后调用 calculatePrime2MeanFields。下面一一来看这两个函数，注意到 addMeanSqrToPrime2Mean 和 calculatePrime2MeanFields 两个函数其实分别是通过调用 addMeanSqrToPrime2MeanType 和 calculatePrime2MeanFieldType 来其作用的，所以这里只看真正执行计算任务的addMeanSqrToPrime2MeanType 和 calculatePrime2MeanFieldType 函数。

addMeanSqrToPrime2MeanType

template<class Type1, class Type2>
void Foam::fieldAverage::addMeanSqrToPrime2MeanType(const label fieldI) const
{
    const word& fieldName = faItems_[fieldI].fieldName();

    if (obr_.foundObject<Type1>(fieldName))
    {
        const Type1& meanField =
            obr_.lookupObject<Type1>(faItems_[fieldI].meanFieldName());

        Type2& prime2MeanField = const_cast<Type2&>
        (
            obr_.lookupObject<Type2>(faItems_[fieldI].prime2MeanFieldName())
        );

        prime2MeanField += sqr(meanField);
    }
}

这个函数很简单，就是直接往上一步的 prime2MeanField 加上上一步的 meanField 。

calculatePrime2MeanFieldType

template<class Type1, class Type2>
void Foam::fieldAverage::calculatePrime2MeanFieldType(const label fieldI) const
{
    const word& fieldName = faItems_[fieldI].fieldName();

    if (obr_.foundObject<Type1>(fieldName))
    {
        const Type1& baseField = obr_.lookupObject<Type1>(fieldName);
        const Type1& meanField =
            obr_.lookupObject<Type1>(faItems_[fieldI].meanFieldName());

        Type2& prime2MeanField = const_cast<Type2&>
        (
            obr_.lookupObject<Type2>(faItems_[fieldI].prime2MeanFieldName())
        );

        scalar dt = obr_.time().deltaTValue();
        scalar Dt = totalTime_[fieldI];

        if (faItems_[fieldI].iterBase())
        {
            dt = 1.0;
            Dt = scalar(totalIter_[fieldI]);
        }

        scalar alpha = (Dt - dt)/Dt;
        scalar beta = dt/Dt;

        if (faItems_[fieldI].window() > 0)
        {
            const scalar w = faItems_[fieldI].window();

            if (Dt - dt >= w)
            {
                alpha = (w - dt)/w;
                beta = dt/w;
            }
        }

        prime2MeanField =
            alpha*prime2MeanField
          + beta*sqr(baseField)
          - sqr(meanField);
    }
}

这个函数里， alpha 与 beta 的定义跟上文计算时均值的函数 calculateMeanFieldType 是一样的，关键在于最后的那一句。
那么，这两个函数结合起来，如何就等价于上文给出的 prime2MeanField 的计算公式了呢？下面做一个简单的分析：
t = 1 时，某个场 $x$ 记为 $x_1$ ， meanField 记为 $m_1$，等于 $x_1$ ， prime2MeanField 记为 $pM_1$，等于 0；

t = 2 时，此刻的 $x$ 记为 $x_2$ ， meanField 记为 $m_2$，等于 $(x_1+x_2)/2$， prime2MeanField 的理论值等于 $\frac{1}{2} [(x_1-m_2)^2+(x_2-m_2)^2]$，而按照代码，此刻的 prime2MeanField （记为$pM_2$）的计算分两步，先是计算当前步的 meanField 之前，执行

$$
prime2MeanField = pM_1 + m_1^2
$$

然后是在计算当前步的 meanField 之后，执行

$$
prime2MeanField = \frac{1}{2}\cdot(prime2MeanField) + \frac{1}{2}\cdot x_2^2 - m_2^2
$$

结合起来，得到

$$
\begin{align}
pM_2 = & \frac{1}{2}\cdot (pM_1+m_1^2) + \frac{1}{2}\cdot x_2^2 - m_2^2 \\
= & \frac{1}{2}\cdot m_1^2 + \frac{1}{2}\cdot x_2^2 - m_2^2 \\
= & \frac{1}{2}\cdot x_1^2 + \frac{1}{2}\cdot x_2^2 - m_2^2
\end{align}
$$

由于 $2\cdot m_2 = x_1 + x_2$ ，于是，上式可以进一步化简，

$$
\begin{align}
pM_2 = & \frac{1}{2}\cdot x_1^2 + \frac{1}{2}\cdot x_2^2 - 2\cdot m_2^2 + m_2^2 \\
= & \frac{1}{2}\cdot x_1^2 + \frac{1}{2}\cdot x_2^2 -(x_1+x_2)\cdot m_2 + m_2^2 \\
= & \frac{1}{2}\cdot (x_1^2 -2\cdot x_1^2 + m_2^2) + \frac{1}{2}\cdot (x_2^2 -2\cdot x_2^2 + m_2^2) = \frac{1}{2}\cdot (x_1-m_2)^2 + \frac{1}{2}\cdot (x_2-m_2)^2
\end{align}
$$

这就与理论值刚好对上了；

t = 3时，也可以类推下去，这里简略写一下：
$$
\begin{align}
pM_3 = & \frac{2}{3}\cdot (pM_2 + m_2^2) + \frac{1}{3}\cdot x_3^2 - m_3^2 \\
= & \frac{1}{3}\cdot (x_1^2 + x_2^2 + x_3^2) - m_3^2 \\
= & \frac{1}{3}\cdot (x_1^2 + x_2^2 + x_3^2) - \frac{2}{3}\cdot (x_1+x_2+x_3)\cdot m_3 + m_3^2 \\
= & \frac{1}{3}\cdot [(x_1-m_3)^2 + (x_2-m_3)^2 + (x_3-m_3)^2]
\end{align}
$$

同样是跟理论值是一样的。

Giskard's CFD Learning Tricks

CFD and Scientific Computing

fieldAverage 使用说明

首要问题：`fieldAverage`对哪个时间区间进行时均计算？

小结

首要问题：fieldAverage对哪个时间区间进行时均计算？

小结

首要问题：`fieldAverage`对哪个时间区间进行时均计算？