Conditions for a well-formed Octtree-ish mesh

[kweide]

In his last response under #2140, @WeiqunZhang wrote:

Blocking factor is usually better to be a power of 2. For proper nesting, it needs to be greater than (ceil(nghost,ref_ratio)+ninterp)*ref_ratio, where ninterp (typically 1) depends on the scheme used for coarse/fine interpolation. Yes this is because only coarse valid cells are used in interpolation.

(I assume the two-argument function ceil stands for "round up to multiple of 2nd arg" - please correct if wrong.)

So the statement of the requirement is

blocking_factor > (ceil(nghost,ref_ratio)+ninterp)*ref_ratio

1. Is this (still) strictly required? Or is a looser condition sufficient - for example, replacing with >= instead of >? (Some of the discussion in Chombo-style "incrementing" FluxRegisters? #2140 suggests that >= could be sufficient.)

2. If we make some common assumptions, is the answer still the same?

   1. ref_ratio == 2
   2. We are only interested in cell-centered data
   3. blocking_factor is a multiple of ref_ratio

3. To make fully octtree-ish assumptions and other assumption that apply in Flash-X, additionally assume
   5. nghost is even
   6. blocking_factor = 2*max_grid_size and max_grid_size is also even
   7. ref_ratio and nghost are independent of spatial direction (for 2D or 3D setups)
   The statement of requirement then becomes

      max_grid_size > nghost + ninterp

Is this the right condition, or is this sufficient but not necessary?

I tend to think of this as a condition on nghost, i.e.,

   nghost < max_grid_size - ninterp

For comparison, in Flash-X we typically require both

   NGUARD >= ceil( max(minLayers, ceil(minLayers/2)+ninterp)  , 2 )
   NGUARD <= max_grid_size/2

where minLayer is the number of ghost cell layers in which correctly filled data is required by some physics algorithm and max_grid_size stands for a hardwired octtree block size.

[kweide]

As an added note, this question is of practical importance. I am trying to figure out whether in a given test configuration we (Flash-X team) have been violating some requirement.

In the problematic configuration we have NDIM=2, NGUARD = 4, max_grid_size = 8 (both directions), ninterp = 2.
We use a PPM Hydro algorithm with flattening, where the flattening can have the effect that actually 4, rather than 3, layers of ghost cells influence the solution. This is in RZ geometry, which may be another source of difficulty...

Results of the test simulation changed noticeably when libamrex was upgraded from version 22.10 to 23.09.

[kweide]

I've been trying to avoid it, more specifically trying to avoid building multiple versions of libamrex. But yes, that would obviously help narrowing it down!

[kweide]

I did the bisection. It turns out that the "culprit" is 6eb91be.

Since "This fixes a very old bug", I shall consider the new results more reliable.
I don't quite understand why this bug fix affected only a couple of tests out of a much larger set.

[WeiqunZhang]

This only affected when there are fine level boxes next to non-periodic boundaries. It's also possible that the bogus values set in the vector for BC type did not affect much because it went to the else part of some if test that called Flash-X boundary filling functions that did not use those BC types set in AMReX.

[kweide]

The situation in our test fits your description: yes, there as a f/c boundary that runs into (i.e., is locally perpendicular to) on non-periodic domain boundary; and this is close to where solution differences were most visible.

Would the if test that you mention be in the Flash-X code or in AMReX code? Could you perhaps give a pointer to where I can find this kind of test?

[WeiqunZhang]

I just looked more carefully into this. When we need to fill a fine cell from spatial interpolation and it is near the physical boundary, we may need to use the physical boundary conditions stored in coarse level ghost cells to compute slopes. The BC type bug does not affect how ghost cells outside physical boundary are filled, because it calls a user provided function that does not use the BC type we set. For Flash-X, this function is gr_fillPhysicalBC. What the BC type bug affects is how slopes are computed. See

amrex/Src/AmrCore/AMReX_MFInterp_C.H

Line 13 in 8755645

Real mf_compute_slopes_x (int i, int j, int k, Array4<Real const> const& u, int nu,

for an example. You may wonder whether we treat ext_dir and hoextrap differently. This is because historically BoxLib codes assumed that the boundary location of ext_dir and hoextrap is exactly on the face, not half cell outside. So this is actually a known defect for codes that assume the location is at half cell outside. Maybe we should introduce a new type called ext_dir_cc. So I assume some Flash-X tests have ext_dir BC. If so, in the interpolation function, it was treated as being at cell centered before that bug fix and now as being on the domain face.

[WeiqunZhang]

Sorry for the typo. ceil(nghost,ref_ratio) should read ceil(nghost/ref_ratio). Yes, >= should be sufficient. So the correct condition should be blocking_factor >= (ceil(nghost/ref_ratio)+ninterp)*ref_ratio. If nghost is even and ref_ratio is 2, that becomes blocking_factor >= nghost + ninterp*2. It can also be written as blocking_factor/2 > nghost/2 + ninterp. The way to think about it can like this. If we look at cells on each level, it's either all cells in a block are refined or none, where the block size is blocking_factor. So the smallest gap between level $\ell$/($\ell + 1$) boundary and $\ell$/($\ell -1$) boundary has blocking_factor/2 level $\ell$ cells. If level $\ell +1$ has nghost ghost cells, the underlying region on level $\ell$ is nghost/2 cells wide. If the interpolation requires ninterp on each side, we require blocking_factor/2 >= nghost/2 + ninterp.

[kweide]

Thank you for the correction and the explanation!

Part of my difficulty in understanding was at first the language: what you call grid, we call block; what you call block, we don't really have a term for. Also, where you write "either all cells in a block are refined or none", it makes more sense to me when I change "are refine" to "are covered by cells of at least one higher level of refinement".

So I think I basically understand now. One of these days I should write down what exactly the conditions are for Flash-X-with-PARAMESH and compare then to those for Flash-X-with-AMReX...

A couple more things...

The condition

   blocking_factor/2 >= ceil(nghost/2) + ninterp

comes from considerations for ghost cell filling at f/c boundaries (filling fine GC's based [at least partially] on coarse cells).
But isn't there a similar condition from similar conditions GC filling at grid boundaries where there is no refine change?
That condition should be

   blocking_factor/2 >= nghost

So that overall the following should apply (still assuming ref_ratio=2 but not necessarily requiring that nghost be even):

   blocking_factor/2 >= max( nghost, ceil(nghost/2) + ninterp )

Do you agree?

[WeiqunZhang]

But isn't there a similar condition from similar conditions GC filling at grid boundaries where there is no refine change?

In that case, we fill as much as possible with the same refinement level data first and then fill the rest with interpolation from the coarse level data.

[kweide]

... and then fill the rest with interpolation from the coarse level data.

I see. Does that also work when "the rest" comes not from coarse level data, but from a boundary condition at the boundary of the simulation domain?

[WeiqunZhang]

Yes, in that case, the physbc function will be called to fill ghost cells outside the domain.