早在2006年以前，计算机性能的提升主要依靠提高CPU主频来实现。随着CPU主频不断提高，其发热量也持续增加。当主频高达3.8GHZ的Pentium4处理器问世之后，继续通过提升主频来获得性能改善已变得十分有限，而由此引发的CPU发热问题却越来越难以控制，系统稳定性也开始下降。“高频高能”的发展路径走到了尽头，各主要CPU制造商纷纷转向多核处理器的研发，即在一个物理CPU内部集成多个计算核心，通过并行处理方式来提升计算机的整体性能。自此，计算机发展迈入了一个新阶段。现今，随着深度学习对大数据计算的要求，GPU横空出世，强势进入高性能，并行计算的战场。

本文主要以并行计算为切入点，概括性的介绍高性能计算（从硬件到软件），并最后结合相关实例，比较不同编程架构的优缺点

一、多核，多进程和多线程与GPU

首先我们需要理解多核，多进程和多线程的区别（粒度上依次减小），此外还有GPU与并行的关系:

1首先理解一下什么是进程与线程（概念很重要）

1.进程是操作系统分配资源的最小单元, 线程是操作系统调度的最小单元。(一个是供分配，一个是供调度)
2.一个应用程序至少包括1个进程，而1个进程包括1个或多个线程，线程的尺度更小。
3.每个进程在执行过程中拥有独立的内存单元，而一个进程的多个线程在执行过程中共享内存。

举个例子：
①一座工厂（类似CPU），假设电力有限，只能供一个车间，即只能运行一个任务，里面有许多的车间（类似进程），执行单个的任务，这时，就是每次只能单个车间运行，若是另一个车间想工作，则当前车间得休息。
②而多核CPU，就像多个工厂，它可以同时让多个车间（类似多个进程）一起工作，当然，是工作在不同工厂里，也就是运行在不同CPU核上。
③而每个车间里有很多的工人，这就类似是线程，一个车间可以有多个工人，也就是一个进程可有多个线程。(引自(8 封私信 / 62 条消息) 一文弄懂多进程与多线程 - 知乎)

而GPU,简单理解，我们就可以看成有超级多核的CPU.

具体参考如下：

(8 封私信 / 62 条消息) 一文弄懂多进程与多线程 - 知乎

(8 封私信 / 64 条消息) 有人能说一下GPU并行计算和CPU多线程计算有什么区别？ - 知乎

第五篇：浅谈CPU 并行编程和 GPU 并行编程的区别 - 穆晨 - 博客园

2.具体来看gpu和cuda的硬件架构和运行逻辑

要深刻理解kernel，必须要对kernel的线程层次结构有一个清晰的认识。首先GPU上很多并行化的轻量级线程。kernel在device上执行时实际上是启动很多线程，一个kernel所启动的所有线程称为一个网格（grid），同一个网格上的线程共享相同的全局内存空间，grid是线程结构的第一层次，而网格又可以分为很多线程块（block），一个线程块里面包含很多线程，这是第二个层次。

此外这里简单介绍一下CUDA的内存模型，如下图所示。可以看到，每个线程有自己的私有本地内存（Local Memory），而每个线程块有包含共享内存（Shared Memory）,可以被线程块中所有线程共享，其生命周期与线程块一致。此外，所有的线程都可以访问全局内存（Global Memory）。还可以访问一些只读内存块：常量内存（Constant Memory）和纹理内存（Texture Memory）。内存结构涉及到程序优化，这里不深入探讨它们。

还有重要一点，你需要对GPU的硬件实现有一个基本的认识。上面说到了kernel的线程组织层次，那么一个kernel实际上会启动很多线程，这些线程是逻辑上并行的，但是在物理层却并不一定。这其实和CPU的多线程有类似之处，多线程如果没有多核支持，在物理层也是无法实现并行的。但是好在GPU存在很多CUDA核心，充分利用CUDA核心可以充分发挥GPU的并行计算能力。GPU硬件的一个核心组件是SM，前面已经说过，SM是英文名是 Streaming Multiprocessor，翻译过来就是流式多处理器。SM的核心组件包括CUDA核心，共享内存，寄存器等，SM可以并发地执行数百个线程，并发能力就取决于SM所拥有的资源数。当一个kernel被执行时，它的gird中的线程块被分配到SM上，一个线程块只能在一个SM上被调度。SM一般可以调度多个线程块，这要看SM本身的能力。那么有可能一个kernel的各个线程块被分配多个SM，所以grid只是逻辑层，而SM才是执行的物理层。SM采用的是SIMT (Single-Instruction, Multiple-Thread，单指令多线程)架构，基本的执行单元是线程束（warps)，线程束包含32个线程，这些线程同时执行相同的指令，但是每个线程都包含自己的指令地址计数器和寄存器状态，也有自己独立的执行路径。所以尽管线程束中的线程同时从同一程序地址执行，但是可能具有不同的行为，比如遇到了分支结构，一些线程可能进入这个分支，但是另外一些有可能不执行，它们只能死等，因为GPU规定线程束中所有线程在同一周期执行相同的指令，线程束分化会导致性能下降。当线程块被划分到某个SM上时，它将进一步划分为多个线程束，因为这才是SM的基本执行单元，但是一个SM同时并发的线程束数是有限的。这是因为资源限制，SM要为每个线程块分配共享内存，而也要为每个线程束中的线程分配独立的寄存器。所以SM的配置会影响其所支持的线程块和线程束并发数量。总之，就是网格和线程块只是逻辑划分，一个kernel的所有线程其实在物理层是不一定同时并发的。所以kernel的grid和block的配置不同，性能会出现差异，这点是要特别注意的。还有，由于SM的基本执行单元是包含32个线程的线程束，所以block大小一般要设置为32的倍数。

CUDA编程的逻辑层和物理层

引用于(7 封私信) CUDA编程入门极简教程 - 知乎！！！！！

3.什么是MPI和OpenMP

OpenMP、MPI、CUDA总结_mpi与cuda-CSDN博客

总结如下：

特性	OpenMP	MPI	CUDA
全称	Open Multi-Processing	Message Passing Interface	Compute Unified Device Architecture
主要硬件目标	多核 CPU（共享内存）	多节点集群（分布式内存）	GPU（异构加速器）
内存模型	共享内存	分布式内存	分层内存（主机+设备）
编程方式	基于线程的并行化（pragma 指令）	基于进程的通信（显式消息传递）	基于核函数（kernel）的 GPU 并行计算
典型语言接口	C/C++、Fortran	C/C++、Fortran、Python（mpi4py）	C/C++、Python（PyCUDA、Numba）、Fortran
并行粒度	细粒度（线程级）	粗粒度（进程级）	极细粒度（线程块级、线程级）

这个类比对照这个blog----(8 封私信 / 12 条消息) 一文弄懂多进程与多线程 - 知乎

类比	含义
OpenMP 像“多线程团队协作”	一台机器上多个线程协作完成任务
MPI 像“多公司合作项目”	各个节点独立运行、通过信件通信
CUDA 像“工厂流水线”	GPU 大量工人同时执行同类操作

问题类型	适合 OpenMP	适合 MPI	适合 CUDA
小规模多核任务	✅ 非常适合	❌ 通信开销太大	⚠️ 数据量需足够大
大规模分布式仿真	⚠️ 内存受限	✅ 最佳选择	⚠️ 可作加速器辅助
向量/矩阵计算	✅ 容易实现	✅ 分布式矩阵分块	✅ 极高性能
深度学习	⚠️ 一般	✅ 分布式训练	✅ 核心加速技术
图像/视频处理	✅ 并行循环	⚠️ 较麻烦	✅ 极高吞吐量
物理仿真/CFD	✅ 小规模	✅ 大规模	✅ 混合并行最优
快速原型开发	✅ 简单上手	❌ 编程复杂	⚠️ 学习曲线陡峭

二、地震波场正演实践

1.公开代码

https://github.com/larsgeb/psvWave

A Python/C++ package for 2D P-SV wave propagation using finite differences and OpenMP. This package was written to facilitate high-throughput numerical wave simulations for Monte Carlo simulation in Seismology. It uses the velocity-stress formulation on a staggered grid from Virieux's classical 1986 paper. For compilation we require only OpenMP and the git subrepos (header-only): Eigen and inih, however installation can be easily done through pip. Used as a PDE-simulation code for this publication.

https://github.com/ar4/deepwave?utm_source=chatgpt.com

Deepwave provides wave propagation modules for PyTorch, for applications such as seismic imaging/inversion. You can use it to perform forward modelling and backpropagation, so it can simulate wave propagation to generate synthetic data, invert for the scattering potential (RTM/LSRTM), other model parameters (FWI), initial wavefields, or source wavelets. You can use it to integrate wave propagation into a larger chain of operations with end-to-end forward and backpropagation. Deepwave enables you to easily experiment with your own objective functions or functions that generate the inputs to the propagator, letting PyTorch's automatic differentiation do the hard work of calculating how to backpropagate through them.

这是

2.简单示例（openMP，cuda）

2.1 openMP

/*
the task ,written for HPC homework, all is about how to FDFD.

Be carefull,l don't want to make it for real job, just for eduation.so after coding, you need to 
know the main progress of FDFD, and get dispersion curve through different transformation
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stddef.h> 
#include <omp.h>
#include <time.h> 
//acutally it's not the esscential,because math.h will tell gcc .
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

/* Indexing helper: flattened 2D -> 1D
   Use Nx x Nz grid for stresses at integer nodes.
   Velocities are stored on same arrays but logically staggered via offsets when differencing.
   [i=0,k=0], [i=0,k=1], ..., [i=0,k=Nz-1],  // 第一行(x=0)
    [i=1,k=0], [i=1,k=1], ..., [i=1,k=Nz-1],  // 第二行(x=1)
    ...
    [i=Nx-1,k=0], ..., [i=Nx-1,k=Nz-1]        // 最后一行(x=Nx-1)
*/
#define IDX(i,k,Nz) ((size_t)(i)*(Nz) + (k))



/* Helpers to clamp */
static inline int min_i(int a, int b){ return a<b? a:b; }
static inline int max_i(int a, int b){ return a>b? a:b; }

/* Write a single-channel PGM (8-bit) from float array (size Nx*Nz) */
void write_pgm(const char *name, float *data, int Nx, int Nz){
    // find min/max
    float mn = data[0], mx = data[0];
    for(int i=0;i<Nx*Nz;i++){ if(data[i]<mn) mn=data[i]; if(data[i]>mx) mx=data[i]; }
    float scale = (mx>mn) ? 255.0f/(mx-mn) : 1.0f;
    FILE *f = fopen(name,"wb");
    if(!f){ perror("fopen"); return; }
    fprintf(f,"P5\n%d %d\n255\n", Nx, Nz); 
    unsigned char *buf = (unsigned char*)malloc(Nx*Nz);
    for(int k=0;k<Nz;k++){        // 外层循环：z方向（图像行）
        for(int i=0;i<Nx;i++){    // 内层循环：x方向（图像列）
            float v = data[IDX(i,k,Nz)];                    // 获取原始数据值
            unsigned char c = (unsigned char)((v-mn)*scale); // 映射到0-255
            buf[k*Nx + i] = c;                              // 存储到缓冲区
        }
    }
    fwrite(buf,1,Nx*Nz,f);
    free(buf);
    fclose(f);
}

/* Ricker wavelet */
float ricker(float t, float f0){
    float x = M_PI*f0*t;
    float x2 = x*x;
    return (1.0f - 2.0f*x2) * expf(-x2);
}

int main(){
    // 添加计时开始
    clock_t start_time = clock();
    double clock_time = omp_get_wtime();


    // OpenMP 设置
    int nthreads = omp_get_max_threads(); // 你也可以设置为 omp_get_max_threads()
    omp_set_num_threads(nthreads);
    printf("Running with %d OpenMP threads\n", nthreads);

    /* Simulation parameters */
    const int Nx = 301*5;      // grid points in x (i)
    const int Nz = 301*5;      // grid points in z (k)
    const float dx = 5.0f;   // m
    const float dz = 5.0f;   // m
    const int nt = 1000;     // time steps
    const float dt = 0.001f; // s (choose to satisfy CFL)

    /* Material: homogeneous model */
    const float rho0 = 2000.0f; // kg/m3
    const float vp0 = 3000.0f;  // P-wave speed m/s
    const float vs0 = 1500.0f;  // S-wave speed m/s

    /* Derived Lame parameters */
    const float mu0 = rho0 * vs0 * vs0;             // shear modulus
    const float lambda0 = rho0 * vp0*vp0 - 2.0f*mu0; // lambda

    /* Source */
    const int isrc = Nx/2;      // source location (stress grid)
    const int ksrc = Nz/2;
    const float f0 = 20.0f;     // dominant frequency
    const float t0 = 1.0f/f0;   // delay
    const float src_amp = 1e9f; // strength (stress units)

    /* Receivers: inline along surface */
    const int nrec = Nx;
    float *rec_time = (float*)malloc(sizeof(float)*nt*nrec);
    if(!rec_time){ perror("malloc rec_time"); return 1; }
    memset(rec_time, 0, sizeof(float)*nt*nrec);

    /* Fixed depth profile data at specific times */
    const int k_fixed_depth = Nz/2;  // 固定深度
    const int num_time_slices = 10;  // 保存10个时间切片
    const int time_slice_interval = nt / num_time_slices; // 时间切片间隔
    float **depth_time_slices = (float**)malloc(num_time_slices * sizeof(float*));
    if(!depth_time_slices){ perror("malloc depth_time_slices"); return 1; }
    for(int i=0; i<num_time_slices; i++){
        depth_time_slices[i] = (float*)malloc(Nx * sizeof(float));
        if(!depth_time_slices[i]){ perror("malloc depth_time_slices[i]"); return 1; }
    }
    int *slice_times = (int*)malloc(num_time_slices * sizeof(int)); // 保存对应的时间步
    if(!slice_times){ perror("malloc slice_times"); return 1; }

    /* Snapshot interval */
    const int SNAPSHOT_INTERVAL = 200; // steps

    /* Damping sponge parameters */
    const int npad = 40; // width of absorbing sponge
    const float damp_max = 0.015f; // damping coefficient magnitude

    /* Allocate fields */
    size_t nxy = (size_t)Nx * Nz;
    float *sxx = (float*)calloc(nxy, sizeof(float));
    float *szz = (float*)calloc(nxy, sizeof(float));
    float *sxz = (float*)calloc(nxy, sizeof(float));
    float *vx  = (float*)calloc(nxy, sizeof(float));
    float *vz  = (float*)calloc(nxy, sizeof(float));
    float *rho = (float*)malloc(sizeof(float)*nxy);
    float *lambda = (float*)malloc(sizeof(float)*nxy);
    float *mu = (float*)malloc(sizeof(float)*nxy);
    float *damp = (float*)calloc(nxy, sizeof(float));
    if(!sxx || !szz || !sxz || !vx || !vz || !rho || !lambda || !mu || !damp){ perror("malloc"); return 1; }

    /* Initialize model */
    for(int i=0;i<Nx;i++){
        for(int k=0;k<Nz;k++){
            size_t id = IDX(i,k,Nz);
            rho[id] = rho0;
            mu[id] = mu0;
            lambda[id] = lambda0;
        }
    }

    /* Precompute damping profile */
    for(int i=0;i<Nx;i++){
        for(int k=0;k<Nz;k++){
            float cx = 0.0f, cz = 0.0f;
            if(i < npad) { float xi = (npad - i) / (float)npad; cx = xi*xi; }
            if(i >= Nx-npad) { float xi = (i - (Nx-npad-1)) / (float)npad; if(xi<0) xi=0; cx = cx + xi*xi; }
            if(k < npad) { float zk = (npad - k) / (float)npad; cz = zk*zk; }
            if(k >= Nz-npad) { float zk = (k - (Nz-npad-1)) / (float)npad; if(zk<0) zk=0; cz = cz + zk*zk; }
            float val = damp_max * (cx + cz);
            damp[IDX(i,k,Nz)] = val;
        }
    }

    /* FD coefficients for 4th-order centered first derivative */
    const float c1 = 8.0f/12.0f; // 2/3
    const float c2 = 1.0f/12.0f; // 1/12

    /* Time stepping loop */
    for(int it=0; it<nt; it++){
        float time = it * dt;

        // Save depth profile at specific times
        if(it % time_slice_interval == 0){
            int slice_idx = it / time_slice_interval;
            if(slice_idx < num_time_slices){
                slice_times[slice_idx] = it;
                for(int i=0; i<Nx; i++){
                    depth_time_slices[slice_idx][i] = vz[IDX(i,k_fixed_depth,Nz)];
                }
                printf("Saved depth profile at time step %d (t=%.3f s)\n", it, time);
            }
        }
        #pragma omp parallel for collapse(2)
        // update velocities from stresses
        for(int i=2; i < Nx-2; i++){
            for(int k=2; k < Nz-2; k++){
                size_t id = IDX(i,k,Nz);
                float dsxx_dx = ( -sxx[IDX(i+2,k,Nz)] + 8.0f*sxx[IDX(i+1,k,Nz)] - 8.0f*sxx[IDX(i-1,k,Nz)] + sxx[IDX(i-2,k,Nz)] ) / (12.0f*dx);
                float dsxz_dz = ( -sxz[IDX(i,k+2,Nz)] + 8.0f*sxz[IDX(i,k+1,Nz)] - 8.0f*sxz[IDX(i,k-1,Nz)] + sxz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
                float invrho = 1.0f / rho[id];
                vx[id] += dt * invrho * (dsxx_dx + dsxz_dz);
            }
        }
        #pragma omp parallel for collapse(2)
        // vz update
        for(int i=2; i < Nx-2; i++){
            for(int k=2; k < Nz-2; k++){
                size_t id = IDX(i,k,Nz);
                float dsxz_dx = ( -sxz[IDX(i+2,k,Nz)] + 8.0f*sxz[IDX(i+1,k,Nz)] - 8.0f*sxz[IDX(i-1,k,Nz)] + sxz[IDX(i-2,k,Nz)] ) / (12.0f*dx);
                float dszz_dz = ( -szz[IDX(i,k+2,Nz)] + 8.0f*szz[IDX(i,k+1,Nz)] - 8.0f*szz[IDX(i,k-1,Nz)] + szz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
                float invrho = 1.0f / rho[id];
                vz[id] += dt * invrho * (dsxz_dx + dszz_dz);
            }
        }
        #pragma omp parallel for collapse(2)
        // update stresses from velocities
        for(int i=2;i<Nx-2;i++){
            for(int k=2;k<Nz-2;k++){
                size_t id = IDX(i,k,Nz);
                float dvx_dx = ( -vx[IDX(i+2,k,Nz)] + 8.0f*vx[IDX(i+1,k,Nz)] - 8.0f*vx[IDX(i-1,k,Nz)] + vx[IDX(i-2,k,Nz)] ) / (12.0f*dx);
                float dvz_dz = ( -vz[IDX(i,k+2,Nz)] + 8.0f*vz[IDX(i,k+1,Nz)] - 8.0f*vz[IDX(i,k-1,Nz)] + vz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
                float dvx_dz = ( -vx[IDX(i,k+2,Nz)] + 8.0f*vx[IDX(i,k+1,Nz)] - 8.0f*vx[IDX(i,k-1,Nz)] + vx[IDX(i,k-2,Nz)] ) / (12.0f*dz);
                float dvz_dx = ( -vz[IDX(i+2,k,Nz)] + 8.0f*vz[IDX(i+1,k,Nz)] - 8.0f*vz[IDX(i-1,k,Nz)] + vz[IDX(i-2,k,Nz)] ) / (12.0f*dx);

                float lam = lambda[id];
                float muval = mu[id];
                sxx[id] += dt * ( (lam + 2.0f*muval) * dvx_dx + lam * dvz_dz );
                szz[id] += dt * ( (lam + 2.0f*muval) * dvz_dz + lam * dvx_dx );
                sxz[id] += dt * ( muval * (dvx_dz + dvz_dx) );
            }
        }

        // add source
        {
            float src = src_amp * ricker(time - t0, f0);
            int i = isrc;
            int k = ksrc;
            size_t id = IDX(i,k,Nz);
            szz[id] += dt * src * 0.5f;
            sxx[id] += dt * src * 0.5f;
        }

        // record receivers at surface
        int krec = 2;
        for(int ir=0; ir<nrec; ir++){
            rec_time[it*nrec + ir] = vz[IDX(ir,krec,Nz)];
        }

        // snapshots
        if (it % SNAPSHOT_INTERVAL == 0) {
            float *snap_vz = (float*)malloc(sizeof(float) * nxy);
            if (snap_vz) {
                for (size_t ii = 0; ii < nxy; ii++) snap_vz[ii] = vz[ii];
                char fname_vz[256];
                snprintf(fname_vz, sizeof(fname_vz), "snapvz_t%04d.pgm", it);
                write_pgm(fname_vz, snap_vz, Nx, Nz);
                free(snap_vz);
                printf("Wrote %s (t=%.3f s)\n", fname_vz, it * dt);
            }
        }

        if(it % 100 == 0) printf("Step %d / %d (t=%.3f s)\n", it, nt, time);
    }

    // Save depth profile time slices
    FILE *fdepth = fopen("depth_profile_time_slices.txt","w");
    if(fdepth){
        fprintf(fdepth, "# Depth profile at z=%d (%.1f m) at different times\n", k_fixed_depth, k_fixed_depth*dz);
        fprintf(fdepth, "# Columns: x_position(m)");
        for(int i=0; i<num_time_slices; i++){
            fprintf(fdepth, " time_%.3fs", slice_times[i]*dt);
        }
        fprintf(fdepth, "\n");
        
        for(int i=0; i<Nx; i++){
            fprintf(fdepth, "%.1f", i*dx);
            for(int t=0; t<num_time_slices; t++){
                fprintf(fdepth, " %.6e", depth_time_slices[t][i]);
            }
            fprintf(fdepth, "\n");
        }
        fclose(fdepth);
        printf("Saved depth_profile_time_slices.txt\n");
    }

    // Save seismograms
    FILE *frec = fopen("seismograms.txt","w");
    if(frec){
        for(int it=0; it<nt; it++){
            for(int ir=0; ir<nrec; ir++){
                fprintf(frec, "%.6e%c", rec_time[it*nrec + ir], (ir==nrec-1)? '\n' : ' ');
            }
        }
        fclose(frec);
        printf("Saved seismograms.txt\n");
    }

    // Free memory
    free(sxx); free(szz); free(sxz); free(vx); free(vz);
    free(rho); free(lambda); free(mu); free(damp); free(rec_time);
    for(int i=0; i<num_time_slices; i++){
        free(depth_time_slices[i]);
    }
    free(depth_time_slices);
    free(slice_times);


     // 在程序结束前添加计时结束
    clock_t end_time = clock();
    double total_time = (double)(end_time - start_time) / CLOCKS_PER_SEC;
    printf("Total execution time: %.2f seconds\n", total_time);

    double t1 = omp_get_wtime();
    printf("Wall-clock time: %.2f seconds\n", t1 - clock_time);

    return 0;
}

关键解析：

pragma omp parallel for collapse(2)：

collapse(2)：将接下来的2层嵌套循环"折叠"成一个大的循环空间
本质上与#pragma omp parallel for相同，但更智能

特性	`#pragma omp parallel for`	`#pragma omp parallel for collapse(n)`
目标循环	仅下一层循环	紧接的 `n` 层嵌套循环
并行粒度	较粗（外层循环）	更细（合并后的循环空间）
负载均衡	当外层循环次数少时可能较差	通常更好，因为有更多任务可分配
适用场景	外层循环次数多，且内层计算量均匀	外层循环次数少，或总迭代空间大且需要更好负载均衡时
数据局部性	可能更好（一个线程处理连续的内存块）	可能稍差（线程可能处理不连续的内存位置）

注意（细节）：

1.gcc -fopenmp helloworld_openmp.c -o helloworld_openmp -lm 要加入数学库

2.#include <time.h> #include <math.h>这些库不要少了

2.2cuda

/***********************************************************************
CUDA accelerated version of 2D Elastic FDFD/FDTD
Author: Based on your CPU version
Purpose: educational, not production
************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <stddef.h>
#include <time.h>
#include <cuda_runtime.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define IDX(i,k,Nz) ((size_t)(i)*(Nz) + (k))

__device__ inline float ricker_gpu(float t, float f0){
    float x = M_PI*f0*t;
    float x2 = x*x;
    return (1.0f - 2.0f*x2) * expf(-x2);
}

/* ============ CUDA Kernels ============ */
__global__ void update_vx_kernel(
    float *vx, const float *sxx, const float *sxz, const float *rho,
    int Nx, int Nz, float dx, float dz, float dt)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int k = blockIdx.y * blockDim.y + threadIdx.y;
    if(i < 2 || i >= Nx-2 || k < 2 || k >= Nz-2) return;
    size_t id = IDX(i,k,Nz);

    float dsxx_dx = ( -sxx[IDX(i+2,k,Nz)] + 8.0f*sxx[IDX(i+1,k,Nz)]
                    - 8.0f*sxx[IDX(i-1,k,Nz)] + sxx[IDX(i-2,k,Nz)] ) / (12.0f*dx);
    float dsxz_dz = ( -sxz[IDX(i,k+2,Nz)] + 8.0f*sxz[IDX(i,k+1,Nz)]
                    - 8.0f*sxz[IDX(i,k-1,Nz)] + sxz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
    float invrho = 1.0f / rho[id];
    vx[id] += dt * invrho * (dsxx_dx + dsxz_dz);
}

__global__ void update_vz_kernel(
    float *vz, const float *sxz, const float *szz, const float *rho,
    int Nx, int Nz, float dx, float dz, float dt)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int k = blockIdx.y * blockDim.y + threadIdx.y;
    if(i < 2 || i >= Nx-2 || k < 2 || k >= Nz-2) return;
    size_t id = IDX(i,k,Nz);

    float dsxz_dx = ( -sxz[IDX(i+2,k,Nz)] + 8.0f*sxz[IDX(i+1,k,Nz)]
                    - 8.0f*sxz[IDX(i-1,k,Nz)] + sxz[IDX(i-2,k,Nz)] ) / (12.0f*dx);
    float dszz_dz = ( -szz[IDX(i,k+2,Nz)] + 8.0f*szz[IDX(i,k+1,Nz)]
                    - 8.0f*szz[IDX(i,k-1,Nz)] + szz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
    float invrho = 1.0f / rho[id];
    vz[id] += dt * invrho * (dsxz_dx + dszz_dz);
}

__global__ void update_stress_kernel(
    float *sxx, float *szz, float *sxz,
    const float *vx, const float *vz,
    const float *lambda, const float *mu,
    int Nx, int Nz, float dx, float dz, float dt)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int k = blockIdx.y * blockDim.y + threadIdx.y;
    if(i < 2 || i >= Nx-2 || k < 2 || k >= Nz-2) return;
    size_t id = IDX(i,k,Nz);

    float dvx_dx = ( -vx[IDX(i+2,k,Nz)] + 8.0f*vx[IDX(i+1,k,Nz)]
                   - 8.0f*vx[IDX(i-1,k,Nz)] + vx[IDX(i-2,k,Nz)] ) / (12.0f*dx);
    float dvz_dz = ( -vz[IDX(i,k+2,Nz)] + 8.0f*vz[IDX(i,k+1,Nz)]
                   - 8.0f*vz[IDX(i,k-1,Nz)] + vz[IDX(i,k-2,Nz)] ) / (12.0f*dz);
    float dvx_dz = ( -vx[IDX(i,k+2,Nz)] + 8.0f*vx[IDX(i,k+1,Nz)]
                   - 8.0f*vx[IDX(i,k-1,Nz)] + vx[IDX(i,k-2,Nz)] ) / (12.0f*dz);
    float dvz_dx = ( -vz[IDX(i+2,k,Nz)] + 8.0f*vz[IDX(i+1,k,Nz)]
                   - 8.0f*vz[IDX(i-1,k,Nz)] + vz[IDX(i-2,k,Nz)] ) / (12.0f*dx);

    float lam = lambda[id];
    float muval = mu[id];
    sxx[id] += dt * ( (lam + 2.0f*muval) * dvx_dx + lam * dvz_dz );
    szz[id] += dt * ( (lam + 2.0f*muval) * dvz_dz + lam * dvx_dx );
    sxz[id] += dt * ( muval * (dvx_dz + dvz_dx) );
}

__global__ void add_source_kernel(float *sxx, float *szz, int isrc, int ksrc, float src, float dt, int Nz){
    size_t id = IDX(isrc, ksrc, Nz);
    szz[id] += dt * src * 0.5f;
    sxx[id] += dt * src * 0.5f;
}

/* ============ CPU helpers ============ */
void write_pgm(const char *name, float *data, int Nx, int Nz){
    float mn = data[0], mx = data[0];
    for(int i=0;i<Nx*Nz;i++){ if(data[i]<mn) mn=data[i]; if(data[i]>mx) mx=data[i]; }
    float scale = (mx>mn)? 255.0f/(mx-mn) : 1.0f;
    FILE *f = fopen(name,"wb");
    fprintf(f,"P5\n%d %d\n255\n",Nx,Nz);
    unsigned char *buf = (unsigned char*)malloc(Nx*Nz);
    for(int k=0;k<Nz;k++){
        for(int i=0;i<Nx;i++){
            float v=data[IDX(i,k,Nz)];
            unsigned char c=(unsigned char)((v-mn)*scale);
            buf[k*Nx+i]=c;
        }
    }
    fwrite(buf,1,Nx*Nz,f);
    fclose(f);
    free(buf);
}

/* ============ Main ============ */
int main(){
    clock_t start_time = clock();
    const int Nx=301*25, Nz=301*25;
    const float dx=5.0f, dz=5.0f, dt=0.001f;
    const int nt=1000;

    const float rho0=2000.0f, vp0=3000.0f, vs0=1500.0f;
    const float mu0=rho0*vs0*vs0;
    const float lambda0=rho0*vp0*vp0-2.0f*mu0;

    const int isrc=Nx/2, ksrc=Nz/2;
    const float f0=20.0f, t0=1.0f/f0, src_amp=1e9f;

    size_t nxy=(size_t)Nx*Nz;
    float *sxx,*szz,*sxz,*vx,*vz,*rho,*lambda,*mu;
    cudaMallocManaged(&sxx,nxy*sizeof(float));
    cudaMallocManaged(&szz,nxy*sizeof(float));
    cudaMallocManaged(&sxz,nxy*sizeof(float));
    cudaMallocManaged(&vx ,nxy*sizeof(float));
    cudaMallocManaged(&vz ,nxy*sizeof(float));
    cudaMallocManaged(&rho,nxy*sizeof(float));
    cudaMallocManaged(&lambda,nxy*sizeof(float));
    cudaMallocManaged(&mu,nxy*sizeof(float));
    for(size_t id=0;id<nxy;id++){ sxx[id]=szz[id]=sxz[id]=vx[id]=vz[id]=0.0f;
        rho[id]=rho0; mu[id]=mu0; lambda[id]=lambda0; }

    dim3 block(16,16);
    dim3 grid((Nx+block.x-1)/block.x,(Nz+block.y-1)/block.y);

    for(int it=0; it<nt; it++){
        float time=it*dt;
        update_vx_kernel<<<grid,block>>>(vx,sxx,sxz,rho,Nx,Nz,dx,dz,dt);
        update_vz_kernel<<<grid,block>>>(vz,sxz,szz,rho,Nx,Nz,dx,dz,dt);
        update_stress_kernel<<<grid,block>>>(sxx,szz,sxz,vx,vz,lambda,mu,Nx,Nz,dx,dz,dt);

        float src = src_amp * ((1.0f - 2.0f*powf(M_PI*f0*(time-t0),2))*expf(-powf(M_PI*f0*(time-t0),2)));
        add_source_kernel<<<1,1>>>(sxx,szz,isrc,ksrc,src,dt,Nz);

        cudaDeviceSynchronize();

        if(it%100==0) printf("Step %d/%d\n",it,nt);
        if(it%200==0){
            char name[64];
            snprintf(name,sizeof(name),"snap_t%04d.pgm",it);
            write_pgm(name,vz,Nx,Nz);
        }
    }

    cudaFree(sxx); cudaFree(szz); cudaFree(sxz);
    cudaFree(vx); cudaFree(vz);
    cudaFree(rho); cudaFree(lambda); cudaFree(mu);
    double elapsed = (double)(clock()-start_time)/CLOCKS_PER_SEC;
    printf("Total GPU runtime: %.2f s\n",elapsed);
    return 0;
}

结论

以下是原代码，用8个核，28个核（all），cuda的用时。其中total time =core_num*real_time,可以清楚的看到随core_num的增加，任务堵塞引起效率降低。

Total execution time: 73.25 seconds（oringnal）

8-core
Total execution time: 172.88 seconds
Wall-clock time: 22.01 seconds
28 core
Total execution time: 728.32 seconds
Wall-clock time: 27.19 seconds

Total GPU runtime: 83.46 s
cuda

some tools code

for f in snap*.pgm; do convert "$f" "${f%.pgm}.png"; done

总结

提示：这里对文章进行总结：
例如：以上就是今天要讲的内容，本文仅仅简单介绍了pandas的使用，而pandas提供了大量能使我们快速便捷地处理数据的函数和方法。