- 追加された行はこの色です。
- 削除された行はこの色です。
[2018/08 started by Maruoka and Maeda]
*目次 [#t1296280]
#contents
*MEEP [#g2bb328a]
**MEEPとは [#fbfd5425]
時間領域差分法（FDTD法）による電磁気場のシミュレーションを行うソフトです。
オープンソースで開発されており、フリーで使用できます。
Simpetusによってオープンソースで開発されており、フリーで使用できます。
**PythonでMEEP [#fc543fa9]
Linux環境でPythonを用いてMEEPを使えるようすることを目標とします。~
MEEPはLinuxかMacでしか動作しないため、Windowsを使用する場合はWSL(Windows Subsystem for Linux)をインストールしてください。
やり方は[[こちら:http://www.atmarkit.co.jp/ait/articles/1608/08/news039.html]]。~
WSLのインストール方法は[[こちら:http://www.atmarkit.co.jp/ait/articles/1608/08/news039.html]]。~
WindowsでLinuxが使えるようになったらLinuxと同じ手順で環境が構築できます。~
***インストール [#p4cd7dd9]
まずはAnacondaを用いてPythonの環境を構築します。([[参考:https://qiita.com/t2y/items/2a3eb58103e85d8064b6]])~
[[Anacondaのサイト:https://www.anaconda.com/download/]]にてPython3.x系のインストーラーをダウンロードしてください。~
インストーラーはBourne Shell(sh)のシェルスクリプトで書かれているので以下のように実行してください。
$ bash [ダウンロードフォルダ]/Anaconda3-5.2.0-Linux-x86_64.sh
次にPyMeepのインストールをします。([[参考:https://meep.readthedocs.io/en/latest/Installation/]])~
以下のコマンドを実行してください。
$ conda create -n mp -c chogan -c conda-forge pymeep
これで必要なセットアップは終わりました。
これでPyMeepのインストールは終わりました。
***実行 [#padb0b5f]
以下のコマンドでPyMeepの仮想環境をアクティベートしてください。
$ source activate mp
"-bash: activate: No such file or directory" とエラーが出たときは
"-bash: activate: No such file or directory" とエラーが出るときは
export PATH=~/anaconda3/bin:$PATH
を実行してください。
flexやtubeでは
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であるにも関わらず、shのシェルスクリプトをsourceしているのが原因です。
tcshからはtcshスクリプトしかsource出来ません。
bash
と打って、ログインシェルをbashに変更してから実行してください。
実行はhogehoge.pyの部分を任意のファイルに置き換えてください。
(mp) $ python hogehoge.py
サンプルソースはGithub等で入手できるので[[このソースコード:https://github.com/stevengj/meep/blob/master/python/examples/wvg-src.py]]で試してみてください。~
無効化するには以下のコマンドを実行してください。~
PyMeepの仮想環境を無効化するには以下のコマンドを実行してください。~
(mp) $ source deactivate
**MEEPの並列化 [#k1f4cae9]
***インストール([https://meep.readthedocs.io/en/latest/Installation/#mpi 参考]) [#v7b0b7df]
並列処理をサポートしたPyMeepをインストールするには上記の作業をすべて行ったうえで次のコマンドを実行してください。
$ conda create -n pmp -c chogan -c conda-forge pymeep-parallel
***実行 [#j3e26633]
$ source activate pmp
にて並列処理をサポートしたPyMeepの仮想環境をアクティベートしたのち、
$ mpirun -np 4 python <script_name>.py
を実行してください。
" -np 4 "とすると4プロセッサで並列処理されます。
オプションの" -np 4 "は4プロセッサで並列処理することを指定しています。
並列に計算するプロセッサの数は適切な値に変更してください。
***別の並列化 [#e415ae6a]
独立なプログラムを並列に計算することで計算時間を削減することもできます。（Embarassingly parallelといいます）。
meep <file-name>.ctl &
とすると、バックグランドでプログラムが実行されます。
とすると、バックグランドでプログラムを実行できます。
可変パラメータを外部から渡すことで、多数のプログラムを同時に動かせます。
例えば、次のようなシェルスクリプトを使用しています。
入射角を0から85度まで5度ずつ変えたプログラムを同時に走らせています。
入射角を0から85度まで5度ずつ変えた18個のプログラムを同時に走らせています。
for((i=0;i<90;i+=5))
do
meep incident-angle=${i} Main.ctl &
done
**MEEPで用いられるMaxwell方程式 [#t352b35e]
MEEPで用いられるMaxwell方程式は
++単位長さaをUserが設定する。
++真空の誘電率ε_0が1
++真空の透磁率μ_0が1
など通常のMaxwell方程式と異なります。
MEEPで用いられるMaxwell方程式は次のような式変形の最後の形です。
（[https://meep.readthedocs.io/en/latest/Introduction/ ドキュメント参考]）
通常の以下のMaxwell方程式について、
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq1.png
以下のように変形します。
次のように変形します。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq2.png
以下のようにx',y',z',t',ω',E',H',σ'を定めます。
aは計算する系で使用する単位長さです。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq3.png
このとき、Maxwell方程式は次のようになります。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/maxwell_eq4.png
この形のMaxwell方程式を使っています。
この形のMaxwell方程式をMEEPで使っています。
このとき、真空の誘電率、透磁率はどちらも1です。
また、単位長さはaです。
単位時間に真空中の光は単位長さaだけ進みます。
周波数fの真空中の光の波長は1/fです。
x',y',z',t',ω',σ'はすべて無次元です。
**HDF5の編集 [#ic8d7598]
-python向けにh5pyというライブラリが提供されています。
MEEPのシミュレーション結果はHDF5ファイルによって出力されます。
PythonではHDF5ファイルを編集するためのh5pyというライブラリが提供されています。
以下のコマンドでインストールできます。
apt-get install python-h5py
bash on windowsではpip install h5pyによってインストールはできたものの、次のエラーが出てきて実行できませんでした。
bash on windowsではpip install h5pyによってインストールできますが、次のエラーにより実行できません。
ImportError: No module named h5py
一つ目のコマンドでインストールするようにしてください。
**左円偏光 [#u95702e4]
次のコードで左回りの円偏光になります。
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)))
)
)
)
**右円偏光 [#u95702e4]
次のコードで右回りの円偏光になります。
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)))
)
)
)
**45度直線偏光 [#u95702e4]
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)))
)
)
)
**斜め入射 [#e6b41572]
斜め入射のシミュレーションの方法は3つあります。
***全反射直前のTM波(ブロッホ周期境界条件) [#s8084013]
(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のプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ez.gif
Exのプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ex.gif
***光源を2つ並べる。 [#d178c899]
[https://github.com/stevengj/meep/blob/master/scheme/examples/pw-source.ctl コード]
Ezのプロット
Ezのプロット (dfを0.1に変更しています）
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-ez.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/ez.gif
光源が一つしかない場合のEzのプロット
光源が一つしかない場合のEzのプロット（dfを0.1に変更しています）
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-single-ez.gif
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/ez-single.gif
***光源を斜めにたくさん並べる。 [#t6b81813]
[http://kenji-ishikawa.cocolog-nifty.com/plasma/2010/01/post-ea91.html 動画とコード]
**電荷分布を見たい [#l0ec2ad8]
MEEPのbuilt-in関数には電荷分布を求める機能がないので、自作する必要があります。
MEEPには電荷分布を計算する機能がないため、自作する必要があります。
ガウスの法則を差分化すればよいです。
前述の全反射直前のTM波(ブロッホ周期境界条件)のシミュレーションについて例を示します。
電荷分布(+電場ベクトル)のプロット
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.gif
[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/charge.py コード]
**PMLをある方向にだけ設置する。（Python） [#uf2f674c]
PMLをある方向にだけ設置するにはオプションdirectionを設定すればよいです。
例えば、Y方向だけに設置するには以下のようにすればよいです。
例えば、Y方向だけに設置するには次のようにすればよいです。
pml_layers = [mp.PML(thickness=1.0, direction=mp.Y)]
PMLが設置されていない方向は固定端の全反射が起きます。
真ん中にwaveguideを置いて、Y方向にだけPMLを設置したときの電場は次のようになります。
X方向にPMLがないので全反射が起きて定在波が生じています。
X方向ではPMLがないため全反射が起きて定在波が生じています。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/pml-ydirection.gif
コードは次の通りです。
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
)
**物性値のライブラリ [#ue28f25c]
***Schemeの場合 [#z7a9ca6a]
物性値のライブラリが公開されています。
[https://github.com/NanoComp/meep/blob/master/scheme/materials.scm materials-library.scm]をダウンロードしてください。
ライブラリを使用するには、プログラムの最初に
(include "/path/to/materials-library.scm")
を追加してライブラリをロードしてください。
"/path/to/materials-library.scm"は自分が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/ を参照してください。
初期設定ではMEEPの単位長さが1μmであるものとしてmaterials.scmが設計されています。
materials.scmはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、特別に指定する必要があります。
例えば、0.1µmを単位長さにしたいときは
(set! um-scale (* um-scale 0.1))
としてください。
***Pythonの場合 [#d7e5fefa]
物性値のライブラリmaterials.pyをロードするのに特別な操作をする必要はありません。
初期設定ではMEEPの単位長さが1μmであるものとしてmaterials.pyが設計されています。
materials.pyはMEEPの単位長さが1μmであるものとして設計されています。
異なる単位長さを使用したい場合、materials.pyを直接書き換える必要があります。
um_scaleが宣言されている行以降でum_scaleを目的の単位長さに書き換えてください。
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のアドレスは/home/students/maru/anaconda3/envs/pmp/lib/python3.6/site-packages/meep/materials.pyです。
このmaterials.pyの書き換えは(pmp)の環境でだけ有効であることに注意してください。
このmaterials.pyの書き換えは(pmp)の仮想環境でだけ有効であることに注意してください。
((mp)でmaterials.pyの単位長さを変える場合は(mp)のmaterials.pyを書き換えてください）
((mp)の仮想環境でmaterials.pyの単位長さを変える場合は(mp)の仮想環境のmaterials.pyを書き換えてください）
**ユーザー定義のドルーデ伝導 [#df7ff749]
ドルーデ伝導で表される次のような比誘電率を用いる方法を説明する。
（ε_∞は周波数が十分高いところでの比誘電率（定数）。f_pはプラズマ周波数）
http://flex.phys.tohoku.ac.jp/~maru/drive-open/drude_epsilon.png
MEEPのドルーデ伝導はMEEPの単位系のf_p', γ', f'で表されているので、SI単位系のf_p, γ, fとの関係を求める必要がある。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/drude_epsilon_meepunit.png
この式から、SI単位系のf_p, γ とMEEPの単位系のf_p', γ'は次の関係にある。
http://flex.phys.tohoku.ac.jp/~maru/drive-open/relation_meepunit_vs_siunit.PNG
上の関係式から求めたf_p', γ'とε_∞について、
user_defined_material = mp.Medium(epsilon=ε_∞, E_susceptibilities=[mp.DrudeSusceptibility(frequency=f_p', gamma=γ', sigma=1)])
とすると、そのf_p', γ', ε_∞で表されるドルーデ伝導を持つ物質が定義できる。
とすると、そのf_p', γ', ε_∞で表されるドルーデ伝導を持つ物質が定義できます。
この物質を使用する際は、geometryのオプション materials = user_defined_material を設定する必要がある。
この物質を使用する際は、geometryのオプション materials = user_defined_material を設定する必要があります。
**hdf5ファイルの取り扱い [#de839d03]
pythonによるhdf5ファイルの取り扱いを説明します。
**HDF5ファイルの取り扱い [#de839d03]
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]))
で分かります。
**電磁場のh5ファイルをmp4に変換するためのスクリプト [#x4773146]
以下のスクリプトを.bash_profileに貼り付けて使用すると、h5ファイルを動画に自動で変換できて便利です。
以下のスクリプトを.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/Infが区切り回数です。
と表示される数字のうち、1000が区切り回数です。
**印加電場と実際の電場の時間変化を折れ線グラフでプロット [#q9e259c0]
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'][:]
**楕円率をヒストグラムで表したい [#g4d319c5]
Ellは2次元配列でxyに応じて楕円率が入っていると想定~
ヒストグラムの正規化はググると「normed=True」が出てくるがうまくいかなかった
[https://stackoverflow.com/questions/3866520/plotting-histograms-whose-bar-heights-sum-to-1-in-matplotlib/16399202 stackoverflow]が役に立った~
[https://jb102.blogspot.com/2017/10/22-histogram.html 日本語解説ブログ]~
ecのキーワードでヒストグラムの縦線の色を指定できる([https://stackoverflow.com/questions/42542252/cannot-get-histogram-to-show-separated-bins-with-vertical-lines 参考])
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.)
に修正。
**Q&A [#y2d911b6]
***一つのファイル内で複数回シミュレーションをすると計算結果がおかしくなる[#b9c73d94]
-1つのファイル内で複数回シミュレーションをするときは、各シミュレーションを終了するたびに(reset-meep) を呼び出さなければなりません。
1つのファイル内で複数回シミュレーションをするときは、各シミュレーションを終了するたびに(reset-meep) を呼び出さなければなりません。
***（入射波＋反射波）のフラックスが入射波のみの時より、小さくなる [#ua9ad9fe]
--フラックスは座標軸の正の方向を正に取っているためです（ポインティングベクトルで定義されている）。正の方向に入射波が進行し、負の方向に反射波が進行するとき、入射波だけの時よりフラックスは小さくなります。
フラックスは座標軸の正の方向を正に取っているためです（ポインティングベクトルで定義されている）。正の方向に入射波が進行し、負の方向に反射波が進行するとき、入射波だけの時よりフラックスは小さくなります。
***dielectricよりmediumでindexを指定した方が計算が早い（計算結果は同じなのに） [#xe0c2da7]
-理由は不明です。
理由は不明です。
***resolutionをsetしても変化しない [#g280b7dd]
-resolutionはset! geometry-latticeの前に置かないとデフォルト値10になってしまうことがあります。
resolutionはset! geometry-latticeの前に置かないとデフォルト値10になってしまうことがあります。
***変数vを使っていないのにunbounded variable vというエラーが出る。 [#ubd3a6e9]
-(define lambda 5)とするとunbounded variable vとしてエラーが出る。schemeではlambdaを用いて関数を定義するからです（予約語になっています）。
(define lambda 5)とするとunbounded variable vとしてエラーが出る。Schemeではlambdaを用いて関数を定義するからです（予約語になっています）。
***meep: Could not determine normal direction for given grid_volume. [#ta6b98ee]
-fluxを3次元で設定していることが原因です。fluxは2次元の面または1次元の点で設定しなければなりません。
fluxを3次元で設定していることが原因です。fluxは2次元の面または1次元の点で設定しなければなりません。
***実行時にGUILE_WARN_DEPRECATEDという警告が出る。 [#ma3f1cad]
-放っておいても問題ありませんが、.loginに以下を記述すると、非推奨(deprecated)な機能を使用しているというエラーを握りつぶせます。
放っておいても問題ありませんが、.loginに以下を記述すると、非推奨(deprecated)な機能を使用しているというエラーを握りつぶせます。
setenv GUILE_WARN_DEPRECATED no
***unbound variable: とエラーが出て、しかも : 以降が空白になっている。 [#i0e856b1]
「unbound variable: 」とエラーが出てるにも関わらず、unbound variableの名前が表示されないときは、全角空白文字がプログラムのどこかに含まれています。
全角空白を削除すればよいです。
***コマンドラインから変数を渡しているのに無視される [#dbd0a076]
meep hoehoge.ctl var=1
とすると、var=1は無視されます。
meep var=1 hogehoge.ctl
としてください。
***直方体を斜めに配置したい。 [#x6b589ab]
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について、厚さをそれぞれ別々にしたい [#z61ce405]
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))))
としてください。
***負の透過率、反射率、吸収率を示す [#j49576fc]
負の透過率、反射率、吸収率は電磁波が発散しきっていないのにシミュレーションを中断したことが原因のことが多いです。
(run-sources+ T
(stop-when-fields-decayed dt Ex (vector3 x y z) 1e-3))
とすると、T秒間シミュレーションをとりあえず回した後に、(x,y,z)におけるdt秒間のExの最大値が、全ての時間を通しての最大値と比較して1e-3倍以下になるまで回してくれます。
しかしこの関数は落とし穴があって、非常に長時間共振するような構造の場合は、まだ、電磁波の発振が続いているのにシミュレーションを中断する場合があります。
電磁波が発散しきっているかどうかの判断は、h5topngなどで電磁波を画像化してを直接確認するしかありません。
***伝導度の理論値を複数のローレンツの重ね合わせで表したい [#p43befa3]
下のリンクの方法で取り組んでいました。
理論値を複数のローレンツの重ね合わせで合わせるのはよくなさそうです。
https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg05147.html
***インターフェイスはpython, schemeどちらを使うべきか [#m1cc6868]
マニュアルのトップページにschemeがobsoleteでpythonに変えろと書いてあったことと、メーリングリストで（バグはまだ多いですが）pythonに機能が追加されていることもあって、先週からpythonに移行のために動作検証しているところです。
***インターフェイスはPython, Schemeどちらを使うべきか [#m1cc6868]
マニュアルのトップページにSchemeがobsoleteでPythonに変えろと書いてあったことと、メーリングリストで（バグはまだ多いですが）Pythonに機能が追加されていることもあって、先週からPythonに移行のために動作検証しているところです。
***PMLの厚さはどのくらいにすべきか？ [#l8406a8a]
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を参照してください。
***tube,flexにおいてmaterials.pyが使えない。[#m4142e79]
***tubeとflexにおいてmaterials.pyが使えない。[#m4142e79]
高専では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
***error on line 497 of h5file.cpp: error opening HDF5 output file [#s9480bc6]
以下のようなエラーが出た場合は、計算を行っているディレクトリおよびファイルの書き込み・変更を許可してください。
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
***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) [#zcd4847a]
flexのGLIBCのバージョンが古いので実行できていません。
flexではGLIBCのバージョンが2.14よりも古いためPyMeepを実行できません。
t70に移動して実行するとよいです。
tube70に移動してPyMeepを実行してください。
***Internal MPI error! Cannot read from remote process (解決) [#i054b206]
(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
を実行させてから
(pmp)$mpirun -np N python hogehoge.py
で並列計算してくれました。
***静電場はどうシミュレーションすればいいか? [#m05eea2e]
波源の周波数を十分低くします。
さらに波源の立ち上げを緩やかにします。
波源の立ち上げを緩やかにするにはContinuousSrcのwidthというパラメータを変更すればよいです。
***ブロッホの周期的境界条件の下でシミュレーションした結果をh5topngすると"h5topng error: data can have at most two dimensions (try specifying a slice)"というエラーが出る。 [#h2f15931]
ブロッホの周期的境界条件の下でシミュレーションした結果を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"を設定してください。
**参考になるサイト [#f39b395a]
-[https://core.ac.uk/download/pdf/38256115.pdf ターゲット特性による超高強度レーザー生成高速電子の制御に関する研究(大阪大学リポジトリ]
-[http://kenji-ishikawa.cocolog-nifty.com/plasma/ プラズマナノテクノロジー(ブログ)]
*コードの解析 [#c4799b48]
git clone https://github.com/stevengj/meep
./autogen.sh
make
とするとコードの完全体が得られます。
[https://github.com/stevengj/libctl/ libctl]というschemeの自作ライブラリが読み込まれています。
[https://github.com/stevengj/libctl/ libctl]というSchemeの自作ライブラリが読み込まれています。
**meep.scm.in [#o9d56706]
-meep.scm.inのdefine群が読み込まれています。
-ctlファイルで使用する関数の定義が書かれています。
Autotools で処理されるべきファイルに .in という拡張子が使われます。
つまり、その段階ではまだ Scheme のコードのテンプレートであり、処理されることで Scheme のコードとして妥当なものになります。
ビルドプロセスによりますが、通常は configure が済んだ段階では処理済みのもの (この場合は meep.scm) が出来ていると思います。
もちろん make が済んだ時点なら確実に出来ているでしょう。
***クラスの定義 [#qa52258e]
(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)
**meep.hpp [#ybf06d67]
- シミュレーション本体(C++)のすべてのクラスの前方宣言が含まれています。
*他の電磁界シミュレーションソフト [#qf1fc505]
株式会社EEMが無料で公開しているOpenFDTDとOpenMOMは[[OpenFDTD DLサイト:http://www.e-em.co.jp/OpenFDTD/index.html]]、[[OpenMOM DLサイト:http://www.e-em.co.jp/OpenMOM/index.html]]からダウンロードすることができます。
使い方については、DLサイトならびに、研究室が所蔵している[[RFワールド No.39 2017年 08月号 :https://www.amazon.co.jp/RF%E3%83%AF%E3%83%BC%E3%83%AB%E3%83%89-No-39-2017%E5%B9%B4-08-%E6%9C%88%E5%8F%B7/dp/B07233YBBQ]]に記載されています。
書籍の方は現在（2018/6/12)丸岡が借りています。
![[PukiWiki]](image/c60.png "[PukiWiki]")
