File: backwardDdtScheme.c
Fun
File: backwardDdtScheme.c
Function: template<class> tmp<geometricfield<type,> > backwardDdtScheme<type>::fvcDdt ( const volScalarField& rho, const GeometricField<type,>& vf ) { ...... rDeltaT.value()*( coefft*vf.boundaryField()  (coefft0*rho.oldTime().boundaryField() *vf.oldTime().boundaryField()  coefft00*rho.oldTime().oldTime().boundaryField() *vf.oldTime().oldTime().boundaryField() ) ) ...... coefft*vf.boundaryField() where is rho.boundaryField() 
Yes you are quite right rho.bo
Yes you are quite right rho.boundaryField() should be in coefft*vf.boundaryField() , i.e.
coefft*rho.boundaryField()*vf.boundaryField() Thanks for reporting it Henry 
Dr Henry Weller,
Thanks for
Dr Henry Weller,
Thanks for your answer. I am checking the source code. I am happy that I will find a bug in the foam. Maybe the foam will be bugfree in the future. Another thing. Could you kindly explain something I don't understand. See my post, http://www.cfdonline.com/cgibin/Op...cus/discus.cgi Liu Huafei 
Dr Henry Weller,
Thanks for y
Dr Henry Weller,
Thanks for your answer. Could you kindly explain something I don't understand. post, Some questions about the implementation of the convection scheme  By Liu Huafei on Sunday, January 20, 2008  05:18 pm: Edit Post It seems that something wrong with the implementation of the convection scheme. I guess the face value is based on the following formula Φf =ΦC +(xf – xC) *ψ(rf) * (dΦ/ dx)f Or Φf =ΦC +(xf – xC) / (xD– xC) *ψ(rf) * (ΦDΦC) f ,the face C,central node D downwind node ψ(rf) the flux limiter rf is based on the following formula rf=(ΦCΦU)/( ΦDΦC) = 2*grad(ΦC)*(xDxC) /( ΦDΦC) 1 so for the faceflux mf mf>0 Φf =ΦP+fx*ψ(rf) * (ΦNΦP) =(1fx*ψ(rf)) ΦP +f x*ψ(rf)* ΦN mf<0 Φf =ΦN+(1fx)*ψ(rf) * (ΦPΦN) =(1fx) *ψ(rf)* ΦP +[1(1fx) *ψ(rf)]* * ΦN fx is the geometric interpolation factor, defined by fx=(xNxf)/(xNxP) P is owner cell N is neighbour cell but in the limitedSurfaceInterpolationScheme.H and limitedSurfaceInterpolationScheme.c template<class> tmp<surfacescalarfield> limitedSurfaceInterpolationScheme<type>::weights ( const GeometricField<type,>& phi ) const { const fvMesh& mesh = this>mesh(); tmp<surfacescalarfield> tWeights(this>limiter(phi)); this isψ(rf)!!!!!!! surfaceScalarField& Weights = tWeights(); const surfaceScalarField& CDweights = mesh.surfaceInterpolation::weights(); fx !!!! scalarField& pWeights = Weights.internalField(); forAll(pWeights, face){ pWeights[face] = pWeights[face]*CDweights[face] + (1.0  pWeights[face])*pos(faceFlux_[face]); It is right?????? } Maybe openfoam is based on other formula. Could you kindly explain clearly it? 
I don't see the problem. The
I don't see the problem. The only thing in the code which may be confusing is the reuse of Weights, it is first used to store the limiter field and then becomes the weights field because after each faceweight is calculated the corresponding limiter value is no longer needed. I can add a comment on this if you think it would be useful.
Henry 
Thanks for your answer.
Thanks for your answer.

If convection scheme is done b
If convection scheme is done by such way,some contradiction arises.
LDUMatrix contains the coefficients, lowerPtr,upperPtr,diagPtr ,and the address reference, uPtr and lPtr. The upperPtr stores the contribution of the neighbour cell(N) to the owner cell(P) and the lowerPtr store the contribution of the owner(P) to the neighbour(N). For upwind scheme, Φf =max(mf,0)/mf * ΦP +min(mf,0)/mf * ΦN, so in upwind.h the weight factor is calculated by tmp<surfacescalarfield> weights() const { return pos(this>faceFlux_); }. It is OK. From this code,I guess the weight is corresponding to the following code fvm.lower() = weights.internalField()*faceFlux.internalField(); fvm.upper() = fvm.lower() + faceFlux.internalField(); fvm.negSumDiag(); fvm.lower is the contribution of the owner(P) to the neighbour(N),fvm.upper the contribution of the neighbour(N) to the owner(P). But for the high order convection scheme, mf>0 Φf =ΦP+fx*ψ(rf) * (ΦNΦP) = (1fx*ψ(rf)) ΦP +f x*ψ(rf)* ΦN mf<0 Φf =ΦN+(1fx)*ψ(rf) * (ΦPΦN) = (1fx) *ψ(rf)* ΦP +[1(1fx) *ψ(rf)]* ΦN the weight factor is calculated by pWeights[face] = pWeights[face]*CDweights[face] + (1.0  pWeights[face])*pos(faceFlux_[face]); It is right? It is corresponding to the fvm.upper rather than fvm.lower,why? fvm.lower is not the contribution of P to N? Another question: Diag[i] should be substract the contribution from all its neighbouring cell. In fvMatrix, the implementation of function negSumDiag() is as follows negSumDiag: for (register label face=0; face<l.size(); face++) { Diag[l[face]] = Lower[face]; // AP[IP]=AP[IP]AL[IFF] is correct?????? Diag[u[face]] = Upper[face]; // AP[IN]=AP[IN]AU[IFF] is correct?????? } why not the following code? for (register label face=0; face<l.size(); face++) { Diag[l[face]] = Upper[face]; (Upper is the contribution of N to P) Diag[u[face]] = Lower[face]; (Lower is the contribution of P to N) } Liu Huafei 
All times are GMT 4. The time now is 04:39. 