Factors for the cubic inerpolation scheme in OpenFOAM

Hoshang · March 10, 2017, 11:36

Dear all,

I have a question regarding the used factors for the cubic interpolation scheme in OpenFOAM. Below you can see a part of the code for the cubic interpolation scheme:

Code:

virtual bool corrected() const
{
    return true;
}

//- Return the explicit correction to the face-interpolate
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
correction
(
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
{
       const fvMesh& mesh = this->mesh();

       // calculate the appropriate interpolation factors
       const surfaceScalarField& lambda = mesh.weights();

       const surfaceScalarField kSc
       (
                 lambda*(scalar(1) - lambda*(scalar(3) - scalar(2)*lambda))    
       );

       const surfaceScalarField kVecP(sqr(scalar(1) - lambda)*lambda);

       const surfaceScalarField kVecN(sqr(lambda)*(lambda - scalar(1)));

       tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsfCorr
       (
             new GeometricField<Type, fvsPatchField, surfaceMesh>
             (
                    IOobject
                    (
                        "cubic::correction(" + vf.name() +')',
                        mesh.time().timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::NO_WRITE,
                        false
                    ),
                    surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc)
             )
       );
       GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr.ref();

       for (direction cmpt=0; cmpt<pTraits<Type>::nComponents;cmpt++)
       {
            sfCorr.replace
            (
                 cmpt,
                 sfCorr.component(cmpt)
             + (
                     surfaceInterpolationScheme
                     <
                          typename outerProduct
                          <
                                 vector,
                                 typename pTraits<Type>::cmptType
                          >::type
                      >::interpolate
                      (
                             fv::gaussGrad
                             <typename pTraits<Type>::cmptType>(mesh)
                            .grad(vf.component(cmpt)),
                             kVecP,
                             kVecN
                       ) & mesh.Sf()
                  )/mesh.magSf()/mesh.surfaceInterpolation::deltaCoeffs()
              );
         }

The cubic interpolation is done using the linear interpolated term plus a correction term which is defined above in the code. As I understand from the Code the correction uses the factor kSc for the owner point value and -kSc for the neighbor point value. In addition, the gradient of the owner point is added with the factor kVecP and the gradient of the neighbor point with the factor kVecN. With this we get for the value at the surface the following correlation:

phi_E=lambda*phi_P+(1-lambda)*phi_N+kSc*phi_P-kSc*phi_N+kVecP*gradPhi_P+kVecN*gradphi_N

Here phi_P and phi_N are the values at the owner and neighbor point and gradphi_P and gradphi_N the gradients at the mentioned points.
But by the use of a normal cubic polynomial function (f(x)=ax³+bx²+cx+d) for the interpolation of the surface value with given values and gradients for the owner and neighbor points as conditions I get another correlation with the defined factors in the above code:

phi_E=lambda*phi_P+(1-lambda)*phi_N-kSc*phi_P+kSc*phi_N-kVecN*gradPhi_P-kVecP*gradphi_N

The differences are the opposite signs and the exchange of the factors.
I have tried lot of other numerical interpolation ways to come on the same correlation as used in OpenFOAM but i wasnt sucessful till yet.
Does anyone knows how one gets the correlation which is used in OpenFOAM respectively which interpolation technique gives this correlation with the given values and gradients at the owner and neighbor point?

Thanks in advance for answers and hopefully someone can give me an advise regarding this problem.

sahas · March 28, 2017, 13:16

I have no deep look into the problem but I am also interesting in that. When lambda is equal to 0.5 both formulas give the same result. Could you post your derivation of the formula please? Of course if it does not cause you any inconvenience

chengyu · August 2, 2017, 08:26

Hey, guys. I have found the same problem as Hoshang met. I post my derivation of cubic interpolation here, if someone can find any errors, plz feel free to point out.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Given the value of the variable and its gradient at point P and N, assume a cubic polynomial variation in the section PN, derive the variable value at the face f.

$\phi(x) = ax^3+bx^2+cx+d$
$\phi(x=0)=d=\phi_P$
$\phi'(x=0)=c=\nabla \phi_P$
$\phi(x=1)=a+b+c+d=\phi_N$
$\phi'(x=1)=3a+2b+c=\nabla \phi_N$

thus a,b,c,d can be solved,

$a=2\phi_P - 2\phi_N + \nabla \phi_P + \nabla \phi_N$
$b= -3\phi_P + 3\phi_N - 2\nabla \phi_P - \nabla \phi_N$
$c=\nabla \phi_P,\quad d=\phi_P$

then the variation of $\phi$ in PN can be expressed as

$\phi(x) = (2\phi_P - 2\phi_N + \nabla \phi_P + \nabla \phi_N)x^3+ (-3\phi_P + 3\phi_N - 2\nabla \phi_P - \nabla \phi_N)x^2 + \nabla \phi_P * x + \phi_P$
$\phi(x) = \phi_P(2x^3-3x^2+1) + \phi_N(-2x^3+3x^2) + \nabla \phi_P(x^3-2x^2+x) + \nabla \phi_N(x^3-x^2)$
$\phi_f = \phi(1-\lambda) = \lambda \phi_P + (1-\lambda)\phi_N + \lambda (1-\lambda)(1-2\lambda)(\phi_N -\phi_P)$ $- \lambda^2(\lambda-1) \nabla \phi_P - \lambda (1-\lambda)^2 \nabla \phi_N$
$\phi_f = \phi_{f,linear} + kSc*\phi_N - kSc*\phi_P - kVecN*\nabla \phi_P - kVecP* \nabla \phi_N$

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The last equation is what exactly Hoshang gave, which is not the same implemented in OpenFOAM.

Best wishes
Yu Cheng

sahas · August 2, 2017, 10:54

Hello, Yu Cheng!
I have checked your derivation and found it to be correct, except one misprint in the formula for $\phi'(x=1)$ : it should be equal to $\nabla \phi_N$ .

I think one should report the bug (https://bugs.openfoam.org/bug_report_page.php).

chengyu · August 2, 2017, 22:03

Yes, you are right, I have changed the original formula, and I think I should better report this to the OpenFOAM official. Thank you.

Best wishes.

Yu Cheng

sahas · August 4, 2017, 12:37

Bug report: https://bugs.openfoam.org/view.php?id=2650

chengyu · August 6, 2017, 04:25

Hi Alexander,
Thank you for the report.

Best wishes,

Yu Cheng

lac · September 15, 2017, 12:05

Hello everyone!

This is my first post in this awsome forum. I hope this question would not be considered as off-topic, as it is related to the cubic scheme.
I started to use openFoam recently, and I'm sure, it is possible to find the answer from the code, but I can't.
Therefore I'd gratefull, if anyone could tell me how the gradients are calculated for the cubic scheme. I mean, are the schemes used from fvSchemes, or are they 'hard wired' in the code for this interpolation scheme.

Thanks in advance.

sahas · September 15, 2017, 15:20

Quote:

Originally Posted by lac

Therefore I'd gratefull, if anyone could tell me how the gradients are calculated for the cubic scheme. I mean, are the schemes used from fvSchemes, or are they 'hard wired' in the code for this interpolation scheme.

Hello!
The scheme for gradient evaluation is hard-coded for cubic interpolation and it is Gauss.

lac · September 16, 2017, 04:58

Quote:

Originally Posted by sahas

Hello!
The scheme for gradient evaluation is hard-coded for cubic interpolation and it is Gauss.

Thank you. I guess it is linear then?

What do you think, is it possible to change it to fourth for instance?
Is it enough to modify the lines:
#include "gaussGrad.H" to #include "fourthGrad.H"
and
>::interpolate ( fv::gaussGrad
to
>::interpolate ( fv::fourthGrad

?

Or use some limiting somehow?

sahas · September 16, 2017, 17:54

Quote:

Originally Posted by lac

Thank you. I guess it is linear then?

What do you think, is it possible to change it to fourth for instance?
Is it enough to modify the lines:
#include "gaussGrad.H" to #include "fourthGrad.H"
and
>::interpolate ( fv::gaussGrad
to
>::interpolate ( fv::fourthGrad

?

Or use some limiting somehow?

I think it is possible and quite simple. I tried to change the scheme to leastSquares and I had no problems.

lac · September 17, 2017, 07:33

Quote:

Originally Posted by sahas

I think it is possible and quite simple. I tried to change the scheme to leastSquares and I had no problems.

I also tried the leastSquares and it was easy. But the fourth is not working and I don't know how to limit the gradients.

sahas · September 18, 2017, 02:04

Quote:

Originally Posted by lac

I also tried the leastSquares and it was easy. But the fourth is not working and I don't know how to limit the gradients.

I don't look it in deep but I think that you may create your own scheme cubicGrad with runtime-selectable scheme for gradient as it is done for linearUpwind scheme.

lac · September 20, 2017, 04:08

Quote:

Originally Posted by sahas

I don't look it in deep but I think that you may create your own scheme cubicGrad with runtime-selectable scheme for gradient as it is done for linearUpwind scheme.

It's a good idea, although I don't know how to do it, yet. Anyway, thanks.

Tobi · September 22, 2017, 05:19

Hi Sahas,

I read your bug report and the comments. Interesting fact. Right now I could imagine that if you have shock-wave propagation, the results might be different based on the gradients of neighbor/owner which are way much more pronounced in such cases than normal ones. Which cases did you check out for your report?

lac · September 22, 2017, 05:51

Quote:

Originally Posted by Tobi

Hi Sahas,

I read your bug report and the comments. Interesting fact. Right now I could imagine that if you have shock-wave propagation, the results might be different based on the gradients of neighbor/owner which are way much more pronounced in such cases than normal ones. Which cases did you check out for your report?

I have also seen that report. However, I'm not completely sure about my derivation, but if you change how λ is defined above, from f(x=1-λ) to f(x=λ), the same coefficients as implemented in openFoam currently can be derived.

sahas · September 22, 2017, 06:52

Quote:

Originally Posted by Tobi

Hi Sahas,

I read your bug report and the comments. Interesting fact. Right now I could imagine that if you have shock-wave propagation, the results might be different based on the gradients of neighbor/owner which are way much more pronounced in such cases than normal ones. Which cases did you check out for your report?

Hi, Tobi,

I have tried a simple unsteady case of 1D flow with sinusoidal varying velocity at inlet, using uniform and non-uniform grids. The difference between realizations of the cubic scheme is observed only on non-uniform grid and it is quite small (switching to linear or linearUpwind scheme has greater effect).
As for shock-wave case, due to the fact that cubic scheme is central scheme it would be hard to get non-oscillating solution for this case. And again, differences in realizations can be only seen on non-uniform grids, where other numerical errors begin to play the role.
To summarize, the current implementation gives small errors in a final solution. Perhaps, the effect is strong in the some cases of DNS of turbulence. But DNS is not very suitable to be a simple test case

bany · March 30, 2024, 08:55

Quote:

Originally Posted by Hoshang

Dear all,

I have a question regarding the used factors for the cubic interpolation scheme in OpenFOAM. Below you can see a part of the code for the cubic interpolation scheme:

Code:

virtual bool corrected() const
{
    return true;
}

//- Return the explicit correction to the face-interpolate
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
correction
(
     const GeometricField<Type, fvPatchField, volMesh>& vf
 ) const
{
       const fvMesh& mesh = this->mesh();

       // calculate the appropriate interpolation factors
       const surfaceScalarField& lambda = mesh.weights();

       const surfaceScalarField kSc
       (
                 lambda*(scalar(1) - lambda*(scalar(3) - scalar(2)*lambda))    
       );

       const surfaceScalarField kVecP(sqr(scalar(1) - lambda)*lambda);

       const surfaceScalarField kVecN(sqr(lambda)*(lambda - scalar(1)));

       tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsfCorr
       (
             new GeometricField<Type, fvsPatchField, surfaceMesh>
             (
                    IOobject
                    (
                        "cubic::correction(" + vf.name() +')',
                        mesh.time().timeName(),
                        mesh,
                        IOobject::NO_READ,
                        IOobject::NO_WRITE,
                        false
                    ),
                    surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc)
             )
       );
       GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr.ref();

       for (direction cmpt=0; cmpt<pTraits<Type>::nComponents;cmpt++)
       {
            sfCorr.replace
            (
                 cmpt,
                 sfCorr.component(cmpt)
             + (
                     surfaceInterpolationScheme
                     <
                          typename outerProduct
                          <
                                 vector,
                                 typename pTraits<Type>::cmptType
                          >::type
                      >::interpolate
                      (
                             fv::gaussGrad
                             <typename pTraits<Type>::cmptType>(mesh)
                            .grad(vf.component(cmpt)),
                             kVecP,
                             kVecN
                       ) & mesh.Sf()
                  )/mesh.magSf()/mesh.surfaceInterpolation::deltaCoeffs()
              );
         }

The cubic interpolation is done using the linear interpolated term plus a correction term which is defined above in the code. As I understand from the Code the correction uses the factor kSc for the owner point value and -kSc for the neighbor point value. In addition, the gradient of the owner point is added with the factor kVecP and the gradient of the neighbor point with the factor kVecN. With this we get for the value at the surface the following correlation:

phi_E=lambda*phi_P+(1-lambda)*phi_N+kSc*phi_P-kSc*phi_N+kVecP*gradPhi_P+kVecN*gradphi_N

Here phi_P and phi_N are the values at the owner and neighbor point and gradphi_P and gradphi_N the gradients at the mentioned points.
But by the use of a normal cubic polynomial function (f(x)=ax³+bx²+cx+d) for the interpolation of the surface value with given values and gradients for the owner and neighbor points as conditions I get another correlation with the defined factors in the above code:

phi_E=lambda*phi_P+(1-lambda)*phi_N-kSc*phi_P+kSc*phi_N-kVecN*gradPhi_P-kVecP*gradphi_N

The differences are the opposite signs and the exchange of the factors.
I have tried lot of other numerical interpolation ways to come on the same correlation as used in OpenFOAM but i wasnt sucessful till yet.
Does anyone knows how one gets the correlation which is used in OpenFOAM respectively which interpolation technique gives this correlation with the given values and gradients at the owner and neighbor point?

Thanks in advance for answers and hopefully someone can give me an advise regarding this problem.

Hi, HoShang. I have a question how "the linear interpolated term plus a correction term" defined in code? I did not find the corresponding code.
Thus, I can only turn the code in cubic.H to the following form:
phi_E=kSc*phi_P-kSc*phi_N+kVecP*gradPhi_P+kVecN*gradphi_N

Any advices are appreciated.

Kaley2 · May 14, 2024, 04:33

I haven't delved deeply into https://greenleafguru.com the problem, but I'm also interested in it. When lambda equals 0.5, both formulas yield the same result. Could you please share https://herbalhighsociety.com your derivation of the formula? Of course, only if it's not inconvenient for you.

March 10, 2017, 11:36	Factors for the cubic inerpolation scheme in OpenFOAM	#1
Hoshang New Member Join Date: Mar 2017 Posts: 1 Rep Power: 0	Dear all, I have a question regarding the used factors for the cubic interpolation scheme in OpenFOAM. Below you can see a part of the code for the cubic interpolation scheme: Code: virtual bool corrected() const { return true; } //- Return the explicit correction to the face-interpolate virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> correction ( const GeometricField<Type, fvPatchField, volMesh>& vf ) const { const fvMesh& mesh = this->mesh(); // calculate the appropriate interpolation factors const surfaceScalarField& lambda = mesh.weights(); const surfaceScalarField kSc ( lambda(scalar(1) - lambda(scalar(3) - scalar(2)lambda)) ); const surfaceScalarField kVecP(sqr(scalar(1) - lambda)lambda); const surfaceScalarField kVecN(sqr(lambda)(lambda - scalar(1))); tmp<GeometricField<Type, fvsPatchField, surfaceMesh>> tsfCorr ( new GeometricField<Type, fvsPatchField, surfaceMesh> ( IOobject ( "cubic::correction(" + vf.name() +')', mesh.time().timeName(), mesh, IOobject::NO_READ, IOobject::NO_WRITE, false ), surfaceInterpolationScheme<Type>::interpolate(vf, kSc, -kSc) ) ); GeometricField<Type, fvsPatchField, surfaceMesh>& sfCorr = tsfCorr.ref(); for (direction cmpt=0; cmpt<pTraits<Type>::nComponents;cmpt++) { sfCorr.replace ( cmpt, sfCorr.component(cmpt) + ( surfaceInterpolationScheme < typename outerProduct < vector, typename pTraits<Type>::cmptType >::type >::interpolate ( fv::gaussGrad <typename pTraits<Type>::cmptType>(mesh) .grad(vf.component(cmpt)), kVecP, kVecN ) & mesh.Sf() )/mesh.magSf()/mesh.surfaceInterpolation::deltaCoeffs() ); } The cubic interpolation is done using the linear interpolated term plus a correction term which is defined above in the code. As I understand from the Code the correction uses the factor kSc for the owner point value and -kSc for the neighbor point value. In addition, the gradient of the owner point is added with the factor kVecP and the gradient of the neighbor point with the factor kVecN. With this we get for the value at the surface the following correlation: phi_E=lambdaphi_P+(1-lambda)phi_N+kScphi_P-kScphi_N+kVecPgradPhi_P+kVecNgradphi_N Here phi_P and phi_N are the values at the owner and neighbor point and gradphi_P and gradphi_N the gradients at the mentioned points. But by the use of a normal cubic polynomial function (f(x)=ax³+bx²+cx+d) for the interpolation of the surface value with given values and gradients for the owner and neighbor points as conditions I get another correlation with the defined factors in the above code: phi_E=lambdaphi_P+(1-lambda)phi_N-kScphi_P+kScphi_N-kVecNgradPhi_P-kVecP*gradphi_N The differences are the opposite signs and the exchange of the factors. I have tried lot of other numerical interpolation ways to come on the same correlation as used in OpenFOAM but i wasnt sucessful till yet. Does anyone knows how one gets the correlation which is used in OpenFOAM respectively which interpolation technique gives this correlation with the given values and gradients at the owner and neighbor point? Thanks in advance for answers and hopefully someone can give me an advise regarding this problem.

August 2, 2017, 08:26	Cubic scheme derivation	#3
chengyu New Member Yu Cheng Join Date: Aug 2014 Posts: 15 Rep Power: 12	Hey, guys. I have found the same problem as Hoshang met. I post my derivation of cubic interpolation here, if someone can find any errors, plz feel free to point out. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Given the value of the variable and its gradient at point P and N, assume a cubic polynomial variation in the section PN, derive the variable value at the face f. $\phi(x) = ax^3+bx^2+cx+d$ $\phi(x=0)=d=\phi_P$ $\phi'(x=0)=c=\nabla \phi_P$ $\phi(x=1)=a+b+c+d=\phi_N$ $\phi'(x=1)=3a+2b+c=\nabla \phi_N$ thus a,b,c,d can be solved, $a=2\phi_P - 2\phi_N + \nabla \phi_P + \nabla \phi_N$ $b= -3\phi_P + 3\phi_N - 2\nabla \phi_P - \nabla \phi_N$ $c=\nabla \phi_P,\quad d=\phi_P$ then the variation of $\phi$ in PN can be expressed as $\phi(x) = (2\phi_P - 2\phi_N + \nabla \phi_P + \nabla \phi_N)x^3+ (-3\phi_P + 3\phi_N - 2\nabla \phi_P - \nabla \phi_N)x^2 + \nabla \phi_P * x + \phi_P$ $\phi(x) = \phi_P(2x^3-3x^2+1) + \phi_N(-2x^3+3x^2) + \nabla \phi_P(x^3-2x^2+x) + \nabla \phi_N(x^3-x^2)$ $\phi_f = \phi(1-\lambda) = \lambda \phi_P + (1-\lambda)\phi_N + \lambda (1-\lambda)(1-2\lambda)(\phi_N -\phi_P)$ $- \lambda^2(\lambda-1) \nabla \phi_P - \lambda (1-\lambda)^2 \nabla \phi_N$ $\phi_f = \phi_{f,linear} + kSc\phi_N - kSc\phi_P - kVecN\nabla \phi_P - kVecP \nabla \phi_N$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The last equation is what exactly Hoshang gave, which is not the same implemented in OpenFOAM. Best wishes Yu Cheng Last edited by chengyu; August 2, 2017 at 22:00. Reason: correct misprint term pointed out by Alexander

August 2, 2017, 22:03	Change the mis-spelling	#5
chengyu New Member Yu Cheng Join Date: Aug 2014 Posts: 15 Rep Power: 12	Yes, you are right, I have changed the original formula, and I think I should better report this to the OpenFOAM official. Thank you. Best wishes. Yu Cheng

September 15, 2017, 12:05	Gradients in cubic	#8
lac New Member Join Date: Apr 2016 Posts: 12 Rep Power: 10	Hello everyone! This is my first post in this awsome forum. I hope this question would not be considered as off-topic, as it is related to the cubic scheme. I started to use openFoam recently, and I'm sure, it is possible to find the answer from the code, but I can't. Therefore I'd gratefull, if anyone could tell me how the gradients are calculated for the cubic scheme. I mean, are the schemes used from fvSchemes, or are they 'hard wired' in the code for this interpolation scheme. Thanks in advance.

September 22, 2017, 05:19		#15
Tobi Super Moderator Tobias Holzmann Join Date: Oct 2010 Location: Bad Wörishofen Posts: 2,711 Blog Entries: 6 Rep Power: 52	Hi Sahas, I read your bug report and the comments. Interesting fact. Right now I could imagine that if you have shock-wave propagation, the results might be different based on the gradients of neighbor/owner which are way much more pronounced in such cases than normal ones. Which cases did you check out for your report? __________________ Keep foaming, Tobias Holzmann

March 28, 2017, 13:16		#2
sahas Member Alexander Join Date: Mar 2009 Posts: 49 Rep Power: 17	I have no deep look into the problem but I am also interesting in that. When lambda is equal to 0.5 both formulas give the same result. Could you post your derivation of the formula please? Of course if it does not cause you any inconvenience

August 2, 2017, 10:54		#4
sahas Member Alexander Join Date: Mar 2009 Posts: 49 Rep Power: 17	Hello, Yu Cheng! I have checked your derivation and found it to be correct, except one misprint in the formula for $\phi'(x=1)$ : it should be equal to $\nabla \phi_N$ . I think one should report the bug (https://bugs.openfoam.org/bug_report_page.php).

August 4, 2017, 12:37		#6
sahas Member Alexander Join Date: Mar 2009 Posts: 49 Rep Power: 17	Bug report: https://bugs.openfoam.org/view.php?id=2650

August 6, 2017, 04:25		#7
chengyu New Member Yu Cheng Join Date: Aug 2014 Posts: 15 Rep Power: 12	Hi Alexander, Thank you for the report. Best wishes, Yu Cheng

May 14, 2024, 04:33		#19
Kaley2 New Member Lavinia Cummerata Join Date: Feb 2024 Posts: 4 Rep Power: 2	I haven't delved deeply into https://greenleafguru.com the problem, but I'm also interested in it. When lambda equals 0.5, both formulas yield the same result. Could you please share https://herbalhighsociety.com your derivation of the formula? Of course, only if it's not inconvenient for you. Last edited by Kaley2; May 21, 2024 at 02:53.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Getting Started with OpenFOAM	wyldckat	OpenFOAM	26	June 21, 2024 07:54
Map of the OpenFOAM Forum - Understanding where to post your questions!	wyldckat	OpenFOAM	10	September 2, 2021 06:29
Radiation With View Factors in OPENFOAM	MadiS	OpenFOAM	1	August 23, 2012 12:21
Second order upwinding scheme in OpenFoam	subash	OpenFOAM Running, Solving & CFD	4	June 20, 2012 19:28
Surface interpolation scheme with tensors as weighting factors	jutta	OpenFOAM Running, Solving & CFD	1	August 27, 2007 12:00