![[PukiWiki]](image/c60.png "[PukiWiki]")
[2018/08 started by Maruoka and Maeda]
MEEPに関すること分からないことをmaru@flex.phys.tohoku.ac.jpに送っていただければ答える努力をします。居室（842号室）に来てもらってもOKです。
時間領域差分法（FDTD法）による電磁気場のシミュレーションを行うソフトです。
Simpetusによってオープンソースで開発されており、フリーで使用できます。
Linux環境でPythonを用いてMEEPを使えるようすることを目標とします。
MEEPはLinuxかMacでしか動作しないため、Windowsを使用する場合はWSL(Windows Subsystem for Linux)をインストールしてください。
WSLのインストール方法はこちら。
WindowsでLinuxが使えるようになったらLinuxと同じ手順で環境が構築できます。
まずはAnacondaを用いてPythonの環境を構築します。(参考)
AnacondaのサイトにてPython3.x系のインストーラーをダウンロードしてください。
インストーラーはBourne Shell(sh)のシェルスクリプトで書かれているので以下のように実行してください。
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_64.sh
次にPyMeepのインストールをします。(参考)
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これでPyMeepのインストールは終わりました。
以下のコマンドでPyMeepの仮想環境をアクティベートしてください。
$ source activate mp
"-bash: activate: No such file or directory" とエラーが出るときは
export PATH=~/anaconda3/bin:$PATH
を実行してください。 flexやtube（齋藤グループのワークステーションの名前）ではPyMeepの仮想環境をアクティベートしようとすると
#!/bin/sh _CONDA_ROOT="/home/students/maru/anaconda3" \. "$_CONDA_ROOT/etc/profile.d/conda.sh" || return $? _conda_activate "$@"
というエラーが出ます。
これはログイン先のシェルがtcshであるにも関わらず、shのシェルスクリプトをsourceしているのが原因です。
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してから実行してください。
実行はhogehoge.pyの部分を任意のファイルに置き換えてください。
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるのでこのソースコードで試してみてください。
PyMeepの仮想環境を無効化するには以下のコマンドを実行してください。
(mp) $ source deactivate
並列処理をサポートしたPyMeepをインストールするには上記の作業をすべて行ったうえで次のコマンドを実行してください。
$ conda create -n pmp -c chogan -c conda-forge pymeep-parallel
$ source activate pmp
にて並列処理をサポートしたPyMeepの仮想環境をアクティベートしたのち、
$ mpirun -np 4 python <script_name>.py
を実行してください。
オプションの" -np 4 "は4プロセッサで並列処理することを指定しています。
並列に計算するプロセッサの数は適切な値に変更してください。
独立なプログラムを並列に計算することで計算時間を削減することもできます。（Embarassingly parallelといいます）。
meep <file-name>.ctl &
とすると、バックグランドでプログラムを実行できます。
例えば、次のようなシェルスクリプトを使用しています。
入射角を0から85度まで5度ずつ変えた18個のプログラムを同時に走らせています。
for((i=0;i<90;i+=5))
do
meep incident-angle=${i} Main.ctl &
done
MEEPで用いられるMaxwell方程式は
など通常のMaxwell方程式と異なります。
MEEPで用いられるMaxwell方程式は次のような式変形の最後の形です。
（ドキュメント参考）
通常の以下のMaxwell方程式について、
次のように変形します。
以下のようにx',y',z',t',ω',E',H',σ'を定めます。
aは計算する系で使用する単位長さです。
このとき、Maxwell方程式は次のようになります。
この形のMaxwell方程式をMEEPで使っています。
このとき、真空の誘電率、透磁率はどちらも1です。
また、単位長さはaです。
単位時間に真空中の光は単位長さaだけ進みます。
周波数fの真空中の光の波長は1/fです。
x',y',z',t',ω',σ'はすべて無次元です。
MEEPのシミュレーション結果はHDF5ファイルによって出力されます。
PythonではHDF5ファイルを編集するためのh5pyというライブラリが提供されています。
以下のコマンドでインストールできます。
apt-get install python-h5py
bash on windowsではpip install h5pyによってインストールできますが、次のエラーにより実行できません。
ImportError: No module named h5py
一つ目のコマンドでインストールするようにしてください。
次のコードで左回りの円偏光になります。
amplitudeによって、ソースの位相を変化することができます。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
)
)
次のコードで右回りの円偏光になります。
amplitudeによって、ソースの位相を変化することができます。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
)
)
x方向に直線偏光した平面波とy方向に直線偏光した平面波を 中心(0,0,8)、大きさ50×50×0の平面から、xy平面について45度方向((x,y)=(1,1)方向)に直線偏光した平面波を出力するコード。
(set! sources
(list
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ex)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
(make source
(src (make continuous-src (frequency (/ 1 7))))
(component Ey)
(center 0 0 8)
(size 50 50 0)
(amplitude (exp (* 0+1i 1.570796327)))
)
)
)
斜め入射のシミュレーションの方法は3つあります。
(set! geometry-lattice (make lattice (size 10 no-size 10)))
(set-param! resolution 10)
(define (pw-amp k x0) (lambda (x)
(exp (* 0+1i (vector3-dot k (vector3+ x x0))))))
(define-param theta-deg 48)
(define theta-rad (deg->rad theta-deg))
(define epsilon1 2.25)
(define epsilon2 1.25)
(define-param fcen 0.8) ; pulse center frequency
(define-param df 0.02) ; turn-on bandwidth
(define-param kdir (vector3 (sin theta-rad) 0 (cos theta-rad))) ; direction of k (length is irrelevant)
(define k (vector3-scale (* 2 pi fcen (sqrt epsilon1))
(unit-vector3 kdir))) ; k with correct length
(set! pml-layers (list (make pml (thickness 2)(direction Z))))
(set! k-point (vector3* (* fcen (sqrt epsilon1) ) (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! geometry
(list
(make block
(center 0 0 -2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon1)))
)
(make block
(center 0 0 2.5)
(size 10 infinity 5)
(material (make dielectric(epsilon epsilon2)))
)
)
)
(set! sources
(list
(make source
(src (make continuous-src (frequency fcen) (fwidth df)))
(component Ez) (center 0 0 -3) (size 10 0 0)
(amp-func (pw-amp k (vector3 0 0 -3))))))
(define-param T 40) ; run time
(run-until T
(to-appended "ez" (at-every 0.1 output-efield-z))
(to-appended "ex" (at-every 0.1 output-efield-x))
)
Ezのプロット
Exのプロット
Ezのプロット (dfを0.1に変更しています）
光源が一つしかない場合のEzのプロット（dfを0.1に変更しています）
; sz : computational domain size in z direction
; theta-rad : incident angle in radian
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param no-obstacle? false)
(define-param theta-deg 0)
(define theta-rad (deg->rad theta-deg))
(define fcen 3)
(define df 0.01)
(define epsilon1 1)
(define epsilon2 9)
(define nfreq 1)
(set! dimensions (if (= theta-rad 0) 1 3))
(set! k-point (vector3* fcen (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex)
(center 0 0 -1))))
(if (not no-obstacle?)
(set! geometry
(list
(make block
(center 0 0 -0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilon1)))))
(make block
(center 0 0 0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilon2))))))))
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 0.5)))
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -0.5))))
(run-sources+ 50
(stop-when-fields-decayed 50 Ex (vector3 0 0 0) 1e-9))
(display-fluxes trans refl)
; sz : computational domain size in z direction
; theta-rad : incident angle in radian
; fcen : pulse center frequency (500 nm)
; df : pulse width in frequenhcy
; nfreq : number of frequencies at which to compute flux
; unit length: 100 micro meter
(set-param! resolution 200)
(define-param sz 3)
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param no-obstacle? false)
(define-param theta-deg 0)
(define theta-rad (deg->rad theta-deg))
(define fcen 3)
(define df 0.01)
(define epsilonA 1)
(define epsilonB 9)
(define nfreq 1)
(set! k-point (vector3* fcen (vector3 (sin theta-rad) 0 (cos theta-rad))))
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ey)
(center 0 0 -1)
)))
(if (not no-obstacle?)
(set! geometry
(list
(make block
(center 0 0 -0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilonA))))
)
(make block
(center 0 0 0.5)
(size infinity infinity 1)
(material (make medium(index (sqrt epsilonB))))
)
)
)
)
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 0.5)
)
)
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -0.5)
)
)
)
(run-sources+ 50
(stop-when-fields-decayed 50 Ey (vector3 0 0 0) 1e-9)
;(at-beginning output-epsilon)
;(to-appended "ex" (at-every 1 output-efield-x))
)
(display-fluxes trans refl)
アンドープのグラフェンを誘電体で挟むと特殊な条件で吸収率が2%から50%に上昇します。
(S. A. Nulli, M. S. Ukhtary, R. Saito, Appl. Phys. Lett. 112, 073101-1-4, (2018))。
MEEPで再現できます。
; a : unit length in the program
; e : elementary charge
; hbar : reduced planck constant
; epsilon_vaccum : permittivity of vaccum
; sigma : the conductivity of graphene in the defenition of 2D material
; sigma_d : the conductivity of graphene in the defenition of MEEP
; d : the thickness of graphene in the SI unit
; dg : the thickness of graphene in this program unit
; fcen : pulse center frequency (500 nm)
; df : pulse width in frequency
; nfreq : number of frequencies at which to compute flux
; isEmpty : Whether graphene exist or not.If false, we calculate in vacuum
; sz : the size of computatinal domain in Z direction in this program unit
; Z_0 : the impedance of the vaccum
; s : the number of repetitions of the AB layers before mirroring.
; alpha : When the ratio of index of dielectric A and B is alpha, the absorbance is up to 50%
; da : thickness of dielectric A
; db : thickness of dielectric B
(define a 1e-9)
(define e (* 1.60217662 1e-19))
(define c 299792458)
(define hbar (* 1.05457180013 1e-34))
(define epsilon_vacuum (* 8.854187817 1e-12))
(define pi (* 4 (atan 1)))
(define d (* 0.3 1e-9))
(define dg 0.3)
(define sigma (/ (* e e) (* 4 hbar)))
(define sigma_d (/ (* sigma a) (* c epsilon_vacuum d)))
(define fcen 0.002)
(define df 0.0001)
(define nfreq 100)
(define s 2)
(define Z_0 376.730313461)
(define alpha (expt
(/ 2 (* Z_0 sigma))
(/ (* 2 s))))
(define epsilon_a (* alpha alpha))
(define epsilon_b (* 1 1))
(define da (/ (* 4 fcen (sqrt epsilon_a) )))
(define db (/ (* 4 fcen (sqrt epsilon_b) )))
(set-param! resolution 10)
(define sz (* 2 400))
(set! geometry-lattice (make lattice(size no-size no-size sz)))
(define-param isEmpty "DO NOT SET VALUE DIRECTLY")
(set! dimensions 1)
(if (string? isEmpty) ((print "Set the parameter \"isEmpty\" true or false from command line.\nDo not use the parameter isEmpty as it is.\n")(exit)))
(if (not isEmpty)
(set! geometry
(list
(make block
(center 0 0 0)
(size infinity infinity dg)
(material (make medium (epsilon 1)(D-conductivity sigma_d)))
)
(make block
(center 0 0 (+ (/ dg 2) (/ db 2)))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (+ (/ dg 2) db (/ da 2)))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (+ (/ dg 2) db da (/ db 2)))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (+ (/ dg 2) db da db (/ da 2)))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (- (+ (/ dg 2) (/ db 2))) )
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db (/ da 2))))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db da (/ db 2))))
(size infinity infinity db)
(material (make dielectric(epsilon epsilon_b)))
)
(make block
(center 0 0 (- (+ (/ dg 2) db da db (/ da 2))))
(size infinity infinity da)
(material (make dielectric(epsilon epsilon_a)))
)
)
)
)
(set! sources (list
(make source
(src (make gaussian-src (frequency fcen) (fwidth df)))
(component Ex)
(center 0 0 350)
)))
(set! pml-layers (list (make pml (thickness 1))))
(define trans
(add-flux fcen df nfreq
(make flux-region
(center 0 0 -330)
)
)
)
(define refl
(add-flux fcen df nfreq
(make flux-region
(center 0 0 330)
)
)
)
(run-sources+ 150000
(stop-when-fields-decayed 50 Ex (vector3 0 0 -330)
1e-3
)
(at-beginning output-epsilon)
(to-appended "ex" (at-every 10 output-efield-x))
)
(display-fluxes trans refl)
全反射における電場ベクトルの分布
シミュレーションのためのコード
(define-param dx 1); cell size in x axis (define-param dy 2); cell size along y axis (define-param fp 1.5) (define-param tpml 0.1); PML thickness (define-param deg 48) (define-param theta (* pi 0.88889)); incident angle of the plane wave (set! theta (* deg (/ pi 180))) (set! resolution 40);define the resolution for simulation domain ;define the lattice (set! geometry-lattice (make lattice (size (+ dx (* 2 tpml)) (+ dy (* 2 tpml)) no-size))) ;define PML thickness (set! pml-layers (list (make pml (thickness tpml)))) ;define the blocks (set! geometry (list (make block (center 0 (* 0.25 dy))(size infinity (+ tpml (* 0.5 dy)) infinity) (material (make dielectric (epsilon 1)))) (make block (center 0 (* -0.25 dy))(size infinity (+ tpml (* 0.5 dy)) infinity) (material (make dielectric (epsilon 2.25)))) )) (set! pml-layers (list (make pml (thickness tpml)))) ;(+ tpml (* 0.5 dy)) ;define the wave vector (define kx (* fp (sqrt 2.25) (sin theta))) ;give the amplitude function (define (f_amp p) (exp (* 0+2i pi kx (vector3-x p)))) (set! k-point (vector3 kx 0 0)) ;define the simulation domain to be complex field (set! force-complex-fields? true) ;define the Gaussian beam (set! sources (list (make source (src (make continuous-src (frequency fp))) (component Ex) (component Ey) (center 0 (* dy -0.50)) (size (+ dx (* 2 tpml)) 0 ) (amp-func f_amp)))) (set! pml-layers (list (make pml (thickness tpml) (direction Y)))) ;extract the data to .png file (run-until 40 (at-beginning output-epsilon) ;(at-end output-efield-y) (to-appended "ey"(at-every 0.05 (in-volume (volume (center 0 0 0)(size dx dy 0))output-efield-y))) (to-appended "ex"(at-every 0.05 (in-volume (volume (center 0 0 0)(size dx dy 0))output-efield-x))) )
電場ベクトルをプロットするには、HDF5ファイルから直接プロットするコードを自分で書く必要があります。
pythonのmatplotlibのquiver関数を用いるといいと思います。
MEEPには電荷分布を計算する機能がないため、自作する必要があります。
ガウスの法則を差分化すればよいです。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレーションについて例を示します。
電荷分布(+電場ベクトル)のプロット
誘電体界面で誘電分極が起きていることと波源で電荷の移動があることが分かります。
PMLをある方向にだけ設置するにはオプションdirectionを設定すればよいです。
例えば、Y方向だけに設置するには次のようにすればよいです。
pml_layers = [mp.PML(thickness=1.0, direction=mp.Y)]
PMLが設置されていない方向は固定端の全反射が起きます。
真ん中にwaveguideを置いて、Y方向にだけPMLを設置したときの電場は次のようになります。
X方向ではPMLがないため全反射が起きて定在波が生じています。
コードは次の通りです。
import meep as mp
cell = mp.Vector3(16,8,0)
geometry = [mp.Block(mp.Vector3(1e20,1,1e20),
center=mp.Vector3(0,0),
material=mp.Medium(epsilon=9))]
sources = [mp.Source(mp.ContinuousSource(frequency=0.1),
component=mp.Ez,
center=mp.Vector3(-7,0))]
pml_layers = [mp.PML(thickness=1.0, direction=mp.Y)]
resolution = 10
sim = mp.Simulation(cell_size=cell,
boundary_layers=pml_layers,
geometry=geometry,
sources=sources,
resolution=resolution)
dT = 0.2
sim.run(
mp.at_beginning(mp.output_epsilon),
mp.to_appended("ez",
mp.at_every(dT, mp.output_efield_z)),
until=100
)
物性値のライブラリが公開されています。
materials-library.scmをダウンロードしてください。
ライブラリを使用するには、プログラムの最初に
(include "/path/to/materials-library.scm")
を追加してライブラリをロードしてください。
"/path/to/materials-library.scm"は自分がmaterials-library.scmを設置した場所のアドレスに書き換えてください。
丸岡がダウンロードしたmaterials.scmを使用する場合は次のパスを使用してください。
(include "/home/students/maru/meep/scheme/materials.scm")
ただし最新のmaterials.scmではない可能性があることに注意してください。
このライブラリはtube61, 62, 63でしか動きません（flex, tube60では動きません）。
詳しい使い方は https://meep.readthedocs.io/en/latest/Materials/ を参照してください。
materials.scmはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、特別に指定する必要があります。
例えば、0.1µmを単位長さにしたいときは
(set! um-scale (* um-scale 0.1))
としてください。
物性値のライブラリmaterials.pyをロードするのに特別な操作をする必要はありません。
materials.pyはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、materials.pyを直接書き換える必要があります。
materials.pyのum_scaleが宣言されている行以降でum_scaleを目的の単位長さに書き換えてください。
例えば、0.1µmを単位長さにしたいときは
um_scale=um_scale*0.1
としてください。
materials.pyの場所は次のようにして調べてください。
(pmp) maru@tube60:~$ python Python 3.6.6 | packaged by conda-forge | (default, Jul 26 2018, 09:53:17) [GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux Type "help", "copyright", "credits" or "license" for more information. >>> import meep Using MPI version 3.1, 1 processes >>> meep.__file__ '/home/students/maru/anaconda3/envs/pmp/lib/python3.6/site-packages/meep/__init__.py'
この場合、materials.pyのアドレスは/home/students/maru/anaconda3/envs/pmp/lib/python3.6/site-packages/meep/materials.pyです。
このmaterials.pyの書き換えは(pmp)の仮想環境でだけ有効であることに注意してください。
((mp)の仮想環境でmaterials.pyの単位長さを変える場合は(mp)の仮想環境のmaterials.pyを書き換えてください）
ドルーデ伝導で表される次のような比誘電率を用いる方法を説明する。
（ε_∞は周波数が十分高いところでの比誘電率（定数）。f_pはプラズマ周波数）
MEEPのドルーデ伝導はMEEPの単位系のf_p', γ', f'で表されているので、SI単位系のf_p, γ, fとの関係を求める必要がある。
この式から、SI単位系のf_p, γ とMEEPの単位系のf_p', γ'は次の関係にある。
上の関係式から求めたf_p', γ'とε_∞について、
user_defined_material = mp.Medium(epsilon=ε_∞, E_susceptibilities=[mp.DrudeSusceptibility(frequency=f_p', gamma=γ', sigma=1)])
とすると、そのf_p', γ', ε_∞で表されるドルーデ伝導を持つ物質が定義できます。
この物質を使用する際は、geometryのオプション materials = user_defined_material を設定する必要があります。
PythonによるHDF5ファイルの取り扱いを説明します。
import h5py #以降のプログラムのpukiwikiでの掲載ではこの行は省略する
f = h5py.File("file-name.h5", "r") #rとすると読み込み専用で開く
print(f)
とすると、<HDF5 file "file-name" (mode r)>と表示されます。
print(list(f).keys())
とすると、 h5fileのkeyが表示される。
例えば、['ex']と表示されたします。
Ex = f['ex'][:]
とすると、キー'ex'に紐づいた配列をExに取り出せます。
注意として、f['ex'][i]のようにして、何回もアクセスしないでください。
このようなアクセスは時間がかかります。
print(Ex) #print(h5file['ex'].value)も同様の動作をする。
とすると、Exの配列の中身が表示されます。
要素が多すぎるときは・・・と表示されて省略されます。
Exが二次元配列でEx[i][j]と参照するとしたとき、iを第一インデックス、jを第二インデックスと呼ぶことにします。
第一インデックスの長さ(iの動く長さは)
print(len(Ex))
で分かります。
第二インデックスの長さは
print(len(Ex[0]))
で分かります。
以下のスクリプトを.bash_profileに貼り付けて使用すると、h5ファイルをmp4ファイルに自動で変換できて便利です。
function h5mp4(){
mkdir $1
cp $1.h5 $1
cd $1
h5topng $1.h5 -Z -c dkbluered -R -t 0:$2
convert $1*png $1.gif
ffmpeg -i $1.gif -qscale:v 1 $1.mp4
}
例えば時間の区切り回数が1000回あるようなhogehoge.h5を動画にしたいときは
h5mp4 hogehoge 999
とするとhogehoge.gifとhogehoge.mp4とhogehoge-<数字>.pngが自動で生成されます。
時間の区切り回数が知りたいときは
h5ls hogehoge.h5
としてください。
ez Dataset{100, 200, 1000/Inf}
と表示される数字のうち、1000が区切り回数です。
import h5py
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from PIL import Image
import moviepy.editor as edit
from mpl_toolkits.axes_grid1 import make_axes_locatable
plt.rcParams["font.size"] = 18
#----------input----------------
surfaceExFilename ="After-ex_+0.5.h5"
surfaceEyFilename = "After-ey_+0.5.h5"
surfacePxFilename = "Pre-ex_+0.5.h5"
surfacePyFilename = "Pre-ey_+0.5.h5"
savename = "current-line-withPreE-takingMiddleCurrent.png"
startTime = 0
endTime = 1000
waveLength = 6
timeResolution = 10
xid=16 #array[xid][yid][:]の電場をプロット
yid=46
#--------------------------------
def plotf(arr,label_):
time=np.array([])
amplitude=np.array([])
for i in range(startTime, endTime):
x=i * 1./(waveLength * timeResolution)
y=arr[xid][yid][i]
time=np.append(time,[x])
amplitude=np.append(amplitude,[y])
plt.plot(time,amplitude,label=label_)
plt.legend(bbox_to_anchor=(1.05, 1), loc="upper left")
h5file_A_x = h5py.File(surfaceExFilename, "r")
h5file_A_y = h5py.File(surfaceEyFilename, "r")
h5file_P_x = h5py.File(surfacePxFilename, "r")
h5file_P_y = h5py.File(surfacePyFilename, "r")
Ax = h5file_A_x['ex'][:]
Ay = h5file_A_y['ey'][:]
Px = h5file_P_x['ex'][:]
Py = h5file_P_y['ey'][:]
Ellは2次元配列でxyに応じて楕円率が入っていると想定
ヒストグラムの正規化はググると「normed=True」が出てくるがうまくいかなかった
stackoverflowが役に立った
日本語解説ブログ
ecのキーワードでヒストグラムの縦線の色を指定できる(参考)
Ell_ = Ell.reshape(-1,)#plt.histは1次元配列を対象とするのでEllを1次元化
weights_ = np.ones(len(Ell_))/float(len(Ell_))#ヒストグラムを正規化するための処理
fig = plt.figure()
ax = fig.add_subplot(111)
#binsはヒストグラムの縦棒の数
ax.hist(Ell_, weights=weights_, bins=21, range=(-1,0,1.0),ec='black')
ax.set_xlim(-1.0,1.0)
ax.set_ylim(0.0,0.5)
fig.savefig("histgram.png")
array.reshape(-1)の挙動
>>> a = np.linspace(1,8,num=8)
>>>
>>> a
array([1., 2., 3., 4., 5., 6., 7., 8.])
>>> b = a.reshape((2,4))
>>> b
array([[1., 2., 3., 4.],
[5., 6., 7., 8.]])
>>> c = a.reshape((2,2,2))
>>> c
array([[[1., 2.],
[3., 4.]],
[[5., 6.],
[7., 8.]]])
>>> d = c.reshape(-1,)
>>> d
array([1., 2., 3., 4., 5., 6., 7., 8.])
代わりにweightsオプションを指定すればこの問題は解決できる。
これはdataの1つの値の重みを指定するパラメータで, デフォルトでは実質的に weights = np.ones(len(data)) を指定していることになる。
いま総和を1にしたいのだから, ひとつのデータの重みは1/len(data)です。
ndarrayの場合は配列 / scalarが可能。
>>> x = [1,2,3] >>> x [1, 2, 3] >>> x /2 Traceback (most recent call last): File "<stdin>", line 1, in <module> TypeError: unsupported operand type(s) for /: 'list' and 'int' >>> y = np.array([1,2,3]) >>> y / 2 array([0.5, 1. , 1.5])
楕円率のヒストグラムの作成の際の hist()のキーワードに
range=(-1,1)
を入れなないと空気を読んで，最大値と最小値の間をbinsの数だけ分割するようです。
ax.set_xlim(-1., 0, 1., 0)
だと，-1から0までにしかプロットされないので
ax.set_xlim(-1.,1.)
に修正。
1つのファイル内で複数回シミュレーションをするときは、各シミュレーションを終了するたびに(reset-meep) を呼び出さなければなりません。
フラックスは座標軸の正の方向を正に取っているためです（ポインティングベクトルで定義されている）。正の方向に入射波が進行し、負の方向に反射波が進行するとき、入射波だけの時よりフラックスは小さくなります。
理由は不明です。
resolutionはset! geometry-latticeの前に置かないとデフォルト値10になってしまうことがあります。
(define lambda 5)とするとunbounded variable vとしてエラーが出る。Schemeではlambdaを用いて関数を定義するからです（予約語になっています）。
fluxを3次元で設定していることが原因です。fluxは2次元の面または1次元の点で設定しなければなりません。
放っておいても問題ありませんが、.loginに以下を記述すると、非推奨(deprecated)な機能を使用しているというエラーを握りつぶせます。
setenv GUILE_WARN_DEPRECATED no
「unbound variable: 」とエラーが出てるにも関わらず、unbound variableの名前が表示されないときは、全角空白文字がプログラムのどこかに含まれています。
全角空白を削除すればよいです。
meep hoehoge.ctl var=1
とすると、var=1は無視されます。
meep var=1 hogehoge.ctl
としてください。
e1,e2,e3[vector3]により、方向を定めることができます。
例えば、xy平面内でx軸からΘだけ傾けて置きたいならば、
(make block (size a b c) (e1 (vector3 (cos theta) (sin theta) 0)) (e2 (vector3 (- (sin theta)) (cos theta) 0) (e3 (vector3 0 0 1)) )
X, Y, Z方向をそれぞれ厚さa, b, cにしたいなら、
(set! pml-layers (list (make pml (thickness a)(direction X)))) (set! pml-layers (list (make pml (thickness b)(direction Y)))) (set! pml-layers (list (make pml (thickness c)(direction Z))))
としてください。
負の透過率、反射率、吸収率は電磁波がPMLの外に発散しきっていないのにシミュレーションを中断したことが原因のことが多いです。
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x y z) 1e-3))
とすると、T秒間シミュレーションをとりあえず回した後に、(x,y,z)におけるdt秒間のExの最大値が、全ての時間を通しての最大値と比較して1e-3倍以下になるまで回してくれます。
しかしこの関数は落とし穴があって、非常に長時間共振するような構造の場合は、まだ、電磁波の発振が続いているのにシミュレーションを中断する場合があります。
電磁波が発散しきっているかどうかの判断は、h5topngなどで電磁波を画像化してを直接確認するしかありません。
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくなさそうです。
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg05147.html
マニュアルのトップページにSchemeがobsoleteでPythonに変えろと書いてあったことと、メーリングリストで（バグはまだ多いですが）Pythonに機能が追加されていることもあって、先週からPythonに移行のために動作検証しているところです。
PMLに垂直に入射するときは半波長。 以下、MEEP作成者のPMLに関するドキュメン(http://math.mit.edu/~stevenj/18.369/pml.pdf)からの引用
With PML, however, the constant factor is very good to start with, so experience shows that a simple quadratic or cubic turn-on of the PML absorption usually produces negligible refections for a PML layer of only half a wavelength or thinner [1,21]. (Increasing the resolution also increases the resolution alos increase the effectiveness of the PML, because it approaches the exact wave equation.)
斜め入射のときはPMLの吸収率が落ちるので、厚みを増す必要があります。
反射率は入射角をXとすると、ある定数δを用いて e^{-2*k*d*cos(X)*δ} と書けます。
dはPMLの厚さ、kは波数。
詳細はhttps://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=477075を参照してください。
高専ではmaterials.pyが使用できるが、tube,flexでは使用できていません(2019/01/20:使用できるようになりました)。
次のようなエラーが出ます。
tube60:~$ python Python 3.6.7 | packaged by conda-forge | (default, Nov 21 2018, 03:09:43) [GCC 7.3.0] on linux Type "help", "copyright", "credits" or "license" for more information. >>> import meep as mp >>> from meep.materials import Au Traceback (most recent call last): File "<stdin>", line 1, in <module> ModuleNotFoundError: No module named 'meep.materials'
高専の環境では~/anaconda3/envs/mp/lib/python3.6/site-packages/meep/にmaterials.pyが存在するが、flex,tubeには存在せずこのことが原因の可能性があります。
flex,tubeの~/anaconda3/pkgs/pymeep-1.5-py36_nomkl_2が高専では~/anaconda3/pkgs/pymeep-1.6-py36_nomkl_2になっており、バージョンを1.5から1.6にすることで使用できるようになる可能性があります。
ネット上からmaterias.pyをダウンロードしてホームディレクトリに置いたのちに実行すると、実行はできるものの次のようなワーニングが出てくる。
>>> from materials import Au /home/guest/maeda/.conda/envs/mp/lib/python3.6/site-packages/meep/geom.py:175: RuntimeWarning: Epsilon < 1 may require adjusting the Courant parameter. See the 'Numerical Stability' entry under the 'Materials' section of the documentation warnings.warn(eps_warning, RuntimeWarning)
高専ではワーニングなしに実行できます。
conda list
とすると、
pymeep 1.5
となり、PyMeepのバージョンが1.5であることが分かりました。
オプションをつけてPyMeepを1.5から1.6にupdateすると実行できるようになりました。
(mp) maru@flex:~$ conda update -c chogan -c conda-forge pymeep
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/envs/mp
added / updated specs:
- pymeep
The following packages will be downloaded:
package | build
---------------------------|-----------------
libctl-4.1.4 | 1 97 KB chogan
libgdsii-0.2.dev | 0 800 KB chogan
python-3.6.6 | hd21baee_1003 29.0 MB conda-forge
openssl-1.0.2p | h14c3975_1002 3.1 MB conda-forge
mpb-1.7.0 | nomkl_1 74 KB chogan
pymeep-1.7.0 | py36_nomkl_1 3.3 MB chogan
openblas-0.3.5 | h9ac9557_1000 15.8 MB conda-forge
nomkl-3.0 | 0 48 KB
ca-certificates-2018.11.29 | ha4d7672_0 143 KB conda-forge
certifi-2018.11.29 | py36_1000 145 KB conda-forge
------------------------------------------------------------
Total: 52.5 MB
The following NEW packages will be INSTALLED:
libgdsii: 0.2.dev-0 chogan
nomkl: 3.0-0
The following packages will be UPDATED:
ca-certificates: 2018.03.07-0 --> 2018.11.29-ha4d7672_0 conda-forge
certifi: 2018.11.29-py36_0 --> 2018.11.29-py36_1000 conda-forge
libctl: 4.1.2-0 chogan --> 4.1.4-1 chogan
mpb: 1.6.2-4 chogan --> 1.7.0-nomkl_1 chogan [nomkl]
openblas: 0.2.20-8 conda-forge --> 0.3.5-h9ac9557_1000 conda-forge
pymeep: 1.5-py36_nomkl_2 chogan [nomkl] --> 1.7.0-py36_nomkl_1 chogan [nomkl]
The following packages will be DOWNGRADED:
openssl: 1.1.1a-h7b6447c_0 --> 1.0.2p-h14c3975_1002 conda-forge
python: 3.6.8-h0371630_0 --> 3.6.6-hd21baee_1003 conda-forge
Proceed ([y]/n)? y
Downloading and Extracting Packages
libctl-4.1.4 | 97 KB | ########################################################################################### | 100%
libgdsii-0.2.dev | 800 KB | ########################################################################################### | 100%
python-3.6.6 | 29.0 MB | ########################################################################################### | 100%
openssl-1.0.2p | 3.1 MB | ########################################################################################### | 100%
mpb-1.7.0 | 74 KB | ########################################################################################### | 100%
pymeep-1.7.0 | 3.3 MB | ########################################################################################### | 100%
openblas-0.3.5 | 15.8 MB | ########################################################################################### | 100%
nomkl-3.0 | 48 KB | ########################################################################################### | 100%
ca-certificates-2018 | 143 KB | ########################################################################################### | 100%
certifi-2018.11.29 | 145 KB | ########################################################################################### | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
pymeep-parallelにおいてmaterials.pyを使用したいときは、次のように1.7.0-py36_nomklh39e3cac_1にアップデートしてください。
(pmp) maru@tube60:~$ conda update -c chogan -c conda-forge pymeep-parallel
Solving environment: done
## Package Plan ##
environment location: /home/students/maru/anaconda3/envs/pmp
added / updated specs:
- pymeep-parallel
The following packages will be downloaded:
package | build
---------------------------|-----------------
libgcc-ng-7.3.0 | hdf63c60_0 6.1 MB conda-forge
pymeep-parallel-1.7.0 |py36_nomklh39e3cac_1 3.3 MB chogan
------------------------------------------------------------
Total: 9.4 MB
The following NEW packages will be INSTALLED:
libgdsii: 0.2.dev-0 chogan
nomkl: 3.0-0
The following packages will be UPDATED:
ca-certificates: 2018.8.24-ha4d7672_0 conda-forge --> 2018.11.29-ha4d7672_0 conda-forge
certifi: 2018.8.24-py36_1 conda-forge --> 2018.11.29-py36_1000 conda-forge
libctl: 4.1.2-0 chogan --> 4.1.4-1 chogan
libgcc-ng: 7.2.0-hdf63c60_3 conda-forge --> 7.3.0-hdf63c60_0 conda-forge
mpb: 1.6.2-4 chogan --> 1.7.0-nomkl_1 chogan [nomkl]
openssl: 1.0.2p-h470a237_0 conda-forge --> 1.0.2p-h14c3975_1002 conda-forge
pymeep-parallel: 1.5-py36_nomklh39e3cac_4 chogan [nomkl] --> 1.7.0-py36_nomklh39e3cac_1 chogan [nomkl]
Proceed ([y]/n)? y
Downloading and Extracting Packages
libgcc-ng-7.3.0 | 6.1 MB |
################################################################################################################################################################################################## | 100%
pymeep-parallel-1.7. | 3.3 MB |
##################################################################################################################################################### ############################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done
以下のようなエラーが出た場合は、計算を行っているディレクトリおよびファイルの書き込み・変更を許可してください。
Working in 3D dimensions.
Computational cell is 20 x 20 x 20 with resolution 10
block, center = (0,0,0)
size (1,3,1)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (1,1,1)
time for set_epsilon = 11.7491 s
-----------
creating output file "./hogehoge.h5"...
on time step 1 (time=0.05), 14.4708 s/step
HDF5-DIAG: Error detected in HDF5 (1.10.2) MPI-process 0:
#000: H5F.c line 445 in H5Fcreate(): unable to create file
major: File accessibilty
minor: Unable to open file
meep: error on line 497 of h5file.cpp: error opening HDF5 output file
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1
flexではGLIBCのバージョンが2.14よりも古いためPyMeepを実行できません。
tube70に移動してPyMeepを実行してください。
(pmp) maeda@tube60:/abinitio2/maeda/pymeep/Rectangle/drude/Omote$ mpirun -np 5 python after-Rec-drude-Au.py 45
Working in 3D dimensions.
Computational cell is 20 x 20 x 20 with resolution 10
block, center = (0,0,0)
size (1,3,1)
axes (1,0,0), (0,1,0), (0,0,1)
dielectric constant epsilon diagonal = (1,1,1)
time for set_epsilon = 5.05642 s
lorentzian susceptibility: frequency=3.47141, gamma=2.01155
lorentzian susceptibility: frequency=2.39466, gamma=0.701702
lorentzian susceptibility: frequency=0.66944, gamma=0.278261
lorentzian susceptibility: frequency=0.33472, gamma=0.19438
drude susceptibility: frequency=1e-10, gamma=0.0427474
-----------
creating output file "./after-Rec-drude-Au-ex-surface-45deg.h5"...
creating output file "./after-Rec-drude-Au-ey-surface-45deg.h5"...
creating output file "./after-Rec-drude-Au-ez-inner-45deg.h5"...
creating output file "./after-Rec-drude-Au-ez-outer-45deg.h5"...
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=92, req_array=0x9adaad0, status_array=0xc9e3670) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(623).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=85, req_array=0xa328890, status_array=0xd21c7f0) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(658).......: fail failed
MPID_nem_handle_pkt(1439).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
Fatal error in PMPI_Waitall: Other MPI error, error stack:
PMPI_Waitall(405)...............: MPI_Waitall(count=85, req_array=0x9c2b8f0, status_array=0xcc31ed0) failed
MPIR_Waitall_impl(221)..........: fail failed
PMPIDI_CH3I_Progress(658).......: fail failed
MPID_nem_handle_pkt(1439).......: fail failed
pkt_RTS_handler(317)............: fail failed
do_cts(662).....................: fail failed
MPID_nem_lmt_dcp_start_recv(302): fail failed
dcp_recv(165)...................: Internal MPI error! Cannot read from remote process
Two workarounds have been identified for this issue:
1) Enable ptrace for non-root users with:
echo 0 | sudo tee /proc/sys/kernel/yama/ptrace_scope
2) Or, use:
I_MPI_SHM_LMT=shm
解決策 エラーメッセージの2つめの方法で解決しました。
export I_MPI_SHM_LMT=shm
を実行するとエラーは出なくなりました。
波源の周波数を十分低くします。
さらに波源の立ち上げを緩やかにします。
波源の立ち上げを緩やかにするにはContinuousSrcのwidthというパラメータを変更すればよいです。
ブロッホの周期的境界条件の下でシミュレーションした結果をh5topngすると
h5topng error: data can have at most two dimensions (try specifying a slice)
というエラーが出ます。
これはブロッホの周期的境界条件の下でシミュレーションすると計算結果に実部と虚部の両方が含まれるためです。
オプション"-d"によって出力するデータを指定してください。
出力できるデータはh5lsで調べることができます。
h5ls output.h5
ez.i Dataset {101, 101, 133/Inf}
ez.r Dataset {101, 101, 133/Inf}
となっているときは、電場のz成分の虚部ez.iと実部ez.rのどちらかを選択できます。
電場のz成分の実部ez.rを出力するときはh5opngにオプション"-d ez.r"を設定してください。
共鳴周波数付近では、波源の立ち上げの電磁場に対する応答のノイズがなかなかPMLの外に放出されません。
ノイズを防ぐには波源の立ち上げを緩やかにしてください。例えばPyMeepを使用している場合はmp.ContinuousSourceのオプションwidthを設定してください。
git clone https://github.com/stevengj/meep ./autogen.sh make
とするとコードの完全体が得られます。
libctlというSchemeの自作ライブラリが読み込まれています。
Autotools で処理されるべきファイルに .in という拡張子が使われます。
つまり、その段階ではまだ Scheme のコードのテンプレートであり、処理されることで Scheme のコードとして妥当なものになります。
ビルドプロセスによりますが、通常は configure が済んだ段階では処理済みのもの (この場合は meep.scm) が出来ていると思います。
もちろん make が済んだ時点なら確実に出来ているでしょう。
(define-class A B (define-property C D E))とすると、Aという名前のクラスが定義されます。
Bはオーバーライドするクラスで、BのプロパティをAは引き継ぎます。
Bをno-parentにすると、オーバーライドは起きません。
Aに対してCというプロパティが設定され、その初期値はDとなります。
またCの型はEとなります。
以下にmeep.scm.inに含まれるクラスの定義の例を示します。
(define-class susceptibility no-parent (define-property sigma-offdiag (vector3 0 0 0) 'vector3) (define-property sigma-diag no-default 'vector3)) (define-class lorentzian-susceptibility susceptibility (define-property frequency no-default 'number)
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMはOpenFDTD DLサイト、OpenMOM DLサイトからダウンロードすることができます。
使い方はダウンロードサイトならびに研究室が所蔵しているRFワールド No.39 2017年 08月号 に記載されています。
書籍は丸岡が借りています（2018/6/12)。
