Question about setting physical boundary value while doing the initial projection

[ruohai0925]

Hello all,

I have a question about the function doMLMGNodalProjection. In this line, https://github.com/AMReX-Codes/IAMR/blob/29f74453b9ef0bb4c69db43354887a840cbc6fa4/Source/Projection.cpp#L2467 , it uses the function set_boundary_velocity. But it shows that this function sets velocity in ghost cells to zero except for inflow. I am confused about why it should be like this?

Here are my thoughts. To do initial cell centered velocity projection, we need to fill the ghost cells. But the ghost cells should be set to reflect the physical boundary condition (if only coarsest levels is considered), not just simply set as 0. So I do not understand the above function.

Could someone please give me some comments and tips?

Thanks so much.

Ruohai

[ajnonaka]

@WeiqunZhang knows more details but he explained this to me once. Right now the way the stencils are set up, the non-inflow ghost cells need to contain zeros. It is currently on the MLMG solver TODO list to eliminate the need for a user to call set_boundary_velocity and just handle these issues internal to the solver. I can't be of much more help than this except to just use this function for now. If you want a cleaner version that should work for incompressible flow I have a simpler version in my MAESTROeX github repo in MaestroNodalProj.cpp that you can take. The main difference is that I took out some code that dealt with strongly non-zero "S" in the div(u)=S constraint at outflow --- stuff you won't need for incompressible flow anyway.

[ruohai0925]

@ajnonaka Thanks for your explanation. I am not sure if I understand it very clearly, here is my thought.

So you mean that we can now just use the set_boundary_velocity. I think I can use it, but I still do not know why the ghost cells should be 0 except the inflow. For example, if it is the periodic physical boundary condition, the ghost cells of velocity should reflect this boundary condition. (I understand that the physical boundary condition belongs to one of the internal boundary conditions, but after using the set_boundary_velocity function, I can not get the correct value of the ghost cells.)

When you say "stencils", I do not quite understand from a user's perspective. Do you mean according to different stencils, should users set up the physical bc in a different way? Like the codes here: https://github.com/AMReX-Codes/IAMR/blob/29f74453b9ef0bb4c69db43354887a840cbc6fa4/Exec/run2d/PROB_2D.F#L1564

I will read your Maestro code. It is really clear to divide the projection process into different situations and assign different Multifabs for coefficients. Thanks for sharing.

Ruohai

[WeiqunZhang]

Ghost cell velocity is set to zero for convenience. This is a detail that should have been hidden from the user. For nodal projection with homogeneous Neumann boundary, div u on the boundary is supposed to be computed with only interior velocity (see Almgren et al. for operator D).

[ajnonaka]

After you create MLNodeLaplacian, at some point you need to call setDomainBC to tell it what kind of physical domain boundary conditions you have (both IAMR and MAESTROeX do this). This internally tells the MLNodeLaplacian object and any solvers associated with this about any modified stencils it needs to use near domain boundaries. The FORT_DENFILL routine is ultimately called via FillPatch routines when you simply need to fill in ghost cells for whatever reason. So this doesn't have anything to do with the MLNodeLaplacian or solvers.

[ruohai0925]

@WeiqunZhang Thanks for your explanation. For the doMLMGNodalProjection subroutine, there are many possibilities for the inputs parameter vel. If it is the initial velocity projection step, the input vel is U_old. If it is the initial sync projection step, the input vel is (U*-U_old)/dt. So does it mean no matter which situation appears, we all should use set_boundary_velocity?

Note the ghost cell value between U_old and (U*-U_old)/dt are different because of different physical boundary conditions. But if the rhs:=div (vel) does not use the ghost cell value, that will be fine for the subroutine.

I am writing the cell-centered pressure codes, so maybe I should consider the ghost cell value for calculating the div (vel).

Jordan

[WeiqunZhang]

yes, because the set_boundary_velocity function only writes zero. so it doesn't matter what physical unit is.