diff --git a/src/cl/lbvh.cl b/src/cl/lbvh.cl index e0423d86..1c6cafe9 100644 --- a/src/cl/lbvh.cl +++ b/src/cl/lbvh.cl @@ -137,9 +137,8 @@ morton_t zOrder(float fx, float fy, int i){ // return 0; } - // TODO - - return 0; + morton_t morton_code = spreadBits(x) | (spreadBits(y) << 1); + return (morton_code << 32) | i; } __kernel void generateMortonCodes(__global const float *pxs, __global const float *pys, @@ -191,25 +190,143 @@ void __kernel merge(__global const morton_t *as, __global morton_t *as_sorted, u int findSplit(__global const morton_t *codes, int i_begin, int i_end, int bit_index) { - // TODO + if (getBit(codes[i_begin], bit_index) == getBit(codes[i_end-1], bit_index)) { + return -1; + } + + int l = i_begin, r = i_end; + while (r - l > 1) { + int m = (l + r) / 2; + int a = getBit(codes[m], bit_index); + if (a == 0) { + l = m; + } else { + r = m; + } + } + return r; } void findRegion(int *i_begin, int *i_end, int *bit_index, __global const morton_t *codes, int N, int i_node) { - // TODO + int dir = 0; + int i_bit = NBITS-1; + for (; i_bit >= 0; --i_bit) { + int a = getBit(codes[i_node - 1], i_bit); + int b = getBit(codes[i_node], i_bit); + int c = getBit(codes[i_node + 1], i_bit); + if (a == 0 && b == 0 && c == 1) { + dir = -1; + break; + } + if (a == 0 && b == 1 && c == 1) { + dir = 1; + break; + } + } + + int K = NBITS - i_bit; + morton_t pref0 = getBits(codes[i_node], i_bit, K); + + // граница зоны ответственности - момент, когда префикс перестает совпадать + int i_node_end = -1; + + if (dir == 1) { + int l = i_node, r = N; + while (r - l > 1) { + int m = (l + r) / 2; + if (getBits(codes[m], i_bit, K) == pref0) { + l = m; + } else { + r = m; + } + } + i_node_end = l; + } else { + int l = -1, r = i_node; + while (r - l > 1) { + int m = (l + r) / 2; + if (getBits(codes[m], i_bit, K) == pref0) { + r = m; + } else { + l = m; + } + } + i_node_end = r; + } + + *bit_index = i_bit - 1; + + if (dir > 0) { + *i_begin = i_node; + *i_end = i_node_end + 1; + } else { + *i_begin = i_node_end; + *i_end = i_node + 1; + } } void initLBVHNode(__global struct Node *nodes, int i_node, __global const morton_t *codes, int N, __global const float *pxs, __global const float *pys, __global const float *mxs) { - // TODO + clear(&nodes[i_node].bbox); + nodes[i_node].mass = 0; + nodes[i_node].cmsx = 0; + nodes[i_node].cmsy = 0; + + if (i_node >= N-1) { + nodes[i_node].child_left = -1; + nodes[i_node].child_right = -1; + int i_point = i_node - (N-1); + + int i = getIndex(codes[i_point]); + float center_mass_x = pxs[i], center_mass_y = pys[i]; + float mass = mxs[i]; + + growPoint(&nodes[i_node].bbox, center_mass_x, center_mass_y); + nodes[i_node].cmsx = center_mass_x; + nodes[i_node].cmsy = center_mass_y; + nodes[i_node].mass = mass; + + return; + } + + // инициализируем внутреннюю ноду + + int i_begin = 0, i_end = N, bit_index = NBITS-1; + // если рассматриваем не корень, то нужно найти зону ответственности ноды и самый старший бит, с которого надо начинать поиск разреза + if (i_node) { + findRegion(&i_begin, &i_end, &bit_index, codes, N, i_node); + } + + bool found = false; + for (int i_bit = bit_index; i_bit >= 0; --i_bit) { + int split = findSplit(codes, i_begin, i_end, i_bit); + if (split < 0) continue; + + if (split < 1) { + printf("112974317235116477\n"); + } + + nodes[i_node].child_left = split - 1 + (split - 1 == i_begin ? N - 1 : 0); + nodes[i_node].child_right = split + (split == i_end - 1 ? N - 1 : 0); + + found = true; + break; + } + + if (!found) { + printf("380527406311728313\n"); + } } __kernel void buidLBVH(__global const float *pxs, __global const float *pys, __global const float *mxs, __global const morton_t *codes, __global struct Node *nodes, int N) { - // TODO + int i = get_global_id(0); + if (i > 2 * N - 1) return; + initLBVHNode(nodes, i, codes, N, pxs, pys, mxs); } void initFlag(__global int *flags, int i_node, __global const struct Node *nodes, int level) @@ -299,7 +416,55 @@ bool barnesHutCondition(float x, float y, __global const struct Node *node) void calculateForce(float x0, float y0, float m0, __global const struct Node *nodes, __global float *force_x, __global float *force_y) { - // TODO + int stack[2 * NBITS_PER_DIM]; + stack[0] = 0; + int stack_size = 1; + float fx = 0.0, fy = 0.0; + while (stack_size) { + const __global struct Node *node = &nodes[stack[--stack_size]]; + + if (isLeaf(node)) { + continue; + } + + { + const __global struct Node *left = &nodes[node->child_left]; + const __global struct Node *right = &nodes[node->child_right]; + if (contains(&left->bbox, x0, y0) && contains(&right->bbox, x0, y0)) { + continue; + } + } + + for (int child_ind = 0; child_ind < 2; child_ind++) { + int i_child = child_ind == 0 ? node->child_left : node->child_right; + const __global struct Node *child = &nodes[i_child]; + if (!contains(&child->bbox, x0, y0) && barnesHutCondition(x0, y0, child)) { + float x1 = child->cmsx; + float y1 = child->cmsy; + float m1 = child->mass; + + float dx = x1 - x0; + float dy = y1 - y0; + float dr2 = max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + fx += m1 * ex * dr2_inv * GRAVITATIONAL_FORCE; + fy += m1 * ey * dr2_inv * GRAVITATIONAL_FORCE; + } else { + stack[stack_size++] = i_child; + if (stack_size >= 2 * NBITS_PER_DIM) { + printf("496146404636797101\n"); + } + } + } + } + *force_x = fx; + *force_y = fy; } __kernel void calculateForces( @@ -311,7 +476,11 @@ __kernel void calculateForces( int N, int t) { - // TODO + int i = get_global_id(0); + if (i > N) + return; + + calculateForce(pxs[i], pys[i], mxs[i], nodes, &dvx2d[t * N + i], &dvy2d[t * N + i]); } __kernel void integrate( @@ -325,7 +494,7 @@ __kernel void integrate( { unsigned int i = get_global_id(0); - if (i >= N) + if (i > N) return; __global float * dvx = dvx2d + t * N; diff --git a/src/cl/nbody.cl b/src/cl/nbody.cl index c70482cc..e6a06f0c 100644 --- a/src/cl/nbody.cl +++ b/src/cl/nbody.cl @@ -26,7 +26,31 @@ __kernel void nbody_calculate_force_global( float y0 = pys[i]; float m0 = mxs[i]; - // TODO + for (int j = 0; j < N; ++j) { + if (j == i) { + continue; + } + + float x1 = pxs[j]; + float y1 = pys[j]; + float m1 = mxs[j]; + + float dx = x1 - x0; + float dy = y1 - y0; + float dr2 = max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + float fx = ex * dr2_inv * GRAVITATIONAL_FORCE; + float fy = ey * dr2_inv * GRAVITATIONAL_FORCE; + + dvx[i] += m1 * fx; + dvy[i] += m1 * fy; + } } __kernel void nbody_integrate( diff --git a/src/main_lbvh.cpp b/src/main_lbvh.cpp index ff486364..f7eb6255 100644 --- a/src/main_lbvh.cpp +++ b/src/main_lbvh.cpp @@ -20,7 +20,7 @@ #define OPENCL_DEVICE_INDEX 0 // TODO включить чтобы начали запускаться тесты -#define ENABLE_TESTING 0 +#define ENABLE_TESTING 1 // имеет смысл отключать при оффлайн симуляции больших N, но в итоговом решении стоит оставить #define EVALUATE_PRECISION 1 @@ -34,8 +34,8 @@ // TODO на сервер лучше коммитить самую простую конфигурацию. Замеры по времени получатся нерелевантные, но зато быстрее отработает CI // TODO локально интересны замеры на самой сложной версии, которую получится дождаться #define NBODY_INITIAL_STATE_COMPLEXITY 0 -//#define NBODY_INITIAL_STATE_COMPLEXITY 1 -//#define NBODY_INITIAL_STATE_COMPLEXITY 2 +// #define NBODY_INITIAL_STATE_COMPLEXITY 1 +// #define NBODY_INITIAL_STATE_COMPLEXITY 2 // использовать lbvh для построения начального состояния. Нужно на очень больших N (>1000000) #define ENABLE_LBVH_STATE_INITIALIZATION 0 @@ -195,11 +195,9 @@ morton_t zOrder(const Point &coord, int i){ int x = coord.x; int y = coord.y; - throw std::runtime_error("not implemented"); -// morton_t morton_code = TODO -// -// // augmentation -// return (morton_code << 32) | i; + morton_t morton_code = spreadBits(x) | (spreadBits(y) << 1); + // augmentation + return (morton_code << 32) | i; } #pragma pack (push, 1) @@ -368,12 +366,10 @@ void calculateForce(float x0, float y0, float m0, const std::vector &nodes // и не спускаться внутрь, если точка запроса не пересекает ноду, а заменить на взаимодействие с ее центром масс int stack[2 * NBITS_PER_DIM]; - int stack_size = 0; - // TODO кладем корень на стек - throw std::runtime_error("not implemented"); - /* while (stack_size) { - // TODO берем ноду со стека - throw std::runtime_error("not implemented"); + stack[0] = 0; + int stack_size = 1; + while (stack_size) { + auto node = nodes[stack[--stack_size]]; if (node.isLeaf()) { continue; @@ -400,17 +396,30 @@ void calculateForce(float x0, float y0, float m0, const std::vector &nodes // Но, с точки зрения физики, замена гравитационного влияния всех точек в регионе на взаимодействие с суммарной массой в центре масс - это точное решение только в однородном поле (например, на поверхности земли) // У нас поле неоднородное, и такая замена - лишь приближение. Чтобы оно было достаточно точным, будем спускаться внутрь ноды, пока она не станет похожа на точечное тело (маленький размер ее ббокса относительно нашего расстояния до центра масс ноды) if (!child.bbox.contains(x0, y0) && barnesHutCondition(x0, y0, child)) { - // TODO посчитать взаимодействие точки с центром масс ноды - throw std::runtime_error("not implemented"); + float x1 = child.cmsx; + float y1 = child.cmsy; + float m1 = child.mass; + + float dx = x1 - x0; + float dy = y1 - y0; + float dr2 = std::max(100.f, dx * dx + dy * dy); + + float dr2_inv = 1.f / dr2; + float dr_inv = std::sqrt(dr2_inv); + + float ex = dx * dr_inv; + float ey = dy * dr_inv; + + *force_x += m1 * ex * dr2_inv * GRAVITATIONAL_FORCE; + *force_y += m1 * ey * dr2_inv * GRAVITATIONAL_FORCE; } else { - // TODO кладем ребенка на стек - throw std::runtime_error("not implemented"); + stack[stack_size++] = i_child; if (stack_size >= 2 * NBITS_PER_DIM) { throw std::runtime_error("0420392384283"); } } } - }*/ + } } void integrate(int i, std::vector &pxs, std::vector &pys, std::vector &vxs, std::vector &vys, float *dvx, float *dvy, int coord_shift) @@ -975,8 +984,17 @@ int findSplit(const std::vector &codes, int i_begin, int i_end, int bi // } // } - // TODO бинпоиск для нахождения разбиения области ответственности ноды - throw std::runtime_error("not implemented"); + int l = i_begin, r = i_end; + while (r - l > 1) { + int m = (l + r) / 2; + int a = getBit(codes[m], bit_index); + if (a == 0) { + l = m; + } else { + r = m; + } + } + return r; // избыточно, так как на входе в функцию проверили, что ответ существует, но приятно иметь sanity-check на случай если набагали throw std::runtime_error("4932492039458209485"); @@ -1029,8 +1047,17 @@ void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector= 0; --i_bit) { - // TODO найти dir и значащий бит - throw std::runtime_error("not implemented"); + int a = getBit(codes[i_node - 1], i_bit); + int b = getBit(codes[i_node], i_bit); + int c = getBit(codes[i_node + 1], i_bit); + if (a == 0 && b == 0 && c == 1) { + dir = -1; + break; + } + if (a == 0 && b == 1 && c == 1) { + dir = 1; + break; + } } if (dir == 0) { @@ -1057,8 +1084,29 @@ void findRegion(int *i_begin, int *i_end, int *bit_index, const std::vector 1) { + int m = (l + r) / 2; + if (getBits(codes[m], i_bit, K) == pref0) { + l = m; + } else { + r = m; + } + } + i_node_end = l; + } else { + int l = -1, r = i_node; + while (r - l > 1) { + int m = (l + r) / 2; + if (getBits(codes[m], i_bit, K) == pref0) { + r = m; + } else { + l = m; + } + } + i_node_end = r; + } *bit_index = i_bit - 1; @@ -1109,28 +1157,20 @@ void initLBVHNode(std::vector &nodes, int i_node, const std::vector= 0; --i_bit) { - /* - int split = TODO + int split = findSplit(codes, i_begin, i_end, i_bit); if (split < 0) continue; if (split < 1) { throw std::runtime_error("043204230042342"); } - */ - throw std::runtime_error("not implemented"); - - - // TODO проинициализировать nodes[i_node].child_left, nodes[i_node].child_right на основе i_begin, i_end, split - // не забудьте на N-1 сдвинуть индексы, указывающие на листья - - throw std::runtime_error("not implemented"); + nodes[i_node].child_left = split - 1 + (split - 1 == i_begin ? N - 1 : 0); + nodes[i_node].child_right = split + (split == i_end - 1 ? N - 1 : 0); found = true; break; @@ -1240,11 +1280,10 @@ void buildBBoxes(std::vector &nodes, std::vector &flags, int N, bool int n_updated = 0; #pragma omp parallel for if(use_omp) reduction(+:n_updated) for (int i_node = 0; i_node < N-1; ++i_node) { - // TODO если находимся на нужном уровне (нужный flag), проинициализируем ббокс и центр масс ноды -// if (TODO) { -// TODO -// ++n_updated; -// } + if (flags[i_node] == level) { + growNode(nodes[i_node], nodes); + ++n_updated; + } }