Discuss on Neumann boundary treatment

[frankyugit]

Dear guys

Wish to discuss whether Neumann boundary is treated explicitly or implicitly
In source code "updateBoundaryGhosts.inl"
switch (flow->boundaries[XMINUS][0].type)
{
case DIRICHLET:
qx[j][-1] = flow->boundaries[XMINUS][0].valuemesh->dy[j];
break;
case CONVECTIVE:
beta = flow->boundaries[XMINUS][0].value * dt / mesh->dx[0];
qx[j][-1] = (1-beta)qx[j][-1] + betaqx[j][0];
if (timeStep == startStep)
qx[j][-1] = qx[j][0];
break;
case NEUMANN:
qx[j][-1] = qx[j][0];
break;
default:
break;
}
Looks like Neumann boundary is treated explicitly like Directly boundary
However in source "generateBC1.inl"
switch (flow->boundaries[XMINUS][0].type)
{
case CONVECTIVE:
case DIRICHLET:
bc1x[j][0] += coeffMinusqx[j][-1]/dy[j];
break;
default:
break;
}
Here, treatment of Neumann boundary is missing, which may only be the reason of implicit treatment of Neumann boundary, is that a small bug, or there's some other design to deal with it

[piyueh]

Hi @frankyugit

Neumann BC is still under development. Since we haven't had the need of Neumann BC in our applications and research so far, so its development is in slow progress. If you are going to implement it and willing to make a pull request, we will appreciate it!

Or you can also wait for the next release of PetIBM. We can have that ready in the next release.

Any comment @mesnardo ?

[mesnardo]

As @piyueh mentioned, Neumann boundary conditions for the velocity field are not properly implemented in PetIBM right now and should not be used.
I believe issue #8 was created for that reason.

So far, when dealing with external flow around bodies, we have used a convective condition at the outlet and did not feel the need to implement Neumann conditions.

@frankyugit, I have also seen many examples using a zero-gradient condition at the outlet for the flow around a cylinder.
In the past, I have been using OpenFOAM to solve the flow around a bluff body at Reynolds numbers of O(10^3) and found it was difficult to properly advect the wake vortices outside the domain with a simple zero-gradient outlet condition for the velocity.
This is mostly why we have implemented a convective boundary condition in PetIBM and not yet Neumann condition.

This feature should definitely be available in the next release.

[frankyugit]

@piyueh @mesnardo
Guys, thx for the reply and detailed answer. Now I know Neumann BC status in PetIBM.
Also wish to discuss a further issue:
Based on the code in "updateBoundaryGhosts.inl", it looks like PetIBM is trying to treat Mixed(Convective) BC explicitly. As for every time step, the ghost point value is explictly updated based on the gradient value and same ghost point value last step. But based on my understanding, explicit treatment of Neumann BC or Neumann part in Mixed BC may bring in possible numerical instability.
So any consideration for this design?

[piyueh]

This design is for the symmetricity of the coefficient matrix (DBnG) in Poisson equation. In other words, this is to make sure the gradient (G) and divergence (D) operators have the relationship G = - D^T. If we use implicit treatment, there will be additional modifications to the G and D operators due to BCs, and then the symmetricity may be broken. This is also caused by the fact that we don't use any pressure BCs and pressure ghost points.

Convective boundary condition is indeed a complex topic. But for our application (incompressible and low Reynolds number), we found this simple treatment already gives us satisfying results.