8489阅读
33回复

[RFEDA原创]Python的实现

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

20楼发表于: 2010-07-30 17:01:55

fd2d_3.3.c平面波源：

from pylab import *
from enthought.mayavi import mlab

NSTEPS         = 115; NPML        = 8
ddx = 0.01; dt = ddx/6.0e8; epsz = 8.8e-12
T   = 0;    t0 = 20.0;     spread = 8.0

IE = 60;   ic = IE/2; ia = 7; ib = IE-ia-1
JE = 60;   jc = JE/2; ja = 7; jb = JE-ja-1

dz = zeros((IE,JE)); hx = zeros((IE,JE)); ihx = zeros((IE,JE))
ez = zeros((IE,JE)); hy = zeros((IE,JE)); ihy = zeros((IE,JE))
ga = ones((IE,JE))
gi2 = ones(IE); gi3 = ones(IE); fi1= zeros(IE); fi2 = ones(IE); fi3 = ones(IE)
gj2 = ones(JE); gj3 = ones(JE); fj1= zeros(JE); fj2 = ones(JE); fj3 = ones(JE)

ez_inc = zeros(JE); ez_inc_low_m1 = 0.0; ez_inc_high_m1 = 0.0
hx_inc = zeros(JE); ez_inc_low_m2 = 0.0; ez_inc_high_m2 = 0.0

xnum = NPML-arange(NPML)

xn = 0.33*(xnum/float(NPML))**3
gi2[0:NPML] =      1.0/(1.0+xn); gi2[-1:-1-NPML:-1] =      1.0/(1.0+xn)
gi3[0:NPML] = (1.0-xn)/(1.0+xn); gi3[-1:-1-NPML:-1] = (1.0-xn)/(1.0+xn)
gj2[0:NPML] =      1.0/(1.0+xn); gj2[-1:-1-NPML:-1] =      1.0/(1.0+xn)
gj3[0:NPML] = (1.0-xn)/(1.0+xn); gj3[-1:-1-NPML:-1] = (1.0-xn)/(1.0+xn)

xn = 0.25*((xnum-0.5)/float(NPML))**3
fi1[0:NPML] =               xn; fi1[-2:-2-NPML:-1] =               xn
fi2[0:NPML] =      1.0/(1.0+xn); fi2[-2:-2-NPML:-1] =      1.0/(1.0+xn)
fi3[0:NPML] = (1.0-xn)/(1.0+xn); fi3[-2:-2-NPML:-1] = (1.0-xn)/(1.0+xn)
fj1[0:NPML] =               xn; fj1[-2:-2-NPML:-1] =               xn
fj2[0:NPML] =      1.0/(1.0+xn); fj2[-2:-2-NPML:-1] =      1.0/(1.0+xn)
fj3[0:NPML] = (1.0-xn)/(1.0+xn); fj3[-2:-2-NPML:-1] = (1.0-xn)/(1.0+xn)

for n in range(NSTEPS):
    T += 1
    ez_inc[1:] += 0.5*(hx_inc[:-1]-hx_inc[1:])

    ez_inc[0]     = ez_inc_low_m2; ez_inc[JE-1]   = ez_inc_high_m2
    ez_inc_low_m2 = ez_inc_low_m1; ez_inc_high_m2 = ez_inc_high_m1
    ez_inc_low_m1 = ez_inc[1];     ez_inc_high_m1 = ez_inc[JE-2]

    for i in range(1,IE):
        for j in range(1,JE):
            dz[j] = gi3*gj3[j]*dz[j]+gi2*gj2[j]*0.5*(hy[j]-hy[i-1][j]-hx[j]+hx[j-1])

    pulse = exp(-0.5*((t0-T)/spread)**2); ez_inc[3] = pulse;

    for i in range(ia,ib+1):
        dz[ja] += 0.5*hx_inc[ja-1]
        dz[jb] += -0.5*hx_inc[jb]

    ez = ga*dz

    for i in range(IE-1): ez[0] = 0.0; ez[-1] = 0.0
    for j in range(JE-1): ez[0][j] = 0.0; ez[-1][j] = 0.0

    hx_inc[:-1] += 0.5*(ez_inc[:-1]-ez_inc[1:])

    for i in range(IE):
        for j in range(JE-1):
            curl_e   = ez[j]-ez[j+1]; ihx[j] += fi1*curl_e
            hx[j] = fj3[j]*hx[j]+fj2[j]*0.5*(curl_e+ihx[j])

    for i in range(ia,ib+1):
        hx[ja-1] += 0.5*ez_inc[ja]
        hx[jb]   += -0.5*ez_inc[jb]

    for i in range(IE-1):
        for j in range(JE-1):
            curl_e   = ez[i+1][j]-ez[j]; ihy[j] += fj1[j]*curl_e
            hy[j] = fi3*hy[j]+fi2*0.5*(curl_e+ihy[j])

    for j in range(ja,jb+1):
        hy[ia-1][j] += -0.5*ez_inc[j]
        hy[ib][j]   += 0.5*ez_inc[j]

x, y = np.ogrid[0:IE,0:JE]
pl = mlab.surf(x, y, ez, warp_scale="auto")

坚持到底就是胜利！

图片:fd2d_3.3.JPG

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

21楼发表于: 2010-08-02 18:01:51

fd2d_3.4.c中间存在圆柱介质：

from pylab import *
from enthought.mayavi import mlab

NSTEPS         = 100; NPML        = 8;
ddx = 0.01; dt = ddx/6.0e8; epsz = 8.8e-12
T   = 0;    t0 = 25.0;     spread = 8.0

IE = 60;   ic = IE/2; ia = 7; ib = IE-ia-1
JE = 60;   jc = JE/2; ja = 7; jb = JE-ja-1

dz = zeros((IE,JE)); hx = zeros((IE,JE)); ihx = zeros((IE,JE))
ez = zeros((IE,JE)); hy = zeros((IE,JE)); ihy = zeros((IE,JE))
ga = ones((IE,JE)); gb = zeros((IE,JE)); iz = zeros((IE,JE))
gi2 = ones(IE); gi3 = ones(IE); fi1= zeros(IE); fi2 = ones(IE); fi3 = ones(IE)
gj2 = ones(JE); gj3 = ones(JE); fj1= zeros(JE); fj2 = ones(JE); fj3 = ones(JE)

ez_inc = zeros(JE); ez_inc_low_m1 = 0.0; ez_inc_high_m1 = 0.0
hx_inc = zeros(JE); ez_inc_low_m2 = 0.0; ez_inc_high_m2 = 0.0

xnum = NPML-arange(NPML)

xn = 0.33*(xnum/float(NPML))**3
gi2[0:NPML] =      1.0/(1.0+xn); gi2[-1:-1-NPML:-1] =      1.0/(1.0+xn)
gi3[0:NPML] = (1.0-xn)/(1.0+xn); gi3[-1:-1-NPML:-1] = (1.0-xn)/(1.0+xn)
gj2[0:NPML] =      1.0/(1.0+xn); gj2[-1:-1-NPML:-1] =      1.0/(1.0+xn)
gj3[0:NPML] = (1.0-xn)/(1.0+xn); gj3[-1:-1-NPML:-1] = (1.0-xn)/(1.0+xn)

xn = 0.25*((xnum-0.5)/float(NPML))**3
fi1[0:NPML] =               xn; fi1[-2:-2-NPML:-1] =               xn
fi2[0:NPML] =      1.0/(1.0+xn); fi2[-2:-2-NPML:-1] =      1.0/(1.0+xn)
fi3[0:NPML] = (1.0-xn)/(1.0+xn); fi3[-2:-2-NPML:-1] = (1.0-xn)/(1.0+xn)
fj1[0:NPML] =               xn; fj1[-2:-2-NPML:-1] =               xn
fj2[0:NPML] =      1.0/(1.0+xn); fj2[-2:-2-NPML:-1] =      1.0/(1.0+xn)
fj3[0:NPML] = (1.0-xn)/(1.0+xn); fj3[-2:-2-NPML:-1] = (1.0-xn)/(1.0+xn)

radius = 10; epsilon = 30.0; sigma = 0.3

for j in range(ja,jb):
    for i in range(ia,ib):
        xdist = (ic-i)
        ydist = (jc-j)
        dist = sqrt(xdist**2+ydist**2)
        if dist<=radius:
            ga[j] = 1.0/(epsilon+sigma*dt/epsz)
            gb[j] = sigma*dt/epsz

for n in range(NSTEPS+1):
    T += 1
    ez_inc[1:] += 0.5*(hx_inc[:-1]-hx_inc[1:])

    ez_inc[0]     = ez_inc_low_m2; ez_inc[JE-1]   = ez_inc_high_m2
    ez_inc_low_m2 = ez_inc_low_m1; ez_inc_high_m2 = ez_inc_high_m1
    ez_inc_low_m1 = ez_inc[1];     ez_inc_high_m1 = ez_inc[JE-2]

    for i in range(1,IE):
        for j in range(1,JE):
            dz[j] = gi3*gj3[j]*dz[j]+gi2*gj2[j]*0.5*(hy[j]-hy[i-1][j]-hx[j]+hx[j-1])

    pulse = exp(-0.5*((t0-T)/spread)**2); ez_inc[3] = pulse;

    for i in range(ia,ib+1):
        dz[ja] += 0.5*hx_inc[ja-1]
        dz[jb] += -0.5*hx_inc[jb]

    ez = ga*(dz-iz)
    iz += gb*ez

    for i in range(IE-1): ez[0] = 0.0; ez[-1] = 0.0
    for j in range(JE-1): ez[0][j] = 0.0; ez[-1][j] = 0.0

    hx_inc[:-1] += 0.5*(ez_inc[:-1]-ez_inc[1:])

    for i in range(IE):
        for j in range(JE-1):
            curl_e   = ez[j]-ez[j+1]; ihx[j] += fi1*curl_e
            hx[j] = fj3[j]*hx[j]+fj2[j]*0.5*(curl_e+ihx[j])

    for i in range(ia,ib+1):
        hx[ja-1] += 0.5*ez_inc[ja]
        hx[jb]   += -0.5*ez_inc[jb]

    for i in range(IE-1):
        for j in range(JE-1):
            curl_e   = ez[i+1][j]-ez[j]; ihy[j] += fj1[j]*curl_e
            hy[j] = fi3*hy[j]+fi2*0.5*(curl_e+ihy[j])

    for j in range(ja,jb+1):
        hy[ia-1][j] += -0.5*ez_inc[j]
        hy[ib][j]   += 0.5*ez_inc[j]

x, y = np.ogrid[0:IE,0:JE]
pl = mlab.surf(x, y, ez, warp_scale="auto")

图片:fd2d_3.4.JPG

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

22楼发表于: 2010-08-04 13:29:44

希望楼主把该贴子整理成pdf文档，后放出来让版主加技术分或者威望
谢谢！

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

23楼发表于: 2010-08-04 13:39:09

我给大家整理了一个pdf，请下载！
帖子整理

python-fdtd.pdf (326 K) 下载次数:12

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

24楼发表于: 2010-08-05 17:27:10

3Q，我很大程度是对FDTD好奇，写代码也没有认真检查过。

关于python调c，有很多种方式，我觉着最好的方式是把c函数编译成dll，然后在python导入ctypes，通过ctypes提供的方法直接调c函数。

3维的fdtd不加边界条件很不鲁棒，数组下标稍微有所变化，结果就相差十万八千里。
第四章三维FDTD第一题fd3d_4.1.c，偶极天线，未加边界条件：
from pylab import *
from enthought.mayavi import mlab

nsteps = 40
IE = 40; ic = IE/2
JE = 40; jc = JE/2
KE = 40; kc = KE/2
gax = ones((IE,JE,KE)); dx = zeros((IE,JE,KE)); ex = zeros((IE,JE,KE)); hx = zeros((IE,JE,KE))
gay = ones((IE,JE,KE)); dy = zeros((IE,JE,KE)); ey = zeros((IE,JE,KE)); hy = zeros((IE,JE,KE))
gaz = ones((IE,JE,KE)); dz = zeros((IE,JE,KE)); ez = zeros((IE,JE,KE)); hz = zeros((IE,JE,KE))
ddx = 0.01; dt = ddx/6e8; epsz = 8.8e-12

for k in range(11,30):
    gaz[ic][jc][k] = 0.0
gaz[ic][jc][kc] = 1.0

t0 = 20.0; spread = 6.0; T = 0

for n in range(1,nsteps+1):
    T += 1;
    for k in range(1,KE):
        for j in range(1,JE):
            for i in range(1,IE):
                dx[j][k] += 0.5*(hz[j][k]-hz[j-1][k]-hy[j][k]+hy[j][k-1])
                dy[j][k] += 0.5*(hx[j][k]-hx[j][k-1]-hz[j][k]+hz[i-1][j][k])
                dz[j][k] += 0.5*(hy[j][k]-hy[i-1][j][k]-hx[j][k]+hx[j-1][k])
    pulse = exp(-0.5*((t0-T)/spread)**2);
    dz[ic][jc][kc] = pulse;
    for k in range(1,KE-1):
        for j in range(1,JE-1):
            for i in range(1,IE-1):
                ex[j][k]  = gax[j][k]*dx[j][k]
                ey[j][k]  = gay[j][k]*dy[j][k]
                ez[j][k]  = gaz[j][k]*dz[j][k]
    for k in range(1,KE-1):
        for j in range(1,JE-1):
            for i in range(1,IE):
                hx[j][k] += 0.5*(ey[j][k+1]-ey[j][k]-ez[j+1][k]+ez[j][k])
    for k in range(1,KE-1):
        for j in range(1,JE):
            for i in range(1,IE-1):
                hy[j][k] += 0.5*(ez[i+1][j][k]-ez[j][k]-ex[j][k+1]+ex[j][k])
    for k in range(1,KE):
        for j in range(1,JE-1):
            for i in range(1,IE-1):
                hz[j][k] += 0.5*(ex[j+1][k]-ex[j][k]-ey[i+1][j][k]+ey[j][k])
    print T
x, y = np.ogrid[0:IE,0:JE]
ez0 = zeros((IE,JE))
for i in range(IE):
    for j in range(JE):
        ez0[j] = ez[j][kc]
pl = mlab.surf(x, y, ez0, warp_scale="auto")

图片:fd3d_4.1.JPG

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

25楼发表于: 2010-08-08 09:54:00

谢谢楼主，我已经将我的FDTD程序利用python做了一个简单界面，但是有一个问题想问一下：

mathplotlib中的画图的函数中，最终显示是用show()，但是该函数仅仅能够显示出一副图片，只有关闭该显示界面后，下面的代码才能够继续进行，如果想显示波传播那样的问题，就不很多方便。
我目前是利用savefig将需要显示的图片保存下来，但是这样的方式不直接。

有没有什么好的办法能够让python将显示的图像显示在自己定义的界面上的某个区域中？而且还可以实时刷新。
谢谢

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

26楼发表于: 2010-08-08 12:54:24

可以试一下mayavi。比matplotlib要简单好用，图形也更美观，可参考上面的np.ogrid，mlab.surf，可以实时刷新，也可以嵌入。附件是一篇介绍。

附件:

full_text.pdf (767 K) 下载次数:5

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

27楼发表于: 2010-08-08 14:32:11

谢谢楼主，我试试看看，

还有一个问题想问一下楼主，就是C++与python之间的数组如何传递，我目前的方式是做一个函数get_list(int i,int j)将数组的结果通过下标i和j返回

由于不是每个步骤都需要返回结果，只有在需要显示的时间步上才导出，因此现在还不觉得有多慢

我希望有一种方式将numpy中的数组（一维或二维数组）从C++中返回给py或者由py指派给C++，不知道怎样解决

其它的数据结构的对应问题已经解决了，我用的是boost.python

离线samuelp

UID ：63797

注册：2010-07-23
登录：2010-10-19
发帖：28
等级：仿真二级

28楼发表于: 2010-08-09 10:33:28

numpy可以通过ctypes调c，如果调c++可以用一个c函数转一下，这样c/c++的代码是可以原封不动的。
http://hyry.dip.jp/pydoc/ctypes_numpy.html

boost不了解，可能要改c++的代码吧：
01.#include <boost/python.hpp>
02.#include "boost/python/extract.hpp"
03.#include "boost/python/numeric.hpp"
04.#include <iostream>
05.
06.using namespace boost::python;
07.
08.// Functions to demonstrate extraction
09.void setArray(boost::python::numeric::array data) {
10.
11.    // Access a built-in type (an array)
12.    boost::python::numeric::array a = data;
13.    // Need to <extract> array elements because their type is unknown
14.    std::cout << "First array item: " << extract<int>(a[0]) << std::endl;
15.}
16.
17.// Expose classes and methods to Python
18.BOOST_PYTHON_MODULE(TestNumPy) {
19.    boost::python::numeric::array::set_module_and_type("numpy", "ndarray");
20.
21.    def("setArray", &setArray);
22.
23.}

离线caocheng82

UID ：10116

注册：2008-03-28
登录：2025-05-26
发帖：697
等级：积极交流六级

29楼发表于: 2010-08-09 19:59:13

你这种方式好像只能setArray
我有时想要的是getArray

发帖回复