[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のインストール方法はこちら。
WSL上ではLinuxと同じ手順で環境が構築できます。
まずはAnacondaを用いてPythonの実行環境を構築します。(参考)
AnacondaのサイトにてPython3.x系のインストーラーをダウンロードしてください。
インストーラーはBourne Shell(sh)のシェルスクリプトで書かれているので以下のように実行してください。
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_64.sh
次のエラーが出るときは、~/anaconda3/binにパスを通してください。
bash: conda: command not found
次のコマンドを用いると、~/anaconda3/binにパスが通ります。
setenv PATH $PATH\:<path to anaconda3/bin>
相対パスではなく絶対パスを使用することに注意してください。
以下のコマンドを実行してください。
$ 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に変更してからPyMeepの仮想環境をアクティベートしてください。
実行はhogehoge.pyの部分を任意のファイルに置き換えてください。
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるのでこのソースコードで試してみてください。
PyMeepの仮想環境を無効化するには以下のコマンドを実行してください。
(mp) $ source deactivate
並列処理の可能なPyMeepをインストールするにはAnacondaのインストールを行ったうえで次のコマンドを実行してください。
$ 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 &
とすると、バックグランドでプログラムを実行できます。
例えば、次のようなシェルスクリプトを使用しています。
for((i=0;i<90;i+=5)) do meep incident-angle=${i} Main.ctl & done
入射角を0から85度まで5度ずつ変えた18個のプログラムを同時に走らせています。
MEEPで用いられるMaxwell方程式は
など通常のMaxwell方程式と異なります。
MEEPで用いられるMaxwell方程式は次のような式変形の最後の形です。
(ドキュメント参考)
通常の以下の形のMaxwell方程式を考えます。
次のように変形します。
aはシミュレーションで使用する単位長さです。
cは光速です。
以下のようにx',y',z',t',ω',E',H',σ'を定めます。
このとき、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
"apt-get install python-h5py"によってインストールするようにしてください。
次のコードで左回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx方向に比べてπ/2だけ進めています。
(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))) ) ) )
次のコードで右回りの円偏光になります。
sourcesクラスのamplitudeによって、y方向の電磁場の位相をx方向に比べてπ/2だけ遅らせています。
(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))) ) ) )
次のコードで45度方向((x,y)=(1,1)方向)に直線偏光した平面波になります。
x方向に直線偏光した平面波とy方向に直線偏光した平面波を重ねることで45度方向の直線偏光にしています。
(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) ) ) )
斜め入射のシミュレーションの方法は3つあります。
ブロッホの周期的境界条件を用いて電磁波の斜め入射をシミュレーションすることができます。
次の図のように画面に平行な方向にxz平面を取ります。
右方向にx軸、縦方向にz軸、画面垂直方向にy軸を取ります。
z<0には比誘電率2.25、z>0には比誘電率1.25の誘電体を敷き詰めます。
誘電体界面に入射角48度で電磁場が入射したときのシミュレーションを考えます。
全反射が入射角arcsin(sqrt(1.25/2.25))=48.2度で起きるので入射角48度は全反射直前の角度です。
z軸に垂直な境界にはPMLを、x軸に垂直な境界にはブロッホの周期的境界条件を用います。
計算領域の大きさはx,y,z方向それぞれについて10,0,10です。
電磁場はy方向依存性を持たないので、y方向の計算領域の大きさは0となっています。
入射角をΘとすると、波数ベクトルkはk=k_0(sinΘ,0,cosΘ)です。
波源はx軸に平行な線として定義するので、位相のx依存性としてexp(i(x k_0 sinΘ))を持つ必要があります。
Schemeではsourceクラスのオプションであるamp-funcを設定することで位相のx座標依存性を設定できます。
また、計算する領域内(今回の場合は大きさ10×0×10の長方形)を単位胞として見た場合の基本並進ベクトルをRとすると、
x軸と垂直な計算領域の境界の点rについてE(r)exp(i k・R) = E(r+R)が成り立ちます。
磁場についても同様です。
点rと点r+Rはともに境界上の点になっています。E(r)は点rでの電場ベクトルです。
このような境界条件はブロッホの周期的境界条件を用いることで実現できます。
Schemeの場合はk-point変数を設定することでブロッホの周期的境界条件を設定することができます。
(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のプロット
屈折角が90度近いので、透過光のExの大きさは小さくなっています。
入射角45度の場合に限ってのみ、2つの波原で斜め入射の平面波が作れます。
計算領域がxz平面上に乗っている場合を例として考えます。
境界条件は四方全てPMLとします。
左側の境界と下側の境界にそれぞれ線状の波原を設置します。
波原の位相のx,z依存性は波数ベクトルk=(k_0sinθ, 0, k_0cosθ)であることを考えるとそれぞれexp(ixk_0sinθ),exp(ixk_0cosθ)です。
Schemeではsourceクラスのオプションであるamp-funcを設定することで位相のx,z座標依存性を設定できます。
Ezのプロット (dfを0.1に変更しています)
光源が一つしかない場合のEzのプロット(dfを0.1に変更しています)
このプロットから一様な振幅の大きさを持つ斜め入射の電磁場を作るには光源が二つ必要なことが分かります。
位相の揃った点光源を一直線に斜めに並べると、光源の並んだ方向と垂直な方向に進行する平面波が作れます。
直線状には点光源の強度がガウス分布になるように並べています。
光源の強度をすべて同じにすると干渉によって強度分布が不均一になります。
半値幅を調節することでビームの横幅を変更できます。
(define-param pi 3.14159265) (define-param dpml 3.0) (define-param len (+ 20 dpml)) (define-param wid (+ 20 dpml)) (set! geometry-lattice (make lattice (size len wid no-size))) (define-param beam-waist 2.5) ; beam sigma (gaussian beam width) (define-param rotation-angle (* (/ 22.5 360) 2 pi)) ; Degrees if you like them (define-param source-points 60) ; should be a big number (define-param source-size (* 4 beam-waist)) ; should be bigger than beam-waist (define-param src_list (list )) (do (( r_0 (/ source-size -2) (+ r_0 (/ source-size (- source-points 1))) )) ((> r_0 (/ source-size 2))) ;for r_0 = -source-size/2 : source- size/source-points : source-size/2 (set! src_list (append src_list (list (make source (src (make gaussian-src (wavelength 1) (width 3) )) (amplitude (exp (- 0 (/ (* r_0 r_0) (* 2 beam-waist beam-waist))))) (component Ez) (center (* r_0 (sin rotation-angle)) (* r_0 (cos rotation-angle))) )) )) ) (set! sources src_list) (set! pml-layers (list (make pml (thickness dpml) ))) (set! resolution 10) (use-output-directory) (run-until (* 2 len) (to-appended "ez" (at-every 0.1 output-efield-z)) )
plot of Ez
TM表面プラズモンを励起するには、入射光の波数の界面に平行な成分および周波数が、励起されるTM表面プラズモンの波長および周波数と一致する必要があります。
この励起状態を満たすためにOtto geometryという構造が使われることがあります。
import meep as mp import math import sys resolution = 100 # 系の大きさ size_x,size_y,size_z sx = 10 sy = 10 sz = 0 cell = mp.Vector3(sx, sy, 0) # PML層の厚さ(MEEP単位)depth_pml dpml = 1 pml_layers = [mp.PML(thickness=dpml)] # 何秒(MEEP単位)おきにデータを取るか dT = 0.2 sources = [] for i in range(100): x_=-4 + 4 * i / 100. y_=0 + 4 * i / 100. r=math.sqrt(abs(x_+2)**2+abs(y_-2)**2) coe=math.exp(-r*r) sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1), center=mp.Vector3(x=x_, y=y_, z=0), component=mp.Ex, amplitude=coe )) sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1), center=mp.Vector3(x=x_, y=y_, z=0), component=mp.Ey, amplitude=coe )) material_=mp.Medium( epsilon=1, E_susceptibilities = [mp.DrudeSusceptibility(frequency=math.sqrt(3), gamma=0, sigma=1)], ) geometry = [ mp.Block(size=mp.Vector3(x=20, y=5, z=mp.inf), center=mp.Vector3(x=0, y=2.6, z=0), material=mp.Medium(epsilon=4) ), mp.Block(size=mp.Vector3(x=20, y=5, z=mp.inf), center=mp.Vector3(x=0, y=-2.6, z=0), material=material_ ) ] sim = mp.Simulation(cell_size=cell, boundary_layers=pml_layers, geometry=geometry, sources=sources, resolution=resolution ) sim.run( mp.to_appended("ey", mp.at_every(dT, mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx, y=sy, z=0)), mp.output_efield_y))), until=20 )
屈折率nはn=√mu√eps(mu,epsはそれぞれ比透磁率と比誘電率)と書けるので、muとepsがともに負のときは屈折率nが負になります。
このとき波の位相速度は通常の位相速度と進行方向が逆向きになります。
シミュレーションの際に、負の屈折率を持つ媒質をPMLと重ねないように注意してください。
PMLと負の屈折率を持つ媒質を重ねると、PML中で電磁波は吸収されずに増大します。
import meep as mp import math import sys resolution = 20 # 系の大きさ size_x,size_y,size_z sx = 20 sy = 15 sz = 0 cell = mp.Vector3(sx, sy, 0) # PML層の厚さ(MEEP単位)depth_pml dpml = 1 pml_layers = [mp.PML(thickness=dpml)] # 何秒(MEEP単位)おきにデータを取るか dT = 0.5 sources = [] for i in range(100): x_=-4 + 4 * i / 100. y_=0 + 4 * i / 100. r=math.sqrt(abs(x_+2)**2+abs(y_-2)**2) coe=math.exp(-r*r) sources.append(mp.Source(src=mp.ContinuousSource(wavelength=1), center=mp.Vector3(x=x_, y=y_, z=0), component=mp.Ez, amplitude=coe )) material_=mp.Medium( epsilon=1, E_susceptibilities = [mp.DrudeSusceptibility(frequency=2, gamma=0.00145, sigma=1)], mu=1, H_susceptibilities = [mp.DrudeSusceptibility(frequency=2, gamma=0.00145, sigma=1)], ) geometry = [ mp.Block(size=mp.Vector3(x=15, y=2, z=mp.inf), center=mp.Vector3(x=0, y=-3, z=0), material=material_ ) ] sim = mp.Simulation(cell_size=cell, boundary_layers=pml_layers, geometry=geometry, sources=sources, resolution=resolution ) sim.run( mp.to_appended("ez", mp.at_every(dT, mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx, y=sy, z=0)), mp.output_efield_z))), until=30 )
ある物体に対する反射率、吸収率、透過率を計算する方法は2つあります。
1つはガウシアンソースを使ったときの、ポインティングベクトルを物体のあるときとないとき両方について、物体の前後に相当する位置でポインティングベクトルを積分して計算する方法です。
例えば透過率ならば透過光のポインティングベクトルの積分/入射光のポインティングベクトルの積分を求めれば良いです。
もうひとつの方法はContinuous Sourceを使ったときの透過光、反射光の振幅の大きさと入射光の振幅の大きさを比較することです。
ここでは一つ目の方法を紹介します。
; 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))。
(J. Blumberg, M. S. Ukhtary, R. Saito, Phys. Rev. Applied 10, 064015-1-14, (2018))。
MEEPで再現できます。
グラフェンを二種類の誘電体A, Bで挟んで例えばABABGBABAのような並びにします。
ここでGはグラフェンの層です。
それぞれの誘電体の厚みは垂直入射したときの光路長が1/4波長になるようにします。
誘電体Aの層の数を2nとします。
このとき、誘電体AとBの誘電率の比がε_A/ε_B={2/(Zσ)}^{1/n}のときグラフェン近傍の電場が増強されて光の吸収率が50%に達します。
ここで、Zは真空のインピーダンス、σはグラフェンの電気伝導率です。
シミュレーションではグラフェンの伝導率σはσ=πe^2/(2h)としています(eは電荷素量、hはプランク定数)。
平面波が誘電体とグラフェンの層に対して垂直入射するとき電磁場の値は層と平行な方向に対して依存性を持ちません。
従ってシミュレーションは一次元で行うことができます(マクスウェル方程式に含まれる微分のうち、層と平行な方向の微分を0に置換するとマクスウェル方程式に含まれる位置の変数が一つになります)。
; 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)
電場分布は次のようになります(パルス幅や誘電率、入射光の波長を上のコードから変更しています。nは屈折率です)。
全反射における電場ベクトルの分布
シミュレーションのためのコード
(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 )
MEEPでは後に置いた物体が先に置いた物体を完全に上書きするため、S字や螺旋を小数のパーツで作ることは難しいです。
しかし、たくさんのディスクを配置することでS字や螺旋を形成することができます。
import meep as mp import math resolution = 10 # 系の大きさ size_x,size_y,size_z sx = 20 sy = 20 sz = 20 cell = mp.Vector3(sx, sy, sz) # PML層の厚さ(MEEP単位)depth_pml dpml = 3 pml_layers = [mp.PML(dpml)] # 何秒(MEEP単位)おきにデータを取るか dT = 0.1 sources = [] geometry_list=[] for i in range(1000): theta=i/1000.*2*math.pi*3 radius_=i/200. geometry_list.append(mp.Cylinder(radius=0.5,height=0.35,center=mp.Vector3(x=radius_*math.cos(theta),y=radius_*math.sin(theta),z=0),material=mp.Medium(epsilon=12))) geometry=geometry_list sim = mp.Simulation(cell_size=cell, boundary_layers=pml_layers, geometry=geometry, sources=sources, resolution=resolution) sim.run( mp.at_beginning(mp.in_volume(mp.Volume(center=mp.Vector3(0, 0, 0), size=mp.Vector3(x=sx - 2 * dpml, y=sy - 2 * dpml, z=0)), mp.output_epsilon)), until=1 )
z=0平面における誘電率のプロットは次の通りです(PMLは除いた範囲をプロットしています)。
物性値のライブラリが公開されています。
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ではない可能性があることに注意してください。
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'][:]
MEEPでは近傍界のフラックスデータから遠方界のフラックスデータを構成することができます。
近傍界を基に遠方界のフラックスを得る利点は小さい系で計算させるためシミュレーションが速く終わることです。
MEEPに遠方界を計算させるための条件
アンテナは計算領域の中心に配置してください。
フラックス面はPMLより内側に設置してください。
パルス波源を使う場合は、電磁波がPMLの外に完全に発散するまでシミュレーションを続ける必要があることに注意してください。
以下に点光源を疑似的ダイポールアンテナと見立てた場合の遠方界のプロットを目的としたシミュレーションのコードを示します。
該当ドキュメント
注意<<<
2019/02/14現在,ドキュメントに書かれているコードではプロットがうまくいかないので,作成者のgitに古いバージョンが挙がっているのでそちらを使った方が安全。
このシミュレーションでは円で遠方界を計算します。
計算領域はxy平面の二次元です。
nearfield_boxに遠方場に必要なデータが入っているので,
get_farfieldで引数にnearfield_boxを引数にとり,遠方界の座標をVector3()で与えます。
プロット用のコード
実行コマンド群(上からシミュレーション、データの整形、プロット)
$ python antenna-radiation.py | tee source_Jz_farfields.out $ grep farfield: source_Jz_farfields.out | cut -d , -f2- > source_Jz_farfields.dat $ python stable-antenna-radiation-plot.py
出力結果
この結果は"antenna-radiation.py"のsrc_cmptにmp.Ex,mp.Ey,mp.Ezのどれかを入れれば
上から順にこれらの画像を作成することができます。
それぞれxy平面をz軸方向から見ている図に相当しています。
ダイポールの振動方向と垂直な方向への放射が最も強いことが分かります。
放射の強さの最大値が1となるように規格化しています。
1つのファイル内で複数回シミュレーションをするときは、各シミュレーションを終了するたびに(reset-meep) を呼び出さなければなりません。
(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方向のPMLをそれぞれ厚さ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の外に発散しきっていないのにシミュレーションを中断したことが原因のことが多いです。
特にstop-when-fields-decayed関数を使用しているときは注意してください。
(run-sources+ T (stop-when-fields-decayed dt Ex (vector3 x y z) 1e-3))
とすると、T秒間シミュレーションした後に、(x,y,z)におけるdt秒間毎のExの最大値が、全ての時間を通しての最大値と比較して10の-3乗倍以下になるまでシミュレーションを続けてくれます。
しかしこの関数は、長時間共鳴するような構造の場合は、まだ、電磁波の発振が続いているのにシミュレーションを中断する場合があります。
電磁波が発散しきっているかどうかの判断は、h5topngなどで電磁波を画像化して直接確認するしかありません。
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくなさそうです。
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg05147.html
マニュアルのトップページにSchemeがobsoleteでPythonに変えろと書いてあったことと、メーリングリストで(バグはまだ多いですが)Pythonに機能が追加されていることもあって、先週からPythonに移行のために動作検証しているところです。
電磁波がPMLに垂直入射するときは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は波数。
詳細はUPMLの論文を参照してください。
PyMeepのバージョンが1.5だと次のエラーによりmaterials.pyが実行できません。
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'
PyMeepのバージョンは
conda list
から確認できます。
次のようにオプションをつけてPyMeepを1.5から1.6にアップデートするとmaterials.pyが実行できるようになります。
(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を実行できません。
ImportError: /lib/x86_64-linux-gnu/libc.so.6: version `GLIBC_2.14' not found (required by /home/students/maru/anaconda3/envs/mp/lib/python3.6/site-packages/meep/_meep.so)
tube70に移動してPyMeepを実行してください。
仙台高専からflexにsshしてmpiを使用すると次のエラーが出ます。
(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
を実行するとmpiが実行できるようになります。
波源の周波数を十分低くします。
さらに波源の立ち上げを緩やかにします。
波源の立ち上げを緩やかにするにはContinuousSrcのwidthというパラメータを変更してください。
Schemeを使用している場合、ブロッホの周期的境界条件の下でシミュレーションした結果を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を設定してください。
PMLと物体は1/4波長、離してください。
PMLと物体が近接していると、PMLに対する電磁波の入射角が大きくなるためPMLの吸収率が減少してしまいます。
The proof of the divergence free nature of the Yee grid is derived by in Taflove's book on FDTD.By "satisfying natural boundary conditions" I mean that if you have a material interface meandering through your grid, you will not have to do anything special to simulate it accurately. If it were not for the Yee grid, you would observe crazy things like 800% reflection from a surface for some polarizations.The only remaining problem really is the staircasing of a curved surface due to the Cartesian grid. There are fixes for this like dielectric smoothing, unstructured grid, etc.
https://www.youtube.com/watch?v=hv5lIx4u8mY&t=566s
fluxを取る面と波源を重ねないでください。
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMはOpenFDTD DLサイト、OpenMOM DLサイトからダウンロードすることができます。
使い方はダウンロードサイトならびに研究室が所蔵しているRFワールド No.39 2017年 08月号 に記載されています。
書籍は丸岡が借りています(2018/6/12)。