CUDA Accelerated Heat Conduction

做这个小程序最初的动机是，为了解决Blobby生成的空间函数系统不连续的问题。详细的推导过程可以结合Lipschitz条件进行。根据推测，AfterBurn应该是没有解决这个问题，当局部地区粒子数目很低的时候。这个问题还可以归结为是CurvatureFlow的问题，在这里就不多说了。假设我们需要散射2D空间内的一个标量场，那么非常直接的就是用数值法求解PDE（ODE），无论是Forward Euler还是RK等等都可以做。下面是散射一个标量场计算的结果，dt = 0.001，

测试表明这种需要大量迭代的计算过程可以方便的映射到GPU上去进行计算，由于显存带宽要比系统总线带宽大的多，所以复制状态就会非常的高效率。经过散射后的场将逐渐的填补不连续的地方，当然当迭代次数大到一定程度后就没什么意义了，所以迭代次数是由制作人员控制的，我就没必要插手了。这里是CUDA的核心代码和部分执行代码，其实就是一个简单的中心差分与前向欧拉法罢了。

texture < float ,  2 , cudaReadModeElementType >  Tex2DRef;

__global__  void  HeatCond2D( float *  Data, unsigned  int  Width, unsigned  int  Height)
{
     int  BlockSize  =  blockDim.x * blockDim.y;
     int  GridOffset  =  blockIdx.x * BlockSize;
     int  Idx  =  ( threadIdx.y * blockDim.x  +  threadIdx.x )  +  GridOffset;

     // float2 ST = make_float2( float()+0.5, float()+0.5 );
     float  X  =  Idx  %  Width;
     float  Y  =  ( Idx  -  X )  /  Width;

     // Four corners
    float2 CoordC  =  make_float2(X + 0.5 , Y + 0.5 );

     float  C  =  tex2D( Tex2DRef, CoordC.x, CoordC.y );
     float  L  =  tex2D( Tex2DRef, CoordC.x  -   1.0 , CoordC.y );
     float  R  =  tex2D( Tex2DRef, CoordC.x  +   1.0 , CoordC.y );
     float  U  =  tex2D( Tex2DRef, CoordC.x, CoordC.y  +   1.0  );
     float  D  =  tex2D( Tex2DRef, CoordC.x, CoordC.y  -   1.0  );
    
     float  Laplacian  =   - 4.0 * C  +  L  +  R  +  U  +  D;

    Data[Idx]  =  C  +  Laplacian * 0.01 ;
}