From 1e80a66b50eaaf7f92f44d56f58c0f7bf76d8f6e Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 5 Mar 2017 17:20:38 -0800 Subject: [PATCH] add fdtd and test --- examples/test_fdtd.py | 175 +++++++++++++++++++++++++++++++ fdfd_tools/fdtd.py | 239 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 414 insertions(+) create mode 100644 examples/test_fdtd.py create mode 100644 fdfd_tools/fdtd.py diff --git a/examples/test_fdtd.py b/examples/test_fdtd.py new file mode 100644 index 0000000..1a25be4 --- /dev/null +++ b/examples/test_fdtd.py @@ -0,0 +1,175 @@ +""" +Example code for running an OpenCL FDTD simulation + +See main() for simulation setup. +""" + +import sys +import time + +import numpy +import h5py + +from fdfd_tools import fdtd +from masque import Pattern, shapes +import gridlock +import pcgen + + +def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern: + """ + Generate a masque.Pattern object containing a perturbed L3 cavity. + + :param a: Lattice constant. + :param radius: Hole radius, in units of a (lattice constant). + :param kwargs: Keyword arguments: + hole_dose, trench_dose, hole_layer, trench_layer: Shape properties for Pattern. + Defaults *_dose=1, hole_layer=0, trench_layer=1. + shifts_a, shifts_r: passed to pcgen.l3_shift; specifies lattice constant (1 - + multiplicative factor) and radius (multiplicative factor) for shifting + holes adjacent to the defect (same row). Defaults are 0.15 shift for + first hole, 0.075 shift for third hole, and no radius change. + xy_size: [x, y] number of mirror periods in each direction; total size is + 2 * n + 1 holes in each direction. Default [10, 10]. + perturbed_radius: radius of holes perturbed to form an upwards-driected beam + (multiplicative factor). Default 1.1. + trench width: Width of the undercut trenches. Default 1.2e3. + :return: masque.Pattern object containing the L3 design + """ + + default_args = {'hole_dose': 1, + 'trench_dose': 1, + 'hole_layer': 0, + 'trench_layer': 1, + 'shifts_a': (0.15, 0, 0.075), + 'shifts_r': (1.0, 1.0, 1.0), + 'xy_size': (10, 10), + 'perturbed_radius': 1.1, + 'trench_width': 1.2e3, + } + kwargs = {**default_args, **kwargs} + + xyr = pcgen.l3_shift_perturbed_defect(mirror_dims=kwargs['xy_size'], + perturbed_radius=kwargs['perturbed_radius'], + shifts_a=kwargs['shifts_a'], + shifts_r=kwargs['shifts_r']) + xyr *= a + xyr[:, 2] *= radius + + pat = Pattern() + pat.name = 'L3p-a{:g}r{:g}rp{:g}'.format(a, radius, kwargs['perturbed_radius']) + pat.shapes += [shapes.Circle(radius=r, offset=(x, y), + dose=kwargs['hole_dose'], + layer=kwargs['hole_layer']) + for x, y, r in xyr] + + maxes = numpy.max(numpy.fabs(xyr), axis=0) + pat.shapes += [shapes.Polygon.rectangle( + lx=(2 * maxes[0]), ly=kwargs['trench_width'], + offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2)), + dose=kwargs['trench_dose'], layer=kwargs['trench_layer']) + for s in (-1, 1)] + return pat + + +def main(): + dtype = numpy.float32 + max_t = 8000 # number of timesteps + + dx = 40 # discretization (nm/cell) + pml_thickness = 8 # (number of cells) + + wl = 1550 # Excitation wavelength and fwhm + dwl = 200 + + # Device design parameters + xy_size = numpy.array([10, 10]) + a = 430 + r = 0.285 + th = 170 + + # refractive indices + n_slab = 3.408 # InGaAsP(80, 50) @ 1550nm + n_air = 1.0 # air + + # Half-dimensions of the simulation grid + xy_max = (xy_size + 1) * a * [1, numpy.sqrt(3)/2] + z_max = 1.6 * a + xyz_max = numpy.hstack((xy_max, z_max)) + pml_thickness * dx + + # Coordinates of the edges of the cells. The fdtd package can only do square grids at the moment. + 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 #### + grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3) + grid.draw_slab(surface_normal=gridlock.Direction.z, + center=[0, 0, 0], + thickness=th, + eps=n_slab**2) + mask = perturbed_l3(a, r) + + grid.draw_polygons(surface_normal=gridlock.Direction.z, + center=[0, 0, 0], + thickness=2 * th, + eps=n_air**2, + polygons=mask.as_polygons()) + + print(grid.shape) + # #### Create the simulation grid #### + epsilon = [eps.astype(dtype) for eps in grid.grids] + + dt = .99/numpy.sqrt(3) + e = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)] + h = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)] + + update_e = fdtd.maxwell_e(dt) + update_h = fdtd.maxwell_h(dt) + + # PMLs in every direction + pml_e_funcs = [] + pml_h_funcs = [] + pml_fields = {} + for d in (0, 1, 2): + for p in (-1, 1): + ef, hf, psis = fdtd.cpml(direction=d, polarity=p, dt=dt, epsilon=epsilon, epsilon_eff=n_slab**2, dtype=dtype) + pml_e_funcs.append(ef) + pml_h_funcs.append(hf) + pml_fields.update(psis) + + # Source parameters and function + w = 2 * numpy.pi * dx / wl + fwhm = dwl * w * w / (2 * numpy.pi * dx) + alpha = (fwhm ** 2) / 8 * numpy.log(2) + delay = 7/numpy.sqrt(2 * alpha) + + def field_source(i): + t0 = i * dt - delay + return numpy.sin(w * t0) * numpy.exp(-alpha * t0**2) + + # #### Run a bunch of iterations #### + output_file = h5py.File('simulation_output.h5', 'w') + start = time.perf_counter() + for t in range(max_t): + [f(e, h, epsilon) for f in pml_e_funcs] + update_e(e, h, epsilon) + + e[1][tuple(grid.shape//2)] += field_source(t) + [f(e, h) for f in pml_h_funcs] + update_h(e, h) + + print('iteration {}: average {} iterations per sec'.format(t, (t+1)/(time.perf_counter()-start))) + sys.stdout.flush() + + if t % 20 == 0: + r = sum([(f * f * e).sum() for f, e in zip(e, epsilon)]) + print('E sum', r) + + # Save field slices + if (t % 20 == 0 and (max_t - t <= 1000 or t <= 2000)) or t == max_t-1: + print('saving E-field') + for j, f in enumerate(e): + output_file['/E{}_t{}'.format('xyz'[j], t)] = f[:, :, round(f.shape[2]/2)] + +if __name__ == '__main__': + main() diff --git a/fdfd_tools/fdtd.py b/fdfd_tools/fdtd.py new file mode 100644 index 0000000..c4aa897 --- /dev/null +++ b/fdfd_tools/fdtd.py @@ -0,0 +1,239 @@ +from typing import List, Callable, Tuple, Dict +import numpy + +from . import dx_lists_t, field_t + +__author__ = 'Jan Petykiewicz' + + +functional_matrix = Callable[[field_t], field_t] + + +def curl_h(dxes: dx_lists_t = None) -> functional_matrix: + """ + Curl operator for use with the H field. + + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :return: Function for taking the discretized curl of the H-field, F(H) -> curlH + """ + if dxes: + dxyz_b = numpy.meshgrid(*dxes[1], indexing='ij') + + def dh(f, ax): + return (f - numpy.roll(f, 1, axis=ax)) / dxyz_b[ax] + else: + def dh(f, ax): + return f - numpy.roll(f, 1, axis=ax) + + def ch_fun(h: field_t) -> field_t: + e = [dh(h[2], 1) - dh(h[1], 2), + dh(h[0], 2) - dh(h[2], 0), + dh(h[1], 0) - dh(h[0], 1)] + return e + + return ch_fun + + +def curl_e(dxes: dx_lists_t = None) -> functional_matrix: + """ + Curl operator for use with the E field. + + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header + :return: Function for taking the discretized curl of the E-field, F(E) -> curlE + """ + if dxes is not None: + dxyz_a = numpy.meshgrid(*dxes[0], indexing='ij') + + def de(f, ax): + return (numpy.roll(f, -1, axis=ax) - f) / dxyz_a[ax] + else: + def de(f, ax): + return numpy.roll(f, -1, axis=ax) - f + + def ce_fun(e: field_t) -> field_t: + h = [de(e[2], 1) - de(e[1], 2), + de(e[0], 2) - de(e[2], 0), + de(e[1], 0) - de(e[0], 1)] + return h + + return ce_fun + + +def maxwell_e(dt: float, dxes: dx_lists_t = None) -> functional_matrix: + curl_h_fun = curl_h(dxes) + + def me_fun(e: field_t, h: field_t, epsilon: field_t): + ch = curl_h_fun(h) + for ei, ci, epsi in zip(e, ch, epsilon): + ei += dt * ci / epsi + return e + + return me_fun + + +def maxwell_h(dt: float, dxes: dx_lists_t = None) -> functional_matrix: + curl_e_fun = curl_e(dxes) + + def mh_fun(e: field_t, h: field_t): + ce = curl_e_fun(e) + for hi, ci in zip(h, ce): + hi -= dt * ci + return h + + return mh_fun + + +def conducting_boundary(direction: int, + polarity: int + ) -> Tuple[functional_matrix, functional_matrix]: + dirs = [0, 1, 2] + if direction not in dirs: + raise Exception('Invalid direction: {}'.format(direction)) + dirs.remove(direction) + u, v = dirs + + if polarity < 0: + boundary_slice = [slice(None)] * 3 + shifted1_slice = [slice(None)] * 3 + boundary_slice[direction] = 0 + shifted1_slice[direction] = 1 + + def en(e: field_t): + e[direction][boundary_slice] = 0 + e[u][boundary_slice] = e[u][shifted1_slice] + e[v][boundary_slice] = e[v][shifted1_slice] + return e + + def hn(h: field_t): + h[direction][boundary_slice] = h[direction][shifted1_slice] + h[u][boundary_slice] = 0 + h[v][boundary_slice] = 0 + return h + + return en, hn + + elif polarity > 0: + boundary_slice = [slice(None)] * 3 + shifted1_slice = [slice(None)] * 3 + shifted2_slice = [slice(None)] * 3 + boundary_slice[direction] = -1 + shifted1_slice[direction] = -2 + shifted2_slice[direction] = -3 + + def ep(e: field_t): + e[direction][boundary_slice] = -e[direction][shifted2_slice] + e[direction][shifted1_slice] = 0 + e[u][boundary_slice] = e[u][shifted1_slice] + e[v][boundary_slice] = e[v][shifted1_slice] + return e + + def hp(h: field_t): + h[direction][boundary_slice] = h[direction][shifted1_slice] + h[u][boundary_slice] = -h[u][shifted2_slice] + h[u][shifted1_slice] = 0 + h[v][boundary_slice] = -h[v][shifted2_slice] + h[v][shifted1_slice] = 0 + return h + + return ep, hp + + else: + raise Exception('Bad polarity: {}'.format(polarity)) + + +def cpml(direction:int, + polarity: int, + dt: float, + epsilon: field_t, + thickness: int = 8, + epsilon_eff: float = 1, + dtype: numpy.dtype = numpy.float32, + ) -> Tuple[Callable, Callable, Dict[str, field_t]]: + + if direction not in range(3): + raise Exception('Invalid direction: {}'.format(direction)) + + if polarity not in (-1, 1): + raise Exception('Invalid polarity: {}'.format(polarity)) + + if thickness <= 2: + raise Exception('It would be wise to have a pml with 4+ cells of thickness') + + if epsilon_eff <= 0: + raise Exception('epsilon_eff must be positive') + + m = (3.5, 1) + sigma_max = 0.8 * (m[0] + 1) / numpy.sqrt(epsilon_eff) + alpha_max = 0 # TODO: Decide what to do about non-zero alpha + transverse = numpy.delete(range(3), direction) + u, v = transverse + + xe = numpy.arange(1, thickness+1, dtype=float) + xh = numpy.arange(1, thickness+1, dtype=float) + if polarity > 0: + xe -= 0.5 + elif polarity < 0: + xh -= 0.5 + xe = xe[::-1] + xh = xh[::-1] + else: + raise Exception('Bad polarity!') + + expand_slice = [None] * 3 + expand_slice[direction] = slice(None) + + def par(x): + sigma = ((x / thickness) ** m[0]) * sigma_max + alpha = ((1 - x / thickness) ** m[1]) * alpha_max + p0 = numpy.exp(-(sigma + alpha) * dt) + p1 = sigma / (sigma + alpha) * (p0 - 1) + return p0[expand_slice], p1[expand_slice] + + p0e, p1e = par(xe) + p0h, p1h = par(xh) + + region = [slice(None)] * 3 + if polarity < 0: + region[direction] = slice(None, thickness) + elif polarity > 0: + region[direction] = slice(-thickness, None) + else: + raise Exception('Bad polarity!') + + if direction == 1: + se = 1 + else: + se = -1 + + # TODO check if epsilon is uniform? + shape = list(epsilon[0].shape) + shape[direction] = thickness + psi_e = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)] + psi_h = [numpy.zeros(shape, dtype=dtype), numpy.zeros(shape, dtype=dtype)] + + fields = { + 'psi_e_u': psi_e[0], + 'psi_e_v': psi_e[1], + 'psi_h_u': psi_h[0], + 'psi_h_v': psi_h[1], + } + + def pml_e(e: field_t, h: field_t, epsilon: field_t) -> Tuple[field_t, field_t]: + psi_e[0] *= p0e + psi_e[0] += p1e * (h[v][region] - numpy.roll(h[v], 1, axis=direction)[region]) + psi_e[1] *= p0e + psi_e[1] += p1e * (h[u][region] - numpy.roll(h[u], 1, axis=direction)[region]) + e[u][region] += se * dt * psi_e[0] / epsilon[u][region] + e[v][region] -= se * dt * psi_e[1] / epsilon[v][region] + return e, h + + def pml_h(e: field_t, h: field_t) -> Tuple[field_t, field_t]: + psi_h[0] *= p0h + psi_h[0] += p1h * (numpy.roll(e[v], -1, axis=direction)[region] - e[v][region]) + psi_h[1] *= p0h + psi_h[1] += p1h * (numpy.roll(e[u], -1, axis=direction)[region] - e[u][region]) + h[u][region] -= se * dt * psi_h[0] + h[v][region] += se * dt * psi_h[1] + return e, h + + return pml_e, pml_h, fields