Skip to content

Commit

Permalink
done
Browse files Browse the repository at this point in the history
  • Loading branch information
AvvAL committed Dec 10, 2023
1 parent 930b2f9 commit 89df99b
Show file tree
Hide file tree
Showing 4 changed files with 503 additions and 287 deletions.
1 change: 1 addition & 0 deletions src/cl/clion_defines.cl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ gentype clamp (gentype x, float minval, float maxval);
gentype degrees (gentype radians);
gentype max (gentype x, gentype y);
gentype min (gentype x, gentype y);
gentype sqrt (gentype x);
gentype mix (gentype x, gentype y, gentype a);
gentype radians (gentype degrees);
gentype sign (gentype x);
Expand Down
193 changes: 182 additions & 11 deletions src/cl/lbvh.cl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -130,16 +130,15 @@ morton_t zOrder(float fx, float fy, int i){

if (x < 0 || x >= (1 << NBITS_PER_DIM)) {
printf("098245490432590890\n");
// return 0;
//return 0;
}
if (y < 0 || y >= (1 << NBITS_PER_DIM)) {
printf("432764328764237823\n");
// return 0;
//return 0;
}
morton_t morton_code = (spreadBits(y) << 1) | spreadBits(x);

// TODO

return 0;
return (morton_code << 32) | i;
}

__kernel void generateMortonCodes(__global const float *pxs, __global const float *pys,
Expand Down Expand Up @@ -191,25 +190,132 @@ 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 start = i_begin;
while (i_begin < i_end) {
int mid = (i_begin + i_end) / 2;
if (getBit(codes[start], bit_index) != getBit(codes[mid], bit_index))
i_end = mid;
else
i_begin = mid + 1;
}
return i_begin;
}

void findRegion(int *i_begin, int *i_end, int *bit_index, __global const morton_t *codes, int N, int i_node)
{
// TODO
if (i_node < 1 || i_node > N - 2) {
printf("842384298293482\n");
return;
}
int dir = 0;
int i_bit = NBITS - 1;
for (; i_bit >= 0; --i_bit) {
if (getBit(codes[i_node - 1], i_bit) == 0 && getBit(codes[i_node + 1], i_bit) == 1) {
dir = getBit(codes[i_node], i_bit) == 0 ? -1 : 1;
break;
}
}

if (dir == 0) {
printf("8923482374983\n");
return;
}

int K = NBITS - i_bit;
morton_t pref0 = getBits(codes[i_node], i_bit, K);


int start = i_node, end = dir == 1 ? N : -1;
while (start + dir != end) {
int mid = start + (end - start) / 2;
if (getBits(codes[mid], i_bit, K) == pref0) {
start = mid;
} else {
end = mid;
}
}
int i_node_end = start;
*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
__global struct Node* node = &nodes[i_node];
clear(&node->bbox);
node->mass = 0;
node->cmsx = 0;
node->cmsy = 0;

if (i_node >= N - 1) {
node->child_left = -1;
node->child_right = -1;
int i_point = i_node - (N - 1);

int idx = getIndex(codes[i_point]);
float center_mass_x = pxs[idx], center_mass_y = pys[idx];
float mass = mxs[idx];

growPoint(&node->bbox, center_mass_x, center_mass_y);
node->cmsx = center_mass_x;
node->cmsy = center_mass_y;
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("043204230042342\n");
return;
}

// проинициализировать nodes[i_node].child_left, nodes[i_node].child_right на основе i_begin, i_end, split
node->child_left = split - i_begin < 2 ? split - 1 + (N - 1) : split - 1;
node->child_right = i_end - split < 2 ? split + (N - 1) : split;
// не забудьте на N-1 сдвинуть индексы, указывающие на листья

found = true;
break;
}

if (!found) {
printf("54356549645\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 tree_size = LBVHSize(N);
int i_node = get_global_id(0);
if (i_node >= tree_size)
return;

initLBVHNode(nodes, i_node, codes, N, pxs, pys, mxs);
}

void initFlag(__global int *flags, int i_node, __global const struct Node *nodes, int level)
Expand Down Expand Up @@ -299,7 +405,66 @@ 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];
int stack_size = 0;
stack[stack_size++] = 0;
while (stack_size) {
int i_node = stack[--stack_size];
__global const struct Node *node = &nodes[i_node];

if (isLeaf(node)) {
continue;
}

// если запрос содержится и а левом и в правом ребенке - то они в одном пикселе
{
__global const struct Node *left = &nodes[node->child_left];
__global const struct Node *right = &nodes[node->child_right];
if (contains(&left->bbox, x0, y0) && contains(&right->bbox, x0, y0)) {
if (!equals(&left->bbox, &right->bbox)) {
printf("42357987645432456547\n");
}
if (!equalsPoint(&left->bbox, x0, y0)) {
printf("5446456456435656\n");
}
continue;
}
}
int idx_children[] = {node->child_left, node->child_right};
for (int i = 0; i < 2; ++i) {
int i_child = idx_children[i];
__global const struct Node *child = &nodes[i_child];
// С точки зрения ббоксов заходить в ребенка, ббокс которого не пересекаем, не нужно (из-за того, что в листьях у нас точки и они не высовываются за свой регион пространства)
// Но, с точки зрения физики, замена гравитационного влияния всех точек в регионе на взаимодействие с суммарной массой в центре масс - это точное решение только в однородном поле (например, на поверхности земли)
// У нас поле неоднородное, и такая замена - лишь приближение. Чтобы оно было достаточно точным, будем спускаться внутрь ноды, пока она не станет похожа на точечное тело (маленький размер ее ббокса относительно нашего расстояния до центра масс ноды)
if (!contains(&child->bbox, x0, y0) && barnesHutCondition(x0, y0, child)) {
float dx = child->cmsx - x0;
float dy = child->cmsy - 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;

*force_x += child->mass * fx;
*force_y += child->mass * fy;
} else {
stack[stack_size++] = i_child;
if (stack_size >= 2 * NBITS_PER_DIM) {
printf("0420392384283\n");
return;
}
}
}
}
}

__kernel void calculateForces(
Expand All @@ -311,7 +476,13 @@ __kernel void calculateForces(
int N,
int t)
{
// TODO
int i = get_global_id(0);
float x0 = pxs[i];
float y0 = pys[i];
float m0 = mxs[i];
__global float *force_x = dvx2d + t * N + i;
__global float *force_y = dvy2d + t * N + i;
calculateForce(x0, y0, m0, nodes, force_x, force_y);
}

__kernel void integrate(
Expand Down
27 changes: 26 additions & 1 deletion src/cl/nbody.cl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,32 @@ __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(
Expand Down
Loading

0 comments on commit 89df99b

Please sign in to comment.