Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Task10 Роман Гостило HSE #329

Closed
wants to merge 4 commits into from
Closed

Task10 Роман Гостило HSE #329

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
aeb2de9
add lbvh (apple silicon does not compile it locally so i will use ci
unicoooorn Jan 8, 2025
184b2d8
uncomment testing impls 0-2
unicoooorn Jan 8, 2025
69938ac
cleanup a bit
unicoooorn Jan 8, 2025
7c950c6
improve gpu naive impl
unicoooorn Jan 8, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
217 changes: 203 additions & 14 deletions src/cl/lbvh.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -126,20 +126,21 @@ morton_t zOrder(float fx, float fy, int i){
int x = fx + 0.5;
int y = fy + 0.5;

// у нас нет эксепшенов, но можно писать коды ошибок просто в консоль, и следить чтобы вывод был пустой

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;
}

// TODO
morton_t spread_bits_x = spreadBits(x);
morton_t spread_bits_y = spreadBits(y);

morton_t morton_code = 2 * spread_bits_x + spread_bits_y;

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

__kernel void generateMortonCodes(__global const float *pxs, __global const float *pys,
Expand Down Expand Up @@ -191,25 +192,138 @@ 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 left = i_begin;
int right = i_end;

while (left < right - 1) {
int mid = (left + right) / 2;

if (getBit(codes[mid], bit_index)) {
right = mid;
} else {
left = mid;
}
}

return right;
}

void findRegion(int *i_begin, int *i_end, int *bit_index, __global const morton_t *codes, int N, int i_node)
{
// TODO
morton_t lhs = codes[i_node - 1];
morton_t rhs = codes[i_node + 1];

morton_t diff = lhs ^ rhs;

*bit_index = log2((double)diff);

int dir = getBit(codes[i_node], *bit_index) ? 1 : -1;

if (dir == 1) {
*i_begin = i_node;
} else {
*i_end = i_node + 1;
}

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

int left = *i_begin - 1;
int right = *i_end;
int split;
while (left < right - 1) {
int mid = (left + right) / 2;

if ((getBits(codes[mid], *bit_index, K) == pref0) == (dir == 1)) {
left = mid;
} else {
right = mid;
}
}
split = right;

if (dir == 1) {
*i_begin = i_node;
*i_end = split;
} else if (dir == -1){
*i_begin = split;
*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];
float 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("043204230042342");
return;
}

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("54356549645");
return;
}
}

__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 gid = get_global_id(0);

if (gid >= tree_size) {
return;
}

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

void initFlag(__global int *flags, int i_node, __global const struct Node *nodes, int level)
Expand All @@ -219,7 +333,7 @@ void initFlag(__global int *flags, int i_node, __global const struct Node *nodes
__global const struct Node *node = &nodes[i_node];
if (isLeaf(node)) {
printf("9423584385834\n");
// return;
return;
}

if (!empty(&node->bbox)) {
Expand Down Expand Up @@ -263,7 +377,7 @@ void growNode(__global struct Node *root, __global struct Node *nodes)

if (root->mass <= 1e-8) {
printf("04230420340322\n");
// return;
return;
}

root->cmsx = (left->cmsx * m0 + right->cmsx * m1) / root->mass;
Expand Down Expand Up @@ -299,7 +413,73 @@ 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;

float t_force_x = 0.0;
float t_force_y = 0.0;

while (stack_size) {
__global const struct Node *node = &nodes[stack[--stack_size]];

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");
return;
}
if (!equalsPoint(&left->bbox, x0, y0)) {
printf("5446456456435656");
return;
}
continue;
}
}

int children[2] = {node->child_left, node->child_right};
for (int i = 0; i < 2; ++i) {
int i_child = children[i];
__global const 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;

float fx = ex * dr2_inv * GRAVITATIONAL_FORCE;
float fy = ey * dr2_inv * GRAVITATIONAL_FORCE;

t_force_x += m1 * fx;
t_force_y += m1 * fy;
} else {
stack[stack_size++] = i_child;

if (stack_size >= 2 * NBITS_PER_DIM) {
printf("0420392384283");
return;
}
}
}
}
*force_x = t_force_x;
*force_y = t_force_y;
}

__kernel void calculateForces(
Expand All @@ -311,7 +491,16 @@ __kernel void calculateForces(
int N,
int t)
{
// TODO
int gid = get_global_id(0);

float x0 = pxs[gid];
float y0 = pys[gid];
float m0 = mxs[gid];

__global float * dvx = dvx2d + t * N;
__global float * dvy = dvy2d + t * N;

calculateForce(x0, y0, m0, nodes, &dvx[gid], &dvy[gid]);
}

__kernel void integrate(
Expand Down
26 changes: 25 additions & 1 deletion src/cl/nbody.cl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 (i == j) {
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