Shamrock 2025.10.0
Astrophysical Code
Loading...
Searching...
No Matches
morton.hpp
Go to the documentation of this file.
1// -------------------------------------------------------//
2//
3// SHAMROCK code for hydrodynamics
4// Copyright (c) 2021-2026 Timothée David--Cléris <tim.shamrock@proton.me>
5// SPDX-License-Identifier: CeCILL Free Software License Agreement v2.1
6// Shamrock is licensed under the CeCILL 2.1 License, see LICENSE for more information
7//
8// -------------------------------------------------------//
9
10#pragma once
11
19#include "bmi.hpp"
20#include "shambackends/math.hpp"
21#include "shambackends/sycl.hpp"
22#include "shambackends/vec.hpp"
23#include "shammath/CoordRangeTransform.hpp"
24#include <type_traits>
25
26namespace shamrock::sfc {
27
28 template<class morton_t>
29 struct MortonInfo {
30 static constexpr morton_t err_code;
31 };
32
33 template<>
34 struct MortonInfo<u32> {
35 static constexpr u32 err_code = 4294967295U;
36 };
37
38 template<>
39 struct MortonInfo<u64> {
40 static constexpr u64 err_code = 18446744073709551615UL;
41 };
42
43 template<class Umorton, u32 dim>
44 class MortonCodes {};
45
46 template<>
47 class MortonCodes<u32, 3> {
48 public:
49 using int_vec_repr_base = u16;
50 using int_vec_repr = u16_3;
51 static constexpr int_vec_repr_base dimension = 3;
52 static constexpr int_vec_repr_base max_val = 1024 - 1;
53 static constexpr int_vec_repr_base val_count = 1024;
54 static constexpr int_vec_repr_base significant_bits_p_coord = 10;
55 static constexpr int_vec_repr_base significant_bits = dimension * significant_bits_p_coord;
56
57 static constexpr u32 err_code = 4294967295U;
58
59 inline static u32 icoord_to_morton(u32 x, u32 y, u32 z) {
60 u32 xx = bmi::expand_bits<u32, 2>((u32) x);
61 u32 yy = bmi::expand_bits<u32, 2>((u32) y);
62 u32 zz = bmi::expand_bits<u32, 2>((u32) z);
63 return xx * 4 + yy * 2 + zz;
64 }
65
66 inline static bool is_morton_bounding_box(int_vec_repr min, int_vec_repr max) noexcept {
67 return min.x() == 0 && min.y() == 0 && min.z() == 0 && max.x() == max_val
68 && max.y() == max_val && max.z() == max_val;
69 }
70
71 template<class flt>
72 inline static u32 coord_to_morton(flt x, flt y, flt z) {
73
74 constexpr bool ok_type = std::is_same<flt, f32>::value || std::is_same<flt, f64>::value;
75 static_assert(ok_type, "unknown input type");
76
77 if constexpr (std::is_same<flt, f32>::value) {
78
79 x = sycl::fmin(sycl::fmax(x * 1024.F, 0.F), 1024.F - 1.F);
80 y = sycl::fmin(sycl::fmax(y * 1024.F, 0.F), 1024.F - 1.F);
81 z = sycl::fmin(sycl::fmax(z * 1024.F, 0.F), 1024.F - 1.F);
82
83 return icoord_to_morton(x, y, z);
84
85 } else if constexpr (std::is_same<flt, f64>::value) {
86
87 x = sycl::fmin(sycl::fmax(x * 1024., 0.), 1024. - 1.);
88 y = sycl::fmin(sycl::fmax(y * 1024., 0.), 1024. - 1.);
89 z = sycl::fmin(sycl::fmax(z * 1024., 0.), 1024. - 1.);
90
91 return icoord_to_morton(x, y, z);
92 }
93 }
94
95 inline static u16_3 morton_to_icoord(u32 morton) {
96
97 u16_3 pos;
98 pos.s0() = (u16) bmi::contract_bits<u32, 2>((morton & 0x24924924U) >> 2U);
99 pos.s1() = (u16) bmi::contract_bits<u32, 2>((morton & 0x12492492U) >> 1U);
100 pos.s2() = (u16) bmi::contract_bits<u32, 2>((morton & 0x09249249U) >> 0U);
101
102 return pos;
103 }
104
105 inline static u16_3 get_offset(u32 clz_) {
106 u16_3 mx;
107 mx.s0() = 1024U >> ((clz_ - 0) / 3);
108 mx.s1() = 1024U >> ((clz_ - 1) / 3);
109 mx.s2() = 1024U >> ((clz_ - 2) / 3);
110 return mx;
111 }
112 };
113
114 template<>
115 class MortonCodes<u64, 3> {
116 public:
117 using int_vec_repr_base = u32;
118 using int_vec_repr = u32_3;
119 static constexpr int_vec_repr_base dimension = 3;
120 static constexpr int_vec_repr_base max_val = 2097152 - 1;
121 static constexpr int_vec_repr_base val_count = 2097152;
122 static constexpr int_vec_repr_base significant_bits_p_coord = 21;
123 static constexpr int_vec_repr_base significant_bits = dimension * significant_bits_p_coord;
124
125 static constexpr u64 err_code = 18446744073709551615UL;
126
127 inline static u64 icoord_to_morton(u64 x, u64 y, u64 z) {
128 u64 xx = bmi::expand_bits<u64, 2>((u64) x);
129 u64 yy = bmi::expand_bits<u64, 2>((u64) y);
130 u64 zz = bmi::expand_bits<u64, 2>((u64) z);
131 return xx * 4 + yy * 2 + zz;
132 }
133
134 inline static bool is_morton_bounding_box(int_vec_repr min, int_vec_repr max) noexcept {
135 return min.x() == 0 && min.y() == 0 && min.z() == 0 && max.x() == max_val
136 && max.y() == max_val && max.z() == max_val;
137 }
138
139 template<class flt>
140 inline static u64 coord_to_morton(flt x, flt y, flt z) {
141
142 constexpr bool ok_type = std::is_same<flt, f32>::value || std::is_same<flt, f64>::value;
143 static_assert(ok_type, "unknown input type");
144
145 if constexpr (std::is_same<flt, f32>::value) {
146
147 x = sycl::fmin(sycl::fmax(x * 2097152.F, 0.F), 2097152.F - 1.F);
148 y = sycl::fmin(sycl::fmax(y * 2097152.F, 0.F), 2097152.F - 1.F);
149 z = sycl::fmin(sycl::fmax(z * 2097152.F, 0.F), 2097152.F - 1.F);
150
151 return icoord_to_morton(x, y, z);
152
153 } else if constexpr (std::is_same<flt, f64>::value) {
154
155 x = sycl::fmin(sycl::fmax(x * 2097152., 0.), 2097152. - 1.);
156 y = sycl::fmin(sycl::fmax(y * 2097152., 0.), 2097152. - 1.);
157 z = sycl::fmin(sycl::fmax(z * 2097152., 0.), 2097152. - 1.);
158
159 return icoord_to_morton(x, y, z);
160 }
161 }
162
163 inline static int_vec_repr morton_to_icoord(u64 morton) {
164
165 u32_3 pos;
166 pos.x() = bmi::contract_bits<u64, 2>((morton & 0x4924924924924924U) >> 2U);
167 pos.y() = bmi::contract_bits<u64, 2>((morton & 0x2492492492492492U) >> 1U);
168 pos.z() = bmi::contract_bits<u64, 2>((morton & 0x1249249249249249U) >> 0U);
169
170 return pos;
171 }
172
173 inline static int_vec_repr get_offset(u32 clz_) {
174 u32_3 mx;
175 mx.s0() = 2097152U >> ((clz_ + 1) / 3);
176 mx.s1() = 2097152U >> ((clz_ - 0) / 3);
177 mx.s2() = 2097152U >> ((clz_ - 1) / 3);
178 return mx;
179 }
180 };
181
182 template<class morton_t, class _pos_t, u32 dim>
183 class MortonConverter {
184
185 public:
186 using Morton = MortonCodes<morton_t, dim>;
187
188 using pos_t = _pos_t;
189 using coord_t = typename shambase::VectorProperties<pos_t>::component_type;
190 using ipos_t = typename Morton::int_vec_repr;
191 using int_t = typename Morton::int_vec_repr_base;
192
193 using CoordTransform = shammath::CoordRangeTransform<ipos_t, pos_t>;
194
195 private:
196 static constexpr bool implemented_int = std::is_same<pos_t, u32_3>::value
197 || std::is_same<pos_t, u64_3>::value
198 || std::is_same<pos_t, i64_3>::value;
199
200 static constexpr bool implemented_float
201 = std::is_same<pos_t, f32_3>::value || std::is_same<pos_t, f64_3>::value;
202
203 static_assert(implemented_int || implemented_float, "not implemented");
204
205 public:
206 static CoordTransform get_transform(pos_t bounding_box_min, pos_t bounding_box_max) {
207 return CoordTransform(
208 shammath::CoordRange<ipos_t>{
209 {0, 0, 0}, {Morton::val_count, Morton::val_count, Morton::val_count}},
210 shammath::CoordRange<pos_t>{bounding_box_min, bounding_box_max});
211 }
212
213 inline static ipos_t to_morton_grid(pos_t pos, CoordTransform transform) {
214
215 ipos_t unit_coord = transform.reverse_transform(pos);
216
217 constexpr int_t zero = 0;
218
219 unit_coord.x() = sycl::min(sycl::max(unit_coord.x(), zero), Morton::max_val);
220 unit_coord.y() = sycl::min(sycl::max(unit_coord.y(), zero), Morton::max_val);
221 unit_coord.z() = sycl::min(sycl::max(unit_coord.z(), zero), Morton::max_val);
222
223 return unit_coord;
224 }
225
226 inline static pos_t to_real_space(ipos_t pos, CoordTransform transform) {
227
228 return transform.transform(pos);
229 }
230 };
231
232} // namespace shamrock::sfc
233
234//%Impl status : Good
235
236namespace morton_3d {
237
242 template<class morton_repr>
243 struct [[deprecated]] morton_types {
244 using int_vec_repr_base = std::void_t<>;
245 using int_vec_repr = std::void_t<>;
246 };
247
248 template<>
249 struct [[deprecated]] morton_types<u32> {
250 using int_vec_repr_base = u16;
251 using int_vec_repr = u16_3;
252 static constexpr int_vec_repr_base max_val = 1024 - 1;
253
254 // not possible yet in hipsycl
255 // static constexpr int_vec_repr max_vec = int_vec_repr{max_val};
256 };
257
258 template<>
259 struct [[deprecated]] morton_types<u64> {
260 using int_vec_repr_base = u32;
261 using int_vec_repr = u32_3;
262 static constexpr int_vec_repr_base max_val = 2097152 - 1;
263
264 // not possible yet in hipsycl
265 // static constexpr int_vec_repr max_vec = int_vec_repr{max_val};
266 };
267
268 template<class morton_prec, class fp_prec>
269 [[deprecated]] morton_prec coord_to_morton(fp_prec x, fp_prec y, fp_prec z);
270
271 template<class morton_prec>
272 [[deprecated]] typename morton_types<morton_prec>::int_vec_repr morton_to_ipos(
273 morton_prec morton);
274
275 template<class morton_prec>
276 [[deprecated]] typename morton_types<morton_prec>::int_vec_repr get_offset(u32 clz_);
277
278 template<>
279 inline u64 coord_to_morton<u64, f64>(f64 x, f64 y, f64 z) {
280 x = sycl::fmin(sycl::fmax(x * 2097152., 0.), 2097152. - 1.);
281 y = sycl::fmin(sycl::fmax(y * 2097152., 0.), 2097152. - 1.);
282 z = sycl::fmin(sycl::fmax(z * 2097152., 0.), 2097152. - 1.);
283
284 u64 xx = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) x);
285 u64 yy = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) y);
286 u64 zz = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) z);
287 return xx * 4 + yy * 2 + zz;
288 }
289
290 template<>
291 inline u64 coord_to_morton<u64, f32>(f32 x, f32 y, f32 z) {
292 x = sycl::fmin(sycl::fmax(x * 2097152.F, 0.F), 2097152.F - 1.F);
293 y = sycl::fmin(sycl::fmax(y * 2097152.F, 0.F), 2097152.F - 1.F);
294 z = sycl::fmin(sycl::fmax(z * 2097152.F, 0.F), 2097152.F - 1.F);
295
296 u64 xx = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) x);
297 u64 yy = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) y);
298 u64 zz = shamrock::sfc::bmi::expand_bits<u64, 2>((u64) z);
299 return xx * 4 + yy * 2 + zz;
300 }
301
302 template<>
303 inline u32 coord_to_morton<u32, f64>(f64 x, f64 y, f64 z) {
304 x = sycl::fmin(sycl::fmax(x * 1024., 0.), 1024. - 1.);
305 y = sycl::fmin(sycl::fmax(y * 1024., 0.), 1024. - 1.);
306 z = sycl::fmin(sycl::fmax(z * 1024., 0.), 1024. - 1.);
307
308 u32 xx = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) x);
309 u32 yy = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) y);
310 u32 zz = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) z);
311 return xx * 4 + yy * 2 + zz;
312 }
313
314 template<>
315 inline u32 coord_to_morton<u32, f32>(f32 x, f32 y, f32 z) {
316 x = sycl::fmin(sycl::fmax(x * 1024.F, 0.F), 1024.F - 1.F);
317 y = sycl::fmin(sycl::fmax(y * 1024.F, 0.F), 1024.F - 1.F);
318 z = sycl::fmin(sycl::fmax(z * 1024.F, 0.F), 1024.F - 1.F);
319
320 u32 xx = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) x);
321 u32 yy = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) y);
322 u32 zz = shamrock::sfc::bmi::expand_bits<u32, 2>((u32) z);
323 return xx * 4 + yy * 2 + zz;
324 }
325
326 template<>
327 inline u32_3 morton_to_ipos<u64>(u64 morton) {
328
329 u32_3 pos;
330 pos.x() = shamrock::sfc::bmi::contract_bits<u64, 2>((morton & 0x4924924924924924U) >> 2U);
331 pos.y() = shamrock::sfc::bmi::contract_bits<u64, 2>((morton & 0x2492492492492492U) >> 1U);
332 pos.z() = shamrock::sfc::bmi::contract_bits<u64, 2>((morton & 0x1249249249249249U) >> 0U);
333
334 return pos;
335 }
336
337 template<>
338 inline u16_3 morton_to_ipos<u32>(u32 morton) {
339
340 u16_3 pos;
341 pos.s0() = (u16) shamrock::sfc::bmi::contract_bits<u32, 2>((morton & 0x24924924U) >> 2U);
342 pos.s1() = (u16) shamrock::sfc::bmi::contract_bits<u32, 2>((morton & 0x12492492U) >> 1U);
343 pos.s2() = (u16) shamrock::sfc::bmi::contract_bits<u32, 2>((morton & 0x09249249U) >> 0U);
344
345 return pos;
346 }
347
348 template<>
349 inline u32_3 get_offset<u64>(uint clz_) {
350 u32_3 mx;
351 mx.s0() = 2097152U >> ((clz_ + 1) / 3);
352 mx.s1() = 2097152U >> ((clz_ - 0) / 3);
353 mx.s2() = 2097152U >> ((clz_ - 1) / 3);
354 return mx;
355 }
356
357 template<>
358 inline u16_3 get_offset<u32>(uint clz_) {
359 u16_3 mx;
360 mx.s0() = 1024U >> ((clz_ - 0) / 3);
361 mx.s1() = 1024U >> ((clz_ - 1) / 3);
362 mx.s2() = 1024U >> ((clz_ - 2) / 3);
363 return mx;
364 }
365
366} // namespace morton_3d
f64
double f64
Alias for double.
Definition aliases_float.hpp:20
f32
float f32
Alias for float.
Definition aliases_float.hpp:19
u32
std::uint32_t u32
32 bit unsigned integer
Definition aliases_int.hpp:27
u64
std::uint64_t u64
64 bit unsigned integer
Definition aliases_int.hpp:26
u16
std::uint16_t u16
16 bit unsigned integer
Definition aliases_int.hpp:28
bmi.hpp
Bit manipulation instruction implementation for SYCL.
morton_3d::morton_types
Helper struct to get types corresponding to a morton code representation.
Definition morton.hpp:243
Generated on Wed May 13 2026 13:21:49 for Shamrock by doxygen 1.9.8