Bug in fvc::grad()
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:
I have found some similar threads regarding fvc::grad, for instance:
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
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...
My ugly solution...
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:
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.
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.
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.
Thanks. Good to know that. Very nice explanation...
|All times are GMT -4. The time now is 19:57.|