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;
}