Code of the Week #9: tolerance 1e-06

如需留言，请移步 CFD Forum

接 Code of the Week #8: solver PCG，这一次来讲讲容差的设置。这个设置位于 system/fvSolution ，部分如下：

solvers                       // 求解矩阵的算法设置
{
    p                         // 求解变量 p
    {
        solver          PCG;  // 求解的方法，preconditioned conjugate gradient
        preconditioner  DIC;  // 预处理方法，
        tolerance       1e-06;// 绝对误差
        relTol          0.1;  // 相对误差
    }
}
...

在 OpenFOAM 里面，误差 (tolerance) 分为两种，一种是绝对误差 (tolerance)，一种是相对误差 (relTol)。如上的一段示例设置表示：收敛的条件是，绝对残差小于 1e-6，相对残差小于 0.1。以下解释设未知量为 $x$ ，需要求解的方程是 $Ax=b$

绝对残差： 当前步的残差 $\delta^{n}$ ， $\delta^{n} =Ax^{n+1} - b$ ， 归一化因子 $u$ ， $u = \sum \left({|Ax - A\bar{x}}|+ |{b - A\bar{x}} |\right)$ ， $\bar x$ 为 $x$ 的平均值。 归一化后的残差表达式为 $$ \delta^n_u = \frac{1}{u} \sum |b-Ax| $$ 需要注意的是，计算的方法可能跟 solver 有关系，不同的 solver 可能计算有差别。但是总体方法是一致的。 相对残差： 当前步的残差与第一步残差的比值， $\delta_u^n\%=\dfrac{||\delta_u^n||}{||\delta_u^0||}$ 。

通过上述的条件设定，收敛的条件就是 $$ ||\delta^{n}_u|| \leqslant \text{tolerance} $$ 或 $$ ||\delta^{n}_u\%|| \leqslant \text{relTol} $$ 如果 relTol 设置为 0 ，则强制执行 tolerance 设定的条件。

以矩阵求解器 CG (位于 src/lduSolvers/lduSolver/cgSolver.C) 为例，其 solve 函数为（暂不对程序做详细注解）：

Foam::lduSolverPerformance Foam::cgSolver::solve
(
    scalarField& x, 
    const scalarField& b, 
    const direction cmpt
) const
{
    // Prepare solver performance
    lduSolverPerformance solverPerf(typeName, fieldName());

    scalarField wA(x.size());
    scalarField rA(x.size());

    // Calculate initial residual
    matrix_.Amul(wA, x, coupleBouCoeffs_, interfaces_, cmpt);

    // Use rA as scratch space when calculating the normalisation factor
    scalar normFactor = this->normFactor(x, b, wA, rA, cmpt);  // 归一化因子

    if (lduMatrix::debug >= 2)
    {
        Info<< "   Normalisation factor = " << normFactor << endl;
    }

    // Calculate residual
    forAll (rA, i)
    {
        rA[i] = b[i] - wA[i];
    }

    solverPerf.initialResidual() = gSumMag(rA)/normFactor;  // 进行归一操作
    solverPerf.finalResidual() = solverPerf.initialResidual();

    if (!stop(solverPerf))
    {
        scalar rho = matrix_.great_;
        scalar rhoOld = rho;

        scalar alpha, beta, wApA;

        scalarField pA(x.size(), 0);

        do
        {
            rhoOld = rho;

            // Execute preconditioning
            preconPtr_->precondition(wA, rA, cmpt);

            // Update search directions
            rho = gSumProd(wA, rA);

            beta = rho/rhoOld;

            forAll (pA, i)
            {
                pA[i] = wA[i] + beta*pA[i];
            }


            // Update preconditioned residual
            matrix_.Amul(wA, pA, coupleBouCoeffs_, interfaces_, cmpt);

            wApA = gSumProd(wA, pA);


            // Check for singularity
            if (solverPerf.checkSingularity(mag(wApA)/normFactor))
            {
                break;
            }

            // Update preconditioned residual
            matrix_.Amul(wA, pA, coupleBouCoeffs_, interfaces_, cmpt);

            wApA = gSumProd(wA, pA);


            // Check for singularity
            if (solverPerf.checkSingularity(mag(wApA)/normFactor))
            {
                break;
            }

            // Update solution and residual
            alpha = rho/wApA;

            forAll (x, i)
            {
                x[i] += alpha*pA[i];
            }

            forAll (rA, i)
            {
                rA[i] -= alpha*wA[i];
            }

            solverPerf.finalResidual() = gSumMag(rA)/normFactor; // 进行归一操作
            solverPerf.nIterations()++;
        } while (!stop(solverPerf));
    }

    return solverPerf;
}

归一化因子的计算位于 src/foam/matrices/lduMatrix/lduMatrix/lduMatrixSolve.C：

Foam::scalar Foam::lduMatrix::solver::normFactor
(
    const scalarField& x,
    const scalarField& b,
    const scalarField& Ax,
    scalarField& tmpField,
    const direction cmpt
) const
{
    // Calculate A dot reference value of x
//     matrix_.sumA(tmpField, coupleBouCoeffs_, interfaces_);
//     tmpField *= gAverage(x);

    // Calculate normalisation factor using full multiplication
    // with mean value.  HJ, 5/Nov/2007
    scalarField xRef(x.size(), gAverage(x));

    // Eliminated equations are removed from residual normalisation
    if (!matrix_.eliminatedEqns().empty())
    {
        labelList elim = matrix_.eliminatedEqns().toc();

        forAll (elim, elimI)
        {
            // Set the value of xRef to be identical to the value of x
            // to eliminate the residual
            xRef[elim[elimI]] = x[elim[elimI]];
        }
    }

    matrix_.Amul
    (
        tmpField,
        xRef,
        coupleBouCoeffs_,
        interfaces_,
        cmpt
    );

    return gSum(mag(Ax - tmpField) + mag(b - tmpField)) + matrix_.small_;

    // At convergence this simpler method is equivalent to the above
    // return 2*gSumMag(b) + matrix_.small_;
}

如需留言，请移步 CFD Forum