updateCoeffs与evaluate

[麦迪文]

各位同学在自定义边界条件的时候可能会发现有的边界条件用的是updateCoeffs()，另外一些是evaluate(),那么这两个函数到底有什么不同？选取两个典型的边界条件讨论一下:zeroGradient和fixedFluxPressure两种边界条件

使用QT4.8调试icofoam，case为cavity，压力边界条件分别给出fixedFluxPressure和zeroGradient

调试过程中发现求解器在构造矩阵之前调用updateCoeffs，好比这里使用了 void constrainPressure( volScalarField& p, const volVectorField& U, const surfaceScalarField& phiHbyA, const RAUType& rAU); 调用了压力P的Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs(const scalarField& snGradp )

在求解之后，调用evaluate，在这里使用pEqn.solve()函数分别调用了Foam::fixedGradientFvPatchField<Type>::evaluate(const Pstream::commsTypes)和 Foam::zeroGradientFvPatchField<Type>::evaluate，其中Foam::fixedGradientFvPatchField<Type>::evaluate是跟Foam::fixedFluxPressureFvPatchScalarField::updateCoeffs绑定在一起

所以在自定义边界条件时要注意好他们之间的逻辑，至于openfoam为什么要这样设计，答案只能在评论里了，楼主也猜不透。。。。。。实际开此帖的目的是欢迎研究边界条件的同学关注一下我上一个求助贴。。。。。地址(http://www.cfd-china.com/topic/2543/求助-液池热毛细对流算例-已编译好求解器和算例-但就是和文献自编程序不一样)，请踊跃提出自己的看法

[李东岳]

边界条件一般通过updateCoeffs做计算，比如fixedFluxPressure：

if (!this->updated())
    {
        this->updateCoeffs();
    }

    Field<Type>::operator=
    (
        this->patchInternalField() + gradient_/this->patch().deltaCoeffs()
    );

    fvPatchField<Type>::evaluate();

opetaror是关键步骤，关键步骤之前是提升系数。不知道你是凑巧遇到这个fixedFluxPressure还是偶然，这个边界条件在求解pEqn之前需要计算一些东西，比如面法向梯度，于是在求解器里面每个时间布就要通过updateCoeffs()计算面法向梯度。你自己写了个面法向梯度的代码? 我这面见到的是求解器顶层植入的：

// Update the fixedFluxPressure BCs to ensure flux consistency
    setSnGrad<fixedFluxPressureFvPatchScalarField>
    (
        p_rgh.boundaryFieldRef(),
        (
            phiHbyA.boundaryField()
          - (
                alphaf1.boundaryField()*phi1.boundaryField()
              + alphaf2.boundaryField()*phi2.boundaryField()
            )
        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
    );

[麦迪文]

@麦迪文 在 updateCoeffs与evaluate 中说：

updateCoeffs(const scalarField& snGradp )

是的是的，我用的5.0版本，所以它的梯度是下面这样定义的（/***/部分），我没有写面法向梯度代码，只是浅析了下OF的代码，这个的确是在evaluate()之前先用updateCoeffs()修正了下系数，学到很多，谢谢老师

 template<class RhoType, class RAUType, class MRFType>
 void Foam::constrainPressure
 (
     volScalarField& p,
     const RhoType& rho,
     const volVectorField& U,
     const surfaceScalarField& phiHbyA,
     const RAUType& rhorAU,
     const MRFType& MRF
 )
 {
     const fvMesh& mesh = p.mesh();
 
     volScalarField::Boundary& pBf = p.boundaryFieldRef();
 
     const volVectorField::Boundary& UBf = U.boundaryField();
     const surfaceScalarField::Boundary& phiHbyABf =
         phiHbyA.boundaryField();
     const typename RAUType::Boundary& rhorAUBf =
         rhorAU.boundaryField();
     const surfaceVectorField::Boundary& SfBf =
         mesh.Sf().boundaryField();
     const surfaceScalarField::Boundary& magSfBf =
         mesh.magSf().boundaryField();
 
     forAll(pBf, patchi)
     {
         if (isA<fixedFluxPressureFvPatchScalarField>(pBf[patchi]))
         {
             refCast<fixedFluxPressureFvPatchScalarField>
             (
                 pBf[patchi]
             ).updateCoeffs
             (
                 (
/***********************************************************************
                     phiHbyABf[patchi]
                   - rho.boundaryField()[patchi]
                    *MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
/***********************************************************************
                 )
                /(magSfBf[patchi]*rhorAUBf[patchi])
             );
         }
     }
 }
 

[麦迪文]

             (
                 (
/***********************************************************************
                     phiHbyABf[patchi]
                   - rho.boundaryField()[patchi]
                    *MRF.relative(SfBf[patchi] & UBf[patchi], patchi)
                 )
                /(magSfBf[patchi]*rhorAUBf[patchi])
/***********************************************************************
             );

这个才是。。。