Skip to content

Modifying field values for neighboring cells inside a postprocessor

OpenLB – Open Source Lattice Boltzmann Code Forums on OpenLB General Topics Modifying field values for neighboring cells inside a postprocessor

Viewing 15 posts - 1 through 15 (of 40 total)
  • Author
    Posts
  • #9042
    Danial.Khazaeipoul
    Participant

    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

    #9043
    Adrian
    Keymaster

    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”)

    #9044
    Danial.Khazaeipoul
    Participant

    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

    #9045
    Adrian
    Keymaster

    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.

    #9046
    Danial.Khazaeipoul
    Participant

    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);
    }

    #9051
    Danial.Khazaeipoul
    Participant

    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

    #9056
    Adrian
    Keymaster

    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.

    #9059
    Danial.Khazaeipoul
    Participant

    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

    #9061
    Danial.Khazaeipoul
    Participant

    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.

    #9078
    Danial.Khazaeipoul
    Participant

    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?

    #9079
    Adrian
    Keymaster

    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.

    #9095
    Danial.Khazaeipoul
    Participant

    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.

    #9100
    Danial.Khazaeipoul
    Participant

    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.

    #9106
    Danial.Khazaeipoul
    Participant

    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.

    #9107
    Adrian
    Keymaster

    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.

Viewing 15 posts - 1 through 15 (of 40 total)
  • You must be logged in to reply to this topic.