solveSegragated的问题?


  • OpenFOAM讲师

    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一点儿影响没有。

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


  • OpenFOAM讲师

    哦,

     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]);
                }
            }
        }
    }