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.