solveSegragated的问题？

[程迪]

fvMatrix<T>::solveSegragated()：

template<class Type>
Foam::SolverPerformance<Type> Foam::fvMatrix<Type>::solveSegregated
(
    const dictionary& solverControls
)
{
    if (debug)
    {
        Info.masterStream(this->mesh().comm())
            << "fvMatrix<Type>::solveSegregated"
               "(const dictionary& solverControls) : "
               "solving fvMatrix<Type>"
            << endl;
    }

    GeometricField<Type, fvPatchField, volMesh>& psi =
       const_cast<GeometricField<Type, fvPatchField, volMesh>&>(psi_);

    SolverPerformance<Type> solverPerfVec
    (
        "fvMatrix<Type>::solveSegregated",
        psi.name()
    );
//原始的diag存起来。因为diag是lduMatrix的私有变量，类型是scalarField，fvMatrix<vector>的diag也是scalarField。
    scalarField saveDiag(diag());
//复制source，Field<T>的copy constructor语义是复制
    Field<Type> source(source_);

    // At this point include the boundary source from the coupled boundaries.
    // This is corrected for the implict part by updateMatrixInterfaces within
    // the component loop.
//源项中加入边界邻接单元的值×邻接单元的矩阵系数。源项是vector，所以是一次性更新了所有分量的源项。
    addBoundarySource(source);

//validComponents
//Return a labelType of valid component indicators.
//1 : valid (solved) -1 : invalid (not solved)
//对称张量之类的可以少解很多分量，只求解需要求解的分量即可。
    typename Type::labelType validComponents
    (
        psi.mesh().template validComponents<Type>()
    );	

    for (direction cmpt=0; cmpt<Type::nComponents; cmpt++)
    {
        if (validComponents[cmpt] == -1) continue; //求解未求解的分量。

        //状态向量初值的cmpt分量，复制语义
        scalarField psiCmpt(psi.primitiveField().component(cmpt));
       //？？？
       //耦合边界的内部系数internalCoeffs_中的cmpt分量加入到对角系数中。
        addBoundaryDiag(diag(), cmpt); 

        //源项，复制语义。
        scalarField sourceCmpt(source.component(cmpt));

       //耦合边界的外侧单元系数的分量，复制语义。
        FieldField<Field, scalar> bouCoeffsCmpt
        (
            boundaryCoeffs_.component(cmpt) //component返回tmp
        );

        FieldField<Field, scalar> intCoeffsCmpt
        (
            internalCoeffs_.component(cmpt)
        );

        //换个马甲
       // lduInterface负责处理耦合边界
        lduInterfaceFieldPtrsList interfaces =
            psi.boundaryField().scalarInterfaces();

        // Use the initMatrixInterfaces and updateMatrixInterfaces to correct
        // bouCoeffsCmpt for the explicit part of the coupled boundary
        // conditions
        initMatrixInterfaces
        (
            true,
            bouCoeffsCmpt,
            interfaces,
            psiCmpt,
            sourceCmpt, //这是输出，变的还是源项。
            cmpt
        );

        updateMatrixInterfaces
        (
            true,
            bouCoeffsCmpt,
            interfaces,
            psiCmpt,
            sourceCmpt,//这是输出，变的还是源项。
            cmpt
        );

        solverPerformance solverPerf;

        // Solver call
        solverPerf = lduMatrix::solver::New
        (
            psi.name() + pTraits<Type>::componentNames[cmpt],
            *this, //this->upper_, lower_ 其实没有变过，只有diag_变了。
            bouCoeffsCmpt,//耦合边界的边界系数分量没变过
            intCoeffsCmpt, //耦合边界的内部系数分量没变过
            interfaces, //interface还是那些，没有变过
            solverControls //每个分量的控制
        )->solve(psiCmpt, sourceCmpt, cmpt); //初值没变过，源项变了。
// 总的来说，非对角项从来都没变过。对角项和源项变了。
        if (SolverPerformance<Type>::debug)
        {
            solverPerf.print(Info.masterStream(this->mesh().comm()));
        }

        solverPerfVec.replace(cmpt, solverPerf);
        solverPerfVec.solverName() = solverPerf.solverName();
       
//复制回去
        psi.primitiveFieldRef().replace(cmpt, psiCmpt);
        diag() = saveDiag;//恢复diag
    }
    //L,U,耦合边界的内外单元系数从来没变过，只有源项和对角项有调整。
//实际每次解的是 (L+U+D_new_i )*x_i = s_new_i
// s_new_i = s_i + C_[n->i]_i*x0_i 
// s_new_i会加上初始状态向量中，耦合边界外部单元对内部单元的贡献的i分量。
//D_new_i = D_new_[i-1] + C_[i->n]_i 
//D_new_i 每次计算都会叠加一次。
 // 计算前准备源项和对角项时，都是按分量单独准备的，已经求解的变量对需要添加的部分无影响。

    psi.correctBoundaryConditions();//更新边界条件

    psi.mesh().setSolverPerformance(psi.name(), solverPerfVec);

    return solverPerfVec;
}

fvMatrix<T>继承了lduMatrix，lduMatrix是纯scalar的，只存了L,D,U,调用lduMatrix::solve() 的时候要传入内外耦合边界的矩阵系数和初值以及源项才能求解。而内外耦合边界系数是存在fvMatrix<T>下的FieldField<T> internalCoeffs_和boundaryCoeffs_。所以我觉得从矩阵系数的分布来看，其实是没有分量间的耦合的，v的求解应该对u一点儿影响没有。

但是从实现来看那个对角系数其实是在不断变化的，这点就非常诡异了。

[程迪]

哦，

 diag() = saveDiag;//恢复diag

是在括号之前的，所以每次的修改不会叠加。

之前理解错误的一点是所有的边界有关的系数都是放在internalCoeffs和boundaryCoeffs中的，不仅仅是coupled BC。从初始化时候的size可以看出来。

template<class Type>
Foam::fvMatrix<Type>::fvMatrix
(
    const GeometricField<Type, fvPatchField, volMesh>& psi,
    const dimensionSet& ds
)
:
    lduMatrix(psi.mesh()),
    psi_(psi),
    dimensions_(ds),
    source_(psi.size(), Zero),
    internalCoeffs_(psi.mesh().boundary().size()), //和边界的face数大小一样。
    boundaryCoeffs_(psi.mesh().boundary().size()),
    faceFluxCorrectionPtr_(nullptr)
{
//...

不过非耦合的应该是internalCoeffs放比例系数，boundaryCoeffs放源项，而耦合的BC是internalCoeffs放owner的系数，boundaryCoeffs放的另一侧单元的系数。从addBoundarySource()的源代码可以看出。

template<class Type>
void Foam::fvMatrix<Type>::addBoundarySource
(
    Field<Type>& source,
    const bool couples
) const
{
    forAll(psi_.boundaryField(), patchi)
    {
        const fvPatchField<Type>& ptf = psi_.boundaryField()[patchi];
        const Field<Type>& pbc = boundaryCoeffs_[patchi];

        if (!ptf.coupled())//非耦合边界，只有owner，另一侧没有单元，pbc存的是边界源项。
        {
            addToInternalField(lduAddr().patchAddr(patchi), pbc, source);
        }
        else if (couples)//耦合边界，如果couples==true, 另一侧有单元，pbc存的是另一侧单元的系数。
        {
            const tmp<Field<Type>> tpnf = ptf.patchNeighbourField();
            const Field<Type>& pnf = tpnf();

            const labelUList& addr = lduAddr().patchAddr(patchi);

            forAll(addr, facei)
            {
                source[addr[facei]] += cmptMultiply(pbc[facei], pnf[facei]);
            }
        }
    }
}