[JBUNSW]

Dear Foamers,

I am seeing some weird behaviors from fvc::grad(phi) operator under periodic boundary conditions.

To make sure that fvc::grad(phi) is working correctly, I setup a 1D domain with length 2mm divided into 400 equal length cells. so

delX = 2e-3 / 400

phi is a volScalarField obtained solving a simple transport equation. You can set some random numbers to build phi, doesn't matter actually.

Then I printed out the fvc::grad(phi) along with the phi at all 400 cells. Since the problem is 1D we don't have gradients in y or z direction, therefore, the gradient at internal points can be simply calculated by a central difference scheme as

grad(phi_i) = ( phi_i+1 - phi_i-1 ) / (2 * delX)

I copied the data into a spreadsheet and did the above. Below are the results I'm getting:

Cell 	   phi	        mag(grad(phi))	     abs(d(phi)/dx) 	Err %
				 (OpenFOAM)	      (manual calc)
----------------------------------------------------------------------------
1	244676	6.31092e-317			-		     -	
2	246466	6.38511e-317		358100000		#DIV/0!
3	248257	0				358300000		#DIV/0!
4	250049	0				358400000		#DIV/0!
5	251841	3.58454e+08		358400000		0.015
6	253633	3.5851e+08			358500000		0.003
7	255426	3.58566e+08		358600000		-0.010
8	257219	3.58622e+08		358600000		0.006
9	259012	3.58564e+08		358600000		-0.010
10	260805	3.58506e+08		358500000		0.002
11	262597	3.58447e+08		358400000		0.013
12	264389	3.58503e+08		358500000		0.001
13	266182	3.58444e+08		358500000		-0.016
14	267974	3.58271e+08		358300000		-0.008
15	269765	3.58098e+08		358100000		-0.001
16	271555	3.57924e+08			-		     -	
.	.		.					.		     .
.	.		.					.		     .
.	.		.					.		     .
392	259012	3.58564e+08			-		     -	
393	257219	3.58622e+08		358600000		0.006
394	255426	3.58566e+08		358600000		-0.009
395	253633	3.5851e+08			358500000		0.003
396	251841	3.58454e+08		358400000		0.015
397	250049	3.58398e+08		358400000		-0.001
398	248257	3.58228e+08		358300000		-0.020
399	246466	3.58058e+08		358100000		-0.012
400	244676	1.79043e+08			-		    -

I am using periodic boundary condition in OpenFOAM 2.0.x. Cell 1 is residing at x=0 and cell 400 ends at x = 2mm. As you can see for all internal nodes safely away from the boundaries the error is negligible. But when you get close to the boundary nodes the periodic boundary condition is affecting the way gradient of phi is being calculated inside the domain. The interesting thing is that it is only happening in x=0 end of the domain and the other end apparently is fine.

I have found some similar threads regarding fvc::grad, for instance:
Gradient Operator
Or some threads on periodic boundary condition in general, for example:
How to setup cyclic BCs in simpleFOAM
Or threads like this on gradient at boundaries in general (read the last post):
Gradient at the boundary



My questions:

1. Am I right?! Is this a bug?!

2. Any body knows how to fix this?! Any thoughts?! Because I need to use the periodic boundary condition in a 2D turbulent combustion context, so I need a reliable gradient operator, otherwise I have to hardcode my own implementation of gradient operator, which is not an elegant solution.

3. Speaking of manually coding my custom version of gradient operator, does anybody know how the gradient is calculated at the boundaries?!

4. How the periodic boundary condition is implemented in OpenFOAM?! This wikipedia article gives an idea, but I want to know how and where it is implemented in OpenFOAM. Anybody knows a good starting point?! Which class?!



Appreciate your time for reading this and perhaps trying to answer either of the above...

Jalal

[JBUNSW]

Hi,

Just for the convenience of prospective readers facing the same bug in the future, I post the following lines:

Answer to 1 & 2:
The reason for this strange behaviour is still unknown to me. However, I could find a workaround. Previously I was using this piece of code to get the gradient of phi which was dependent on another geometricField, say a:

const volScalarField& phi = a;
const volScalarField& jGrad = mag(fvc::grad(phi));

By changing it to the below, it worked fine. I was not going to waste memory by creating a local volScalarField, however, this ugly solution works:

const volScalarField& phi = a;
const volScalarField jGrad = mag(fvc::grad(phi));

Hope this is rectified in the future releases of OpenFOAM.

Answer to 3 & 4:
For the boundary noeds under PBC scheme, the gradient is divided by two and assigned to each end of the domain, as the boundary nodes are essentially the clone of each other in the big picture. For internal nodes the gradient can be calculated using a central difference scheme, for the boundary nodes, however, one can use the forward or backward difference formulae.

Cheers,
Jalal

[niklas]

Originally Posted by JBUNSW

const volScalarField& jGrad = mag(fvc::grad(phi));

This line of code is a bug.
You are creating a reference to a temporary object and the result can be anything...which is also what you get. The fact that you are getting mostly correct numbers is just 'luck' and the boundaries of the allocated memory used by the mag-grad call has already been partially overwritten.

Look at the code below.

    const volScalarField magGrad = mag(fvc::grad(T));
    const volScalarField& refMagGrad = mag(fvc::grad(T));

    forAll(magGrad, i)
    {
      Info << i << ": " << magGrad[i] << ", " 
	   << refMagGrad[i] << ", " << magGrad[i]-refMagGrad[i] << endl;
    }

run this on a small grid ( I ran it on a 2x2x1) mesh with T the same for all.
This should produce only zeros, but if you run it, you will get something like this

0: 0, 4.52321e-317, -4.52321e-317
1: 0, 9.88131e-324, -9.88131e-324
2: 0, 2.122e-314, -2.122e-314
3: 0, 4.06912e-320, -4.06912e-320
End

Notice how the magGrad values are exactly zero, while the refMagGrad values are something else.
They look close to zero, but thats just coincidence. They are pointing to a memory location that is currently not being used.

Try removing the const statements and you will get a compiler error (as you should) for the reference call.

As for the other questions, I dont know if the problem will go away if you write your code correctly.

[JBUNSW]

Hi Niklas,

Thanks. Good to know that. Very nice explanation...

Jalal