网格数值计算中周期性边界条件的处理

时间：2015-06-20 22:02:21 收藏：0 阅读：739

涉及到网格的数值计算中，边界条件的处理总是一个比较烦人的东西。

一者，本来好好的逐格处理过程，到边界这里总是要中断一下，或者加两个if，或者for循环中要精心处理下标关系，以免混乱，动不动就给你来个数组越界。

二者，并行计算的时候，最好的情况是所有网格统一处理。尤其是在CUDA编程中，代码中出现if很有可能大大降低执行效率。如果只对中间的网格进行并行处理，那么边界点的网格，单独做处理边界的kernel太烦，拷贝回内存处理更是得不偿失。

最近几天正在学习Level Set算法，研究从网上下载的代码的时候，发现了一个对付周期性边界条件的奇技淫巧，动一动脑筋的话，应该也对处理其他情况的边界条件有所帮助。关键部分的代码贴在下面：

       for i=1:N
         ip(i)=i+1;
         im(i)=i-1;
       end
       im(1)=N;
       ip(N)=1;
%
%      begin simulation loop
       for iter=1:nit
           for i=1:N
             for j=1:N
              dmx=(phi(i,j)-phi(im(i),j))/h;                     % x backward difference
              dpx=(phi(ip(i),j)-phi(i,j))/h;                     % x forward difference
              dmy=(phi(i,j)-phi(i,im(j)))/h;                     % y backward difference
              dpy=(phi(i,ip(j))-phi(i,j))/h;                     % y forward difference
              convx=max(u(i,j),0)*dmx+min(u(i,j),0)*dpx;        % if u(i,j)>0, use backward difference dmx, or else use dpx
              convy=max(v(i,j),0)*dmy+min(v(i,j),0)*dpy;
              phin(i,j)=phi(i,j)-(convx+convy)*dt;         % advance by dt
            end
           end 
           phi=phin;                                             % update
%
%         Plotting
          contour(X,Y,phi,[0,0],‘r‘);
          axis([-1 1 -1 1])
          axis(‘square‘)
          pause(.001)
       end

完整版的代码可以参考这个网页中第四部分。：http://www.math.lsa.umich.edu/~psmereka/LEVELSET/levelsetcodes.html

在上面代码中，u、v、phi等都是方形矩阵，具体跟Level Set算法有关。重点看代码中的 1~6行，这里构造了两个一维数组，内容很简单，大体上是按坐标递增的，但是在边界点，变成了另一头的坐标：im(1)=N;ip(N)=1;

这样，在后面的计算中，直接从相应位置取出im或者ip的值，当做坐标使用就可以了，程序自然会去数组的另一头取需要的数值进来代入计算。如果是其他类型的边界条件，应该也有类似的办法进行处理，不过我现在还没有想到……

这样做的坏处是多出了一个一维数组，并且需要多一次取数值的操作。但是这个方法可以做到将运算完全向量化，到底能不能使计算更快还没有测试，但是起码可以让代码看起来更干净。

来自为知笔记(Wiz)

原文：http://www.cnblogs.com/metorm/p/4591009.html