Compare commits
31 Commits
Author | SHA1 | Date |
---|---|---|
Jan Petykiewicz | b703f1ee20 | 4 months ago |
Jan Petykiewicz | 5b1e758c27 | 4 months ago |
Jan Petykiewicz | d0011cb1f9 | 4 months ago |
Jan Petykiewicz | 654f9c165b | 4 months ago |
jan | 689b3176cc | 1 year ago |
Jan Petykiewicz | 77c10feead | 1 year ago |
Jan Petykiewicz | 08a9328390 | 1 year ago |
Jan Petykiewicz | f489eb6e06 | 3 years ago |
Jan Petykiewicz | b2043ce841 | 3 years ago |
Jan Petykiewicz | c2e5c94e14 | 3 years ago |
Jan Petykiewicz | c75a7e67c9 | 3 years ago |
Jan Petykiewicz | 26ae8578c9 | 3 years ago |
Jan Petykiewicz | 4c81f33d05 | 3 years ago |
Jan Petykiewicz | 1403b9e1a3 | 3 years ago |
Jan Petykiewicz | e58efd2127 | 3 years ago |
Jan Petykiewicz | a6384429ee | 3 years ago |
Jan Petykiewicz | 292fde05d5 | 3 years ago |
Jan Petykiewicz | bbd33844b3 | 3 years ago |
Jan Petykiewicz | 2b96e5dd84 | 3 years ago |
Jan Petykiewicz | 922161c0e6 | 3 years ago |
Jan Petykiewicz | da3216133a | 3 years ago |
Jan Petykiewicz | 45e047c50e | 3 years ago |
Jan Petykiewicz | 8afe1d1f26 | 4 years ago |
Jan Petykiewicz | 8cbb0e9864 | 4 years ago |
Jan Petykiewicz | 4e0bd8b3c6 | 4 years ago |
Jan Petykiewicz | a4dd031666 | 5 years ago |
Jan Petykiewicz | 0c7b2fd3a2 | 5 years ago |
Jan Petykiewicz | 23d65d6ebb | 5 years ago |
Jan Petykiewicz | cf0db63a3f | 5 years ago |
Jan Petykiewicz | 89d99187b6 | 5 years ago |
Jan Petykiewicz | 1627d0da05 | 5 years ago |
@ -0,0 +1,29 @@
|
|||||||
|
[flake8]
|
||||||
|
ignore =
|
||||||
|
# E501 line too long
|
||||||
|
E501,
|
||||||
|
# W391 newlines at EOF
|
||||||
|
W391,
|
||||||
|
# E241 multiple spaces after comma
|
||||||
|
E241,
|
||||||
|
# E302 expected 2 newlines
|
||||||
|
E302,
|
||||||
|
# W503 line break before binary operator (to be deprecated)
|
||||||
|
W503,
|
||||||
|
# E265 block comment should start with '# '
|
||||||
|
E265,
|
||||||
|
# E123 closing bracket does not match indentation of opening bracket's line
|
||||||
|
E123,
|
||||||
|
# E124 closing bracket does not match visual indentation
|
||||||
|
E124,
|
||||||
|
# E221 multiple spaces before operator
|
||||||
|
E221,
|
||||||
|
# E201 whitespace after '['
|
||||||
|
E201,
|
||||||
|
# E741 ambiguous variable name 'I'
|
||||||
|
E741,
|
||||||
|
|
||||||
|
|
||||||
|
per-file-ignores =
|
||||||
|
# F401 import without use
|
||||||
|
*/__init__.py: F401,
|
@ -1,9 +1,11 @@
|
|||||||
.idea/
|
.idea/
|
||||||
*.h5
|
*.h5
|
||||||
|
|
||||||
__pycache__
|
__pycache__/
|
||||||
*.py[cod]
|
*.py[cod]
|
||||||
build/
|
build/
|
||||||
dist/
|
dist/
|
||||||
*.egg-info/
|
*.egg-info/
|
||||||
|
|
||||||
|
.pytest_cache/
|
||||||
|
.mypy_cache/
|
||||||
|
@ -0,0 +1 @@
|
|||||||
|
../LICENSE.md
|
@ -0,0 +1 @@
|
|||||||
|
../README.md
|
@ -1,5 +1,5 @@
|
|||||||
from .simulation import Simulation, type_to_C
|
from .simulation import Simulation, type_to_C
|
||||||
|
|
||||||
__author__ = 'Jan Petykiewicz'
|
__author__ = 'Jan Petykiewicz'
|
||||||
|
__version__ = '0.4'
|
||||||
version = '0.4'
|
version = __version__
|
||||||
|
@ -1,36 +0,0 @@
|
|||||||
/*
|
|
||||||
* Update E-field, including any PMLs.
|
|
||||||
*
|
|
||||||
* Template parameters:
|
|
||||||
* common_header: Rendered contents of common.cl
|
|
||||||
* pmls: [('x', 'n'), ('z', 'p'),...] list of pml axes and polarities
|
|
||||||
* pml_thickness: Number of cells (integer)
|
|
||||||
*
|
|
||||||
* OpenCL args:
|
|
||||||
* E, H, dt, S, oS
|
|
||||||
*/
|
|
||||||
|
|
||||||
{{common_header}}
|
|
||||||
|
|
||||||
//////////////////////////////////////////////////////////////////////
|
|
||||||
|
|
||||||
|
|
||||||
/*
|
|
||||||
* Calculate S from oS (pre-calculated components)
|
|
||||||
*/
|
|
||||||
__global ftype *Sx = S + XX;
|
|
||||||
__global ftype *Sy = S + YY;
|
|
||||||
__global ftype *Sz = S + ZZ;
|
|
||||||
|
|
||||||
// Use unscaled S components from H locations
|
|
||||||
__global ftype *oSxy = oS + 0 * field_size;
|
|
||||||
__global ftype *oSyz = oS + 1 * field_size;
|
|
||||||
__global ftype *oSzx = oS + 2 * field_size;
|
|
||||||
__global ftype *oSxz = oS + 3 * field_size;
|
|
||||||
__global ftype *oSyx = oS + 4 * field_size;
|
|
||||||
__global ftype *oSzy = oS + 5 * field_size;
|
|
||||||
|
|
||||||
ftype s_factor = dt * 0.125;
|
|
||||||
Sx[i] = (oSxy[i] + oSxz[i] + oSxy[i + my] + oSxz[i + mz]) * s_factor;
|
|
||||||
Sy[i] = (oSyz[i] + oSyx[i] + oSyz[i + mz] + oSyx[i + mx]) * s_factor;
|
|
||||||
Sz[i] = (oSzx[i] + oSzy[i] + oSzx[i + mx] + oSzy[i + my]) * s_factor;
|
|
@ -0,0 +1,381 @@
|
|||||||
|
import unittest
|
||||||
|
import numpy
|
||||||
|
|
||||||
|
from opencl_fdtd import Simulation
|
||||||
|
from meanas import fdtd
|
||||||
|
|
||||||
|
|
||||||
|
class BasicTests():
|
||||||
|
def test_initial_fields(self):
|
||||||
|
# Make sure initial fields didn't change
|
||||||
|
e0 = self.es[0]
|
||||||
|
h0 = self.hs[0]
|
||||||
|
mask = self.src_mask
|
||||||
|
|
||||||
|
self.assertEqual(e0[mask], self.j_mag / self.epsilon[mask])
|
||||||
|
self.assertFalse(e0[~mask].any())
|
||||||
|
self.assertFalse(h0.any())
|
||||||
|
|
||||||
|
|
||||||
|
def test_initial_energy(self):
|
||||||
|
e0 = self.es[0]
|
||||||
|
h0 = self.hs[0]
|
||||||
|
h1 = self.hs[1]
|
||||||
|
mask = self.src_mask[1]
|
||||||
|
dxes = self.dxes if self.dxes is not None else tuple(tuple(numpy.ones(s) for s in e0.shape[1:]) for _ in range(2))
|
||||||
|
dV = numpy.prod(numpy.meshgrid(*dxes[0], indexing='ij'), axis=0)
|
||||||
|
u0 = self.j_mag * self.j_mag / self.epsilon[self.src_mask] * dV[mask]
|
||||||
|
args = {
|
||||||
|
'dxes': self.dxes,
|
||||||
|
'epsilon': self.epsilon,
|
||||||
|
}
|
||||||
|
|
||||||
|
# Make sure initial energy and E dot J are correct
|
||||||
|
energy0 = fdtd.energy_estep(h0=h0, e1=e0, h2=self.hs[1], **args)
|
||||||
|
e_dot_j_0 = fdtd.delta_energy_j(j0=(e0 - 0) * self.epsilon, e1=e0, dxes=self.dxes)
|
||||||
|
self.assertTrue(numpy.allclose(energy0[mask], u0))
|
||||||
|
self.assertFalse(energy0[~mask].any(), msg=f'{energy0=}')
|
||||||
|
self.assertTrue(numpy.allclose(e_dot_j_0[mask], u0))
|
||||||
|
self.assertFalse(e_dot_j_0[~mask].any(), msg=f'{e_dot_j_0=}')
|
||||||
|
|
||||||
|
|
||||||
|
def test_energy_conservation(self):
|
||||||
|
e0 = self.es[0]
|
||||||
|
u0 = fdtd.delta_energy_j(j0=(e0 - 0) * self.epsilon, e1=e0, dxes=self.dxes).sum()
|
||||||
|
args = {'dxes': self.dxes,
|
||||||
|
'epsilon': self.epsilon}
|
||||||
|
|
||||||
|
for ii in range(1, 8):
|
||||||
|
with self.subTest(i=ii):
|
||||||
|
u_hstep = fdtd.energy_hstep(e0=self.es[ii-1], h1=self.hs[ii], e2=self.es[ii], **args)
|
||||||
|
u_estep = fdtd.energy_estep(h0=self.hs[ii], e1=self.es[ii], h2=self.hs[ii + 1], **args)
|
||||||
|
self.assertTrue(numpy.allclose(u_hstep.sum(), u0), msg=f'u_hstep: {u_hstep.sum()}\n{numpy.moveaxis(u_hstep, -1, 0)}')
|
||||||
|
self.assertTrue(numpy.allclose(u_estep.sum(), u0), msg=f'u_estep: {u_estep.sum()}\n{numpy.moveaxis(u_estep, -1, 0)}')
|
||||||
|
|
||||||
|
|
||||||
|
def test_poynting(self):
|
||||||
|
for ii in range(1, 3):
|
||||||
|
with self.subTest(i=ii):
|
||||||
|
s = fdtd.poynting(e=self.es[ii], h=self.hs[ii+1] + self.hs[ii])
|
||||||
|
sf = numpy.moveaxis(s, -1, 0)
|
||||||
|
ss = numpy.moveaxis(self.ss[ii], -1, 0)
|
||||||
|
self.assertTrue(numpy.allclose(s, self.ss[ii], rtol=1e-4),
|
||||||
|
msg=f'From ExH:\n{sf}\nFrom sim.S:\n{ss}')
|
||||||
|
|
||||||
|
|
||||||
|
def test_poynting_divergence(self):
|
||||||
|
args = {
|
||||||
|
'dxes': self.dxes,
|
||||||
|
'epsilon': self.epsilon,
|
||||||
|
}
|
||||||
|
|
||||||
|
dxes = self.dxes if self.dxes is not None else tuple(tuple(numpy.ones(s) for s in self.epsilon.shape[1:]) for _ in range(2))
|
||||||
|
dV = numpy.prod(numpy.meshgrid(*dxes[0], indexing='ij'), axis=0)
|
||||||
|
|
||||||
|
u_eprev = None
|
||||||
|
for ii in range(1, 8):
|
||||||
|
with self.subTest(i=ii):
|
||||||
|
u_hstep = fdtd.energy_hstep(e0=self.es[ii-1], h1=self.hs[ii], e2=self.es[ii], **args)
|
||||||
|
u_estep = fdtd.energy_estep(h0=self.hs[ii], e1=self.es[ii], h2=self.hs[ii + 1], **args)
|
||||||
|
|
||||||
|
du_half_h2e = u_estep - u_hstep
|
||||||
|
div_s_h2e = self.dt * fdtd.poynting_divergence(e=self.es[ii], h=self.hs[ii], dxes=self.dxes) * dV
|
||||||
|
|
||||||
|
du_half_h2e_f = numpy.moveaxis(du_half_h2e, -1, 0)
|
||||||
|
div_s_h2e_f = -numpy.moveaxis(div_s_h2e, -1, 0)
|
||||||
|
self.assertTrue(numpy.allclose(du_half_h2e, -div_s_h2e, rtol=1e-4),
|
||||||
|
msg=f'du_half_h2e\n{du_half_h2e_f}\ndiv_s_h2e\n{div_s_h2e_f}')
|
||||||
|
|
||||||
|
if u_eprev is None:
|
||||||
|
u_eprev = u_estep
|
||||||
|
continue
|
||||||
|
|
||||||
|
# previous half-step
|
||||||
|
du_half_e2h = u_hstep - u_eprev
|
||||||
|
div_s_e2h = self.dt * fdtd.poynting_divergence(e=self.es[ii-1], h=self.hs[ii], dxes=self.dxes) * dV
|
||||||
|
du_half_e2h_f = numpy.moveaxis(du_half_e2h, -1, 0)
|
||||||
|
div_s_e2h_f = -numpy.moveaxis(div_s_e2h, -1, 0)
|
||||||
|
self.assertTrue(numpy.allclose(du_half_e2h, -div_s_e2h, rtol=1e-4),
|
||||||
|
msg=f'du_half_e2h\n{du_half_e2h_f}\ndiv_s_e2h\n{div_s_e2h_f}')
|
||||||
|
u_eprev = u_estep
|
||||||
|
|
||||||
|
|
||||||
|
def test_poynting_planes(self):
|
||||||
|
args = {
|
||||||
|
'dxes': self.dxes,
|
||||||
|
'epsilon': self.epsilon,
|
||||||
|
}
|
||||||
|
dxes = self.dxes if self.dxes is not None else tuple(tuple(numpy.ones(s) for s in self.epsilon.shape[1:]) for _ in range(2))
|
||||||
|
dV = numpy.prod(numpy.meshgrid(*dxes[0], indexing='ij'), axis=0)
|
||||||
|
|
||||||
|
mx = numpy.roll(self.src_mask, (-1, -1), axis=(0, 1))
|
||||||
|
my = numpy.roll(self.src_mask, -1, axis=2)
|
||||||
|
mz = numpy.roll(self.src_mask, (+1, -1), axis=(0, 3))
|
||||||
|
px = numpy.roll(self.src_mask, -1, axis=0)
|
||||||
|
py = self.src_mask.copy()
|
||||||
|
pz = numpy.roll(self.src_mask, +1, axis=0)
|
||||||
|
|
||||||
|
u_eprev = None
|
||||||
|
for ii in range(1, 8):
|
||||||
|
with self.subTest(i=ii):
|
||||||
|
u_hstep = fdtd.energy_hstep(e0=self.es[ii-1], h1=self.hs[ii], e2=self.es[ii], **args)
|
||||||
|
u_estep = fdtd.energy_estep(h0=self.hs[ii], e1=self.es[ii], h2=self.hs[ii + 1], **args)
|
||||||
|
|
||||||
|
s_h2e = -fdtd.poynting(e=self.es[ii], h=self.hs[ii]) * self.dt
|
||||||
|
s_h2e[0] *= dxes[0][1][None, :, None] * dxes[0][2][None, None, :]
|
||||||
|
s_h2e[1] *= dxes[0][0][:, None, None] * dxes[0][2][None, None, :]
|
||||||
|
s_h2e[2] *= dxes[0][0][:, None, None] * dxes[0][1][None, :, None]
|
||||||
|
planes = [s_h2e[px].sum(), -s_h2e[mx].sum(),
|
||||||
|
s_h2e[py].sum(), -s_h2e[my].sum(),
|
||||||
|
s_h2e[pz].sum(), -s_h2e[mz].sum()]
|
||||||
|
du = (u_estep - u_hstep)[self.src_mask[1]]
|
||||||
|
self.assertTrue(numpy.allclose(sum(planes), (u_estep - u_hstep)[self.src_mask[1]]),
|
||||||
|
msg=f'planes: {planes} (sum: {sum(planes)})\n du:\n {du}')
|
||||||
|
|
||||||
|
if u_eprev is None:
|
||||||
|
u_eprev = u_estep
|
||||||
|
continue
|
||||||
|
|
||||||
|
s_e2h = -fdtd.poynting(e=self.es[ii - 1], h=self.hs[ii]) * self.dt
|
||||||
|
s_e2h[0] *= dxes[0][1][None, :, None] * dxes[0][2][None, None, :]
|
||||||
|
s_e2h[1] *= dxes[0][0][:, None, None] * dxes[0][2][None, None, :]
|
||||||
|
s_e2h[2] *= dxes[0][0][:, None, None] * dxes[0][1][None, :, None]
|
||||||
|
planes = [s_e2h[px].sum(), -s_e2h[mx].sum(),
|
||||||
|
s_e2h[py].sum(), -s_e2h[my].sum(),
|
||||||
|
s_e2h[pz].sum(), -s_e2h[mz].sum()]
|
||||||
|
du = (u_hstep - u_eprev)[self.src_mask[1]]
|
||||||
|
self.assertTrue(numpy.allclose(sum(planes), (u_hstep - u_eprev)[self.src_mask[1]]),
|
||||||
|
msg=f'planes: {du} (sum: {sum(planes)})\n du:\n {du}')
|
||||||
|
|
||||||
|
# previous half-step
|
||||||
|
u_eprev = u_estep
|
||||||
|
|
||||||
|
|
||||||
|
class Basic2DNoDXOnlyVacuum(unittest.TestCase, BasicTests):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 5, 1]
|
||||||
|
self.dt = 0.5
|
||||||
|
self.epsilon = numpy.ones(shape, dtype=float)
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = None
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 2, 0] = True
|
||||||
|
|
||||||
|
e = numpy.zeros_like(self.epsilon)
|
||||||
|
h = numpy.zeros_like(self.epsilon)
|
||||||
|
e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask]
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True)
|
||||||
|
|
||||||
|
for _ in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
class Basic2DUniformDX3(unittest.TestCase, BasicTests):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 1, 5]
|
||||||
|
self.dt = 0.5
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = tuple(tuple(numpy.full(s, 2.0) for s in shape[1:]) for _ in range(2))
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros(shape, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 0, 2] = True
|
||||||
|
|
||||||
|
self.epsilon = numpy.full(shape, 1, dtype=float)
|
||||||
|
self.epsilon[self.src_mask] = 2
|
||||||
|
|
||||||
|
e = numpy.zeros_like(self.epsilon)
|
||||||
|
h = numpy.zeros_like(self.epsilon)
|
||||||
|
e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask]
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True)
|
||||||
|
|
||||||
|
for _ in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
class Basic3DUniformDXOnlyVacuum(unittest.TestCase, BasicTests):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 5, 5]
|
||||||
|
self.dt = 0.5
|
||||||
|
self.epsilon = numpy.ones(shape, dtype=float)
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2))
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 2, 2] = True
|
||||||
|
|
||||||
|
e = numpy.zeros_like(self.epsilon)
|
||||||
|
h = numpy.zeros_like(self.epsilon)
|
||||||
|
e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask]
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True)
|
||||||
|
|
||||||
|
for _ in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
class Basic3DUniformDXUniformN(unittest.TestCase, BasicTests):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 5, 5]
|
||||||
|
self.dt = 0.5
|
||||||
|
self.epsilon = numpy.full(shape, 2.5, dtype=float)
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2))
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros_like(self.epsilon, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 2, 2] = True
|
||||||
|
|
||||||
|
e = numpy.zeros_like(self.epsilon)
|
||||||
|
h = numpy.zeros_like(self.epsilon)
|
||||||
|
e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask]
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True)
|
||||||
|
|
||||||
|
for _ in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
class Basic3DUniformDX(unittest.TestCase, BasicTests):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 5, 5]
|
||||||
|
self.dt = 0.33
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = tuple(tuple(numpy.ones(s) for s in shape[1:]) for _ in range(2))
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros(shape, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 2, 2] = True
|
||||||
|
|
||||||
|
self.epsilon = numpy.full(shape, 1, dtype=float)
|
||||||
|
self.epsilon[self.src_mask] = 2
|
||||||
|
|
||||||
|
e = numpy.zeros_like(self.epsilon)
|
||||||
|
h = numpy.zeros_like(self.epsilon)
|
||||||
|
e[self.src_mask] = self.j_mag / self.epsilon[self.src_mask]
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True)
|
||||||
|
|
||||||
|
for _ in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
class JdotE_3DUniformDX(unittest.TestCase):
|
||||||
|
def setUp(self):
|
||||||
|
shape = [3, 5, 5, 5]
|
||||||
|
self.dt = 0.5
|
||||||
|
self.j_mag = 32
|
||||||
|
self.dxes = tuple(tuple(numpy.full(s, 2.0) for s in shape[1:]) for _ in range(2))
|
||||||
|
|
||||||
|
self.src_mask = numpy.zeros(shape, dtype=bool)
|
||||||
|
self.src_mask[1, 2, 2, 2] = True
|
||||||
|
|
||||||
|
self.epsilon = numpy.full(shape, 4, dtype=float)
|
||||||
|
self.epsilon[self.src_mask] = 2
|
||||||
|
|
||||||
|
e = numpy.random.randint(-128, 128 + 1, size=shape).astype(float)
|
||||||
|
h = numpy.random.randint(-128, 128 + 1, size=shape).astype(float)
|
||||||
|
self.es = [e]
|
||||||
|
self.hs = [h]
|
||||||
|
self.ss = []
|
||||||
|
|
||||||
|
sim = Simulation(epsilon=self.epsilon, pmls=[], dt=self.dt, dxes=self.dxes,
|
||||||
|
initial_fields={'E': e, 'H': h}, do_poynting=True, float_type=numpy.float64)
|
||||||
|
|
||||||
|
eh2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes)
|
||||||
|
eh2e = fdtd.maxwell_e(dt=self.dt, dxes=self.dxes)
|
||||||
|
for ii in range(9):
|
||||||
|
e = e.copy()
|
||||||
|
h = h.copy()
|
||||||
|
sim.update_H([]).wait()
|
||||||
|
sim.update_E([]).wait()
|
||||||
|
|
||||||
|
if ii == 1:
|
||||||
|
nn = numpy.where(self.src_mask.flat)[0][0]
|
||||||
|
self.e_before = sim.E.get().reshape(shape)
|
||||||
|
sim.E[nn] += self.j_mag / self.epsilon[self.src_mask][0]
|
||||||
|
self.e_after = sim.E.get().reshape(shape)
|
||||||
|
self.j_dot_e = self.j_mag * e[self.src_mask]
|
||||||
|
|
||||||
|
self.hs.append(sim.H.get().reshape(shape))
|
||||||
|
self.es.append(sim.E.get().reshape(shape))
|
||||||
|
self.ss.append(sim.S.get().reshape(shape))
|
||||||
|
|
||||||
|
|
||||||
|
def test_j_dot_e(self):
|
||||||
|
h0 = self.hs[2]
|
||||||
|
e0 = self.es[1]
|
||||||
|
e1 = self.es[2]
|
||||||
|
j0 = numpy.zeros_like(e0)
|
||||||
|
j0[self.src_mask] = self.j_mag
|
||||||
|
args = {
|
||||||
|
'dxes': self.dxes,
|
||||||
|
'epsilon': self.epsilon,
|
||||||
|
}
|
||||||
|
e2h = fdtd.maxwell_h(dt=self.dt, dxes=self.dxes)
|
||||||
|
|
||||||
|
#ee = j0 * (2 * e0 - j0)
|
||||||
|
#hh = h0 * e2h(j0, numpy.zeros_like(h0))
|
||||||
|
#u0 = fdtd.dxmul(ee, hh, **args)
|
||||||
|
u0 = fdtd.delta_energy_j(j0=j0, e1=(e0 + e1), dxes=self.dxes)
|
||||||
|
|
||||||
|
#uh = [fdtd.energy_hstep(e0=self.es[ii-1], h1=self.hs[ii], e2=self.es[ii], **args) for ii in range(1, 9)]
|
||||||
|
#ue = [fdtd.energy_estep(h0=self.hs[ii], e1=self.es[ii], h2=self.hs[ii + 1], **args) for ii in range(8)]
|
||||||
|
#uht = [uu.sum() for uu in uh]
|
||||||
|
#uet = [uu.sum() for uu in ue]
|
||||||
|
|
||||||
|
u_hstep = fdtd.energy_hstep(e0=self.es[0], h1=self.hs[1], e2=self.es[1], **args)
|
||||||
|
u_estep = fdtd.energy_estep(h0=self.hs[-2], e1=self.es[-2], h2=self.hs[-1], **args)
|
||||||
|
#breakpoint()
|
||||||
|
du = (u_estep - u_hstep).sum()
|
||||||
|
self.assertTrue(numpy.allclose(u0.sum(), (u_estep - u_hstep).sum()), msg=f'{u0.sum()} != {du}')
|
@ -1,213 +0,0 @@
|
|||||||
"""
|
|
||||||
Routines for creating normalized 2D lattices and common photonic crystal
|
|
||||||
cavity designs.
|
|
||||||
"""
|
|
||||||
|
|
||||||
from typing import List
|
|
||||||
|
|
||||||
import numpy
|
|
||||||
|
|
||||||
|
|
||||||
def triangular_lattice(dims: List[int],
|
|
||||||
asymmetrical: bool=False
|
|
||||||
) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
Return an ndarray of [[x0, y0], [x1, y1], ...] denoting lattice sites for
|
|
||||||
a triangular lattice in 2D. The lattice will be centered around (0, 0),
|
|
||||||
unless asymmetrical=True in which case there will be extra holes in the +x
|
|
||||||
direction.
|
|
||||||
|
|
||||||
:param dims: Number of lattice sites in the [x, y] directions.
|
|
||||||
:param asymmetrical: If true, each row in x will contain the same number of
|
|
||||||
lattice sites. If false, the structure is symmetrical around (0, 0).
|
|
||||||
:return: [[x0, y0], [x1, 1], ...] denoting lattice sites.
|
|
||||||
"""
|
|
||||||
dims = numpy.array(dims, dtype=int)
|
|
||||||
|
|
||||||
if asymmetrical:
|
|
||||||
k = 0
|
|
||||||
else:
|
|
||||||
k = 1
|
|
||||||
|
|
||||||
positions = []
|
|
||||||
ymax = (dims[1] - 1)/2
|
|
||||||
for j in numpy.linspace(-ymax, ymax, dims[0]):
|
|
||||||
j_odd = numpy.floor(j) % 2
|
|
||||||
|
|
||||||
x_offset = j_odd * 0.5
|
|
||||||
y_offset = j * numpy.sqrt(3)/2
|
|
||||||
|
|
||||||
num_x = dims[0] - k * j_odd
|
|
||||||
xmax = (dims[0] - 1)/2
|
|
||||||
xs = numpy.linspace(-xmax, xmax - k * j_odd, num_x) + x_offset
|
|
||||||
ys = numpy.full_like(xs, y_offset)
|
|
||||||
|
|
||||||
positions += [numpy.vstack((xs, ys)).T]
|
|
||||||
|
|
||||||
xy = numpy.vstack(tuple(positions))
|
|
||||||
return xy[xy[:, 0].argsort(), ]
|
|
||||||
|
|
||||||
|
|
||||||
def square_lattice(dims: List[int]) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
Return an ndarray of [[x0, y0], [x1, y1], ...] denoting lattice sites for
|
|
||||||
a square lattice in 2D. The lattice will be centered around (0, 0).
|
|
||||||
|
|
||||||
:param dims: Number of lattice sites in the [x, y] directions.
|
|
||||||
:return: [[x0, y0], [x1, 1], ...] denoting lattice sites.
|
|
||||||
"""
|
|
||||||
xs, ys = numpy.meshgrid(range(dims[0]), range(dims[1]), 'xy')
|
|
||||||
xs -= dims[0]/2
|
|
||||||
ys -= dims[1]/2
|
|
||||||
xy = numpy.vstack((xs.flatten(), ys.flatten())).T
|
|
||||||
return xy[xy[:, 0].argsort(), ]
|
|
||||||
|
|
||||||
# ### Photonic crystal functions ###
|
|
||||||
|
|
||||||
|
|
||||||
def nanobeam_holes(a_defect: float,
|
|
||||||
num_defect_holes: int,
|
|
||||||
num_mirror_holes: int
|
|
||||||
) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
Returns a list of [[x0, r0], [x1, r1], ...] of nanobeam hole positions and radii.
|
|
||||||
Creates a region in which the lattice constant and radius are progressively
|
|
||||||
(linearly) altered over num_defect_holes holes until they reach the value
|
|
||||||
specified by a_defect, then symmetrically returned to a lattice constant and
|
|
||||||
radius of 1, which is repeated num_mirror_holes times on each side.
|
|
||||||
|
|
||||||
:param a_defect: Minimum lattice constant for the defect, as a fraction of the
|
|
||||||
mirror lattice constant (ie., for no defect, a_defect = 1).
|
|
||||||
:param num_defect_holes: How many holes form the defect (per-side)
|
|
||||||
:param num_mirror_holes: How many holes form the mirror (per-side)
|
|
||||||
:return: Ndarray [[x0, r0], [x1, r1], ...] of nanobeam hole positions and radii.
|
|
||||||
"""
|
|
||||||
a_values = numpy.linspace(a_defect, 1, num_defect_holes, endpoint=False)
|
|
||||||
xs = a_values.cumsum() - (a_values[0] / 2) # Later mirroring makes center distance 2x as long
|
|
||||||
mirror_xs = numpy.arange(1, num_mirror_holes + 1) + xs[-1]
|
|
||||||
mirror_rs = numpy.ones_like(mirror_xs)
|
|
||||||
return numpy.vstack((numpy.hstack((-mirror_xs[::-1], -xs[::-1], xs, mirror_xs)),
|
|
||||||
numpy.hstack((mirror_rs[::-1], a_values[::-1], a_values, mirror_rs)))).T
|
|
||||||
|
|
||||||
|
|
||||||
def ln_defect(mirror_dims: List[int], defect_length: int) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
N-hole defect in a triangular lattice.
|
|
||||||
|
|
||||||
:param mirror_dims: [x, y] mirror lengths (number of holes). Total number of holes
|
|
||||||
is 2 * n + 1 in each direction.
|
|
||||||
:param defect_length: Length of defect. Should be an odd number.
|
|
||||||
:return: [[x0, y0], [x1, y1], ...] for all the holes
|
|
||||||
"""
|
|
||||||
if defect_length % 2 != 1:
|
|
||||||
raise Exception('defect_length must be odd!')
|
|
||||||
p = triangular_lattice([2 * d + 1 for d in mirror_dims])
|
|
||||||
half_length = numpy.floor(defect_length / 2)
|
|
||||||
hole_nums = numpy.arange(-half_length, half_length + 1)
|
|
||||||
holes_to_keep = numpy.in1d(p[:, 0], hole_nums, invert=True)
|
|
||||||
return p[numpy.logical_or(holes_to_keep, p[:, 1] != 0), ]
|
|
||||||
|
|
||||||
|
|
||||||
def ln_shift_defect(mirror_dims: List[int],
|
|
||||||
defect_length: int,
|
|
||||||
shifts_a: List[float]=(0.15, 0, 0.075),
|
|
||||||
shifts_r: List[float]=(1, 1, 1)
|
|
||||||
) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
N-hole defect with shifted holes (intended to give the mode a gaussian profile
|
|
||||||
in real- and k-space so as to improve both Q and confinement). Holes along the
|
|
||||||
defect line are shifted and altered according to the shifts_* parameters.
|
|
||||||
|
|
||||||
:param mirror_dims: [x, y] mirror lengths (number of holes). Total number of holes
|
|
||||||
is 2 * n + 1 in each direction.
|
|
||||||
:param defect_length: Length of defect. Should be an odd number.
|
|
||||||
:param shifts_a: Percentage of a to shift (1st, 2nd, 3rd,...) holes along the defect line
|
|
||||||
:param shifts_r: Factor to multiply the radius by. Should match length of shifts_a
|
|
||||||
:return: [[x0, y0, r0], [x1, y1, r1], ...] for all the holes
|
|
||||||
"""
|
|
||||||
if not hasattr(shifts_a, "__len__") and shifts_a is not None:
|
|
||||||
shifts_a = [shifts_a]
|
|
||||||
if not hasattr(shifts_r, "__len__") and shifts_r is not None:
|
|
||||||
shifts_r = [shifts_r]
|
|
||||||
|
|
||||||
xy = ln_defect(mirror_dims, defect_length)
|
|
||||||
|
|
||||||
# Add column for radius
|
|
||||||
xyr = numpy.hstack((xy, numpy.ones((xy.shape[0], 1))))
|
|
||||||
|
|
||||||
# Shift holes
|
|
||||||
# Expand shifts as necessary
|
|
||||||
n_shifted = max(len(shifts_a), len(shifts_r))
|
|
||||||
|
|
||||||
tmp_a = numpy.array(shifts_a)
|
|
||||||
shifts_a = numpy.ones((n_shifted, ))
|
|
||||||
shifts_a[:len(tmp_a)] = tmp_a
|
|
||||||
|
|
||||||
tmp_r = numpy.array(shifts_r)
|
|
||||||
shifts_r = numpy.ones((n_shifted, ))
|
|
||||||
shifts_r[:len(tmp_r)] = tmp_r
|
|
||||||
|
|
||||||
x_removed = numpy.floor(defect_length / 2)
|
|
||||||
|
|
||||||
for ind in range(n_shifted):
|
|
||||||
for sign in (-1, 1):
|
|
||||||
x_val = sign * (x_removed + ind + 1)
|
|
||||||
which = numpy.logical_and(xyr[:, 0] == x_val, xyr[:, 1] == 0)
|
|
||||||
xyr[which, ] = (x_val + numpy.sign(x_val) * shifts_a[ind], 0, shifts_r[ind])
|
|
||||||
|
|
||||||
return xyr
|
|
||||||
|
|
||||||
|
|
||||||
def r6_defect(mirror_dims: List[int]) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
R6 defect in a triangular lattice.
|
|
||||||
|
|
||||||
:param mirror_dims: [x, y] mirror lengths (number of holes). Total number of holes
|
|
||||||
is 2 * n + 1 in each direction.
|
|
||||||
:return: [[x0, y0], [x1, y1], ...] specifying hole centers.
|
|
||||||
"""
|
|
||||||
xy = triangular_lattice([2 * d + 1 for d in mirror_dims])
|
|
||||||
|
|
||||||
rem_holes_plus = numpy.array([[1, 0],
|
|
||||||
[0.5, +numpy.sqrt(3)/2],
|
|
||||||
[0.5, -numpy.sqrt(3)/2]])
|
|
||||||
rem_holes = numpy.vstack((rem_holes_plus, -rem_holes_plus))
|
|
||||||
|
|
||||||
for rem_xy in rem_holes:
|
|
||||||
xy = xy[(xy != rem_xy).any(axis=1), ]
|
|
||||||
|
|
||||||
return xy
|
|
||||||
|
|
||||||
|
|
||||||
def l3_shift_perturbed_defect(mirror_dims: List[int],
|
|
||||||
perturbed_radius: float=1.1,
|
|
||||||
shifts_a: List[float]=(),
|
|
||||||
shifts_r: List[float]=()
|
|
||||||
) -> numpy.ndarray:
|
|
||||||
"""
|
|
||||||
3-hole defect with perturbed hole sizes intended to form an upwards-directed
|
|
||||||
beam. Can also include shifted holes along the defect line, intended
|
|
||||||
to give the mode a more gaussian profile to improve Q.
|
|
||||||
|
|
||||||
:param mirror_dims: [x, y] mirror lengths (number of holes). Total number of holes
|
|
||||||
is 2 * n + 1 in each direction.
|
|
||||||
:param perturbed_radius: Amount to perturb the radius of the holes used for beam-forming
|
|
||||||
:param shifts_a: Percentage of a to shift (1st, 2nd, 3rd,...) holes along the defect line
|
|
||||||
:param shifts_r: Factor to multiply the radius by. Should match length of shifts_a
|
|
||||||
:return: [[x0, y0, r0], [x1, y1, r1], ...] for all the holes
|
|
||||||
"""
|
|
||||||
xyr = ln_shift_defect(mirror_dims, 3, shifts_a, shifts_r)
|
|
||||||
|
|
||||||
abs_x, abs_y = (numpy.fabs(xyr[:, i]) for i in (0, 1))
|
|
||||||
|
|
||||||
# Sorted unique xs and ys
|
|
||||||
# Ignore row y=0 because it might have shifted holes
|
|
||||||
xs = numpy.unique(abs_x[abs_x != 0])
|
|
||||||
ys = numpy.unique(abs_y)
|
|
||||||
|
|
||||||
# which holes should be perturbed? (xs[[3, 7]], ys[1]) and (xs[[2, 6]], ys[2])
|
|
||||||
perturbed_holes = ((xs[a], ys[b]) for a, b in ((3, 1), (7, 1), (2, 2), (6, 2)))
|
|
||||||
for row in xyr:
|
|
||||||
if numpy.fabs(row) in perturbed_holes:
|
|
||||||
row[2] = perturbed_radius
|
|
||||||
return xyr
|
|
@ -0,0 +1,46 @@
|
|||||||
|
[build-system]
|
||||||
|
requires = ["hatchling"]
|
||||||
|
build-backend = "hatchling.build"
|
||||||
|
|
||||||
|
[project]
|
||||||
|
name = "opencl_fdtd"
|
||||||
|
description = "OpenCL FDTD solver"
|
||||||
|
readme = "README.md"
|
||||||
|
license = { file = "LICENSE.md" }
|
||||||
|
authors = [
|
||||||
|
{ name="Jan Petykiewicz", email="jan@mpxd.net" },
|
||||||
|
]
|
||||||
|
homepage = "https://mpxd.net/code/jan/opencl_fdtd"
|
||||||
|
repository = "https://mpxd.net/code/jan/opencl_fdtd"
|
||||||
|
keywords = [
|
||||||
|
"FDTD",
|
||||||
|
"finite",
|
||||||
|
"difference",
|
||||||
|
"time",
|
||||||
|
"domain",
|
||||||
|
"simulation",
|
||||||
|
"optics",
|
||||||
|
"electromagnetic",
|
||||||
|
"dielectric",
|
||||||
|
]
|
||||||
|
classifiers = [
|
||||||
|
"Programming Language :: Python :: 3",
|
||||||
|
"Development Status :: 4 - Beta",
|
||||||
|
"Intended Audience :: Developers",
|
||||||
|
"Intended Audience :: Manufacturing",
|
||||||
|
"Intended Audience :: Science/Research",
|
||||||
|
"License :: OSI Approved :: GNU Affero General Public License v3",
|
||||||
|
"Topic :: Scientific/Engineering",
|
||||||
|
]
|
||||||
|
requires-python = ">=3.8"
|
||||||
|
dynamic = ["version"]
|
||||||
|
dependencies = [
|
||||||
|
"numpy~=1.21",
|
||||||
|
"pyopencl",
|
||||||
|
"jinja2",
|
||||||
|
"meanas>=0.3",
|
||||||
|
]
|
||||||
|
|
||||||
|
[tool.hatch.version]
|
||||||
|
path = "opencl_fdtd/__init__.py"
|
||||||
|
|
@ -1,30 +0,0 @@
|
|||||||
#!/usr/bin/env python3
|
|
||||||
|
|
||||||
from setuptools import setup, find_packages
|
|
||||||
import opencl_fdtd
|
|
||||||
|
|
||||||
with open('README.md', 'r') as f:
|
|
||||||
long_description = f.read()
|
|
||||||
|
|
||||||
setup(name='opencl_fdtd',
|
|
||||||
version=opencl_fdtd.version,
|
|
||||||
description='OpenCL FDTD solver',
|
|
||||||
long_description=long_description,
|
|
||||||
long_description_content_type='text/markdown',
|
|
||||||
author='Jan Petykiewicz',
|
|
||||||
author_email='anewusername@gmail.com',
|
|
||||||
url='https://mpxd.net/code/jan/opencl_fdtd',
|
|
||||||
packages=find_packages(),
|
|
||||||
package_data={
|
|
||||||
'opencl_fdfd': ['kernels/*']
|
|
||||||
},
|
|
||||||
install_requires=[
|
|
||||||
'numpy',
|
|
||||||
'pyopencl',
|
|
||||||
'jinja2',
|
|
||||||
'fdfd_tools>=0.3',
|
|
||||||
],
|
|
||||||
extras_require={
|
|
||||||
},
|
|
||||||
)
|
|
||||||
|
|
Loading…
Reference in New Issue