7.3_热传导模拟
7.3 热传导模拟
物理模拟问题或许是在计算上最具挑战性的问题之一。这类问题通常在计算精度与计算复杂性上存在着某种权衡。近年来，计算机模拟正变得越来越重要，这在很大程度上要归功于计算精度的不断提高，而这正是并行计算革命带来的好处。由于许多物理模拟计算都可以很容易地并行化，因此我们将在这个示例中看到一种非常简单的模拟模型。
7.3.1 简单的传热模型
为了说明一种可以有效使用纹理内存的情况，我们将构造一个简单的二维热传导模拟。首先假设有一个矩形房间，将其分成一个格网。在格网中随机散布一些“热源”，它们有着不同的固定温度。在图7.2中给出了这个房间的示意图。
图7.2 一个带有不同温度“热源”的房间
在给定了矩形格网以及热源分布后，我们可以计算格网中每个单元的温度随时间的变化情况。为了简单，热源单元本身的温度将保持不变。在时间递进的每个步骤中，我们假设热量在某个单元及其邻接单元之间“流动”。如果某个单元的邻接单元的温度更高，那么热量将从邻接单元传导到该单元。相反地，如果某个单元的温度比邻接单元的温度高，那么它将变冷。图7.3给出了这种热量流动示意图。
图7.3 热量从更热的单元传导到更冷的单元
在热传导模型中，我们对单元中新温度的计算方法为，将单元与邻接单元的温差相加起来，然后加上原有温度，计算公式如等式7.1所示。
等式7.1
在上面计算单元温度的等式中，常量k表示模拟过程中热量的流动速率。k值越大，表示系统会更快地达到稳定温度，而k值越小，则温度梯度将存在更长的时间。由于我们只考虑4个邻接单元（上，下，左，右）并且等式中的 和 都是常数，因此把等式7.1展开后将如等式7.2所示。
等式7.2
就像第6章的光线跟踪示例一样，在产品代码中不能使用这个模型（事实上，它与实际的物理热传导情况相差很远）。我们极大地简化了这个模型，从而将重点放到所使用的技术上。因此，让我们来看看如何在GPU上实现等式7.2中的更新。
7.3.2 温度更新的计算
我们稍后将介绍每个步骤的具体细节，首先给出更新流程的基本介绍：
1）给定一个包含初始输入温度的格网，将其中作为热源的单元温度值复制到格网相应的单元中。这将覆盖这些单元之前计算出的温度，因此也就确保了“加热单元将保持恒温”这个条件。这个复制操作是在copy_const_kernel()中执行的。
2）给定一个输入温度格网，根据等式7.2中的更新公式计算输出温度格网。这个更新操作是在blend_kernel()中执行的。
3）将输入温度格网和输出温度格网交换，为下一个步骤的计算做好准备。当模拟下一个时间步时，在步骤2中计算得到的输出温度格网将成为步骤1中的输入温度格网。
在开始模拟之前，我们假设已经获得了一个格网。格网中大多数单元的温度值都是0，但有些单元包含了非0的温度值，这些单元就是拥有固定温度的热源。在模拟过程中，缓冲区中这些常量值不会发生变化，并且在每个时间步中读取。
根据我们对热传导的建模方式，首先获得前一个时间步的输出温度格网并将其作为当前时间步的输入温度格网。然后，根据步骤1，将作为热源的单元的温度值复制到输出格网中，从而覆盖该单元之前计算出的温度值。这么做是因为我们已经假设这些热源单元的温度将保持不变。通过以下核函数将格网中的热源单元复制到输入格网中：
__global__void copy_const_kernel(float *iptr, const float *cptr) { //将threadIdx/BlockIdx映射到像素位置 int $\mathbf{x} =$ threadIdx.x $^+$ blockIdx.x \*blockDim.x; int $\mathrm{y} =$ threadIdx.y $^+$ blockIdx.y \*blockDim.y; int offset $=$ x $^+$ y \*blockDim.x \*gridDim.x; if (cptr[offset] != 0) iptr[offset] $=$ cptr[offset]; }前三行代码看上去非常熟悉。前两行代码将线程的threadIdx和blockIdx转换为x坐标和y坐标。第三行代码计算在输入缓冲区中的线性偏移。在加粗的一行中把cptr[]中的热源温度复制到iptr[]的输入单元中。注意，只有当格网中单元的温度值非0时，才会执行复制操作。我们这么做是为了维持非热源单元在上一个时间步中计算得到的温度值。热源单元在cptr[]中对应的元素为非零值，因此调用这个复制核函数后，热源单元的温度值在不同的时间步中将保持不变。
算法步骤2中包含的计算最多。为了执行这些更新操作，可以在模拟过程中让每个线程都负责计算一个单元。每个线程都将读取对应单元及其邻接单元的温度值，执行前面给出的更新运算，然后用计算得到的新值更新它的温度。这个核函数的大部分代码与之前使用的技术都是类似的。
__global__void blend_kernel(float *outSrc, const float *inSrc) { //将threadIdx/BlockIdx映射到像素位置 int $\mathbf{x} =$ threadIdx.x $^+$ blockIdx.x \*blockDim.x; int $\mathrm{y} =$ threadIdx.y $^+$ blockIdx.y \*blockDim.y; int offset $=$ x $^+$ y \* blockDim.x \*gridDim.x; int left $=$ offset-1; int right $=$ offset+1; if $(\mathbf{x} = = 0)$ left++; if $(\mathbf{x} = =$ DIM-1) right--; int top $=$ offset-DIM; int bottom $=$ offset+DIM; if $(\mathbf{y} = = 0)$ top $= =$ DIM; if $(\mathbf{y} = =$ DIM-1) bottom $= =$ DIM; outSrc[offset] $=$ inSrc[offset] $^+$ SPEED \* (inSrc[top] + inSrc(bottom] $^+$ inSrc[left] $^+$ inSrc(right] - inSrc[offset]*4);
}注意，代码的开头与生成输出图像示例的开头是一样的，只是此时线程不是计算像素的颜色值，而是计算模拟单元中的温度值。首先，代码将threadIdx和blockIdx转换为x，y和offset。现在，你可以在睡觉的时候都能复诵这些代码行（但我们并不真的希望你在睡觉时还想到它们）。
接下来，代码会计算出上，下，左，右四个邻接单元的偏移并读取这些单元的温度。我们需要这些值来计算当前单元中的已更新温度。这里唯一需要注意的地方就是需要调整边界的索引，从而在边缘的单元上不会回卷。最后，在加粗的代码行中执行了等式7.2的更新运算，将原来的温度与该单元温度和邻接单元的温度的温差相加。
7.3.3 模拟过程动态演示
剩下的代码主要是设置好单元，然后显示热量的动画输出。下面给出了这段代码：
include“cuda.h" #include“../common/book.h" #include“../common/cpu anim.h" #define DIM 1024 #define PI 3.1415926535897932f #define MAX_TEMP 1.0f #define MIN_TEMP 0.0001f #define SPEED 0.25f//更新函数中需要的全局变量
struct DataBlock{ unsigned char \*output bitmap; float \*dev_inSrc; float \*dev_outSrc; float \*dev_constSrc; CPUAnimBitmap \* bitmap; JudaEvent_t start,stop; float totalTime; float frames;
}；
void anim_gpu(DataBlock \*d，int ticks）{ HANDLE_ERROR( CUDAEventRecord(d->start，0)）; dim3 blocks(DIM/16,DIM/16); dim3 threads(16,16); CPUAnimBitmap \* bitmap $=$ d-> bitmap;for (int i=0; i<90; i++) { copy_const_kernel<<blocks, threads>>>(d->dev_inSrc, d->dev_constSrc); blend_kernel<<blocks, threads>>>(d->dev_outSrc, d->dev_inSrc); swap(d->dev_inSrc, d->dev_outSrc); } float_to_color<<blocks, threads>>>(d->output_mask, d->dev_inSrc); HANDLE_ERROR(udaMemcpy( bitmap->get_ptr(), d->output_mask, bitmap->image_size(),udaMemcpyDeviceToHost)); HANDLE_ERROR(udaEventRecord(d->stop, 0)); HANDLE_ERROR(udaEventSynchronize(d->stop)); float elapsedTime; HANDLE_ERROR(udaEventElapsedTime(& elapsedTime, d->start, d->stop)); d->totalTime += elapsedTime; ++d->frames; printf("Average Time per frame: %3.1f ms\n", d->totalTime/d->frames); } void anim_exit(DataBlock *d) {udaFree(d->dev_inSrc);udaFree(d->dev_outSrc);udaFree(d->dev_constSrc); HANDLE_ERROR(udaEventDestroy(d->start)); HANDLE_ERROR(udaEventDestroy(d->stop)); }在代码中加入了基于事件的计时功能，这与第6章中光线跟踪示例一样。示例中计时代码的作用与前面的示例相同：由于我们希望提高算法的执行速度，因此在代码中添加了测量性能的机制从而确保成功提升了性能。
动画框架的每一帧都将调用函数anim_gpu()。这个函数的参数是一个指向DataBlock的指针，以及动画已经经历的时间 ticks。在动画示例中，每个线程块包含了256个线程，构成了一个16×16的二维单元。anim_gpu()中for()循环的每次迭代都将计算模拟过程的一个时间步，我们在7.3.2节的开头介绍了这个算法，共包含三个步骤。由于DataBlock包含了表示热源的常量缓冲区，以及在上一个时间步中计算得到的输出值，因此能够表示动画的完整状态，因而
animＧpU(实际上并不需要用到 ticks的值。
值得注意的是，在每帧中将执行90个时间步。这个数值并非有着特殊含义，而是通过实验得出的，这个数值既能避免在每个时间步中都需要复制一张位图图片，又可以避免在每一帧计算过多时间步（这将导致不稳定的动画）。如果你更关心每个模拟步骤的输出结果，而不是实时地播放动画结果，那么可以将这个值修改为在每帧中只计算一个时间步。
在计算了90个时间步后，anim_gpu()就已经准备好将当前动画的位图帧复制回CPU。由于for()循环会交换输入和输出，因此我们首先交换输入缓冲区和输出缓冲区，这样在输出缓冲区中实际上将包含第90个时间步的输出。我们通过核函数float_to_color()将温度转换为颜色，然后将结果图像通过cudaMemcpy()复制回CPU，并将复制的方向指定为cudaMemcpyDeviceToHost。最后，为了准备下一系列时间步，我们将输出缓冲区交换到输入缓冲区以便作为下一个时间步的输入。
int main(void){ DataBlock data; CPUAnimBitmap bitmap(DIM，DIM,&data); data.bitmap $=$ &bitmap; data.totalTime $= 0$ · dataFrames $= 0$ ： HANDLE_ERROR(udaEventCreate(&data.start）); HANDLE_ERROR(udaEventCreate(&data.stop）); HANDLE_ERROR(udaMalloc( $\mathrm{(void**)}$ )&data.output.bitmap, bitmap(image_size()）); //假设float类型的大小为4个字符（即rgbba) HANDLE_ERROR(udaMalloc( $\mathrm{(void**)}$ )&data.dev_inSrc, bitmap(image_size()））； HANDLE_ERROR(udaMalloc( $\mathrm{(void**)}$ )&data.dev_outSrc, bitmap(image_size()））； HANDLE_ERROR(udaMalloc( $\mathrm{(void**)}$ )&data.dev_constSrc, bitmap(image_size()））; float \*temp $=$ (float\*)malloc(bitmap(image_size())）; for（int $\mathrm{i} = 0$ ；i<DIM\*DIM；i++){ temp[i] $= 0$ int $\mathbf{x} = \mathbf{i}\%$ DIM; int $\mathrm{y} = \mathrm{i}/$ DIM; if((x>300)&& $(\mathrm{x} < 600)$ && $(\mathrm{y} > 310)$ && $(\mathrm{y} < 601)$ temp[i] $=$ MAX_TEMP; } temp[DM\*100+100] $=$ (MAX_TEMP + MIN_TEMP)/2;temp[DIM*700+100] = MIN_TEMP;
temp[DIM*300+300] = MIN_TEMP;
temp[DIM*200+700] = MIN_TEMP;
for (int y=800; y<900; y++) {
for (int x=400; x<500; x++) {
temp[x+y*DIM] = MIN_TEMP;
}
}
HANDLE_ERROR(udaMemcpy(data.dev_constSrc, temp, bitmap(image_size(),udaMemcpyHostToDevice));
for (int y=800; y<DIM; y++) {
for (int x=0; x<200; x++) {
temp[x+y*DIM] = MAX_TEMP;
}
}
HANDLE_ERROR(udaMemcpy(data.dev_inSrc, temp, bitmap(image_size(),udaMemcpyHostToDevice));
free(temp);
heap志愿和exit((void (*)(void*, int))anim_gpu, (void (*)(void*))anim_exit);图7.4给出了输出图像的一个示例。注意在这张图像中，某些“热源”看上去像一些像素大小的岛屿，它们打乱了温度分布的连续性。
图7.4 热传导模拟动画的截屏
7.3.4 使用纹理内存
在温度更新计算的内存访问模式中存在着巨大的内存空间局部性。前面曾解释过，这种访问模式可以通过GPU纹理内存来加速。接下来，我们来学习如何使用纹理内存。
首先，需要将输入的数据声明为texture类型的引用。我们使用浮点类型纹理的引用，因为温度数值是浮点类型的。
// 这些变量将位于GPU上
texture<float> texConstSrc;
texture<float> texIn;
texture<float> texOut;下一个需要注意的问题是，在为这三个缓冲区分配了GPU内存后，需要通过CUDABindTexture()将这些变量绑定到内存缓冲区。这相当于告诉CUDA运行时两件事情：
我们希望将指定的缓冲区作为纹理来使用。
我们希望将纹理引用作为纹理的“名字”。
在热传导模拟中分配了这三个内存后，需要将这三个内存绑定到之前声明的纹理引用（texConstSrc、texIn、texOut）。
HANDLE_ERROR(udaMalloc(void\*\*)&data.dev_inSrc, imageSize）);
HANDLE_ERROR(udaMalloc(void\*\*)&data.dev_outSrc, imageSize）);
HANDLE_ERROR(udaMalloc(void\*\*)&data.dev_constSrc, imageSize）);
HANDLE_ERROR(udaBindTexture(NULL，texConstSrc, data.dev_constSrc, imageSize）);
HANDLE_ERROR(udaBindTexture(NULL，texIn, data.dev_inSrc, imageSize）);
HANDLE_ERROR(udaBindTexture(NULL，texOut, data.dev_outSrc, imageSize）);此时，纹理变量已经设置好了，现在可以启动核函数。然而，当读取核函数中的纹理时，需要通过特殊的函数来告诉GPU将读取请求转发到纹理内存而不是标准的全局内存。因此，当读取内存时不再使用方括号从缓冲区中读取，而是将blend_kernel()改为使用tex1Dfetch()。
此外，在全局内存和纹理内存的使用上还存在另一个差异。虽然tex1Dfetch()看上去像一个
函数，但它其实是一个编译器内置函数（Intrinsic）。由于纹理引用必须声明为文件作用域内的全局变量，因此我们不再将输入缓冲区和输出缓冲区作为参数传递给blend_kernel()，因为编译器需要在编译时知道texlDfetch()应该对哪些纹理采样。我们不是像前面那样传递指向输入缓冲区和输出缓冲区的指针，而是将一个布尔标志dstOut传递给blend_kernel()，这个标志会告诉我们使用哪个缓冲区作为输入，以及哪个缓冲区作为输出。下面是对blend_kernel()的修改：
__global__ void blend_kernel(float *dst, bool dstOut) {
// 将threadIdx/BlockIdx映射到像素位置
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
int left = offset - 1;
int right = offset + 1;
if (x == 0) left++;
if (x == DIM-1) right--;
int top = offset - DIM;
int bottom = offset + DIM;
if (y == 0) top += DIM;
if (y == DIM-1) bottom -= DIM;
float t, l, c, r, b;
if (dstOut) {
t = tex1Dfetch(txIn, top);
l = tex1Dfetch(txIn, left);
c = tex1Dfetch(txIn, offset);
r = tex1Dfetch(txIn, right);
b = tex1Dfetch(txIn, bottom);
} else {
t = tex1Dfetch(txOut, top);
l = tex1Dfetch(txOut, left);
c = tex1Dfetch(txOut, offset);
r = tex1Dfetch(txOut, right);
b = tex1Dfetch(txOut, bottom);
}
dst[offset] = c + SPEED * (t + b + r + l - 4 * c);由于核函数copy_const_kernel()将读取包含热源位置和温度的缓冲区，因此同样需要修改为从纹理内存而不是从全局内存中读取：
__global__ void copy_const_kernel(float *iptr) {
// 将threadIdx/BlockIdx映射到像素位置
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
float c = tex1DfetchTEXConstSrc, offset);
if (c != 0)
iptr[offset] = c;
}由于blend_kernel()的函数原型被修改为接收一个标志，并且这个标志表示在输入缓冲区与输出缓冲区之间的切换，因此需要对anim_gpio()函数进行相应的修改。现在，不是交换缓冲区，而是在每组调用之后通过设置dstOut = !dstOut来进行切换：
void anim_gpu(DataBlock \*d,int ticks）{
HANDLE_ERROR( CUDAEventRecord( d->start,0)）;
dim3 blocks(DIM/16,DIM/16);
dim3 threads(16,16);
CPUAnimBitmap \* bitmap $=$ d-> bitmap;
//由于tex是全局并且有界的，因此我们必须通过一个标志来
//选择每次迭代中哪个是输入/输出
volatile bool dstOut $=$ true;
for (int i=0；i<90；i++) { float \*in，\*out; if (dstOut){ in $=$ d->dev_inSrc; out $=$ d->dev_outSrc; } else { out $=$ d->dev_inSrc; in $=$ d->dev_outSrc; } copy_const_kernel<<blocks,threads>>>(in); blend_kernel<<blocks,threads>>>(out,dstOut); dstOut $=$ !dstOut;
} float_to_color<<blocks,threads>>>(d->output bitmap, d->dev_inSrc);
HANDLE_ERROR( CUDAMemcpy( bitmap->get_ptr(), d->output bitmap, bitmap->image_size(), CUDAMempyDeviceToHost )）;HANDLE_ERROR(udaEventRecord(d->stop, 0));
HANDLE_ERROR(udaEventSynchronize(d->stop));
float elapsedTime;
HANDLE_ERROR(udaEventElapsedTime(& elapsedTime, d->start, d->stop));
d->totalTime += elapsedTime;
++d->frames;
printf("Average Time per frame:%3.1f ms\n", d->totalTime/d->frames);对热传导函数的最后一个修改就是在应用程序运行结束后的清理工作。不仅要释放全局缓冲区，还需要清除与纹理的绑定：
// 释放在GPU上分配的内存
void anim_exit(DataBlock *d) {
canadaUnbindTexture(texIn);
canadaUnbindTextureTEXOut);
canadaUnbindTextureTEXConstSrc);
canadaFree(d->dev_inSrc);
canadaFree(d->dev_outSrc);
canadaFree(d->dev_constSrc);
HANDLE_ERROR(cudaEventDestroy(d->start));
HANDLE_ERROR(cudaEventDestroy(d->stop));
}7.3.5 使用二维纹理内存
在本书的开头，我们提到了在某些二维问题中使用二维的线程块和线程格是非常有用的。对于纹理内存来说同样如此。在许多情况下，二维内存空间是非常有用的，它们对于熟悉标准C中多维数组的程序员来说都不会感到陌生。我们来看看如何将热传导应用程序修改为使用二维纹理。
首先，需要修改纹理引用的声明。如果没有具体说明，默认的纹理引用都是一维的，因此我们增加了代表维数的参数2，这表示声明的是一个二维纹理引用。
texture<float,2> texConstSrc;
texture<float,2> texIn;
texture<float,2> texOut;二维纹理将简化blend_kernel()方法的实现。虽然我们需要将tex1Dfetch()调用修改为tex2D()调用，但却不再需要使用通过线性化offset变量以计算top、left、right和bottom等偏移。当使用二维纹理时，可以直接通过x和y来访问纹理。
而且，当使用tex2D()时，我们不再需要担心发生溢出问题。如果x或y小于0，那么tex2D()将返回0处的值。同理，如果某个值大于宽度，那么tex2D()将返回位于宽度处的值。注意，在我们的应用程序刚好需要这种行为，但在其他应用程序中可能需要其他的行为。
这些简化带来的好处之一就是，核函数的代码将变得更加简单。
__global__void blend_kernel(float *dst, bool dstOut) {
//将threadIdx/BlockIdx映射到像素位置
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
float t, l, c, r, b;
if (dstOut) {
t = tex2DTEXIn,x,y-1);
l = tex2DTEXIn,x-1,y);
c = tex2DTEXIn,x,y);
r = tex2DTEXIn,x+1,y);
b = tex2DTEXIn,x,y+1);
} else {
t = tex2DTEXOut,x,y-1);
l = tex2DTEXOut,x-1,y);
c = tex2DTEXOut,x,y);
r = tex2DTEXOut,x+1,y);
b = tex2DTEXOut,x,y+1);
}
dst[offset] = c + SPEED * (t + b + r + 1 - 4 * c);
}由于所有之前对tex1Dfetch()的调用都需要修改为对tex2D()的调用，因此我们需要在copy_const_kernel()中进行相应的修改。与核函数blend_kernel()类似的是，我们不再需要通过offset来访问纹理，而只需使用x和y来访问热源的常量数量：
__global__void copy_const_kernel(float \*iptr）{ //将threadIdx/BlockIdx映射到像素位置 int $x =$ threadIdx.x $^+$ blockIdx.x \*blockDim.x; int $y =$ threadIdx.y $^+$ blockIdx.y \*blockDim.y; int offset $= x + y$ \*blockDim.x \*gridDim.x; floatc $\equiv$ tex2D(texConstSrc,x,y); if（c！ $= 0$ ） iptr[offset] $= c$在使用一维纹理的热传导模拟版本中还剩下一些最后的修改，它们与之前的修改是类似的。在main()中需要对纹理绑定调用进行修改，并告诉运行时：缓冲区将被视为二维纹理而不是一维纹理。
HANDLE_ERROR(udaMalloc( (void**)&data.dev_inSrc, imageSize) );
HANDLE_ERROR(udaMalloc( (void**)&data.dev_outSrc, imageSize) );
HANDLE_ERROR(udaMalloc( (void**)&data.dev_constSrc, imageSize) );
cudaChannelFormatDesc desc =udaCreateChannelDesc<float>();
HANDLE_ERROR(udaBindTexture2D(NULL, texConstSrc, data.dev_constSrc, desc, DIM, DIM, sizeof(float)*DIM) );
HANDLE_ERROR(udaBindTexture2D(NULL, texIn, data.dev_inSrc, desc, DIM, DIM, sizeof(float)*DIM) );
HANDLE_ERROR(udaBindTexture2D(NULL, texOut, data.dev_outSrc, desc, DIM, DIM, sizeof(float)*DIM) );与不使用纹理内存和使用一维纹理内存相同的是，都需要为输入数组事先分配存储空间。接下来的与一维纹理示例有所不同，因为当绑定二维纹理时，CUDA运行时要求提供一个CUDAChannelFormatDesc。在上面的代码中包含了一个对通道格式描述符（Channel Format Descriptor）的声明。在这里可以使用默认的参数，并且只要指定需要的是一个浮点描述符。然后，我们通过CUDABindTexture2D()，纹理的维数（DIM×DIM）以及通道格式描述符（desc）将这三个输入缓冲区绑定为二维纹理。main()函数的其他部分保持不变。
int main(void) {
DataBlock data;
CPUAnimBitmap bitmap(DIM, DIM, &data);
data.bitmap = &bitmap;
data.totalTime = 0;
dataFrames = 0;
HANDLE_ERROR(cudaEventCreate(&data.start));
HANDLE_ERROR(cudaEventCreate(&data.stop));int imageSize $=$ bitmap(image_size();
HANDLE_ERROR(udaMalloc( (void\*\*)&data.output bitmap, imageSize）);
//假设float类型的大小为4个字符（即rgba)
HANDLE_ERROR(udaMalloc( (void\*\*)&data.dev_inSrc, imageSize）);
HANDLE_ERROR(udaMalloc( (void\*\*)&data.dev_outSrc, imageSize）);
HANDLE_ERROR(udaMalloc( (void\*\*)&data.dev_constSrc, imageSize）);
cudaChannelFormatDesc desc $=$ CUDACreateChannelDesc(float $\gimel$ ）;
HANDLE_ERROR(udaBindTexture2D(NULL，texConstSrc, data.dev_constSrc, desc，DIM，DIM, sizeof(float)*DIM）);
HANDLE_ERROR(udaBindTexture2D(NULL，texIn, data.dev_inSrc, desc，DIM，DIM, sizeof(float)\*DIM）;
HANDLE_ERROR(udaBindTexture2D(NULL，texOut, data.dev_outSrc, desc，DIM，DIM, sizeof(float)\*DIM）;//初始化常量数据
float \*temp $=$ (float\*)malloc( imageSize);
for (int i=0; i<DIM\*DIM; i++) { temp[i] = 0; int x = i & DIM; int y = i / DIM; if ((x>300) && (x<600) && (y>310) && (y<601)) temp[i] = MAX_TEMP;
}
temp[DM\*100+100] = (MAX_TEMP + MIN_TEMP)/2;
temp[DM\*700+100] = MIN_TEMP;
temp[DM\*300+300] = MIN_TEMP;
temp[DM\*200+700] = MIN_TEMP;
for (int y=800; y<900; y++) { for (int x=400; x<500; x++) { temp[x+y\*DIM] = MIN_TEMP;}
}
HANDLE_ERROR(udaMemcpy( data.dev_constSrc，temp, imageSize,udaMemcpyHostToDevice）);
//初始化输入数据
for（int $\mathrm{y = 800}$ ：y<DIM；y++) { for（int $\mathbf{x} = 0$ ：x<200； $x + +$ ）{ temp[x+y\*DIM] $\equiv$ MAX_TEMP; }
}
HANDLE_ERROR(udaMemcpy( data.dev_inSrc，temp, imageSize,udaMemcpyHostToDevice））; free（temp）;
bitmap.anim_and_exit( (void $(^{*})$ (void\*,int))anim_gpu, (void $(^{*})$ (void\*))anim_exit）;虽然我们需要通过不同的函数来告诉运行时绑定一维纹理还是绑定二维纹理，但可以通过同一个函数来取消纹理的绑定，即cudaUnbindTexture()。正因为如此，执行释放操作的函数可以保持不变。
// 释放在GPU上分配的内存
void anim_exit(DataBlock \*d){ CUDAUnbindTexture( texIn ); CUDAUnbindTexture( texOut ); CUDAUnbindTexture( texConstSrc ); CUDAFree( d->dev_inSrc ); CUDAFree( d->dev_outSrc ); CUDAFree( d->dev_constSrc ); HANDLE_ERROR( CUDAEventDestroy( d->start ) ); HANDLE_ERROR( CUDAEventDestroy( d->stop ) );
}无论使用二维纹理还是一维纹理，热传导模拟应用程序的性能基本相同。因此，基于性能来选择使用一维纹理还是二维纹理可能没有太大的意义。对于这个应用程序来说，当使用二维纹理时，代码会更简单一些，因为模拟的问题刚好是二维的。但通常来说，如果问题不是二维的，那么究竟选择一维纹理还是二维纹理要视具体情况而定。