Boundary Conditions Handling Performance

[ToguyC]

Hello !

Thank you again in advance for taking the time to answer.

I am currently rewriting a piece of code I had where boundary conditions where handled manually by having a MultiFab dedicated to flag cells as BOUNCE_BACK or BULK (fluid) and do computation accordingly.

However, I find this method not in line with the AMReX philosophy, and so I am trying to change in order to use BCRec and boundary functor.

The model I am doing is a simple lid-driven cavity for an LBM solver in D3Q19 scheme.

I am getting the correct result, that is not the problem, but the performances dropped significantly by doing so.
Is there a better way of doing boundary conditions ?

Here is the pieces of code that I have to handle boundary conditions:

AMREX_FORCE_INLINE AMREX_GPU_DEVICE Real feq(int l, Real rho, Array1D<Real, 0, AMREX_SPACEDIM> u) {
    Real ck_u = c[l][0] * u(0) + c[l][1] * u(1) + c[l][2] * u(2);
    Real usqr = 1.5 * (u(0) * u(0) + u(1) * u(1) + u(2) * u(2));
    return rho * w[l] * (1. + 3. * ck_u + 4.5 * ck_u * ck_u - usqr);
}

// Custom boundary condition functor for lid-driven cavity
struct CavityBCFill {
    Real ulb;  // Lid velocity
    
    void operator()(const Box& bx, FArrayBox& dest, const int dcomp, const int numcomp,
                    Geometry const& geom, const Real /*time*/, 
                    const Vector<BCRec>& /*bcr*/, const int /*bcomp*/, const int /*orig_comp*/) const {
        const auto& domain = geom.Domain();
        const int* domlo = domain.loVect();
        const int* domhi = domain.hiVect();
        
        auto data = dest.array();
        
        // Define the "inside" box (valid domain cells only)
        const IntVect liv(AMREX_D_DECL(domlo[0], domlo[1], domlo[2]));
        const IntVect uiv(AMREX_D_DECL(domhi[0], domhi[1], domhi[2]));
        const Box inside(liv, uiv);
        
        // Process ghost cells - loop over each cell and population direction
        ParallelFor(bx, [=, this] AMREX_GPU_DEVICE(int i, int j, int k) {
            IntVect iv(i, j, k);
            
            // Only process if this is actually a ghost cell
            if (inside.contains(iv)) return;
            
            // Loop over all populations
            for (int q = dcomp; q < dcomp + numcomp; ++q) {
                // Check if this population direction points INTO the domain
                IntVect neighbor(iv[0] + c[q][0], iv[1] + c[q][1], iv[2] + c[q][2]);
                
                if (inside.contains(neighbor)) {
                    // This population crosses the boundary - apply BC
                    
                    // Determine which wall we're on
                    bool at_low_wall[AMREX_SPACEDIM], at_high_wall[AMREX_SPACEDIM];
                    for (int idir = 0; idir < AMREX_SPACEDIM; ++idir) {
                        at_low_wall[idir] = (iv[idir] < domlo[idir]);
                        at_high_wall[idir] = (iv[idir] > domhi[idir]);
                    }
                    
                    // Top wall (j = domhi[1]) is the moving lid
                    if (at_high_wall[1]) {
                        // Moving wall BC: set to equilibrium with lid velocity
                        Real rho = 1.0;
                        Array1D<Real, 0, AMREX_SPACEDIM> vel_bc{ulb, 0.0, 0.0};
                        data(iv, q) = feq(q, rho, vel_bc);
                    } else {
                        // Stationary wall BC: bounce-back
                        data(iv, q) = data(neighbor, opp[q]);
                    }
                }
            }
        });
    }
};

int main(int argc, char *argv[]) {
    // ...
    {
        // ...
        // Boundary Conditions - set to ext_dir for all boundaries
        Vector<BCRec> pop_bcs(POPULATION);

        for (int idim = 0; idim < AMREX_SPACEDIM; ++idim) {
            for (auto &bc : pop_bcs) {
                bc.setLo(idim, BCType::ext_dir);  // External Dirichlet BC
                bc.setHi(idim, BCType::ext_dir);  // External Dirichlet BC
            }
        }

        // Create PhysBCFunct with the custom boundary condition functor
        CavityBCFill cavity_bc{params.ulb};
        PhysBCFunct<CavityBCFill> physbcf(geometry, pop_bcs, cavity_bc);
        
        // ...
        
        for (int step = 1; step < max_iter + 1; ++step) {
            // ...
            physbcf(fin, 0, POPULATION, fin.nGrowVect(), time, 0);
            // ...
        }
    }
    
    Finalize();
    return 0;

[AlexanderSinn]

How are you defining the global variables c, w and opp? This could lead to performance issues if they are placed in the wrong type of memory.

[AlexanderSinn]

I guess it is obvious, but the ParallelFor is iterating over all the cells of the domain but only does something at the edge. Instead, up to six ParallelFors could be used to only iterate over the cells outside the domain. You can do const Box inside = geom.Domain(); directly by the way.