From fb3bef23bfbccd3c6a103163e8a5309633cbad63 Mon Sep 17 00:00:00 2001 From: jan Date: Wed, 10 Dec 2025 21:14:34 -0800 Subject: [PATCH] [examples/fdfd] split fdfd example into two files --- examples/fdfd0.py | 103 ++++++++++++++++++++++ examples/{fdfd.py => fdfd1.py} | 153 +++++++-------------------------- 2 files changed, 135 insertions(+), 121 deletions(-) create mode 100644 examples/fdfd0.py rename examples/{fdfd.py => fdfd1.py} (54%) diff --git a/examples/fdfd0.py b/examples/fdfd0.py new file mode 100644 index 0000000..a6227c5 --- /dev/null +++ b/examples/fdfd0.py @@ -0,0 +1,103 @@ +import numpy +from numpy.linalg import norm +from matplotlib import pyplot, colors +import logging + +import meanas +from meanas import fdtd +from meanas.fdmath import vec, unvec +from meanas.fdfd import waveguide_3d, functional, scpml, operators +from meanas.fdfd.solvers import generic as generic_solver + +import gridlock + + +logging.basicConfig(level=logging.DEBUG) +logging.getLogger('matplotlib').setLevel(logging.WARNING) + +__author__ = 'Jan Petykiewicz' + + +def pcolor(ax, v) -> None: + mappable = ax.pcolor(v, cmap='seismic', norm=colors.CenteredNorm()) + ax.axis('equal') + ax.get_figure().colorbar(mappable) + + +def test0(solver=generic_solver): + 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 #### + grid = gridlock.Grid(edge_coords) + epsilon = grid.allocate(n_air**2, dtype=numpy.float32) + grid.draw_cylinder( + epsilon, + h = dict(axis='z', center=center[2], span=th), + radius = max(radii), + center2d = center[:2], + foreground = n_ring ** 2, + num_points = 24, + ) + grid.draw_cylinder( + epsilon, + h = dict(axis='z', center=center[2], span=th * 1.1), + radius = min(radii), + center2d = center[:2], + foreground = n_air ** 2, + num_points = 24, + ) + + dxes = [grid.dxyz, grid.autoshifted_dxyz()] + for a in (0, 1, 2): + for p in (-1, 1): + dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, + thickness=pml_thickness) + + J = [numpy.zeros_like(epsilon[0], dtype=complex) for _ in range(3)] + J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1 + + + # + # Solve! + # + sim_args = dict( + omega = omega, + dxes = dxes, + epsilon = vec(epsilon), + ) + x = solver(J=vec(J), **sim_args) + + A = operators.e_full(omega, dxes, vec(epsilon)).tocsr() + b = -1j * omega * vec(J) + print('Norm of the residual is ', norm(A @ x - b) / norm(b)) + + E = unvec(x, grid.shape) + + # + # Plot results + # + grid.visualize_slice(E.real, plane=dict(z=0), which_shifts=1, pcolormesh_args=dict(norm=colors.CenteredNorm(), cmap='bwr')) + + +if __name__ == '__main__': + test0() diff --git a/examples/fdfd.py b/examples/fdfd1.py similarity index 54% rename from examples/fdfd.py rename to examples/fdfd1.py index e768ba7..5596639 100644 --- a/examples/fdfd.py +++ b/examples/fdfd1.py @@ -1,6 +1,8 @@ import importlib import numpy from numpy.linalg import norm +from matplotlib import pyplot, colors +import logging import meanas from meanas import fdtd @@ -10,9 +12,6 @@ from meanas.fdfd.solvers import generic as generic_solver import gridlock -from matplotlib import pyplot - -import logging logging.basicConfig(level=logging.DEBUG) logging.getLogger('matplotlib').setLevel(logging.WARNING) @@ -20,86 +19,6 @@ logging.getLogger('matplotlib').setLevel(logging.WARNING) __author__ = 'Jan Petykiewicz' -def test0(solver=generic_solver): - 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 #### - grid = gridlock.Grid(edge_coords) - epsilon = grid.allocate(n_air**2, dtype=numpy.float32) - grid.draw_cylinder( - epsilon, - surface_normal=2, - center=center, - radius=max(radii), - thickness=th, - foreground=n_ring**2, - num_points=24, - ) - grid.draw_cylinder( - epsilon, - surface_normal=2, - center=center, - radius=min(radii), - thickness=th*1.1, - foreground=n_air ** 2, - num_points=24, - ) - - dxes = [grid.dxyz, grid.autoshifted_dxyz()] - for a in (0, 1, 2): - for p in (-1, 1): - dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, - thickness=pml_thickness) - - J = [numpy.zeros_like(epsilon[0], dtype=complex) for _ in range(3)] - J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1 - - - # - # Solve! - # - sim_args = { - 'omega': omega, - 'dxes': dxes, - 'epsilon': vec(epsilon), - } - x = solver(J=vec(J), **sim_args) - - A = operators.e_full(omega, dxes, vec(epsilon)).tocsr() - b = -1j * omega * vec(J) - print('Norm of the residual is ', norm(A @ x - b)) - - E = unvec(x, grid.shape) - - # - # Plot results - # - pyplot.figure() - pyplot.pcolor(numpy.real(E[1][:, :, grid.shape[2]//2]), cmap='seismic') - pyplot.axis('equal') - pyplot.show() - - def test1(solver=generic_solver): dx = 40 # discretization (nm/cell) pml_thickness = 10 # (number of cells) @@ -126,7 +45,7 @@ def test1(solver=generic_solver): # #### Create the grid and draw the device #### grid = gridlock.Grid(edge_coords) epsilon = grid.allocate(n_air**2, dtype=numpy.float32) - grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], foreground=n_wg**2) + grid.draw_cuboid(epsilon, x=dict(center=0, span=8e3), y=dict(center=0, span=w), z=dict(center=0, span=th), foreground=n_wg**2) dxes = [grid.dxyz, grid.autoshifted_dxyz()] for a in (0, 1, 2): @@ -160,17 +79,9 @@ def test1(solver=generic_solver): # grid.draw_cuboid(pmcg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) # grid.visualize_isosurface(pmcg) - def pcolor(v) -> None: - vmax = numpy.max(numpy.abs(v)) - pyplot.pcolor(v, cmap='seismic', vmin=-vmax, vmax=vmax) - pyplot.axis('equal') - pyplot.colorbar() - - 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) + grid.visualize_slice(J.imag, plane=dict(y=6*dx), which_shifts=1, pcolormesh_args=dict(norm=colors.CenteredNorm(), cmap='bwr')) + fig, ax = pyplot.subplots() + ax.pcolormesh((numpy.abs(J).sum(axis=2).sum(axis=0) > 0).astype(float).T, cmap='hot') pyplot.show(block=True) # @@ -196,16 +107,14 @@ def test1(solver=generic_solver): # Plot results # center = grid.pos2ind([0, 0, 0], None).astype(int) - pyplot.figure() - pyplot.subplot(2, 2, 1) - pcolor(numpy.real(E[1][center[0], :, :]).T) - pyplot.subplot(2, 2, 2) - pyplot.plot(numpy.log10(numpy.abs(E[1][:, center[1], center[2]]) + 1e-10)) - pyplot.grid(alpha=0.6) - pyplot.ylabel('log10 of field') - pyplot.subplot(2, 2, 3) - pcolor(numpy.real(E[1][:, :, center[2]]).T) - pyplot.subplot(2, 2, 4) + fig, axes = pyplot.subplots(2, 2) + grid.visualize_slice(E.real, plane=dict(x=0), which_shifts=1, ax=axes[0, 0], finalize=False, pcolormesh_args=dict(norm=colors.CenteredNorm(), cmap='bwr')) + grid.visualize_slice(E.real, plane=dict(z=0), which_shifts=1, ax=axes[0, 1], finalize=False, pcolormesh_args=dict(norm=colors.CenteredNorm(), cmap='bwr')) +# pcolor(axes[0, 0], numpy.real(E[1][center[0], :, :]).T) +# pcolor(axes[0, 1], numpy.real(E[1][:, :, center[2]]).T) + axes[1, 0].plot(numpy.log10(numpy.abs(E[1][:, center[1], center[2]]) + 1e-10)) + axes[1, 0].grid(alpha=0.6) + axes[1, 0].set_ylabel('log10 of field') def poyntings(E): H = functional.e2h(omega, dxes)(E) @@ -219,34 +128,35 @@ def test1(solver=generic_solver): return s0, s1, s2 s0x, s1x, s2x = poyntings(E) - 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) - pyplot.legend() - pyplot.show() + ax = axes[1, 1] + ax.plot(s0x[0].sum(axis=2).sum(axis=1), label='s0', marker='.') + ax.plot(s1x[0].sum(axis=2).sum(axis=1), label='s1', marker='.') + ax.plot(s2x[0].sum(axis=2).sum(axis=1), label='s2', marker='.') + ax.plot(E[1][:, center[1], center[2]].real.T, label='Ey', marker='x') + ax.grid(alpha=0.6) + ax.legend() + + p_in = (-E * J.conj()).sum() / 2 * (dx * dx * dx) + print(f'{p_in=}') q = [] for i in range(-5, 30): e_ovl_rolled = numpy.roll(e_overlap, i, axis=1) - q += [numpy.abs(vec(E) @ vec(e_ovl_rolled).conj())] - pyplot.figure() - pyplot.plot(q, marker='.') - pyplot.grid(alpha=0.6) - pyplot.title('Overlap with mode') - pyplot.show() + q += [numpy.abs(vec(E).conj() @ vec(e_ovl_rolled))] + fig, ax = pyplot.subplots() + ax.plot(q, marker='.') + ax.grid(alpha=0.6) + ax.set_title('Overlap with mode') print('Average overlap with mode:', sum(q[8:32])/len(q[8:32])) + pyplot.show(block=True) + def module_available(name): return importlib.util.find_spec(name) is not None if __name__ == '__main__': - #test0() -# test1() - if module_available('opencl_fdfd'): from opencl_fdfd import cg_solver as opencl_solver test1(opencl_solver) @@ -257,3 +167,4 @@ if __name__ == '__main__': # test1(magma_solver) else: test1() +