|
[Sponsors] |
Simplifying Hilbert curves for domain decomposition |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
August 21, 2022, 07:46 |
Simplifying Hilbert curves for domain decomposition
|
#1 |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
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! |
|
August 21, 2022, 13:05 |
|
#2 |
Senior Member
|
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 |
|
August 21, 2022, 15:11 |
|
#3 |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
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. |
|
August 21, 2022, 15:34 |
|
#4 | |
Super Moderator
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49 |
Quote:
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. |
||
August 21, 2022, 21:42 |
|
#5 | |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
Quote:
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. |
||
August 22, 2022, 08:32 |
|
#6 |
Super Moderator
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49 |
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. |
|
September 22, 2022, 21:05 |
|
#7 | |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
Quote:
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! |
||
September 22, 2022, 21:21 |
|
#8 | |
Senior Member
|
Quote:
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 |
||
September 22, 2022, 23:24 |
|
#9 | |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
Quote:
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. |
||
September 23, 2022, 04:57 |
|
#10 | |
Super Moderator
Alex
Join Date: Jun 2012
Location: Germany
Posts: 3,427
Rep Power: 49 |
Quote:
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. |
||
September 23, 2022, 05:26 |
|
#11 |
Senior Member
Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,285
Rep Power: 34 |
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?
|
|
September 23, 2022, 05:59 |
|
#12 | |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
Quote:
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. |
||
September 23, 2022, 11:24 |
|
#13 | |
Senior Member
|
Quote:
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) |
||
September 23, 2022, 11:40 |
|
#14 | |
Senior Member
Sayan Bhattacharjee
Join Date: Mar 2020
Posts: 495
Rep Power: 8 |
Quote:
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. |
||
September 23, 2022, 11:45 |
|
#15 |
Senior Member
|
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.
|
|
|
|
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 |