※ifortで通ったので修正しました。
I/Oはpythonで書くけど、科学技術計算の部分はfortranに任して高速化してしまおうという話です。以下のページをほぼそのまま使用しました。
PythonからFortranのサブルーチンを呼ぶ。 - いぐにさんにっき
windows 環境で fortran と C++ を連携させる(intel fortranの場合) | Moonmile Solutions Blog
ちなみにnumbaによる高速化も試しましたが、今回のプログラムでは全然早くなりませんでした。どうもリスト内包標記が駄目っぽいです。それだとpythonの良さがなくなってしまうので使えないですね。
速度比較
- モデルは例のsectionQ2Dモデルです。
収束計算はガウス・ザイデルです。(これをADIとかに変えれば良いという話もありますが。)
とりあえず速度の比較です。
- python:231.958 sec
- gfortran:1.222 sec
- ifort:0.706 sec
衝撃の速さです。なお、イテレーションの回数は同じです。
つまずいた点
計算条件等はpythonで設定して、配列のままfortranに渡すことを試みました。環境依存は結構ありそうです。今回は以下の環境です。
- windows7 64bit
- python3.5
- gfortran 4.8.1
- ifort 2013
つまずいた点は以下です。
ソースコード
python版
def sectionq2d3(): un = np.array([[ini_velocity for k in range(nz)] for j in range(ny)]) unnew = np.zeros((ny,nz)) def NumericalSimulation(): for j in range(0,ny): for k in range(0,nz): if j == 0 : # Wall Boundary apy = (nyut2[k]+nyut2[k])*dz**2 ayp = dz**2*nyut2[k]*un[j+1,k] aym = 0 elif j == ny-1 : # Wall Boundary apy = (nyut2[k]+nyut2[k])*dz**2 ayp = 0 aym = dz**2*nyut2[k]*un[j-1,k] else: # Normal apy = (nyut2[k]+nyut2[k])*dz**2 ayp = dz**2*nyut2[k]*un[j+1,k] aym = dz**2*nyut2[k]*un[j-1,k] if k == 0 : # Bed Boundary apz = (0.5*nyut2[k+1]+nyut2[k]+0.5*nyut2[k])*dy**2 azp = dy**2*0.5*(nyut2[k]+nyut2[k+1])*un[j,k+1] azm = 0 elif k == nz-1 : # WaterSurface Boundary: Slip Condition apz = 0.5*(nyut2[k]+nyut2[k-1])*dy**2 azp = 0 azm = dy**2*0.5*(nyut2[k]+nyut2[k-1])*un[j,k-1] else: # Normal apz = (0.5*nyut2[k+1]+nyut2[k]+0.5*nyut2[k-1])*dy**2 azp = dy**2*0.5*(nyut2[k]+nyut2[k+1])*un[j,k+1] azm = dy**2*0.5*(nyut2[k]+nyut2[k-1])*un[j,k-1] un[j,k] = (9.8*ib*dy**2*dz**2+ayp+aym+azp+azm)/(apy+apz) for num in range(1,10000): NumericalSimulation() if np.max(np.abs(unnew-un)) > 0.000001 : unnew[:,:] = un else : print(num) break return un
fortran版
python側でdllを呼ぶ。
from ctypes import * def sectionq2d4(fortType): if fortType == 'gfort' : f = np.ctypeslib.load_library("sq2dFortranGfort.dll",".") f.sectionq2d4_.argtypes = [ np.ctypeslib.ndpointer(dtype=np.float64), np.ctypeslib.ndpointer(dtype=np.float64), POINTER(c_int64),POINTER(c_int64), POINTER(c_double),POINTER(c_double),POINTER(c_double) ] f.sectionq2d4_.restype = c_void_p elif fortType == 'ifort' : f = np.ctypeslib.load_library("sq2dFortranIfort.dll",".") f.SECTIONQ2D4.argtypes = [ np.ctypeslib.ndpointer(dtype=np.float64), np.ctypeslib.ndpointer(dtype=np.float64), POINTER(c_int64),POINTER(c_int64), POINTER(c_double),POINTER(c_double),POINTER(c_double) ] f.SECTIONQ2D4.restype = c_void_p un = np.array([[ini_velocity for k in range(nz)] for j in range(ny)],dtype=np.float64) nyut21 = np.array(nyut2,dtype=np.float64) ny2 = byref(c_int64(ny)) nz2 = byref(c_int64(nz)) dy2 = byref(c_double(dy)) dz2 = byref(c_double(dz)) ib2 = byref(c_double(ib)) if fortType == 'gfort' : f.sectionq2d4_(un,nyut21,ny2,nz2,dy2,dz2,ib2) elif fortType == 'ifort' : f.SECTIONQ2D4(un,nyut21,ny2,nz2,dy2,dz2,ib2) return un
fortranのコードです。
subroutine sectionq2d4(un, nyut2, ny, nz, dy, dz, ib) ! ifortの場合、以下の1文が必要。gfortranの場合は不要。 !DEC$ ATTRIBUTES DLLEXPORT :: sectionq2d4 implicit none integer(8), INTENT(IN) :: ny, nz real(8), INTENT(IN) :: dy, dz, ib real(8), INTENT(IN) :: nyut2(0: nz-1) real(8), INTENT(INOUT) :: un(0: nz-1, 0: ny-1) integer(8) :: num real(8) :: unnew(0: nz-1, 0: ny-1) real(8) :: tmp(0: nz-1, 0: ny-1) integer(8) :: k, j real(8) :: apy,ayp,aym,apz,azp,azm unnew = 0.0 do num = 1, 1000000 call NumericalSimulation if (maxval(abs(unnew-un)) > 0.000001) then unnew = un else print *,num exit end if end do CONTAINS subroutine NumericalSimulation() do j = 0, ny-1 do k = 0, nz-1 if (j == 0) then !Wall Boundary apy = (nyut2(k)+nyut2(k))*dz**2 ayp = dz**2*nyut2(k)*un(k,j+1) aym = 0 else if (j == ny-1) then ! Wall Boundary apy = (nyut2(k)+nyut2(k))*dz**2 ayp = 0 aym = dz**2*nyut2(k)*un(k,j-1) else ! Normal apy = (nyut2(k)+nyut2(k))*dz**2 ayp = dz**2*nyut2(k)*un(k,j+1) aym = dz**2*nyut2(k)*un(k,j-1) end if if (k == 0) then ! Bed Boundary apz = (0.5*nyut2(k+1)+nyut2(k)+0.5*nyut2(k))*dy**2 azp = dy**2*0.5*(nyut2(k)+nyut2(k+1))*un(k+1,j) azm = 0 else if (k == nz-1) then ! WaterSurface Boundary: Slip Condition apz = 0.5*(nyut2(k)+nyut2(k-1))*dy**2 azp = 0 azm = dy**2*0.5*(nyut2(k)+nyut2(k-1))*un(k-1,j) else ! Normal apz = (0.5*nyut2(k+1)+nyut2(k)+0.5*nyut2(k-1))*dy**2 azp = dy**2*0.5*(nyut2(k)+nyut2(k+1))*un(k+1,j) azm = dy**2*0.5*(nyut2(k)+nyut2(k-1))*un(k-1,j) end if un(k,j) = (9.8*ib*dy**2*dz**2+ayp+aym+azp+azm)/(apy+apz) end do end do end subroutine end subroutine
gfortran -shared -o sq2dFortran.dll sq2dFortran.f90 ifort /nologo /dll /iface:cvf -o sq2dFortranIfort.dll sq2dFortran.f90
相当使えそうです。やっぱりfortranはサルでも速いコードがかけるので便利ですね。