Particle kick, boundary condition for phi in non-periodic dimensions

[chongchonghe]

Describe the bug
Here

quokka/src/simulation.hpp

Line 1586 in fab8caa

amrex::GpuBndryFuncFab<setFunctorParticleAccel> boundaryFunctor(setFunctorParticleAccel{});

the boundaryFunctor is no-op by default. Therefore, the non-periodic dimensions of phi_extended is unfilled, causing the code to crash. Is this intentional? is it expected that the user would have to specify the boundary condition for non-periodic dimensions?

@BenWibking Can you add an option for a foextrap B.C. at some specified dimensions? I'll need this for disk simulations with star formation where the z axis is open B.C.

To Reproduce
Steps to reproduce the behavior:

1. Compile this problem '...'
2. Run this problem '...'
3. See error

Additional context
Add any other context about the problem here.

[dosubot]

Yes, the default no-op boundaryFunctor for non-periodic dimensions is intentional—users are expected to implement custom boundary conditions for open boundaries. If you set a boundary to ext_dir (external Dirichlet) and don't provide a specialization of setCustomBoundaryConditions, the ghost cells remain uninitialized, which can cause crashes. This design prevents silent errors and ensures users explicitly handle physical boundaries as needed for their simulations source.

Currently, there is no built-in option for a fixed extrapolation (foextrap) boundary condition. To implement foextrap for a specific dimension (like the z-axis in disk simulations), you should specialize the setCustomBoundaryConditions (or setCustomBoundaryConditionsFaceVar for face-centered variables) static method in your problem generator. This method is called for each ghost cell on boundaries set to ext_dir, and you can fill those cells using your desired extrapolation logic source.

The unified boundary condition interface (quokka::BC and quokka::BCType) lets you specify periodic, Dirichlet, reflecting, and mixed BCs, but foextrap is not included out of the box source. For disk simulations with open z boundaries, you’ll need to implement the extrapolation manually in your boundary condition specialization.

To reply, just mention @dosu.

---

^{How did I do? Good | Irrelevant | Incorrect | Verbose | Hallucination | Report 🐛 | Other}    

[chongchonghe]

@BenWibking Is the hydro equation using 3 ghost cells for phi? Then, can we assume the boundary cells of phi are 'correct' and use them directly in kickParticlesAllLevels?

I think we can simply copy nghost_phi ghost cells (instead of 0) here, right?

quokka/src/simulation.hpp

Line 1576 in fab8caa

amrex::MultiFab::Copy(phi_extended, phi[lev], 0, 0, 1, 0);

[BenWibking]

@BenWibking Is the hydro equation using 3 ghost cells for phi? Then, can we assume the boundary cells of phi are 'correct' and use them directly in kickParticlesAllLevels?

I think we can simply copy nghost_phi ghost cells (instead of 0) here, right?

quokka/src/simulation.hpp

Line 1576 in fab8caa

amrex::MultiFab::Copy(phi_extended, phi[lev], 0, 0, 1, 0);

I don't understand your question. What problem you trying to solve?

[chongchonghe]

Here in KickParticles

quokka/src/simulation.hpp

Line 1586 in fab8caa

amrex::GpuBndryFuncFab<setFunctorParticleAccel> boundaryFunctor(setFunctorParticleAccel{});

the boundaryFunctor is no-op by default. Therefore, only periodic dimensions of phi_extended is filled. If there are non-periodic dimensions, phi_extended will have nan values and acceleration will have nan values as well, causing the code to crash. This is the case in my disk slab simulations where x and y dimensions are periodic and z dimension is outflow.

I understand that the user is expected to specify the boundary condition for phi, by overwriting setFunctorParticleAccel. Is this right?

We can add support for predefined boundary conditions for phi, like foextrap, for an simpler user interface.

[markkrumholz]

We can add support for predefined boundary conditions for phi, like foextrap, for an simpler user interface.

This is mostly right, but the boundary conditions you can select for the potential are a bit different than those for hydro, because the Poisson equation is an elliptic system instead of a hyperbolic one. Thus for example foextrap does not make sense for an elliptic system. The boundary conditions that are physically meaningful are the ones that AMReX supplies as options. In general the boundary conditions on this sort of system specify either the value of the potential on the face (Dirichlet), the derivative on the face (Neumann), or some algebraic combination of the two. Our filler should be based off this choice. For for example for Dirichlet the sensible choice is that the potential in the first ghost zone is such that the average of that and the potential in the first real zone is equal to the Dirichlet value set on the domain edge, for Neumann we should choose the first ghost cell value such that the slope between the final real cell and first ghost cell matches the specified derivative, and so forth.

[chongchonghe]

OK, mathematically that makes sense. However, I don't understand why we can't have foextrap in this case -- it's equivalent to Neumann BC with the derivative equal to the derivative of the first two real cells, i.e. the force is unchanged from the first real cell to the first ghost cell.

The point I want to make is that we can provide pre-defined, simply-to-use BC options so that the user doesn't have to specify the value or first derivative at the boundary, which is usually non-trival, and in practice the user have to extrapolate phi inside setFunctorParticleAccel anyway.

[markkrumholz]

You can't use foextrap as a boundary condition when actually solving the Poisson equation because it makes the problem not well-defined. Because the equation is elliptic, the derivative everywhere depends on the boundary, so if you say that the derivative at the boundary depends on the derivative in the interior, the problem becomes circular.

Now you are certainly right in a mechanical sense that, once you have solved the Poisson equation using some other boundary condition, you could then fill in the ghost zone potential using first order extrapolation. But I think you should not do so, because that choice will not in general be consistent with the choice that was made when the Poisson equation was solved. My proposal wasn't that we shouldn't provide pre-defined BC options that fill the ghost zone, my point is that we in fact already do provide pre-defined BC options when solving the Poisson equation, and we should use this choice to determine how to fill the ghost zones too so that we enforce consistency.

[BenWibking]

Ah, sorry. I forgot I had written that code. I originally wrote it to extrapolate the acceleration at physical boundaries, but that required too much communication when I did the INCITE scaling tests, so I rewrote it to use the potential.

Regarding Mark's point, I agree it would be ideal if we could just set the boundary condition once, but the problem is that amrex::OpenBC requires using exactly 1 ghost cell for phi (I'm not sure whether MLMG has this restriction or not). So, at least when using 3D open boundaries, we always have to perform some extrapolation ourselves when growing the number of ghost cells in the particle kick routine.

[BenWibking]

For the DiskGalaxy problem, I think I forgot about properly resolving this issue because the particles just never get close enough to the boundary that it matters at all.

Given that the only two boundary conditions we support for solving the Poisson equation are open (or setting Dirichlet phi=0 BC as an approximation to this) or periodic, I would be fine if the potential were computed with hoextrap at the non-periodic boundaries. We cannot use foextrap, because then the acceleration would be zero at the boundary, which is wrong. For foextrap, we will get a constant acceleration, which (I think) that is what we want physically.

[BenWibking]

It's important to note that the ghost cell at domain boundaries for the AMReX linear solvers (unlike the hyperbolic solvers) always refers to the value exactly at the boundary, not the value at +0.5*dx away from the boundary.

[BenWibking]

I don't think this causes any mathematical issues about the well-posedness of the solution. We are calculating a series expansion of phi at the boundary and using that to extrapolate the solution, so I think it should be consistent with the original solution up to the order of the expansion (unless there was a discontinuity at the boundary, but I don't think that's physically possible for the gravitational potential).

However, hoextrap does not understand that the 1 ghost cell in phi is exactly at the boundary, so you will have to write your own quadratic extrapolation function to fill the extra ghost cells.

[BenWibking]

Since particles can wander into edges and corners of the ghost region, it is necessary to handle those with 2D and 3D Taylor series expansions.

Alternatively, you could expand the domain size to include the ghost cells needed for particle kicks (for rhs and phi only), do the Poisson solve on the expanded domain, and then you have a fully consistent solution in the ghost cells up to the solver tolerance.