Functors on cells with certain mask value

[shreyasashok]

Hi, I'm trying to assess whether Amrex is appropriate for my use case. I am trying to develop a Lattice-Boltzmann method based code.

For my use case, I would need to be able to easily run functors only on cells with a certain value of a mask. For example, I would be implementing actuator disks, which would encompass <1% of the cells in the domain. I would not want to run a functor on the whole domain just to update these few cells.

Additionally, I want to add some additional data attributes only to these specific cells. Adding the data to the whole MultiFab system would severely impact my memory usage. Is there a good way to go about this? I would want to ensure coalesced memory access wherever possible.

[WeiqunZhang]

What kind of operations do you do on those special cells? Do you need to access regular or special data in other cells?

[WeiqunZhang]

You can probably do something like this. Sorry for not having comments. But I am happy to explain every single line. The code is not tested and is likely to have bugs. But I believe the idea should work.

iMultiFab mask(...);
iMultiFab idx1d(...);
int ntot = 0;
Vector<int> n_special_cells;
for (MFIter mfi(mask); mfi.isValid(); ++mfi)
{
    Box const& b = mfi.validbox();
    int ncells = b.numPts();
    auto const& mfab = mask.array(mfi);
    auto const& ifab = idx1d.aray(mfi);
    int N = Scan::PrefixSum<int>(ncells,
                              [=] AMREX_GPU_DEVICE (int icell) -> int
                              {
                                  GpuArray<int,3> ijk = b.atOffset3d(icell);
                                  int cell_mask = 1; // or 0                                                                                                                 
                                  mfab(ijk[0],ijk[1],ijk[2]) = cell_mask;
                                  return cell_mask;
                              },
                              [=] AMREX_GPU_DEVICE (int icell, int const& s)
                              {
                                  GpuArray<int,3> ijk = b.atOffset3d(icell);
                                  ifab(ijk[0],ijk[1],ijk[2]) = s + ntot;
                              },
                              Scan::Type::exclusive, Scan::retSum);
    n_special_cells.push_back(N);
    ntot += N;
}

TableData<Real,2> spatial_idx({0,0}, {ntot-1,AMREX_SPACEDIM});
auto const& si = spatial_idx.table();
auto const& mma = mask.const_arrays();
auto const& ima = idx1d.const_arrays();
amrex::ParallelFor(mask, [=] AMREX_GPU_DEVICE (int box_no, int i, int j, int k)
{
    if (mma[box_no](i,j,k)) {// special cell                                                                                                                                 
        int m = ima[box_no](i,j,k);
        si(m,0) = box_no; // box #                                                                                                                                           
        si(m,1) = i; // x-direction index                                                                                                                                    
        si(m,2) = j; // y-direction index                                                                                                                                    
        si(m,3) = k; // z-direction index                                                                                                                                    
    }
});

MultiFab regular_data(...);
auto const& rd = regular_data.arrays();

int ncomp; // Assuming the sparse data has ncomp reals.                                                                                                                      
TableData<Real,2> sparse_data({0,0}, {ntot-1, ncomp-1});
auto const& sd = sparse_data.table();
amrex::ParallelFor(ntot, [=] AMREX_GPU_DEVICE (int m)
{
    // sd(m,0:ncomp-1) store the data for special cell # m.                                                                                                                  
    // si(m,0) gives you the box_no                                                                                                                                          
    // si(m,1:3) gives you the 3d spatial index                                                                                                                              
    // If you need access to regular data,                                                                                                                                   
    //     use rd[si(m,0)](si(m,1),si(m,1),si(m,2))                                                                                                                          
});

[shreyasashok]

I think this is along the lines of what I am looking for! I will test this out. Thank you!

[drummerdoc]

Also, is there any real benefit to storing these in the MultiFab (rather than a sparse/distributed set of functor objects)?

[WeiqunZhang]

Accessing neighbor cells is more involved, like rd[si(m,0)](si(m,1),si(m,1),si(m,2)) for regular data and even longer for irregular data. But I suppose they can be wrapped nicely in a class.

Also, communication can be tricky.