Clarifications on ownership of nodal data (via fortran)

[pbrady]

I apologize if these are already answered in the manual by I couldn't find it. As noted in the manual, on nodal data, the boxes associated with an mfiter are not disjoint.

1) Are the boxes returned by tilebox() considered to be valid non-ghost cells even though they overlap? If so, are these overlapping non-ghost cells which are "shared" by different multifabs ever updated via fillboundary calls or is that job of the application developer?

2) Is there an idiomatic way of looping over them such that they appear disjoint? (My main interest here is in interfacing with external solvers which represent the grid as a disjoint union among the ranks)

[asalmgren]

Peter -- have you looked at Figures 4.3 through 4.8 in the User's Guide? Note that the boxes returned by tilebox do not overlap even in the case of nodal data.

[WeiqunZhang]

I think Peter is talking about the boxes returned by tilebox() from different Fabs. Yes, they could overlap in index space, but the memory location in arrays do not.

fillboundary does not modify those overlapping non-ghost cells. If you can give more details about what you want to achieve, maybe there is a way to do this.

As for question 2, MultiFab in C++ has a function called OverlapMask. This or this type of functions might be useful for this problem. If you give more details, that will help us understand.

PS. There are many tilebox types of functions. The Fortran interface is incomplete on this. We will finish the porting as soon as possible.

[pbrady]

Thanks for the responses. For a little context, I'm exploring amrex by implementing a simple incompressible flow solver on a staggered mesh. At the moment I'm restricting the investigation to a single level of refinement. Pressure is cell centered and each velocity component is nodal in its direction. With this setup, I have some initialization code (with added serialization to make sense of the io):

    do i=1,amrex_parallel_nprocs()
       if (i-1.eq.amrex_parallel_myproc()) then
          print *, amrex_parallel_myproc(), "init_pressure"
          call amrex_mfiter_build(mfi, P)
          do while(mfi%next())
             bx = mfi%tilebox()
             call amrex_print(bx)

             dp => P%dataPtr(mfi)
             dp = 0.0d0
          end do
          call amrex_mfiter_destroy(mfi)
       end if
       call MPI_BARRIER(MPI_COMM_WORLD, ierr)
    end do

I have a similar loop for the u velocity and the output for a particular set of parameters is:

$ mpirun -n 4 ./main2d.gnu.DEBUG.MPI.ex inputs 
MPI initialized with 4 MPI processes
           0 init_pressure
((0,0) (31,31) (0,0))
           1 init_pressure
((32,0) (63,31) (0,0))
           2 init_pressure
((0,32) (31,63) (0,0))
           3 init_pressure
((32,32) (63,63) (0,0))
           0 init_u
((0,0) (32,31) (1,0))
           1 init_u
((32,0) (64,31) (1,0))
           2 init_u
((0,32) (32,63) (1,0))
           3 init_u
((32,32) (64,63) (1,0))

In this setup, all u data at i=32 is logically (although not physically) shared by 2 processors due to the overlapping tileboxes (please correct me if I'm wrong). I was hoping to leverage hypre (specifically, the very simple struct interface) as a linear solver for different parts of the algorithm. The initialization phase involves setting up a grid, which for cell centered data is simply:

    call HYPRE_StructGridCreate(comm, n, grid, ierr)

    call amrex_mfiter_build(mfi, phi)
    do while (mfi%next())
       bx = mfi%tilebox()
       call HYPRE_StructGridSetExtents(grid, bx%lo(1:n), bx%hi(1:n), ierr)
    end do
    call amrex_mfiter_destroy(mfi)

    call HYPRE_StructGridAssemble(grid, ierr)

However, this only works if each tilebox specifies a disjoint set of indices. This grid object is then used to construct each part of a distributed Ax=b system. To use the struct interface, I would need some way of specifying the disjoint parts of the indices to initialize the grid object and then a way to communicate the solution to the logically shared non-ghost cells.

Is there a builtin way to accomplish this that hasn't been exposed to fortran yet or would it involve mucking around with amrex internals?

The other option might be to use the sstruct hypre interface which appears to allow some amount of sharing of face data.

[WeiqunZhang]

For single level, we could find some way dealing with this. But it would be much hard for multiple levels.

Yes, Hypre supports types. For example, HYPRE_SSTRUCT_VARIABLE_CELL, HYPRE_SSTRUCT_VARIABLE_XFACE, etc.

[pbrady]

I suspect the single level case is not worth the difficulty. We were planning on migrating to the sstruct interface for amr so it makes sense for me to simply abandon the struct hypre interface and use one that's a more natural fit for amrex. Thanks for clarifying.

[pbrady]

One more clarification question:

Is there a way to communicate the non-ghost overlapping nodal value (on a single level) such that the application can decide who the "true" owner is and have all others see the exact same value rather than having each fab perform an independent update on the logically shared nodal point.

ping @nncarlson

[WeiqunZhang]

I think this can be built the functions we have in AMReX. How do you decide who the "true" owner is? How do you like to pass that information to AMReX?

[pbrady]

One strategy could simply be that the "true" owner is the one with the lowest mpi-rank. I could imagine one could pass this information in as an option when creating the distromap that will then be used to construct multifabs.

[WeiqunZhang]

Suppose a rank is not the true owner of a non-ghost nodal point. Are you trying to avoid writing to that point? We could provide a mask array. Then you can if test the mask to decide if you need to work on a point or not. Note that in a general case a face of a nodal box could have multiple "true" owners. So we cannot provide a logical rectangular region for truly owned points. (We could, but then we would have to break a box into multiple smaller ones.)

Suppose you do update that point and then later a call to a function syncing up the data will overwrite it with the "true" value (or sum up to obtain the true value). Does that work for you?

[pbrady]

In general, I don't see a problem with a non-true owner updating some data which is later overwritten by the true owner. A mask array that reflected the true disjoint ownership of the points should be fine.

Just to make sure I understand your comment not being able to provide logical rectangular regions for truly owned points, is this something that can only happen on refined levels due to the positioning of the patches? I'm thinking of something like:

------|
   2  |---------| 
      |   3     |
      |---------|
      |   1     |
------|---------|  

where the numbers reflect the mpi rank. In this case simply picking the lowest mpi-rank would result in some of the shared faces going to 1 and some going to 2. Is this what you had in mind or is there something deeper?

[WeiqunZhang]

Yes, the figure shows what I mean. This could happen even on a single level run, depending on how the grids are generated.

[WeiqunZhang]

For the mask, do you like a boolean flag indicating the ownership or an int flag of mpi ranks?

[pbrady]

So long as all the communication is handled by amrex then a boolean flag should be fine.

[WeiqunZhang]

Yes, AMReX will handle the communication.

But using mpi rank as a way of deciding the true owner may not work. Two boxes on the same process could also overlap. So maybe we can use the Box index in the global array of boxes. That index is unique and is independent of DistributionMapping.

[drummerdoc]

Note however that we strive to make all operations as independent as possible of the box layout. So, as long as your operations are consistent with that mindset, and that your only using this rule as a way to be concrete, Weiqun's suggestion seems appropriate.

[WeiqunZhang]

Peter, I have implemented some new functions for creating nodal mask and syncing up nodal data.

call mf%override_sync(geom)  ! mf: MultiFab   geom: Geometry

With the call above, the data owner will override others. This is a parallel call containing communication inside. An owner mask is built internally. For nodal data shared by multiple grids, the grid with the lowest index number in the BoxArray is the owner. For periodic boundary, it is possible points in the same grid overlap periodically. In that case, the ownership is determined by lexicographical comparison of points.

One can explicitly build the owner mask provided by amrex.

type(amrex_multifab) :: owner_mask
call amrex_multifab_build_owner_mask(owner_mask, mf, geom)  ! no communication inside

This will build a mask for MultiFab mf. Note that the data type of mask is amrex_multifab and it contains real data. The values should normally be either 1.0 (owner) or 0.0 (non-owner). We did not use logical type here. So in principle, we could have half or one third of ownership.

One could provide their own mask, and call the following to sync up data:

call mf%override_sync(geom, my_mask)

Let us know if this meets your needs.

[pbrady]

Thanks @WeiqunZhang. I'll take a look at this over the next few days.

[nncarlson]

Thanks @WeiqunZhang. I've talked to Peter and this looks like something we can work with for the time being. However using a real type mask array is just plain wrong. I know there issues with Fortran logical type interoperating with C/C++, but is there a good reason for not at least using an integer mask array?

[WeiqunZhang]

Depending on how you view this. If we view them as weights, real makes sense. (And the implicit conversion from double to bool is actually well specified in C/C++.) Anyway, we can certainly change it to integer.

[nncarlson]

That's interesting twist (and I should have read your earlier message more carefully). How does having a nodal value half owned by two fabs work? (scattering from owner to non-owner) Is there a use case for this? Anyway, we would greatly prefer the mask to be an integer -- Fortran has no implicit conversion from real to logical.

[WeiqunZhang]

In the latest development branch, the owner mask is now imultifab. The underlying data type is integer. As expected, zero means false and non-zero means true.

type(amrex_imultifab) :: owner_mask
call amrex_imultifab_build_owner_mask(owner_mask, mf, geom)

[nncarlson]

@WeiqunZhang It looks like these changes are working for us, but we did notice one thing that I'd like you to confirm. @certik reports that It appears that override_sync needs to be called prior to each and every call to fill_boundary. Is that correct? If so, that makes it fairly awkward to use. Would it be possible to have override_sync define how all subsequent fill_boundary calls behave for that multifab?

[certik]

What @nncarlson wrote is indeed what I observed.

[WeiqunZhang]

It is correct that currently you need to call override_sync everytime before calling fill_boundary.

How about adding an optional flag to fill_boundary indicating override_sync should be called? There is a parallel communication call inside override_sync. This will allow us to merge that into the communication perfomed in fill_boundary in the future without affecting the user interface.

[certik]

One complication is that we use a custom mask for the override_sync. Otherwise I think this would work.

[nncarlson]

Dragging around the mask imultifab is what makes the repeated calls to override_sync really awkward. What we are looking for is a single call, like override_sync that says "use this ownership mask for every subsequent call to fill_boundary. If it were just an additional matter of using an optional (scalar) flag to override_sync that would be okay, but I suspect that's not enough with a user-supplied mask.

[drummerdoc]

Could it be that the best way to manage this from your perspective is to make a very singe/thin (user) class that implements the sequence just as you like? Just a thought...

[WeiqunZhang]

Then I guess what you want is something like,

call mf%set_sync_mask(msk)  ! store a copy of this in multifab.  msk can be destroyed after this

If fill_boundary is called after this, it will also fix the synchronization using the mask.

[nncarlson]

@WeiqunZhang yes, this sounds like what I had in mind. However, I don't know the internals of amrex, and this may be at odds with the design and implementation. If that's the case, then it may be better for us to follow @drummerdoc 's advice and write our own layer over amrex that provides the desired behavior.