[ander]

Ok, so you would extend the gradient and flux limiter containers to also contain the ghost cells, but set these values to zero, store boundary conditions defined at the faces in the ghost cells, and then only require one loop over all faces to evaluate the flux balance? I guess this is a tradeoff between memory usage and code complexity then.

I’ve also noticed that I need to need to use function pointers, inherited overloaded functions or something like this when implementing this stuff, when the given setting is read from a config file or similar and not known at compile time. (using C++ in my case, but I guess this need arise regardless of the language).

Not sure how much slower this is compared to versions of the same code where the actual function used (for instance flux scheme or boundary condition) was known at compile time and could be inlined.

[hunt_mat]

Originally Posted by ander

Hi,

I have a question regarding how to implement boundary conditions for an unstructured cell centered fvm for compressible flow.

Say you have a no slip wall boundary, and you have one layer of ghost cells. You then want the flux contribution from the velocity to be zero. In 1D: F = [rho u, rho u + p, E]^T = [0, p, 0]^T

In the book "Computational Fluid Dynamics: Principles and Applications" by Blazek, it is proposed to set u_ghost = - u_1, (assuming u_1 is the cell nearest the wall). Up to this point I see no problem.

However, if a reconstruction procedure is used, what would now happen? The book doesn't consider this as far as I'm aware. If a reconstruction was used, The flux that would be calculated at the boundary would be F(U_L, U_R), where U_L and U_R are reconstructed values at the face. If the constant reconstruction was used, the riemann solver would return a zero velocity flux contribution, but I guess it won't for reconstructed values that doesn't necessarily obey u_R = -u_L.

Should reconstruction be turned off at boundaries?

The way that I've looked at this problem in a very simplistic way is to note that the quantity at a point x_j is approximately the average of its two neighbours. So

f_j=0.5*(f_(j-1)+f_(j+1))

If you have a condition which states that the velocity is zero at a boundary, you have that:

0=0.5*(f_(-1)+f_1) => f_(-1)=-f_1

That should keep the boundary condition into your conservation law.

[sbaffini]

Originally Posted by ander

Ok, so you would extend the gradient and flux limiter containers to also contain the ghost cells, but set these values to zero, store boundary conditions defined at the faces in the ghost cells, and then only require one loop over all faces to evaluate the flux balance? I guess this is a tradeoff between memory usage and code complexity then.

Not that I am actually suggesting this but, yes, it's a possibility. My preferred choice is:

1 - loop over patches with faces sharing the same boundary conditions

2 - in the previous loop, use if statements to detect the patch bc and pass the relative bc function (fp) to another function that actually applies the bc (fb)

3 - fb does a loop on the faces of the current patch and uses the passed in fp to actually apply the bc fluxes

4 - loop over the remaining interior faces and just use the standard flux of the used scheme

The day this is going to be my bottleneck I will reconsider it (not the specific approach I am using, this whole aspect independently from the specific implementation in use). But may I be damned if this will ever happen. There's much more around to keep under control: the actual convective scheme (say, Roe), the linear iterations for the implicit scheme and the gradient computations are typically much more costly

[ander]

This approach seems conceptually quite similar to the approach I have chosen for my hobby code project then, the exception being that I have used ghost cells instead of direct flux computation.

Don't know what language you are using (I assume C++ or Fortran). I use C++ and chose virtual overloaded functions instead of function pointers for the purpose of calculating ghost values from a boundary condition. Not really sure if there is a performance gap between these two. I have used function pointers for the convective scheme though.

My understanding is that it is beneficial to use inlined functions for functions called within a large loop performance wise, but this seems hard to to in general, when various solver settings are not known at compile time.