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

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
归一化因子 uu = \sum \left({|Ax - A\bar{x}}|+ |{b - A\bar{x}} |\right)\bar xx 的平均值。
归一化后的残差表达式为

\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_;
}