CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM

Bug in fvc::grad() and periodic boundary condition

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree2Likes
  • 2 Post By niklas

LinkBack Thread Tools Search this Thread Display Modes
Old   November 14, 2012, 06:33
Default Bug in fvc::grad()
New Member
James Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 27
Rep Power: 14
JBUNSW is on a distinguished road
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...


Last edited by JBUNSW; November 14, 2012 at 07:44.
JBUNSW is offline   Reply With Quote

Old   November 14, 2012, 21:42
Default My ugly solution...
New Member
James Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 27
Rep Power: 14
JBUNSW is on a distinguished road

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.

JBUNSW is offline   Reply With Quote

Old   November 15, 2012, 02:05
Super Moderator
niklas's Avatar
Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 29
niklas will become famous soon enoughniklas will become famous soon enough
Originally Posted by JBUNSW View Post
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

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.
maddalena and New_Old like this.
niklas is offline   Reply With Quote

Old   November 15, 2012, 03:09
Thumbs up
New Member
James Behzadi
Join Date: Oct 2011
Location: Sydney, Australia
Posts: 27
Rep Power: 14
JBUNSW is on a distinguished road
Hi Niklas,

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

JBUNSW is offline   Reply With Quote


cyclic bc, fvc::grad(), gradient at boundary, periodic bc

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On

All times are GMT -4. The time now is 03:49.