meanas/examples/fdfd.py

256 lines
7.9 KiB
Python
Raw Normal View History

import importlib
2016-05-30 22:30:45 -07:00
import numpy
from numpy.linalg import norm
2016-05-30 22:30:45 -07:00
2019-08-04 14:13:51 -07:00
import meanas
2019-11-27 22:59:52 -08:00
from meanas import fdtd
from meanas.fdmath import vec, unvec
from meanas.fdfd import waveguide_3d, functional, scpml, operators
2019-08-04 14:13:51 -07:00
from meanas.fdfd.solvers import generic as generic_solver
2016-05-30 22:30:45 -07:00
import gridlock
from matplotlib import pyplot
2019-11-27 22:59:52 -08:00
import logging
logging.basicConfig(level=logging.DEBUG)
2021-09-05 17:49:08 -07:00
logging.getLogger('matplotlib').setLevel(logging.WARNING)
2016-07-04 16:35:28 -07:00
2016-05-30 22:30:45 -07:00
__author__ = 'Jan Petykiewicz'
def test0(solver=generic_solver):
2016-05-30 22:30:45 -07:00
dx = 50 # discretization (nm/cell)
pml_thickness = 10 # (number of cells)
wl = 1550 # Excitation wavelength
omega = 2 * numpy.pi / wl
# Device design parameters
radii = (1, 0.6)
th = 220
center = [0, 0, 0]
# refractive indices
n_ring = numpy.sqrt(12.6) # ~Si
n_air = 4.0 # air
# Half-dimensions of the simulation grid
xyz_max = numpy.array([1.2, 1.2, 0.3]) * 1000 + pml_thickness * dx
# Coordinates of the edges of the cells.
half_edge_coords = [numpy.arange(dx/2, m + dx, step=dx) for m in xyz_max]
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
# #### Create the grid, mask, and draw the device ####
2022-10-04 12:29:11 -07:00
grid = gridlock.Grid(edge_coords)
epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
grid.draw_cylinder(epsilon,
surface_normal=2,
2016-05-30 22:30:45 -07:00
center=center,
radius=max(radii),
thickness=th,
eps=n_ring**2,
num_points=24)
2022-10-04 12:29:11 -07:00
grid.draw_cylinder(epsilon,
surface_normal=2,
2016-05-30 22:30:45 -07:00
center=center,
radius=min(radii),
thickness=th*1.1,
eps=n_air ** 2,
num_points=24)
2016-07-03 14:23:24 -07:00
dxes = [grid.dxyz, grid.autoshifted_dxyz()]
2016-05-30 22:30:45 -07:00
for a in (0, 1, 2):
for p in (-1, 1):
2019-08-05 00:20:06 -07:00
dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega,
thickness=pml_thickness)
2016-05-30 22:30:45 -07:00
2022-10-04 12:29:11 -07:00
J = [numpy.zeros_like(epsilon[0], dtype=complex) for _ in range(3)]
2019-08-05 00:20:06 -07:00
J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1
2016-05-30 22:30:45 -07:00
2019-09-27 20:48:03 -07:00
'''
Solve!
'''
2019-08-05 00:20:06 -07:00
sim_args = {
'omega': omega,
'dxes': dxes,
2022-10-04 12:29:11 -07:00
'epsilon': vec(epsilon),
2019-08-05 00:20:06 -07:00
}
x = solver(J=vec(J), **sim_args)
2022-10-04 12:29:11 -07:00
A = operators.e_full(omega, dxes, vec(epsilon)).tocsr()
2016-05-30 22:30:45 -07:00
b = -1j * omega * vec(J)
print('Norm of the residual is ', norm(A @ x - b))
2016-05-30 22:30:45 -07:00
E = unvec(x, grid.shape)
'''
Plot results
'''
2016-05-30 22:30:45 -07:00
pyplot.figure()
pyplot.pcolor(numpy.real(E[1][:, :, grid.shape[2]//2]), cmap='seismic')
pyplot.axis('equal')
pyplot.show()
def test1(solver=generic_solver):
2016-05-30 22:30:45 -07:00
dx = 40 # discretization (nm/cell)
pml_thickness = 10 # (number of cells)
wl = 1550 # Excitation wavelength
omega = 2 * numpy.pi / wl
# Device design parameters
w = 600
th = 220
center = [0, 0, 0]
# refractive indices
n_wg = numpy.sqrt(12.6) # ~Si
n_air = 1.0 # air
# Half-dimensions of the simulation grid
xyz_max = numpy.array([0.8, 0.9, 0.6]) * 1000 + (pml_thickness + 2) * dx
# Coordinates of the edges of the cells.
half_edge_coords = [numpy.arange(dx/2, m + dx/2, step=dx) for m in xyz_max]
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
# #### Create the grid and draw the device ####
2022-10-04 12:29:11 -07:00
grid = gridlock.Grid(edge_coords)
epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], eps=n_wg**2)
2016-05-30 22:30:45 -07:00
2016-07-03 14:23:24 -07:00
dxes = [grid.dxyz, grid.autoshifted_dxyz()]
2016-05-30 22:30:45 -07:00
for a in (0, 1, 2):
for p in (-1, 1):
2019-08-04 14:13:51 -07:00
dxes = scpml.stretch_with_scpml(dxes,omega=omega, axis=a, polarity=p,
thickness=pml_thickness)
2016-05-30 22:30:45 -07:00
half_dims = numpy.array([10, 20, 15]) * dx
dims = [-half_dims, half_dims]
dims[1][0] = dims[0][0]
ind_dims = (grid.pos2ind(dims[0], which_shifts=None).astype(int),
grid.pos2ind(dims[1], which_shifts=None).astype(int))
2019-08-26 01:12:36 -07:00
src_axis = 0
2016-05-30 22:30:45 -07:00
wg_args = {
'slices': [slice(i, f+1) for i, f in zip(*ind_dims)],
'dxes': dxes,
2019-08-26 01:12:36 -07:00
'axis': src_axis,
2016-05-30 22:30:45 -07:00
'polarity': +1,
}
2022-10-04 12:29:11 -07:00
wg_results = waveguide_3d.solve_mode(mode_number=0, omega=omega, epsilon=epsilon, **wg_args)
2019-11-27 22:59:52 -08:00
J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'],
2022-10-04 12:29:11 -07:00
omega=omega, epsilon=epsilon, **wg_args)
2019-11-27 22:59:52 -08:00
e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args)
2016-05-30 22:30:45 -07:00
2022-10-04 12:29:11 -07:00
pecg = numpy.zeros_like(epsilon)
# pecg.draw_cuboid(pecg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
# pecg.visualize_isosurface(pecg)
2016-07-03 14:23:51 -07:00
2022-10-04 12:29:11 -07:00
pmcg = numpy.zeros_like(epsilon)
# grid.draw_cuboid(pmcg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
# grid.visualize_isosurface(pmcg)
2016-07-03 16:45:38 -07:00
2022-10-04 12:43:26 -07:00
def pcolor(v) -> None:
2019-08-26 01:12:36 -07:00
vmax = numpy.max(numpy.abs(v))
pyplot.pcolor(v, cmap='seismic', vmin=-vmax, vmax=vmax)
pyplot.axis('equal')
pyplot.colorbar()
2019-08-30 22:03:54 -07:00
ss = (1, slice(None), J.shape[2]//2+6, slice(None))
# pyplot.figure()
# pcolor(J3[ss].T.imag)
# pyplot.figure()
# pcolor((numpy.abs(J3).sum(axis=2).sum(axis=0) > 0).astype(float).T)
pyplot.show(block=True)
2016-07-04 16:35:28 -07:00
'''
Solve!
'''
sim_args = {
'omega': omega,
'dxes': dxes,
2022-10-04 12:29:11 -07:00
'epsilon': vec(epsilon),
'pec': vec(pecg),
'pmc': vec(pmcg),
}
x = solver(J=vec(J), **sim_args)
2016-05-30 22:30:45 -07:00
b = -1j * omega * vec(J)
2019-08-04 14:13:51 -07:00
A = operators.e_full(**sim_args).tocsr()
print('Norm of the residual is ', norm(A @ x - b))
2016-07-04 16:35:28 -07:00
2016-05-30 22:30:45 -07:00
E = unvec(x, grid.shape)
2016-07-04 16:35:28 -07:00
'''
Plot results
'''
2016-05-30 22:30:45 -07:00
center = grid.pos2ind([0, 0, 0], None).astype(int)
pyplot.figure()
pyplot.subplot(2, 2, 1)
2019-08-26 01:12:36 -07:00
pcolor(numpy.real(E[1][center[0], :, :]).T)
2016-05-30 22:30:45 -07:00
pyplot.subplot(2, 2, 2)
pyplot.plot(numpy.log10(numpy.abs(E[1][:, center[1], center[2]]) + 1e-10))
2019-10-27 12:46:12 -07:00
pyplot.grid(alpha=0.6)
pyplot.ylabel('log10 of field')
2016-05-30 22:30:45 -07:00
pyplot.subplot(2, 2, 3)
2019-08-26 01:12:36 -07:00
pcolor(numpy.real(E[1][:, :, center[2]]).T)
2016-05-30 22:30:45 -07:00
pyplot.subplot(2, 2, 4)
def poyntings(E):
2019-08-30 22:03:54 -07:00
H = functional.e2h(omega, dxes)(E)
2019-10-27 12:46:12 -07:00
poynting = fdtd.poynting(e=E, h=H.conj(), dxes=dxes)
2019-08-30 22:03:54 -07:00
cross1 = operators.poynting_e_cross(vec(E), dxes) @ vec(H).conj()
2019-09-27 20:48:03 -07:00
cross2 = operators.poynting_h_cross(vec(H), dxes) @ vec(E).conj() * -1
2019-10-27 12:46:12 -07:00
s1 = 0.5 * unvec(numpy.real(cross1), grid.shape)
s2 = 0.5 * unvec(numpy.real(cross2), grid.shape)
s0 = 0.5 * poynting.real
2019-08-30 22:03:54 -07:00
# s2 = poynting.imag
2019-09-27 20:48:03 -07:00
return s0, s1, s2
2016-05-30 22:30:45 -07:00
2019-09-27 20:48:03 -07:00
s0x, s1x, s2x = poyntings(E)
2019-10-27 12:46:12 -07:00
pyplot.plot(s0x[0].sum(axis=2).sum(axis=1), label='s0', marker='.')
pyplot.plot(s1x[0].sum(axis=2).sum(axis=1), label='s1', marker='.')
pyplot.plot(s2x[0].sum(axis=2).sum(axis=1), label='s2', marker='.')
pyplot.plot(E[1][:, center[1], center[2]].real.T, label='Ey', marker='x')
pyplot.grid(alpha=0.6)
2019-09-27 20:48:03 -07:00
pyplot.legend()
2016-05-30 22:30:45 -07:00
pyplot.show()
q = []
for i in range(-5, 30):
2019-09-05 13:42:39 -07:00
e_ovl_rolled = numpy.roll(e_overlap, i, axis=1)
q += [numpy.abs(vec(E) @ vec(e_ovl_rolled).conj())]
2016-05-30 22:30:45 -07:00
pyplot.figure()
2019-10-27 12:46:12 -07:00
pyplot.plot(q, marker='.')
pyplot.grid(alpha=0.6)
2016-05-30 22:30:45 -07:00
pyplot.title('Overlap with mode')
pyplot.show()
print('Average overlap with mode:', sum(q)/len(q))
def module_available(name):
return importlib.util.find_spec(name) is not None
2016-05-30 22:30:45 -07:00
if __name__ == '__main__':
2019-08-26 01:12:36 -07:00
#test0()
2019-10-27 12:46:12 -07:00
# test1()
if module_available('opencl_fdfd'):
from opencl_fdfd import cg_solver as opencl_solver
test1(opencl_solver)
# from opencl_fdfd.csr import fdfd_cg_solver as opencl_csr_solver
# test1(opencl_csr_solver)
# elif module_available('magma_fdfd'):
# from magma_fdfd import solver as magma_solver
# test1(magma_solver)
else:
test1()