[Started by Maruoka, 2018/10/24]
[Updated by Maruoka, 2018/01/25]
 
*index [#e16fd6bb]
#contents

*How to install [#ldc78cdd]
First of all, let's install Anaconda.
From [[The websit of Anaconda:https://www.anaconda.com/download/]], please download the installer of the series of Python3.x .
Since the installer is of the file of .sh, execute it as below.
 $ bash [download-folder]/Anaconda3-5.2.0-Linux-x86_64.sh
Then, let's install pymeep ([https://meep.readthedocs.io/en/latest/Installation/ ref] )
Execute the command below.
 $ conda create -n mp -c chogan -c conda-forge pymeep
All setup is finished.
*How to execute [#padb0b5f]
Activate the system as below.
 $ source activate mp
Execute the file as below.Please rename hogehoge.py as what you name it.
 (mp) $ python hogehoge.py
You can get sample source from Github.[https://github.com/stevengj/meep/blob/master/python/examples/wvg-src.py example]
If you want to deactivate the system, do as below.
 (mp) $ source deactivate

*parallelization of MEEP [#h659d054]
***installation([https://meep.readthedocs.io/en/latest/Installation/#mpi ref]) [#v7b0b7df]
In order to install pymeep which support parallel calculation, execute the following command.
 $ conda create -n pmp -c chogan -c conda-forge pymeep-parallel
***Execution [#ld3fff3a]
Activate pymeep which support parallel calculation by following command
 $ source activate pmp
Programs run with following command.
 $ mpirun -np 4 python <script_name>
"-np" indicate how many processors will run. In this case, 4 processors run simultaneously.


*How to edit the HDF5 file ? [#qe50ab81]
The h5py package can be used for editing HDF5 file in python.
Install it by the following command.
 apt-get install python-h5py
h5py is already installed on flex and tubes.

On maruoka's WSL,
 pip install h5py
doesn't work saying
 ImportError: No module named h5py


*Useful script for converting hdf5 to mp4 and gif [#b431428c]
It is useful to use script below to automatically converting hdf5 to mp4 by pasting this to .bash_profile
 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
 }

For example, if you want to convert hogehoge.h5 which have 1000 number of time indexes, do as follow.
 h5mp4 hogehoge 999
It output hogehoge.gif and hogehoge.mp4 automatically.

If you want to know the number of time indexes,do as follow.
 h5ls hogehoge.h5
 ez Dataset{100, 200, 1000/Inf}
Here, 1000 is the number of time indexes.

*source for incident light [#r32309cc]
**left-handed circularly polarized light (Scheme)[#u95702e4]
We can modify the phase of source by the 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)))
          )
        )
       )



**right-handed circularly polarized light (Scheme) [#u95702e4]

 (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)))
          )
        )
       )


**the light linearly polarized at x radian (Scheme)[#u95702e4]

 (set! sources
       (list
        (make source
          (src (make continuous-src (frequency (/ 1 7))))
          (component Ex)
          (center 0 0 8)
          (size 50 50 0)
          (amplitude (cos x))
          )
        (make source
          (src (make continuous-src (frequency (/ 1 7))))
          (component Ey)
          (center 0 0 8)
          (size 50 50 0)
          (amplitude (sin x))
          )
        )
       )
**the light at oblique incident angle (Scheme)[#y9a5c654]
There are three ways to implement oblique incident light.

***Bloch's periodic boundary with modifying phase of the source (Scheme)[#s8084013]
TM wave just before total internal reflection.
 (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))
            )

plot of Ez

http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/main-ez.gif




***two sources with modifying phase (Scheme) [#d178c899]
 ; This example creates an approximate Ez-polarized planewave in vacuum
 ; propagating at a 45-degree angle, by using a couple of current sources
 ; with amplitude exp(ikx) corresponding to the desired planewave.
 (define-param s 11) ; the size of the computational cell, not including PML
 (define-param dpml 1) ; thickness of PML layers
 (define sxy (+ s (* 2 dpml))) ; cell size, including PML
 (set! geometry-lattice (make lattice (size sxy sxy no-size)))
 (set! pml-layers (list (make pml (thickness dpml))))
 (set-param! resolution 10)
 ; pw-amp is a function that returns the amplitude exp(ik(x+x0)) at a
 ; given point x.  (We need the x0 because current amplitude functions
 ; in Meep are defined relative to the center of the current source,
 ; whereas we want a fixed origin.)  Actually, it is a function of k
 ; and x0 that returns a function of x ...
 (define (pw-amp k x0) (lambda (x)
   (exp (* 0+1i (vector3-dot k (vector3+ x x0))))))
 (define-param fcen 0.8) ; pulse center frequency
 (define-param df 0.02) ; turn-on bandwidth
 (define-param kdir (vector3 1 1)) ; direction of k (length is irrelevant)
 (define-param n 1) ; refractive index of material containing the source
 (define k (vector3-scale (* 2 pi fcen n)
                          (unit-vector3 kdir))) ; k with correct length
 (set! sources
       (list
        ; left
        (make source
 	 (src (make continuous-src (frequency fcen) (fwidth df)))
 	 (component Ez) (center (* -0.5 s) 0) (size 0 s)
 	 (amp-func (pw-amp k (vector3 (* -0.5 s) 0))))
        ; bottom
        (make source
 	 (src (make continuous-src (frequency fcen) (fwidth df)))
 	 (component Ez) (center 0 (* -0.5 s)) (size s 0)
 	 (amp-func (pw-amp k (vector3 0 (* -0.5 s)))))
        ))
 (define-param T 400) ; run time
 (run-until T (at-end output-efield-z))
[https://github.com/stevengj/meep/blob/master/scheme/examples/pw-source.ctl ref]

plot of Ez

http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-ez.gif

plot of Ez when there is only one source ( why we need two sources )

http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/double-sources/pw-source-single-ez.gif

***array many point-sources  (Scheme)[#bb9694c1]
http://article.gmane.org/gmane.comp.science.electromagnetism.meep.general/3058
 (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

http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/oblique-incident/array.gif

*Plot of Electric field Vector [#b59ce919]

Electric field of Total Internal Reflection

http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/evanescent/evanescent_afterEfield.gif

Code for simulation

 (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)))
 )

If you want to plot electric field vector, you need to plot it from hdf5 file by your code.
The ax.quiver function of matplotlib is good idea to plot electric field vector.


*material library [#ba28d968]
To use material library, add
 (include "/home/students/maru/meep/scheme/materials.scm")
at the beginning of your program.
This library only works at tube61, 62, 63.
For how to use the material library, check https://meep.readthedocs.io/en/latest/Materials/

You need to specify unit to use material library.
The default is 1000 micro meter.
In order to change the unit to  100 micro meter, add 
 (set! um-scale (* um-scale 0.05))

*Q&A [#b8d898d0]
**How thick PML should be ?[#l8406a8a]
In the document of PML by one of the the author of MEEP(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.)

However, when the light coming at oblique incident angle, PML should be thicker.
The reflectivity of PML for incident angle at X can be represented as 
e^{-2*k*d*cos(X)*G}
for some constant G. d is the thickness of PML. k is the wavenumber.(Please check https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=477075 for details).
If X approaches to 90 degree, the efficiency of the absorbance of PML becomes really bad.


**PML shows reflectance. [#o360c5d3]
If you put dielectrics touching to the PML, please confirm you overlap the dielectrics and PML. Otherwise, PML is useless because PML is defined for the boundary of vaccum to absorb the light not for dielectrics.
It is implemented by just multiplying tensors on dielectric constants. So you need to overlap PML and dielectrics.  (Please check https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=477075 for details).

[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/dielectrics_and_pml_touch_each_other.mp4 PML and Dielectrics touch each other but do not overlap]

[http://flex.phys.tohoku.ac.jp/~maru/drive-open/pukiwiki/dielectrics_and_pml_touch_each_other.ctl code]


** ERROR: Unbound variable: E-susceptibilities [#qd0361a8]
Following error occurs, when you try to use material library at FLEX or Tube60. 
 ERROR: Unbound variable: E-susceptibilities
Material library only works at tube61, 62, 63. 
Please go to tube61, 62, 63.

** materials.py doesn't work [#cc53b98c]
If the version is lower than 1.6, materials.py doesn't work.
Check the version of pymeep by
 conda list
If you find
 pymeep           1.5
then update pymeep from 1.5 to 1.6 as follows.

 (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

If you want to use parallel pymeep, then update it also to 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 ? [#af67089a]
If you find error as below, give permission to working directory for writing and edinting.

 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]

The version of GLIBC at flex is old. Go to t70.


Front page   New Page list Search Recent changes   Help   RSS of recent changes