misc example updates

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)

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,8 +48,7 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
         `masque.Pattern` object containing the L3 design
     """
 
-    default_args = {'hole_dose':    1,
+    default_args = {
-                    'trench_dose':  1,
         'hole_layer':   0,
         'trench_layer': 1,
         'shifts_a':     (0.15, 0, 0.075),
@@ -53,38 +59,40 @@ def perturbed_l3(a: float, radius: float, **kwargs) -> Pattern:
         }
     kwargs = {**default_args, **kwargs}
 
-    xyr = pcgen.l3_shift_perturbed_defect(mirror_dims=kwargs['xy_size'],
+    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'])
+        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.name = f'L3p-a{a:g}r{radius:g}rp{kwargs["perturbed_radius"]:g}'
-    pat.shapes += [shapes.Circle(radius=r, offset=(x, y),
+    pat.shapes[(kwargs['hole_layer'], 0)] += [
-                                 dose=kwargs['hole_dose'],
+        Circle(radius=r, offset=(x, y))
-                                 layer=kwargs['hole_layer'])
         for x, y, r in xyr]
 
     maxes = numpy.max(numpy.fabs(xyr), axis=0)
-    pat.shapes += [shapes.Polygon.rectangle(
+    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)),
+            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
+    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)
+    epsilon = grid.allocate(n_air ** 2, dtype=dtype)
-    grid.draw_slab(epsilon,
+    grid.draw_slab(
-                   surface_normal=2,
+        epsilon,
-                   center=[0, 0, 0],
+        slab = dict(axis='z', center=0, span=th),
-                   thickness=th,
+        foreground = n_slab ** 2,
-                   eps=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,
+    print(f'{grid.shape=}')
-                       surface_normal=2,
-                       center=[0, 0, 0],
-                       thickness=2 * th,
-                       eps=n_air**2,
-                       polygons=mask.as_polygons())
 
-    print(grid.shape)
+    dt = dx * 0.99 / numpy.sqrt(3)
-
+    ee = numpy.zeros_like(epsilon, dtype=dtype)
-    dt = .99/numpy.sqrt(3)
+    hh = numpy.zeros_like(epsilon, dtype=dtype)
-    e = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)]
-    h = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)]
 
     dxes = [grid.dxyz, grid.autoshifted_dxyz()]
 
     # PMLs in every direction
-    pml_params = [[cpml_params(axis=dd, polarity=pp, dt=dt,
+    pml_params = [
-                               thickness=pml_thickness, epsilon_eff=1.0**2)
+        [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,
+    update_E, update_H = updates_with_cpml(cpml_params=pml_params, dt=dt, dxes=dxes, epsilon=epsilon)
-                                           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
+    source_phasor, _delay = gaussian_packet(wl=wl, dwl=100, dt=dt, turn_on=1e-5)
-    fwhm = dwl * w * w / (2 * numpy.pi * dx)
+    aa, cc, ss = source_phasor(numpy.arange(max_t))
-    alpha = (fwhm ** 2) / 8 * numpy.log(2)
+    srca_real = aa * cc
-    delay = 7/numpy.sqrt(2 * alpha)
+    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):
+    Jph = numpy.zeros_like(epsilon, dtype=complex)
-        t0 = i * dt - delay
+    Jph[1, *(grid.shape // 2)] = epsilon[1, *(grid.shape // 2)]
-        return numpy.sin(w * t0) * numpy.exp(-alpha * t0**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):
+    for tt in range(max_t):
-        update_E(e, h, epsilon)
+        update_E(ee, hh, epsilon)
 
-        e[1][tuple(grid.shape//2)] += field_source(t)
+        if tt < src_maxt:
-        update_H(e, h)
+            ee[1, *(grid.shape // 2)] -= srca_real[tt]
+        update_H(ee, hh)
 
-        avg_rate = (t + 1)/(time.perf_counter() - start))
+        avg_rate = (tt + 1) / (time.perf_counter() - start)
-        print(f'iteration {t}: average {avg_rate} iterations per sec')
         sys.stdout.flush()
 
-        if t % 20 == 0:
+        if tt % 200 == 0:
-            r = sum([(f * f * e).sum() for f, e in zip(e, epsilon)])
+            print(f'iteration {tt}: average {avg_rate} iterations per sec')
-            print('E sum', r)
+            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:
+        if (tt % 20 == 0 and (max_t - tt <= 1000 or tt <= 2000)) or tt == max_t - 1:
-            print('saving E-field')
+            print(f'saving E-field at iteration {tt}')
-            for j, f in enumerate(e):
+            output_file[f'/E_t{tt}'] = ee[:, :, :, ee.shape[3] // 2]
-                output_file['/E{}_t{}'.format('xyz'[j], t)] = f[:, :, round(f.shape[2]/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()