使用 CUDA 的显式 FDM

Question

我正在为以下代码实现 CUDA。第一个版本是连续写的，第二个版本是用CUDA写的。我确信它在串行版本中的结果。我希望我添加了 CUDA 功能的第二个版本也能给我相同的结果，但内核函数似乎没有做任何事情，它给了我 u 和 v 的初始值。我知道由于缺乏经验，错误可能很明显，但我无法弄清楚。另外，请不要推荐使用 flatten array，因为我很难理解代码中的索引。第一版：

#include <fstream>
#include <iostream>
#include <math.h>
#include <vector>
#include <chrono>
#include <omp.h>
using namespace std;
const int M = 1024; 
const int N = 1024; 
const double A = 1;
const double B = 3;
const double Du = 5 * pow(10, -5);
const double Dv = 5 * pow(10, -6);
const int Max_Itr = 1000;
const double h = 1.0 / static_cast<double>(M - 1);
const double delta_t = 0.0025;
const double s1 = (Du * delta_t) / pow(h, 2);
const double s2 = (Dv * delta_t) / pow(h, 2);
int main() {
    double** u=new double* [M]; 
    double** v=new double* [M];
        for (int i=0; i<M; i++){
            u[i]=new double [N]; 
            v[i]=new double [N]; 
        }
    for (int j = 0; j < M; j++) {
        for (int i = 0; i < N;i++) {
            u[i][j]=0.02;
            v[i][j]=0.02;

        }
    }

    for (int k = 1; k < Max_Itr; k++) {
            
        for (int i = 1; i < N - 1; i++) {
            for (int j = 1; j < M - 1; j++) {

                u[i][j] = ((1 - (4 * s1)) * u[i][j]) + (s1 * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1])) +

                    (A * delta_t) + (delta_t * pow(u[i][j], 2) * v[i][j]) - (delta_t * (B + 1) * u[i][j]);


                v[i][j] = ((1 - (4 * s2)) * v[i][j]) + (s2 * (v[i + 1][j] + v[i - 1][j] + v[i][j + 1] + v[i][j - 1])) + (B * delta_t * u[i][j])
                    - (delta_t * pow(u[i][j], 2) * v[i][j]);

                }

            }


        }
        cout<<"u: "<<u[512][512]<<" v: "<<v[512][512]<<endl; 

    return 0;
}

第二个版本：

#include <fstream>
#include <iostream>
#include <math.h>
#include <vector>
using namespace std;

#define M 1024

#define N 1024


__global__ void my_kernel(double** v, double** u){
    int i= blockIdx.y * blockDim.y + threadIdx.y;
    int j= blockIdx.x * blockDim.x + threadIdx.x;
    double A = 1;
    double B = 3;
    int Max_Itr = 1000;
    double delta_t = 0.0025;
    double Du = 5 * powf(10, -5);
    double Dv = 5 * powf(10, -6);
    double h = 1.0 / (M - 1);
    double s1 = (Du * delta_t) / powf(h, 2);

    double s2 = (Dv * delta_t) / powf(h, 2);
    for (int k = 1; k < Max_Itr; k++) {
        u[i][j] = ((1 - (4 * s1)) 
        * u[i][j]) + (s1 * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1])) +
        (A * delta_t) + (delta_t * pow(u[i][j], 2) * v[i][j]) - (delta_t * (B + 1) * u[i][j]);

        v[i][j] = ((1 - (4 * s2)) 
        * v[i][j]) + (s2 * (v[i + 1][j] + v[i - 1][j] + v[i][j + 1] + v[i][j - 1])) + (B * delta_t * u[i][j])
        - (delta_t * pow(u[i][j], 2) * v[i][j]);

    
       
        __syncthreads();
    }
}

int main() {
    double** u=new double* [M]; 
    double** v=new double* [M];
    for (int i=0; i<M; i++){
    
        u[i]=new double [N]; 
        v[i]=new double [N]; 
    }

    dim3 blocks(32,32);
        dim3 grids(M/32 +1, N/32 + 1);
    for (int j = 0; j < M; j++) {
        for (int i = 0; i < N;i++) {
        u[i][j]=0.02;
        v[i][j]=0.02;
           
        }
    }
     double **u_d, **v_d;
     int d_size = N * M * sizeof(double);

     cudaMalloc(&u_d, d_size);
     cudaMalloc(&v_d, d_size);
     cudaMemcpy(u_d, u, d_size, cudaMemcpyHostToDevice);
         cudaMemcpy(v_d, v, d_size, cudaMemcpyHostToDevice);
     my_kernel<<<grids, blocks>>> (v_d,u_d);
         cudaDeviceSynchronize();
     cudaMemcpy(v, v_d, d_size, cudaMemcpyDeviceToHost);
     cudaMemcpy(u, u_d, d_size, cudaMemcpyDeviceToHost);
     cout<<"u: "<<u[512][512]<<" v: "<<v[512][512]<<endl; 



    return 0;
}

我对第二个版本的期望是：

u: 0.2815 v: 1.7581

Answer 1

您的二维数组 - 在程序的第一个版本中 - 是使用指针数组实现的，每个指针指向一个单独分配的 double 值数组。

在你的第二个版本中，你使用相同的指针到指针到 double 类型，但是 - 你没有为实际数据分配任何 space，只是为了指针数组（而不是将任何数据复制到 GPU - 只是指针；无论如何复制都是无用的，因为它们是指向主机端内存的指针。）

最有可能发生的情况是您的内核试图访问无效地址的内存，并且其执行被中止。

如果你properly check for errors，正如@njuffa 指出的那样，你就会知道发生了什么。

现在，如果您要使用单个数据区域而不是为每个二维一维数组单独分配，则可以避免进行多次内存分配；对于程序的第一个和第二个版本都是如此。那不完全是数组展平。在 this page.

上查看有关如何执行此操作（C 语言风格）的说明

但是请注意，您坚持在内核中执行的双重解引用可能会显着降低速度。

使用 CUDA 的显式 FDM

Explicit FDM with CUDA

cuda

multidimensional-array