What is the best way to create and access cell's flag ?

[ToguyC]

Hello,

I am very new to AMReX, and I am currently stuck on a problem.

I am working on a project where I need to change the computation based on the type of a neighboring cell. I am in a D3Q19 domain (for LBM) and I am currently using a MultiFab called flags that has only 1 component, containing the cell's type.

To clarify, I want to create the LBM roofline model where the lid has a constant fluid velocity. In order to correctly make the streaming step, I need to access all 19 neighbors and apply a specific update computation if the neighbor's flag is set as bounce back.

But I find unexpected behavior when checking the neighbor's flag during computation (e.g.: accessing cell at x=x-1, y=y, z=z-1).

My question is: What is the recommended way of creating and accessing cell's flags ?

[asalmgren]

@ToguyC -- have you created the MultiFab with ghost cells (you'll need to have at least one ghost cell if you want to check nearest neighbors)? Maybe you can share the MFIter / ParallelFor loops where you're checking the value of the cell?

[ToguyC]

Yes I have set 1 ghost cell (as you said necessary for neighbor access).
The lattice I use stores 19*2 components to only use 1 array for old and new array, with the parity value switching back and forth to "swap" between them

The initialization of the MultiFabs is as follows:

MultiFab lattice(box_array, distribution_mapping, POPULATION * 2, 1);
MultiFab flags(box_array, distribution_mapping, 1, 1);

Here is the MFIter / ParallelFor in question:

for (MFIter mfi(lattice); mfi.isValid(); ++mfi) {
	const Box &box = mfi.validbox();
	auto lat = lattice.array(mfi);
	auto _flags = flags.array(mfi);

	ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
	  if (_flags(i, j, k) == static_cast<int>(CellType::bulk)) {
	      auto [rho, u] = macro(lat, i, j, k, parity);
	      Real usqr = 1.5 * (u(0) * u(0) + u(1) * u(1) + u(2) * u(2));

	      for (int l = 0; l < 9; ++l) {
			// clang-format off
			// auto [pop_out, pop_out_opp] = collideBgk(lat, i, j, k, l, rho, u, usqr, params.omega, parity);
			// clang-format on

			int x = i + c_vect(l, 0);
			int y = j + c_vect(l, 1);
			int z = k + c_vect(l, 2);

			// Voluntary assignation for debugging purposes
			if (_flags(x, y, z) == static_cast<int>(CellType::bounce_back)) {
				lat(i, j, k, (1 - parity) * POPULATION + opp_vect(l)) = 0;
			} else {
				lat(x, y, z, (1 - parity) * POPULATION + l) = 1;
			}
			// stream(lat, _flags, usqr, i, j, k, opp_vect(l), pop_out_opp, parity);
	      }

	      // Rest of population (c_l = {0,0,0})
	      for (int l : {9}) {
			Real eq = rho * w_vect(l) * (1. - usqr);
			// clang-format off
			lat(i, j, k, (1 - parity) * POPULATION + l) = (1. - params.omega) * lat(i, j, k, parity * POPULATION + l) + params.omega * eq;
			// clang-format on
	      }
	  }
	});
}

[asalmgren]

What is the unexpected behavior that you see?

[ToguyC]

I have a reference code written in plain C++, and I am comparing the output of the computation of energy of the system after each iteration:

template <class T>
AMREX_FORCE_INLINE AMREX_GPU_DEVICE GpuTuple<Real, Array1D<Real, 0, AMREX_SPACEDIM>> macro(T f, int i, int j,
                                                                             int k, int parity) {
    // clang-format off
    Real X_M1 = f(i, j, k, parity * POPULATION + 0) + f(i, j, k, parity * POPULATION + 3) + f(i, j, k, parity * POPULATION + 4) + f(i, j, k, parity * POPULATION + 5) + f(i, j, k, parity * POPULATION + 6);
    Real X_P1 = f(i, j, k, parity * POPULATION + 10) + f(i, j, k, parity * POPULATION + 13) + f(i, j, k, parity * POPULATION + 14) + f(i, j, k, parity * POPULATION + 15) + f(i, j, k, parity * POPULATION + 16);
    Real X_0  = f(i, j, k, parity * POPULATION + 9) + f(i, j, k, parity * POPULATION + 1) + f(i, j, k, parity * POPULATION + 2) + f(i, j, k, parity * POPULATION + 7) + f(i, j, k, parity * POPULATION + 8) + f(i, j, k, parity * POPULATION + 11) + f(i, j, k, parity * POPULATION + 12) + f(i, j, k, parity * POPULATION + 17) + f(i, j, k, parity * POPULATION + 18);

    Real Y_M1 = f(i, j, k, parity * POPULATION + 1) + f(i, j, k, parity * POPULATION + 3) + f(i, j, k, parity * POPULATION + 7) + f(i, j, k, parity * POPULATION + 8) + f(i, j, k, parity * POPULATION + 14);
    Real Y_P1 = f(i, j, k, parity * POPULATION + 4) + f(i, j, k, parity * POPULATION + 11) + f(i, j, k, parity * POPULATION + 13) + f(i, j, k, parity * POPULATION + 17) + f(i, j, k, parity * POPULATION + 18);

    Real Z_M1 = f(i, j, k, parity * POPULATION + 2) + f(i, j, k, parity * POPULATION + 5) + f(i, j, k, parity * POPULATION + 7) + f(i, j, k, parity * POPULATION + 16) + f(i, j, k, parity * POPULATION + 18);
    Real Z_P1 = f(i, j, k, parity * POPULATION + 6) + f(i, j, k, parity * POPULATION + 8) + f(i, j, k, parity * POPULATION + 12) + f(i, j, k, parity * POPULATION + 15) + f(i, j, k, parity * POPULATION + 17);

    Real rho = X_M1 + X_P1 + X_0;
    Array1D<Real, 0, AMREX_SPACEDIM> u{ (X_P1 - X_M1) / rho, (Y_P1 - Y_M1) / rho, (Z_P1 - Z_M1) / rho };
    // clang-format on

    return {rho, u};
}

Real computeEnergy(MultiFab const &f_in, MultiFab const &flags, int parity) {
    auto const &f_in_arrays = f_in.arrays();
    auto const &flags_arrays = flags.arrays();

    Real energy = ParReduce(
        TypeList<ReduceOpSum>{}, TypeList<Real>{}, f_in, IntVect(0),
        [=] AMREX_GPU_DEVICE(int box_no, int i, int j, int k) noexcept -> Real {
            auto f = f_in_arrays[box_no];
            auto _flags = flags_arrays[box_no];

            if (_flags(i, j, k) == static_cast<int>(CellType::bounce_back)) {
                return 0;
            }

            auto [rho, u] = macro(f, i, j, k, parity);
            Real _energy = u(0) * u(0) + u(1) * u(1) + u(2) * u(2);

            return _energy;
        });

    return 0.5 * energy;
}

As a result, the reference code outputs a value which diverge only by a little amount from the AMReX code.

I mainly think that my lack of understanding the framework is at the root of the problem, thus asking here if I have introduced any obvious mistakes for this kind of computation.

[WeiqunZhang]

Note that ParReduce does not do MPI communication. If you are not doing it somewhere else, you should do it here.
See the end of this. https://amrex-codes.github.io/amrex/docs_html/GPU.html?highlight=parreduce#multifab-reductions
(We need to explicit mention this in Doxygen documentation.)

Instead of MultiFab for flags, you might want to use iMultiFab.

[ToguyC]

is iMultiFab used for integer values instead of real ?

[WeiqunZhang]

Right. iMultiFab's underlying type is int. You could also use your own type, FabArray<BaseFab<CellType>>.

[asalmgren]

I'm not seeing anything obvious. Have you checked that you have correctly set the values of flag in the ghost cells?

[ToguyC]

That might be part of the issue. I have only assigned the values of the flags in a ParallelFor when initializing both the lattice and the flags, but then only called FillBoundary() at each iteration, thinking that it automatically copies over the values in the ghost cells when required.

[asalmgren]

FillBoundary will fill interior ghost cells, i.e. ones that overlap another Box in the BoxArray on which the MultiFab is defined. But it won't fill values outside of that, which are presumably outside the domain (unless you call FillBoundary(geom.periodicity() -- but then it will only fill ghost cells at periodic boundaries). I'm thinking that since you mentioned imposing a lid, you will need the ghost cells outside the physical domain?

[ToguyC]

I am using an array of the same size as the lattice, and apply a fixed formula for the cells located on the edges of the domain. So I don't think I will need such thing.

[WeiqunZhang]

What's the index type of your MultiFabs? Nodal or cell centered?

[ToguyC]

I have kept the default settings for that, so I presume centered.

If you want the whole context, here is the initialization:

        SimulationParameters params;
        readSimulationParameters(&params, true);

        BoxArray box_array; // List of boxes that cover the domain
        Geometry geometry;  // Information like physical domain size,
                            // number of points in the domain, periodicity

        IntVect domain_limit_low(0, 0, 0);
        IntVect domain_limit_high(params.n_cell - 1, params.n_cell - 1, params.n_cell - 1);

        Box domain(domain_limit_low, domain_limit_high);

        // Set the domain and break ti up into chunks no larger than
        // "max_grid_size" along a direction
        box_array.define(domain);
        box_array.maxSize(params.max_grid_size);

        RealBox real_box({0., 0., 0.}, {1., 1., 1.});
        Array<int, AMREX_SPACEDIM> is_periodic{1, 1, 1};

        geometry.define(domain, real_box, CoordSys::cartesian, is_periodic);

        // Extract the delta X, Y, Z of the domain given the parameters
        GpuArray<Real, AMREX_SPACEDIM> dx = geometry.CellSizeArray();

        // How boxes are distributed among MPI processes
        DistributionMapping distribution_mapping(box_array);

        // Would be changed for the f_in and f_out in LBM
        MultiFab lattice(box_array, distribution_mapping, POPULATION * 2, 1);
        MultiFab flags(box_array, distribution_mapping, 1, 1);

PS.: Thank you for your help with this seemingly trivial aspect haha

[ToguyC]

I might be better off rewritting from scratch trying not to include the bug.
If there is no obvious mistake, it can only be induced by a programming error I written somewhere else I guess.

[ToguyC]

Hello again !

I finally found out what the cause was. Whenever I use a max_grid_size different than the number of cell along a direction, I have errors coming in due to the ghost cells update.

I update the value of some neighbors based on some condition, and when the neighbor is a ghost cell, it causes issue.

I do not really know how to handle such case in AMReX though. Is there a way to merge back the ghost cells on anything like that ? Sorry if I missed it in the documentation.

[WeiqunZhang]

What are the possible values of c[l][]?

[ToguyC]

AMREX_GPU_MANAGED IntVect c[POPULATION] = {
        {-1, 0, 0}, { 0,-1, 0}, { 0, 0,-1},
        {-1,-1, 0}, {-1, 1, 0}, {-1, 0,-1},
        {-1, 0, 1}, { 0,-1,-1}, { 0,-1, 1},
        { 0, 0, 0},
        { 1, 0, 0}, { 0, 1, 0}, { 0, 0, 1},
        { 1, 1, 0}, { 1,-1, 0}, { 1, 0, 1},
        { 1, 0,-1}, { 0, 1, 1}, { 0, 1,-1}
    };

[WeiqunZhang]

I am not familiar with the algorithm. Maybe you can do something like

    for (MFIter mfi(lattice); mfi.isValid(); ++mfi) {
        const Box &box = amrex::grow(mfi.validbox(),1);
        auto lat = lattice.array(mfi);
        auto _flags = flags.array(mfi);

        ParallelFor(box, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
            if (_flags(i, j, k) == CellType::bulk) {
                for (int l = 0; l < 1; ++l) {
                    int x = i + c[l][0];
                    int y = j + c[l][1];
                    int z = k + c[l][2];
                    if (lat.contains(x,y,z)) {
                        lat(x, y, z, (1 - parity) * POPULATION + l) = 1;
                    }
                }
            }
        });
    }

Note that the issue you have has nothing to do with AMReX. If you have MPI, you have will the same issue. I am not if my suggestion is the right approach. So maybe you can do some research on how domain decomposition is done in this type of algorithm.

[ToguyC]

Thank you for your help. I narrowed down the source of the problem, so it might be easier to fix now.
Agin, thank you for your help.

[WeiqunZhang]

You are welcome. Basically you need to check box+1 region of flags to update box region of lat. And don't forget about mpi reduce in computeEnergy once you use MPI.