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

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

 November 14, 2012, 07:33 Bug in fvc::grad() #1 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 5 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: Code: ``` 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 Last edited by JBUNSW; November 14, 2012 at 08:44.

 November 14, 2012, 22:42 My ugly solution... #2 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 5 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: Code: ```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: Code: ```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

November 15, 2012, 03:05
#3
Super Moderator

Niklas Nordin
Join Date: Mar 2009
Location: Stockholm, Sweden
Posts: 693
Rep Power: 20
Quote:
 Originally Posted by JBUNSW Code: `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.
Code:
```    const volScalarField magGrad = mag(fvc::grad(T));

{
Info << i << ": " << magGrad[i] << ", "
}```
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.

 November 15, 2012, 04:09 #4 New Member   James Behzadi Join Date: Oct 2011 Location: Sydney, Australia Posts: 27 Rep Power: 5 Hi Niklas, Thanks. Good to know that. Very nice explanation... Jalal