Modifying field values for neighboring cells inside a postprocessor

[Danial.Khazaeipoul]

Dear community,

I would like to be able to modify the field values of neighboring cells inside a post-processor if certain conditions are met. For instance, in one of the stages of the “freeSurfacePostProcessor3D”. As far as I understand, a reference to the “cell” interface is passed as an argument of the template functions as below:

template <typename T, typename DESCRIPTOR>
template <typename CELL, typename PARAMETERS>
void FreeSurfaceMassFlowPostProcessor3D<T, DESCRIPTOR>::apply(CELL& cell, PARAMETERS& vars)
{}

Here, we can modify a field value using the setField method provided by the “Cell” class as follow:

cell.template setField<FreeSurface::MASS>(mass_tmp);

Apart from this, I also need to modify the neighboring values of a user-defined field. I understand that we can access the neighbor cells using the “neighbor()” method. However, I do not understand how to achieve this correctly as the below approach gives a compile time error indicating that I am trying to assign a const value to a reference:

CELL& nbrCell = cell.neighbor(direction);

Would you please provide a solution to this problem or pointing me to the source files needed to be read to understand write-access method provided for such scenarios?

Best regards,
Danial

[Adrian]

This will work if you change it to the correct type returned by CELL::neighbor. Note that this is not necessarily the same type as CELL. I suggest you use auto neighborCell = cell.neighbor(c_i);.

You may also be interested in the core and FAQ sections of the user guide (“how to write your own post processor”)

[Danial.Khazaeipoul]

Thank you for the prompt response, Adrian. Following your suggestion, is then possible to modify the field value of that neighbor cell as below with that same post-processor?

auto nbrCell = cell.neighbor(direction);
nbrCell.template setField<FreeSurface::MASS>(mass_tmp);

By reviewing the cell.h and cell.hh files, it appears to me that the “neighbor()” method returns either of the following types depending on the circumstances:

(1) ConstCell<T,DESCRIPTOR>: Does not provide a setField method
(2) Cell<T,DESCRIPTOR> : Provides a setField method

[Adrian]

Exactly, you can modify the fields of the neighbor cells in this way.

The (Const)Cell implemented in core/cell.h(h) is just one of multiple implementations of the “cell concept” in OpenLB. Specifically, it will not be used for the actual application of post processors on any of our target platforms. Instead, more efficient implementations of the concept are provided there.

The next release will include C++20 concepts declaring explicitly what a cell must provide. In the meantime you can look at the core cell header for the interface but not the concrete types.

[Danial.Khazaeipoul]

Thank you very much for the details, I managed to successfully modify the neighboring cells data using the method you mentioned.

Further in the code, I need to store the neighboring cells of a cell in a container for further processing, if a certain condition is met. I have tried the following approach with no luck. While the it compiles and runs, the value of “FreeSurface::TEST_FIELD” remains unchanged for that specific cell, although “setField” is called. Here is the simplified code to demonstrate the idea:

template <typename CELL, typename T>
void runTest(CELL& cell, const std::uint16_t& seed)
{
std::queue<CELL*> queued_;
queued_.push(&cell);
CELL* iCellPtr = queued_.front();
iCellPtr->template setField<FreeSurface::TEST_FIELD>(seed);
queued_.pop();
}

FYI, the following version works. However, I need a way to store the cells in a container for the algorithm that I am trying to implement.

template <typename CELL, typename T>
void runTest(CELL& cell, const std::uint16_t& seed)
{
cell.template setField<FreeSurface::TEST_FIELD>(seed);
}

[Danial.Khazaeipoul]

Dear Adrian,

Would you please explain why the two following implementations behave differently? In the first implementation, although seed is only a copy of cell, the “TEST_ID” field is modified correctly.

template <typename CELL, typename T>
void test(CELL& cell)
{
auto seed = cell;
seed.template setField<FreeSurface::TEST_ID>(1);
}

On the other hand, the below implementation, although compiles and runs with no issue, but the “TEST_ID” field remains unchanged. I have also tried different methods to be able to store these cell objects in a standard container, e.g., copying, raw pointers, smart pointers, etc, with no luck.

template <typename CELL, typename T>
void test(CELL& cell)
{
std::queue<CELL> queued_;
queued_.push(cell);
auto iCell = queued_.front();
iCell.template setField<FreeSurface::TEST_ID>(1);
}

Best,
Danial

• This reply was modified 1 month, 1 week ago by Danial.Khazaeipoul.

[Adrian]

Your second test implementation just instantiates a queue on the function’s stack, adds an item and throws the entire thing away again on function return. This would be the same for any C++ function written in this way and has nothing todo with OpenLB. If you want to preserve state outside of the function you need to do so explicitly, e.g. by modifying a global variable (not recommended outside of basic testing, see the remainder of this comment).

Updating a global shared structure on each operator application doesn’t mesh well with the parallelization requirements of LBM. Ideally you only update state local to each cell (otherwise you will, at a minium, need to think about possible parallelization conflicts). The restriction of the operator interface to individual cells (and optionally a set of parameters) is by design in other to guide model development towards parallelization friendly and maintainable code.

What exactly do you want to achieve with this list you want to create during post processor application? There may be a better way to do what you want. If not, we can talk about suitable parallelization friendly ways of achieving your requirements.

[Danial.Khazaeipoul]

Thank you Adrian, for the detailed explanation and also the important notes on parallelization aspects.

I am working on a verification case involving a single bubble rising in a liquid column, using free-surface dynamics. The example works so far, with reasonable errors compared to experimental data. However, OpenLB’s free-surface implementation lacks a bubble model that enforces constant bubble pressure, which may explain the discrepancies observed. I am now attempting to implement a bubble model, starting with detecting bubbles in the domain by assigning a unique tag to cells belonging to each bubble. To achieve this, I plan to implement an algorithm like flood-fill to detect connected regions (neighboring gas and interface cells) and store their IDs in a global shared structure. Here’s the pseudo-code outlining the algorithm:

// Start tagging process
int currentTag = 1;

for each cell in cells
{
if ((cell.isGasOrInterface()) && (!cell.hasTag()))
{
// Add the cell to the queue and assign it a unique tag
queue.push(cell);

// Process the queue
while (!queue.empty())
{
// Get the next cell from the queue
Cell currentCell = queue.front();
currentCell.setTag(currentTag);
queue.pop();

// Get the face neighbors of the current cell
for each neighbor in currentCell.getFaceNeighbors()
{
if ((neighbor.isGasOrInterface()) && (!neighbor.hasTag()))
{
// If the neighbor is gas or interface and has no tag, add to queue
queue.push(neighbor);
}
}
}

// Increment the tag for the next group of cells
currentTag++;
}
}

So the purpose of the list is to enable searching for the connected cells, if they are gas or interface.

Best regards,
Danial

• This reply was modified 1 month, 1 week ago by Danial.Khazaeipoul.

[Danial.Khazaeipoul]

I was just wondering that the flood-fill algorithm can be implemented in both recursive or iterative method. Given your explanation, I believe if I try to replace my current iterative approach with the recursive one, it eliminates the need for a queueing container.

It would still be very appreciated if you could explain the best approach to implement the iterative approach which seems to offer better performance while also using less memory.

[Danial.Khazaeipoul]

Dear Adrian,

I have tried to test my flood-fill implementation outside of OpenLB’s framework, normal C++ program running on an arbitrary dataset, and it works as expected.

Implementing algorithms such as flood-fill requires data modification of the neighboring cells of a given cell, and there is no way around this. Is it possible for you to give me a hint on the best approach to implement such algorithms in OpenLB? Or maybe a similar algorithm already implemented in OpenLB which I can go through to get ideas?

[Adrian]

It is perfectly fine to change values in neighboring cells from inside OpenLB operators. However, you are responsible for ensuring that there either A) are no access conflicts due to (in principle) all cells being processed at the same time or B) the access conflicts do not impact the result (note that this is highly probabilistic). This is a general parallel programming issue that you need to think about in any case.

If a (as I suspect you did in c++) sequential implementation is fine from a performance perspective that you can duplicate this in OpenLB as a simple manual loop over the relevant cells.

[Danial.Khazaeipoul]

Yes, I have implemented the stand-alone version using a sequential method, but since the algorithm needs to be executed each time-step, the performance would not be optimal. This is specially true as dealing with bubbles requires an adequate lattice resolution.

I have also found a problem with my previous attempts. As I am new to GPU programming, I learned that you can not use standard C++ containers on device codes. However, the following containers work correctly when compiled for GPU:

1) std::array
2) olb::Vector

Meanwhile, I am still unable to store pointers to a CELL inside a container, which I think would probably work if I compile the code for CPU instead of GPU.

[Danial.Khazaeipoul]

FYI, as also stated in the user-guide, the standard C++ std::vector container does not work with CUDA but it is still used within the “freeSurfaceHelpers.hh” where the calculateCubeOffset function is defined.

template<typename T, typename DESCRIPTOR>
static T calculateCubeOffset(T volume, const Vector<T,DESCRIPTOR::d>& normal) any_platform;

template<typename T, typename DESCRIPTOR>
T calculateCubeOffset(T volume, const Vector<T,DESCRIPTOR::d>& normal) {
std::vector<T> abs_normal(DESCRIPTOR::d, T{0});
for(int i = 0; i < DESCRIPTOR::d; i++){
abs_normal[i] = util::abs(normal[i]);
}

T volume_symmetry = 0.5 – util::abs(volume – 0.5);

std::sort(abs_normal.begin(), abs_normal.end());

if constexpr (DESCRIPTOR::d == 2) {
abs_normal[0] = util::max(normal[0], 1e-5);
} else if (DESCRIPTOR::d == 3){
abs_normal[0] = util::max(normal[0], 1e-12);
abs_normal[1] = util::max(normal[1], 1e-12);
}

T d = offsetHelper<T,DESCRIPTOR>(volume_symmetry, abs_normal);

T sorted_normal_acc = 0;
for(int i = 0; i < DESCRIPTOR::d; i++){
sorted_normal_acc += abs_normal[i];
}

return std::copysign(d – 0.5 * sorted_normal_acc, volume – 0.5);
}

P.S. this is the original code from olb-1.7r0.

[Danial.Khazaeipoul]

Good news for anyone interested, I managed to finally store pointers to “cells” within a container. The critical issue is to avoid using std::vector, queue, or stack when the platform is CUDA. Use either std::array or olb::Vector. The only downside is that with std::array or olb::Vector, you need to set the size of the container beforehand. Here is the sample test code that works on both CPU and GPU.

template <typename CELL, typename V>
void Test(CELL& cell, std::uint32_t& seed)
{
Vector<CELL*, 2> queued = {};
queued[0] = &cell;
auto neighborCell = queued[0]->neighbor(descriptors::c<DESCRIPTOR>(2));
queued[1] = &neighborCell;
queued[1]->template setField<FreeSurface::TEST_ID>(seed);
}

Thank you Adrian for your assistance and clarification.

[Adrian]

Note that this is neither platform transparent nor memory save as the cell classes are only interfaces that are not permanently stored. What you want to do is store the cell IDs. (See the core section of the user guide)

I assume you want to preserve the queue across operator applications? If so you will need to pass a pointer to it as Parameter / use global data and take care of the parallel accesses yourself.

[Author]

Posts