Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    solveSegragated的问题?

    OpenFOAM
    1
    2
    1747
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • 程
      程迪 last edited by

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

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

      github: chengdi123000
      网站:chengdi123000.github.io
      本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

      1 Reply Last reply Reply Quote
      • 程
        程迪 last edited by

        哦,

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

        github: chengdi123000
        网站:chengdi123000.github.io
        本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。

        1 Reply Last reply Reply Quote
        • First post
          Last post

        CFD中文网 | 东岳流体 | 京ICP备15017992号-2
        论坛登录问题反馈可联系 li.dy@dyfluid.com