ADDA + Matplotlib > InitField-Y for various values of Re{m} > v0.1, v0.2 (with 1d interpolation)

動作環境

GeForce GTX 1070 (8GB)
ASRock Z170M Pro4S [Intel Z170chipset]
Ubuntu 16.04 LTS desktop amd64
TensorFlow v1.2.1
cuDNN v5.1 for Linux
CUDA v8.0
Python 3.5.2
IPython 6.0.0 -- An enhanced Interactive Python.
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.4) 5.4.0 20160609
GNU bash, version 4.3.48(1)-release (x86_64-pc-linux-gnu)
ADDA v.1.3b6

This article is related to ADDA (light scattering simulator based on the discrete dipole approximation).

In this articles, I will show a graph plotting of the InitField-Y for various values of Re{m} (real part of refractive indices).

Data preparation

I prepared the following code under ~/ADDA/adda/src/seq/CODE_171028.

loop_rem_171028.py

import numpy as np
import subprocess as sb
import sys

# on Python 3.5.2

# codingrule:PEP8

RUN_PARAM = "-store_int_field -grid 52 -shape chebyshev 0.7 12"

CENTER = 1.45
RANGE = 0.40
mrs = np.linspace(CENTER - RANGE, CENTER + RANGE, 21)

for amr in mrs:
    ami = 1e-5
    cmd = "./adda -m %f %f %s" % (amr, ami, RUN_PARAM)
    print(cmd)
    sb.run(cmd.split(), stdout=sb.DEVNULL)

Under ~/ADDA/adda/src/seq, I executed the following:

$ python3 CODE_171028/loop_rem_171028.py

Then, I moved the run****_chebyshev_* directories into ~/ADDA/adda/src/seq/CALC_171028_HW0001/.

Graph plotting

v0.1

Under ~/ADDA/adda/src/seq/CODE_171028, I made the following Jupyter notebook.

read_intfield_varRem_171028.ipynb

import numpy as np
from scipy import signal, interpolate
from matplotlib import pylab as plt
import glob
from scipy import signal, interpolate

'''
v0.1 Oct. 28, 2017
  - plot graph
  - loop for [IntField-Y] and [log] files
  - add read_center_initfield()
  - add get_real_m()
'''

# codingrule: PEP8


def get_real_m(inpfile):
    IDX_REM = 2  # 'refactive index: 1.05+1e-05'
    with open(inpfile) as logfd:
        lines = logfd.readlines()
        for aline in lines:
            if aline.find('refractive') == -1:
                continue
            # print(aline)
            elem = aline.replace('+', ' ').split(' ')
            return elem[IDX_REM]
        return 0.0


def read_center_initfield(inpfile):
    dat = np.genfromtxt(inpfile, delimiter=' ',
                        skip_header=1)
    ary = []
    for elem in dat:
        ary += [elem]

    vcenter = len(ary) // 2
    hcenter = len(ary[0]) // 2

    return ary[vcenter][hcenter]


listINTF = glob.glob("../CALC_171028_HW0001/*/IntField-Y")
listLOG = glob.glob("../CALC_171028_HW0001/*/log")
xs = []
ys = []
for idx in range(len(listINTF)):
    rem = get_real_m(listLOG[idx])
    intf = read_center_initfield(listINTF[idx])
    xs += [rem]
    ys += [intf]

print(xs)
print(ys)

plt.figure(dpi=200)
plt.plot(xs, ys, "o")
plt.xlabel('real part of refractive index')
plt.ylabel('Init Field-Y')
plt.show()

['1.33', '1.41', '1.25', '1.69', '1.05', '1.85', '1.17', '1.37', '1.49', '1.65', '1.77', '1.53', '1.09', '1.81', '1.13', '1.61', '1.57', '1.29', '1.73', '1.21', '1.45']
[-0.026233794340000002, -0.022115366130000001, -0.023346986680000001, 0.0098408259539999998, -0.0039200259860000001, 0.045906058600000001, -0.016363607209999999, -0.025000819000000001, -0.015049274919999999, 0.0016250158680000001, 0.031175279090000001, -0.01194725152, -0.0078601925140000008, 0.041183378409999999, -0.01215601187, -0.0044200969, -0.0086785273040000006, -0.02554396125, 0.019991008380000001, -0.020176115410000001, -0.018504119140000001]

v0.2 > with interpolation

I performed 1d interpolation using scipy.

read_intfield_varRem_171028.ipynb

import numpy as np
from scipy import signal, interpolate
from matplotlib import pylab as plt
import glob
from scipy import signal, interpolate

'''
v0.2 Oct. 28, 2017
  - interpolate
v0.1 Oct. 28, 2017
  - plot graph
  - loop for [IntField-Y] and [log] files
  - add read_center_initfield()
  - add get_real_m()
'''

# codingrule: PEP8


def get_real_m(inpfile):
    IDX_REM = 2  # 'refactive index: 1.05+1e-05'
    with open(inpfile) as logfd:
        lines = logfd.readlines()
        for aline in lines:
            if aline.find('refractive') == -1:
                continue
            # print(aline)
            elem = aline.replace('+', ' ').split(' ')
            return elem[IDX_REM]
        return 0.0


def read_center_initfield(inpfile):
    dat = np.genfromtxt(inpfile, delimiter=' ',
                        skip_header=1)
    ary = []
    for elem in dat:
        ary += [elem]

    vcenter = len(ary) // 2
    hcenter = len(ary[0]) // 2

    return ary[vcenter][hcenter]


listINTF = glob.glob("../CALC_171028_HW0001/*/IntField-Y")
listLOG = glob.glob("../CALC_171028_HW0001/*/log")
xs = []
ys = []
for idx in range(len(listINTF)):
    rem = get_real_m(listLOG[idx])
    intf = read_center_initfield(listINTF[idx])
    xs += [rem]
    ys += [intf]

print(xs)
print(ys)

# interpolation
xns = np.array(xs, dtype=float)
yns = np.array(ys)
tt = np.linspace(1.05, 1.85, 51)  # 51: arbitrary
f3 = interpolate.interp1d(xns, yns, kind="cubic")
y3 = f3(tt)

#
plt.figure(dpi=200)
plt.plot(xs, ys, "o")
plt.plot(tt, y3, "r")
plt.xlabel('real part of refractive index')
plt.ylabel('Init Field-Y')
plt.show()