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
- This topic has 53 replies, 3 voices, and was last updated 6 days, 2 hours ago by Danial.Khazaeipoul.
-
AuthorPosts
-
October 5, 2024 at 7:08 pm #9322Danial.KhazaeipoulParticipant
Thank you, Mathias, for the explanation. The code I’ve implemented works well in a single partition setup where
localCellID == globalCellID
. However, in configurations with multiple partitions or blocks,localCellID != globalCellID
, leading to non-unique IDs.For local cell IDs, which are of type
CellID
orstd::uint32_t
, there’s an advantage in CCL algorithms as these IDs are unique and remain consistent throughout the simulation. Because of this, I’m looking for a way to generate global cell IDs with similar characteristics. For example, I am thinking about a method like the following which creates a unique global CellID for each cell:std::uint32_t globalCellID = blockID × cellsPerBlock + localCellID
- This reply was modified 4 months, 1 week ago by Danial.Khazaeipoul.
October 5, 2024 at 11:06 pm #9326AdrianKeymasterYes, the cell IDs are block-local. Computing a global ID in the way you posted will work if you use the maximum global block size for the multiplier.
- This reply was modified 4 months, 1 week ago by Adrian.
October 8, 2024 at 8:57 pm #9333Danial.KhazaeipoulParticipantJust to confirm my understanding, by
if you use the maximum global block size
Are you referring to the fact that cuboids may not have an equal number of associated cells, and that using the maximum block size can prevent overlapping global cell IDs?
October 8, 2024 at 10:50 pm #9334AdrianKeymasterYes, the cuboids commonly do not have exactly the same number of cells. You can also see this in the cuboid geometry output.
October 11, 2024 at 7:12 pm #9362Danial.KhazaeipoulParticipantUsing the cell interface, is there any method that can be used to identify the global cuboid number to which the current cell belongs?
October 11, 2024 at 7:16 pm #9363AdrianKeymasterNo, this is not intended as all individual threads of a kernel always only process cells of the same block (all LBM processing is split into per block operators). You should be able to handle the cuboid number / cell ID mapping at the same place where you perform the inter block communication of your algorithm.
However, in principle nothing prevents you from passing the cuboid number as a block-specific parameter to your operator (with application scope PerCellWithParameters)
November 5, 2024 at 8:38 pm #9474Danial.KhazaeipoulParticipantHello Adrian,
Is my understanding correct regarding the use of a block-specific parameter passed to an operator? Suppose I have an offset vector, synchronized across all processors, which records the number of cells in each preceding block, shown below.
struct CELL_COUNT : public descriptors::TYPED_FIELD_BASE<std::size_t,1> { };
std::vector<std::size_t> offset(_sLattice.getCuboidGeometry().getNc(), 0);
The size of the offset vector matches the number of global cuboids. For a domain with 4 blocks, the resulting offset vector might look like {0, 431433, 855297, 1279161}, as an example. Given this setup, will the following method yield the expected results for the
CELL_COUNT
, as shown in the image?for (int iC = 0; iC < _sLattice.getLoadBalancer().size(); ++iC) { auto& block = _sLattice.getBlock(iC); block.template setParameter<FreeSurface::CELL_COUNT>(offset[_sLattice.getLoadBalancer().glob(iC)]); }
Please note that the image is only for demonstration purposes. When an operator processes a cell, will calling the parameter use the correct value of
CELL_COUNT
based on the block in which the cell resides?- This reply was modified 3 months, 1 week ago by Danial.Khazaeipoul.
- This reply was modified 3 months, 1 week ago by Danial.Khazaeipoul.
- This reply was modified 3 months, 1 week ago by Danial.Khazaeipoul.
November 6, 2024 at 10:05 am #9478AdrianKeymasterYes, this should work as you convert the local cuboid index
iC
to the global one for accessing the offset array.November 20, 2024 at 1:02 am #9541Danial.KhazaeipoulParticipantDear Adrian,
I believe I’ve identified the issue with the CCL implementation when MPI is enabled. Thank you for your insights so far in working toward a robust solution. Since you’ve already reviewed the code, you’re familiar with how cell IDs are used to assign bubble IDs and to navigate between cells in each kernel to access their bubble IDs.
I’ll do my best to clearly explain the problem. The CCL algorithm assumes that cell IDs follow a consistent order in each direction within a block, which works fine for a single block. However, when MPI is enabled and the domain is divided into multiple blocks, this consistency breaks down for neighboring cells located in different blocks. As a result, the bubble ID assignment becomes incorrect for blocks other than Block 0. In the worst-case scenario, this can lead to illegal memory access, especially since blocks may have an unequal number of cells, see the below picture.
When the orange cell in Block #2 is being processed, it accesses the bubble ID of the red cell, as the red cell is a valid neighbor. However, since the red cell belongs to a different block, its cell ID does not follow the block-local order of cell IDs in Block #2. When the orange cell adopts the bubble ID of the red cell, it may incorrectly reference a cell within Block #2, as indicated by the blue arrow. This can result in either assigning an incorrect bubble ID to the orange cell or, in the worst-case scenario, causing illegal memory access. This happens because Block #2 may not contain a cell with the cell ID from Block #0.
Is it possible to defer the automatic communication between blocks for a specific global data structure until it is explicitly triggered when needed?
Regards,
Danial- This reply was modified 2 months, 3 weeks ago by Danial.Khazaeipoul.
November 24, 2024 at 1:54 pm #9555AdrianKeymasterSure, the overlap communication only happens at discrete and configurable stages. i.e. you can move all communication related to the CCL implementation into a distinct (or more than one) communicator stage fully within your control. See e.g. User Guide Section 2.2.3.
You can create a new communicator stage for CCL simply by requesting it:
namespace stage { struct CCL { }; } { auto& comm = sLattice.getCommunicator(stage::CCL{}); comm.template requestField<descriptors::POPULATION>(); comm.requestOverlap(1); comm.exchangeRequests(); }
this stage will only be triggered if you explicitly call
sLattice.getCommunicator(stage::CCL{}).communicate()
. It is also possible to request more than one field or other subsets of the overlap than just a constant width if necessary.December 12, 2024 at 10:37 pm #9604Danial.KhazaeipoulParticipantHello Adrian,
Can you confirm that the “apply” operators used in the post-processors operate only on the internal nodes? As I understand it, these operators cannot act directly on overlap nodes. The values of overlap nodes are updated only when communication is triggered using the communicate() method. However, it is still possible to access the values of an overlap node, such as when an internal node has an overlap node as its neighbor.
Is it possible to determine whether a neighbor of an internal node is an overlap node? I have found a temporary workaround to address my issue, as shown in the picture below. However, a cleaner and safer approach would be preferable as I need to avoid using the overlap nodes values until ids are finalized.
After the temporary fix is applied, as you can see, the connected regions on each block have their own unique ID.
December 12, 2024 at 11:22 pm #9605AdrianKeymasterYes, by default all post processors are only applied to the non-overlap nodes and attempts to explicitly set them on overlap nodes are filtered out. In the next release this will change as I also needed more control on this for my FSI implementations. Right now the most straight forward way would be to mark the overlap cells with a dedicated field.
If you need more control on the post processor placement right now I can probably also extract you a patch.
December 12, 2024 at 11:46 pm #9606Danial.KhazaeipoulParticipantThank you for the assistance Adrian.
Right now the most straight forward way would be to mark the overlap cells with a dedicated field.
Are you suggesting that I define a new dedicated field, where for instance overlap cells are assigned a value of “1” and all other cells are assigned a value of “0”? This field could then be used to distinguish overlap cells from regular cells.
If you need more control on the post processor placement right now I can probably also extract you a patch.
That would be greatly appreciated, as it would help me gain a better understanding of the code. Additionally, it could provide me with more control over the process.
December 18, 2024 at 11:29 am #9619AdrianKeymasterAre you suggesting that I define a new dedicated field, where for instance overlap cells are assigned a value of “1” and all other cells are assigned a value of “0”? This field could then be used to distinguish overlap cells from regular cells.
Yes, exactly (this is of course quite ugly but it will work).
If you want to add post processors to the overlap you can drop the
isPadding
check in line 618 ofcore/blockLattice.h
and add them per-cell.December 19, 2024 at 5:58 am #9620Danial.KhazaeipoulParticipantYes, exactly (this is of course quite ugly but it will work).
FYI, I implemented a cleaner method to replace this and now one can check whether a cell is a padding or not directly from the Cell interface. However, this required me to modify the src/core files to introduce a new function inside the BlockStructureD class that retrieves a local LatticR<D> from a linearized CellID.
FYI, I implemented a more streamlined method to replace this, allowing direct checks on whether a cell is padding via the Cell interface. To achieve this, I updated the src/core files to include a new function in the BlockStructureD class, which retrieves a local LatticeR<D> from a linearized CellID.
/// Get D-dimensional LatticeR from 1D cell ID LatticeR<D> getLatticeR(CellID iCell) const { if constexpr (D == 3) { std::int32_t iX = iCell / _projection[0]; std::int32_t remainder = iCell % _projection[0]; std::int32_t iY = remainder / _projection[1]; std::int32_t iZ = remainder % _projection[1]; return LatticeR<D>(iX - _padding, iY - _padding, iZ - _padding); } else { std::int32_t iX = iCell / _projection[0]; std::int32_t iY = iCell % _projection[0]; return LatticeR<D>(iX - _padding, iY - _padding); } };
/// Return whether iCell is valid bool isPadding(CellID iCell) const { return isPadding(getLatticeR(iCell)); };
For more details, please refer to the “feature-CCL” branch in my repository, specifically commits 2ce2d4fc and b28e5305.
-
AuthorPosts
- You must be logged in to reply to this topic.