萘Naie_Char

1
__global__ void computeForcesNaive(const Particle* particles, float* accX, float* accY,
2
                                    int n, float G, float softening)
3
{
4
    int i = blockIdx.x * blockDim.x + threadIdx.x;
5
    if (i >= n) return;
6

7
    float myPosX = particles[i].posX;
8
    float myPosY = particles[i].posY;
9
    float aX = 0.0f, aY = 0.0f;
10

11
    for (int j = 0; j < n; ++j) {
12
        if (j == i) continue;
13

14
        float dx = particles[j].posX - myPosX;
15
        float dy = particles[j].posY - myPosY;
16
        float distSqr = dx * dx + dy * dy + softening * softening;
17

18
        float invDist = rsqrtf(distSqr);   // 快速倒数平方根内建函数
19
        float invDist3 = invDist * invDist * invDist;
20

21
        aX += G * particles[j].mass * dx * invDist3;
22
        aY += G * particles[j].mass * dy * invDist3;
23
    }
24

25
    accX[i] = aX;
26
    accY[i] = aY;
27
}

1
#define TILE_SIZE 256
2

3
struct TileData { float posX, posY, mass; };
4

5
__global__ void computeForcesTiled(const Particle* particles, float* accX, float* accY,
6
                                    int n, float G, float softening)
7
{
8
    __shared__ TileData tile[TILE_SIZE];
9

10
    int i = blockIdx.x * blockDim.x + threadIdx.x;
11

12
    // 注意：这里不能提前 return！即使 i >= n，这个线程仍然要参与
13
    // 后面协作搬运数据到 shared memory 的工作，否则会导致本block内的
14
    // __syncthreads() 死锁——因为 __syncthreads() 要求同一 block 内
15
    // 所有线程都执行到这一行才能继续，提前退出的线程永远不会走到那一行，
16
    // 其余线程会永远卡在那里等
17
    float myPosX = (i < n) ? particles[i].posX : 0.0f;
18
    float myPosY = (i < n) ? particles[i].posY : 0.0f;
19

20
    float aX = 0.0f, aY = 0.0f;
21
    int numTiles = (n + TILE_SIZE - 1) / TILE_SIZE;
22

23
    for (int t = 0; t < numTiles; ++t) {
24
        int loadIdx = t * TILE_SIZE + threadIdx.x;
25

26
        // 协作加载：每个线程只搬1个粒子进shared memory，分摊带宽
27
        if (loadIdx < n) {
28
            tile[threadIdx.x] = { particles[loadIdx].posX,
29
                                   particles[loadIdx].posY,
30
                                   particles[loadIdx].mass };
31
        } else {
32
            tile[threadIdx.x] = { 0.0f, 0.0f, 0.0f }; // 质量为0，不产生引力贡献
33
        }
34

35
        __syncthreads(); // 同步点1：确保整个block都搬完数据，才能开始读
36

37
        for (int k = 0; k < TILE_SIZE; ++k) {
38
            int j = t * TILE_SIZE + k;
39
            if (j == i || j >= n) continue;
40

41
            float dx = tile[k].posX - myPosX;
42
            float dy = tile[k].posY - myPosY;
43
            float distSqr = dx * dx + dy * dy + softening * softening;
44
            float invDist = rsqrtf(distSqr);
45
            float invDist3 = invDist * invDist * invDist;
46

47
            aX += G * tile[k].mass * dx * invDist3;
48
            aY += G * tile[k].mass * dy * invDist3;
49
        }
50

51
        __syncthreads(); // 同步点2：确保整个block都用完数据，才能覆盖写下一tile
52
    }
53

54
    if (i < n) { accX[i] = aX; accY[i] = aY; }
55
}