diff --git a/examples/fdfd1.py b/examples/fdfd1.py index 5596639..f983d37 100644 --- a/examples/fdfd1.py +++ b/examples/fdfd1.py @@ -1,4 +1,5 @@ import importlib +import logging import numpy from numpy.linalg import norm from matplotlib import pyplot, colors @@ -6,12 +7,14 @@ import logging import meanas from meanas import fdtd -from meanas.fdmath import vec, unvec +from meanas.fdmath import vec, unvec, fdfield_t from meanas.fdfd import waveguide_3d, functional, scpml, operators from meanas.fdfd.solvers import generic as generic_solver import gridlock +from matplotlib import pyplot + logging.basicConfig(level=logging.DEBUG) logging.getLogger('matplotlib').setLevel(logging.WARNING) diff --git a/examples/fdtd.py b/examples/fdtd.py index 8378b34..284ce07 100644 --- a/examples/fdtd.py +++ b/examples/fdtd.py @@ -1,18 +1,25 @@ """ -Example code for running an OpenCL FDTD simulation +Example code for running an FDTD simulation See main() for simulation setup. """ import sys import time +import copy import numpy import h5py +from numpy.linalg import norm from meanas import fdtd from meanas.fdtd import cpml_params, updates_with_cpml -from masque import Pattern, shapes +from meanas.fdtd.misc import gaussian_packet + +from meanas.fdfd.operators import e_full +from meanas.fdfd.scpml import stretch_with_scpml +from meanas.fdmath import vec +from masque import Pattern, Circle, Polygon import gridlock import pcgen @@ -41,50 +48,51 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern: `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, - } + default_args = { + '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 = 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 = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}' - pat.shapes += [shapes.Circle(radius=r, offset=(x, y), - dose=kwargs['hole_dose'], - layer=kwargs['hole_layer']) - for x, y, r in xyr] + #pat.name = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}' + pat.shapes[(kwargs['hole_layer'], 0)] += [ + Circle(radius=r, offset=(x, y)) + 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)] + pat.shapes[(kwargs['trench_layer'], 0)] += [ + Polygon.rectangle( + lx=(2 * maxes[0]), ly=kwargs['trench_width'], + offset=(0, s * (maxes[1] + a + kwargs['trench_width'] / 2)) + ) + for s in (-1, 1)] return pat def main(): dtype = numpy.float32 - max_t = 8000 # number of timesteps + max_t = 3600 # number of timesteps dx = 40 # discretization (nm/cell) pml_thickness = 8 # (number of cells) wl = 1550 # Excitation wavelength and fwhm - dwl = 200 + dwl = 100 # Device design parameters xy_size = numpy.array([10, 10]) @@ -107,69 +115,89 @@ def main(): # #### Create the grid, mask, and draw the device #### grid = gridlock.Grid(edge_coords) - epsilon = grid.allocate(n_air**2, dtype=dtype) - grid.draw_slab(epsilon, - surface_normal=2, - center=[0, 0, 0], - thickness=th, - eps=n_slab**2) + epsilon = grid.allocate(n_air ** 2, dtype=dtype) + grid.draw_slab( + epsilon, + slab = dict(axis='z', center=0, span=th), + foreground = n_slab ** 2, + ) + mask = perturbed_l3(a, r) + grid.draw_polygons( + epsilon, + slab = dict(axis='z', center=0, span=2 * th), + foreground = n_air ** 2, + offset2d = (0, 0), + polygons = mask.as_polygons(library=None), + ) - grid.draw_polygons(epsilon, - surface_normal=2, - center=[0, 0, 0], - thickness=2 * th, - eps=n_air**2, - polygons=mask.as_polygons()) + print(f'{grid.shape=}') - print(grid.shape) - - 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)] + dt = dx * 0.99 / numpy.sqrt(3) + ee = numpy.zeros_like(epsilon, dtype=dtype) + hh = numpy.zeros_like(epsilon, dtype=dtype) dxes = [grid.dxyz, grid.autoshifted_dxyz()] # PMLs in every direction - pml_params = [[cpml_params(axis=dd, polarity=pp, dt=dt, - thickness=pml_thickness, epsilon_eff=1.0**2) - for pp in (-1, +1)] - for dd in range(3)] - update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, - dxes=dxes, epsilon=epsilon) + pml_params = [ + [cpml_params(axis=dd, polarity=pp, dt=dt, thickness=pml_thickness, epsilon_eff=n_air ** 2) + for pp in (-1, +1)] + for dd in range(3)] + update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon) + + # sample_interval = numpy.floor(1 / (2 * 1 / wl * dt)).astype(int) + # print(f'Save time interval would be {sample_interval} * dt = {sample_interval * dt:3g}') + # 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) + source_phasor, _delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5) + aa, cc, ss = source_phasor(numpy.arange(max_t)) + srca_real = aa * cc + src_maxt = numpy.argwhere(numpy.diff(aa < 1e-5))[-1] + assert aa[src_maxt - 1] >= 1e-5 + phasor_norm = dt / (aa * cc * cc).sum() - def field_source(i): - t0 = i * dt - delay - return numpy.sin(w * t0) * numpy.exp(-alpha * t0**2) + Jph = numpy.zeros_like(epsilon, dtype=complex) + Jph[1, *(grid.shape // 2)] = epsilon[1, *(grid.shape // 2)] + Eph = numpy.zeros_like(Jph) # #### Run a bunch of iterations #### output_file = h5py.File('simulation_output.h5', 'w') start = time.perf_counter() - for t in range(max_t): - update_E(e, h, epsilon) + for tt in range(max_t): + update_E(ee, hh, epsilon) - e[1][tuple(grid.shape//2)] += field_source(t) - update_H(e, h) + if tt < src_maxt: + ee[1, *(grid.shape // 2)] -= srca_real[tt] + update_H(ee, hh) - avg_rate = (t + 1)/(time.perf_counter() - start)) - print(f'iteration {t}: average {avg_rate} iterations per sec') + avg_rate = (tt + 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) + if tt % 200 == 0: + print(f'iteration {tt}: average {avg_rate} iterations per sec') + E_energy_sum = (ee * ee * epsilon).sum() + print(f'{E_energy_sum=}') # 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 (tt % 20 == 0 and (max_t - tt <= 1000 or tt <= 2000)) or tt == max_t - 1: + print(f'saving E-field at iteration {tt}') + output_file[f'/E_t{tt}'] = ee[:, :, :, ee.shape[3] // 2] + + Eph += (cc[tt] - 1j * ss[tt]) * phasor_norm * ee + + omega = 2 * numpy.pi / wl + Eph *= numpy.exp(-1j * dt / 2 * omega) + b = -1j * omega * Jph + dxes_fdfd = copy.deepcopy(dxes) + for pp in (-1, +1): + for dd in range(3): + stretch_with_scpml(dxes_fdfd, axis=dd, polarity=pp, omega=omega, epsilon_effective=n_air ** 2, thickness=pml_thickness) + A = e_full(omega=omega, dxes=dxes, epsilon=epsilon) + residual = norm(A @ vec(ee) - vec(b)) / norm(vec(b)) + print(f'FDFD residual is {residual}') + if __name__ == '__main__': main()