CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > General Forums > Main CFD Forum

Simplifying Hilbert curves for domain decomposition

Register Blogs Community New Posts Updated Threads Search

Like Tree8Likes
  • 1 Post By sbaffini
  • 1 Post By flotus1
  • 1 Post By flotus1
  • 1 Post By sbaffini
  • 1 Post By aerosayan
  • 1 Post By flotus1
  • 1 Post By sbaffini
  • 1 Post By sbaffini

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 21, 2022, 07:46
Default Simplifying Hilbert curves for domain decomposition
  #1
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Hello everyone,

I implemented Morton ordering (z curve) few years ago for a rudimentary domain partitioning and decomposition. But Morton ordering has a lot of complicated maths which I don't want to include in my work.

I'm thinking of using something like Hilbert curves ( https://en.wikipedia.org/wiki/Hilbert_R-tree ) for doing this, at it seems simpler, but still it's a little bit more complicated than I want it to be right now.

Here's my idea, kindly let me know if it looks good, or if there's some better alternatives which I don't know about :

I'm thinking of doing a quadtree based search for sorting my cells. That is, until each node of the quadtree has at max 2 or 4 cells, we'll keep subdividing the quadtree and sorting the elements.

I'm thinking of sorting the elements in this manner :

For grid cell in each node of our quadtree, we assign a number/marker out of set {1,2,3,4} to them, based on if they're in the 1st,2nd,3rd or 4th "quadrant" of each quadtree node. Then we sort the whole array based on the marker.

In real code, we'll try to use the Hilbert ordering to assign the quadrant marker, so our decomposition will be fractal in nature, but for simplicity, I'm just not describing everything. I'm not sure if we can get a true fractal decomposition.

This will sort the cells in a way that the cells in the same quadrant will be nearby to the other cells in the same quadrant. A data transformation of a list of cells marked like [1,2,4,3,2,1,3,2,4,2] to [1,1,2,2,2,2,3,3,4,4] will happen.

Then, we recursively repeat the same process for each quadrant of the quadtree.

The computation complexity of the algorithm is probably going to be dominated by O(Nlog(N)), i.e the initial sort algorithm for the whole mesh. Then we'll progressively only work on smaller sections of mesh.

Thanks!
Attached Images
File Type: png 106.png (86.5 KB, 22 views)
File Type: gif Figure2_Hilbert.gif (7.6 KB, 16 views)
aerosayan is offline   Reply With Quote

Old   August 21, 2022, 13:05
Default
  #2
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
If I understand what you mean, the SFC ordering is just meant to do that. It basically assigns an index from a scaled set of coordinates.

Given a cuboid that would include your whole domain and with each side split in a power of 2 number of cells, you first identify the index of the cell where your point (e.g., cell center) falls in and then assign it an index representing the order in which the cell would be visited by the particular SFC in use (Morton, Hilbert, etc).

Now, according to the chosen resolution, you might end up having multiple points in the same cell, but that's not a real issue, as long as they are O(1). That is, they being in the same cell means they are close, which is what you need in the end.

The point of using a SFC specific machinery is that it is efficient, you don't actually build a tree to assign the index. For each point you have, basically, an O(logC) cost, where C is linked to the chosen, fixed, resolution.

Of course, to be useful, you finally need to sort the cells according to their index, which typically costs less than assigning the SFC index.

The main difficulty I found in using the Hilbert SFC is that I had to build myself the matrices for the states in order to correct an error in the material I found. That required me some full time dedication, but once you get it, it should be straightforward.

EDIT: I don't think that SFC can get easier than Morton ordering. This seems especially so with this R-tree variant of the Hilbert ordering
aerosayan likes this.
sbaffini is offline   Reply With Quote

Old   August 21, 2022, 15:11
Default
  #3
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
I don't think that SFC can get easier than Morton ordering.
Morton ordering is definitely easy. I was trying to avoid it due to how much complicated bit operations are required, and i don't fully understand everything about it. Hilbert curve has a really good physical interpretation, and seemed easy at the moment.

But morton ordering also is limited by the max number of cells in the x and y direction. We can't morton order the whole mesh at the minimum resolution (i.e min cell size). At least not in one pass.

We would need to recursively morton order different sections of the mesh, so all the different sized cells are well distributed in memory.

You may be right that we won't need a quadtree.

I'm thinking that when we morton order the whole mesh at a coarse resolution, and sort the cells, the fine grained cells (let's say at the compression corner) will be close together in memory, but disorganized locally, and we would need to again morton order and sort those sections locally.

For an int64 morton index, we can create a pretty big grid for morton indexing, and if we can probably index even the finest resolution cells in a few steps/recursions.

Okay, maybe I'll try morton ordering again, but this time, try to make it more flexible and able to automatically and optimally order even the finest cells.
aerosayan is offline   Reply With Quote

Old   August 21, 2022, 15:34
Default
  #4
Super Moderator
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49
flotus1 has a spectacular aura aboutflotus1 has a spectacular aura about
Quote:
But morton ordering also is limited by the max number of cells in the x and y direction. We can't morton order the whole mesh at the minimum resolution (i.e min cell size). At least not in one pass.

We would need to recursively morton order different sections of the mesh, so all the different sized cells are well distributed in memory.
Are you sure that this is a relevant limitation for the grids you are working with? For a 64-bit intermediate z-order, you get 2^21 increments in the Cartesian directions in 3D. I have yet to run into a case where this is not enough.
As for what's more complicated: Without speed improvements, the bit manipulations required to get a z-index are really straightforward. I don't claim to fully understand the math behind the performance optimizations. But honestly: so what. It's not like Hilbert curves are intuitive to implement.
aerosayan likes this.
flotus1 is offline   Reply With Quote

Old   August 21, 2022, 21:42
Default
  #5
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
Are you sure that this is a relevant limitation for the grids you are working with? For a 64-bit intermediate z-order, you get 2^21 increments in the Cartesian directions in 3D. I have yet to run into a case where this is not enough.
As for what's more complicated: Without speed improvements, the bit manipulations required to get a z-index are really straightforward. I don't claim to fully understand the math behind the performance optimizations. But honestly: so what. It's not like Hilbert curves are intuitive to implement.
Could you please explain how you got the 2^21 value?

int64 max is 9223372036854775807 and the square root is 3037000499.98 i.e 3E+09.

You're right that if we use an int64 value, we can probably represent all meshes. I didn't consider that before. We create a Morton order's fake background grid of size 3e+9 in each direction in 2D. I made the wrong decision in the past to use an int32 whose square root of max value is 46340.95.

I didn't consider that before.

Maybe I'll use int64.
aerosayan is offline   Reply With Quote

Old   August 22, 2022, 08:32
Default
  #6
Super Moderator
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49
flotus1 has a spectacular aura aboutflotus1 has a spectacular aura about
IIRC, Morton order transformation works by inserting the bits of the Cartesian integer coordinates into the Z-value. One possible solution In 3D:
1st bit of xpos goes into 1st bit of the z-order
1st bit of ypos goes into 2nd bit of the z-order
1st bit of zpos goes into 3rd bit of the z-order
2nd bit of xpos goes into 4th bit of the z-order
...

After 21 bits in the position values, we have exhausted a 64-bit variable for the z-order.
In 2D, you get even more bits. That's probably where your square root approximation comes from.
I have been using this for 3D for a while now, and haven't even come close to exhausting the value range. Maybe if you are working with non-Cartesian grids, things are different. But yeah, using a 64-Bit variable for computing the intermediate order is kind of necessary for real-world applications in 3D.
aerosayan likes this.
flotus1 is offline   Reply With Quote

Old   September 22, 2022, 21:05
Default
  #7
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
IIRC, Morton order transformation works by inserting the bits of the Cartesian integer coordinates into the Z-value. One possible solution In 3D:
1st bit of xpos goes into 1st bit of the z-order
1st bit of ypos goes into 2nd bit of the z-order
1st bit of zpos goes into 3rd bit of the z-order
2nd bit of xpos goes into 4th bit of the z-order
...

After 21 bits in the position values, we have exhausted a 64-bit variable for the z-order.
In 2D, you get even more bits. That's probably where your square root approximation comes from.
I have been using this for 3D for a while now, and haven't even come close to exhausting the value range. Maybe if you are working with non-Cartesian grids, things are different. But yeah, using a 64-Bit variable for computing the intermediate order is kind of necessary for real-world applications in 3D.
Okay, I tried it, and the maths seems easy. I just want to run it by you first.

We want to find z-index value from a given (x,y) co-ordinate.

In Morton ordering, we create the z-index by interleaving the bits from 2 integers representing the x,y axes, respectively.

If we use a 64 bit z-index number, then our grid will be max(32bit unsigned int) * max(32 bit unsigned int) sized, i.e a (2^32) * (2^32) sized grid, i.e a 4294967296 * 4294967296 grid.


Those are big numbers, and I think should cover ANY grid in use.

Am I correct in my approach?


Thanks!
aerosayan is offline   Reply With Quote

Old   September 22, 2022, 21:21
Default
  #8
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by aerosayan View Post
Okay, I tried it, and the maths seems easy. I just want to run it by you first.

We want to find z-index value from a given (x,y) co-ordinate.

In Morton ordering, we create the z-index by interleaving the bits from 2 integers representing the x,y axes, respectively.

If we use a 64 bit z-index number, then our grid will be max(32bit unsigned int) * max(32 bit unsigned int) sized, i.e a (2^32) * (2^32) sized grid, i.e a 4294967296 * 4294967296 grid.


Those are big numbers, and I think should cover ANY grid in use.

Am I correct in my approach?


Thanks!
Those are indeed big numbers. But something that mprinkey here on this forum made me notice once about all this (which is more relevant in 3D, where you can only pack 3 2**21 indices in a 64 bit integer) is that, then what? If they are not enough, or even if they are but the grid is very clustered (in a non pathological way), you could at most end up having those close points with the same SFC index, and you have to decide how to order them. But, at that point, they are already close, and whatever order you choose is not less legit than the SFC one (if a deeper level was available).

Also note that, even for 2^21 levels, you are already hitting (or close to) the single precision typical of meshes, adding more levels shouldn't really help in any case
aerosayan likes this.
sbaffini is offline   Reply With Quote

Old   September 22, 2022, 23:24
Default
  #9
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
Those are indeed big numbers. But something that mprinkey here on this forum made me notice once about all this (which is more relevant in 3D, where you can only pack 3 2**21 indices in a 64 bit integer) is that, then what? If they are not enough, or even if they are but the grid is very clustered (in a non pathological way), you could at most end up having those close points with the same SFC index, and you have to decide how to order them. But, at that point, they are already close, and whatever order you choose is not less legit than the SFC one (if a deeper level was available).

Also note that, even for 2^21 levels, you are already hitting (or close to) the single precision typical of meshes, adding more levels shouldn't really help in any case

You would recursively apply the same ordering on previously sorted sections of mesh.

Let's say 1000 cells are close in proximity and have the same z-index. Then you could apply the Morton ordering only on those 1000 cells and reorder them.

This will work guaranteed because the reordering is local and even if it might cause the reordered grid to not be a "true" Morton ordered grid, it still works for our purpose of improving locality and thus, cache efficiency while acessing neighbor cells.
sbaffini likes this.
aerosayan is offline   Reply With Quote

Old   September 23, 2022, 04:57
Default
  #10
Super Moderator
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49
flotus1 has a spectacular aura aboutflotus1 has a spectacular aura about
Quote:
Originally Posted by aerosayan View Post
Okay, I tried it, and the maths seems easy. I just want to run it by you first.

We want to find z-index value from a given (x,y) co-ordinate.

In Morton ordering, we create the z-index by interleaving the bits from 2 integers representing the x,y axes, respectively.

If we use a 64 bit z-index number, then our grid will be max(32bit unsigned int) * max(32 bit unsigned int) sized, i.e a (2^32) * (2^32) sized grid, i.e a 4294967296 * 4294967296 grid.


Those are big numbers, and I think should cover ANY grid in use.

Am I correct in my approach?


Thanks!
Seems to check out. The fact that you can use the full range of 32Bit integers for position values in 2D, when using a 64Bit integer for the z-order, is why I had doubts that this approach is limiting.
Sure, you can always think of a mesh topology that breaks the limit with much fewer cells than 2^64. Like e.g. a really thin/long channel. But in those cases, your idea of recursively ordering parts of the mesh will help.
aerosayan likes this.
flotus1 is offline   Reply With Quote

Old   September 23, 2022, 05:26
Default
  #11
Senior Member
 
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,285
Rep Power: 34
arjun will become famous soon enougharjun will become famous soon enough
Can I ask why are you writing this. This i ask because i thought you are doing this because of finding the partition where the cell will belong? Or are you using this to perform paritioning?
arjun is online now   Reply With Quote

Old   September 23, 2022, 05:59
Default
  #12
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by arjun View Post
Can I ask why are you writing this. This i ask because i thought you are doing this because of finding the partition where the cell will belong? Or are you using this to perform paritioning?
Hopefully I understood your question correctly.

It's not only for partitioning.

Reordering the cells using z curve, makes the physically nearby cells to also be near by in memory. The z curve has a very efficient space filling pattern which allows us to map out 2d domain into a 1D linear space. (edit: To correct: Morton curve has a somewhat efficient pattern. Hilbert curve is better for improving locality, but they're harder to implement, and don't have the mathematical assurance of Morton curve)

Cells which are nearby in physical space will also be nearby in memory.

This improves the efficiency of our code in accessing the neighborhood cells.

Same treatment can be applied to nodes.

Partitioning using this, is just a secondary accidental benefit.

Last edited by aerosayan; September 23, 2022 at 08:00.
aerosayan is offline   Reply With Quote

Old   September 23, 2022, 11:24
Default
  #13
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by aerosayan View Post
You would recursively apply the same ordering on previously sorted sections of mesh.

Let's say 1000 cells are close in proximity and have the same z-index. Then you could apply the Morton ordering only on those 1000 cells and reorder them.

This will work guaranteed because the reordering is local and even if it might cause the reordered grid to not be a "true" Morton ordered grid, it still works for our purpose of improving locality and thus, cache efficiency while acessing neighbor cells.
Just out of curiosity, have you ever met a case where this was actually needed? I mean, order 1000s of cells with the same SFC index? Let me highlight again that in 3D this means that, for 2^21 bins (approx 2Mln) covering a domain along each direction, around 1000 points are closer between them than approximately 5e-7 times the domain size.

For example, that would happen with a 1Km resolution on earth on a domain sized 1 Astronomical Unit. Or, on a 100m domain size, with a resolution of 50 micrometers. I mean, this seems far beyond anything approachable today, at least for something that has sense (and increasing resolution so much just locally in a single spot doesn't)
aerosayan likes this.
sbaffini is offline   Reply With Quote

Old   September 23, 2022, 11:40
Default
  #14
Senior Member
 
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8
aerosayan is on a distinguished road
Quote:
Originally Posted by sbaffini View Post
Just out of curiosity, have you ever met a case where this was actually needed? I mean, order 1000s of cells with the same SFC index? Let me highlight again that in 3D this means that, for 2^21 bins (approx 2Mln) covering a domain along each direction, around 1000 points are closer between them than approximately 5e-7 times the domain size.

For example, that would happen with a 1Km resolution on earth on a domain sized 1 Astronomical Unit. Or, on a 100m domain size, with a resolution of 50 micrometers. I mean, this seems far beyond anything approachable today, at least for something that has sense (and increasing resolution so much just locally in a single spot doesn't)
One year ago, for 2D, I would've needed to use recursion, but I stopped the work back then to spend most of my time learning CFD, instead of trying to develop something with limited knowledge.

In 2D, I would've needed to do this, as I was using a 32 bit usinged integer for z-index.

Shifting to 64 bits has been enough to not need recursion anymore.

For 3D, I have another idea, which might work.

In C++, we can define a struct with 3X32bit ints to represent the indices in 3 dimension. Imagine that the 3 ints represents a big 96 bit software defined integer.

96 bit since 3X32 = 96 bit.

Then to calculate the z index, we could do bit manipulation to derive our 96 bit z-index.

Then we can essentially sort our cells based on this 96 bit z-index value.

We would need to create our own < (less than) operator and just pass it to std::sort as a lambda std::sort(cells.begin(), cells.end(), [](cell a, cell b){return a.zindex<b.zindex;}

Then you won't need recursion.

I did mention once that learning bit manipulation is useful.
aerosayan is offline   Reply With Quote

Old   September 23, 2022, 11:45
Default
  #15
Senior Member
 
sbaffini's Avatar
 
Paolo Lampitella
Join Date: Mar 2009
Location: Italy
Posts: 2,190
Blog Entries: 29
Rep Power: 39
sbaffini will become famous soon enoughsbaffini will become famous soon enough
Send a message via Skype™ to sbaffini
Quote:
Originally Posted by aerosayan View Post
In 2D, I would've needed to do this, as I was using a 32 bit usinged integer for z-index.

Shifting to 64 bits has been enough to not need recursion anymore.
Ok, I missed the 32bit part. In my 3D case I was indeed referring to the use of a 64 bit signed integer in Fortran.
aerosayan likes this.
sbaffini is offline   Reply With Quote

Reply


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Domain decomposition of 3D overset grid case tok3rat0r OpenFOAM Running, Solving & CFD 2 February 24, 2023 03:58
FoamExtend4.0 decomposition warnings? EternalSeekerX OpenFOAM Running, Solving & CFD 0 May 23, 2021 08:17
[ICEM] Attached geometry curves niablis ANSYS Meshing & Geometry 1 October 7, 2020 06:59
[snappyHexMesh] How to define to right point for locationInMesh Mirage12 OpenFOAM Meshing & Mesh Conversion 7 March 13, 2016 15:07
meshing of a compound volume in GMSH shawn3531 OpenFOAM 4 March 12, 2015 11:45


All times are GMT -4. The time now is 23:40.