# compressibleFoam ＆求解器编写思路探讨

• 各位：
最近在研究这个求解器：
compressibleFoam，发现里面的求解思路和我之前采用fortran语言是一样的。思路就是按边来循环，分别对internalFace和boundaryFace进行处理，求得单元交界面上的通量。
控制方程如下

其中一些量如下

下面这个是求内部面无粘通量的函数

/// Loop over internal faces and get face flux from L and R cell center
label leftCell , rightCell;
forAll( mesh.owner() , iface ) {
/// Get the left and right cell index
leftCell = mesh.owner()[iface];
rightCell = mesh.neighbour()[iface];
/// Approximate Riemann solver at interface
scalar lambda = (*fluxSolver)( &rho[leftCell]  , &U[leftCell]  , &p[leftCell] , /// Left
&rho[rightCell] , &U[rightCell] , &p[rightCell] ,          /// Right
&massFlux[iface] , &momFlux[iface] , &energyFlux[iface] ,
&nf[iface] );
/// Multiply with face area to get face flux
massFlux[iface] *= mesh.magSf()[iface];
momFlux[iface] *= mesh.magSf()[iface];
energyFlux[iface] *= mesh.magSf()[iface];
localDt[ leftCell ]  += lambda * mesh.magSf()[iface];
localDt[ rightCell ] += lambda * mesh.magSf()[iface];
}


下面这个是求物面边界条件的边界面上的无粘通量的函数

\*---------------------------------------------------------------------------*/
/////////////////////////////////////////////////////////////////////////////////
/// Remember BC's for Gas Dynamics solvers
/// are based on characteristics and cannot be applied
/// to individual variables. Rather it should be applied
/// to the state vector as a whole. OpenFOAM framework
/// does not allow us to do this without breaking the
/// convention.
/////////////////////////////////////////////////////////////////////////////////
/// Loop over boundary patches and determine the
/// type of boundary condition of the patch
/////////////////////////////////////////////////////////////////////////////////
forAll ( mesh.boundaryMesh() , ipatch ) {
scalar bflux[5];
word BCTypePhysical = mesh.boundaryMesh().physicalTypes()[ipatch];
word BCType         = mesh.boundaryMesh().types()[ipatch];
word BCName         = mesh.boundaryMesh().names()[ipatch];
const UList<label> &bfaceCells = mesh.boundaryMesh()[ipatch].faceCells();
/////////////////////////////////////////////////////////////////////////////////
/// Slip Wall BC ///
/////////////////////////////////////////////////////////////////////////////////
if( BCTypePhysical == "slip" || BCTypePhysical == "symmetry" ) {
forAll( bfaceCells  , iface ) {
/// Extrapolate wall pressure
scalar p_e = p[ bfaceCells[iface] ];
vector normal = nf.boundaryField()[ipatch][iface];
scalar face_area = mesh.magSf().boundaryField()[ipatch][iface];
scalar lambda = std::fabs( U[ bfaceCells[iface] ] & normal ) +
std::sqrt ( gama * p_e / rho[ bfaceCells[iface] ] );
momResidue[ bfaceCells[iface] ] += p_e * normal * face_area;
localDt[ bfaceCells[iface] ]    += lambda * face_area;
}
}


之后我们将边界上的通量从owner单元减去，加到nighbour单元去

   /// Clear all residues
massResidue = dimensionedScalar( "", massResidue.dimensions() , 0.0 );
momResidue = dimensionedVector( "" , momResidue.dimensions() , vector( 0.0 , 0.0 , 0.0 ) );
energyResidue = dimensionedScalar( "" , energyResidue.dimensions() ,  0.0 );
/// Loop over each face and add the flux to left and subtract from the right
forAll( mesh.owner() , iface ) {
/// Store left and right cell reference
leftCell = mesh.owner()[iface];
rightCell = mesh.neighbour()[iface];
/// Note that the normal vector to a face will point
/// from the owner cell to the neighbour cell (L to R)
massResidue[leftCell]    += massFlux[iface];
momResidue[leftCell]     += momFlux[iface];
energyResidue[leftCell]  += energyFlux[iface];
/// Subtract from right cell
massResidue[rightCell]   -= massFlux[iface];
momResidue[rightCell]    -= momFlux[iface];
energyResidue[rightCell] -= energyFlux[iface];
}


原来编程就是一直按照这种思路来的，可是转到OpenFOAM 下一时适应不过来。看到官方的求解器，根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗？
大家编写自己的求解器的时候都是怎么实现的？希望看透这一切联系的牛人能给予一点提示和点拨。跪谢！

• 看到官方的求解器，根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗？

粗略看了下，这个求解器完全没用openfoam的思路来，采用的他自己的思路，就像你说的，他在开发这个求解器的时候，可能是有其他的代码（比如同组的Fortran）转过来的。并且完全是面向过程的思想。

如果你要用OpenFOAM求解的思路，还是看OpenFOAM的求解器比较好他这个就是把求解的过程自己实现了，然而没有使用OpenFOAM自带的方法。

他的main函数已经写的很清楚了：

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "setInputValues.H"

#include "createFields.H"

/// Time step loop
/// Posts the non-blocking send/recv of fields
long int iter = 0;
while( runTime.loop() ) {

/// 构造通量
#include "constructFaceFlux.H"

/// 有限体积离散
#include "sumFlux.H"

/// 边界通量修正
#include "boundaryFlux.H"

Info << "Iteration = " << ++iter << "  ";

/// 矩阵计算
#include "stateUpdateLTS.H"

Info << " Max residue = " << rhoResidMax << endl;
/// 输出结果
runTime.write();
}

return 0;
}