jan/meanas
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

update for new Gridlock

Browse Source
This commit is contained in:
Jan Petykiewicz 2022-10-04 12:29:11 -07:00
parent 831c39246a
commit eedcc7b919
4 changed files with 39 additions and 33 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
examples/bloch.py
View File

@ -47,10 +47,10 @@ g = gridlock.Grid([numpy.arange(-x_period/2, x_period/2, dx),
numpy.arange(-1000, 1000, dx), numpy.arange(-1000, 1000, dx),
numpy.arange(-1000, 1000, dx)], numpy.arange(-1000, 1000, dx)],
shifts=numpy.array([[0,0,0]]), shifts=numpy.array([[0,0,0]]),
initial=1.445**2,
periodic=True) periodic=True)
gdata = g.allocate(1.445**2)
g.draw_cuboid([0,0,0], [200e8, 220, 220], eps=3.47**2) g.draw_cuboid(gdata, [0,0,0], [200e8, 220, 220], foreground=3.47**2)
#x_period = y_period = z_period = 13000 #x_period = y_period = z_period = 13000
#g = gridlock.Grid([numpy.arange(3), ]*3, #g = gridlock.Grid([numpy.arange(3), ]*3,
@ -60,9 +60,9 @@ g.draw_cuboid([0,0,0], [200e8, 220, 220], eps=3.47**2)
g2 = g.copy() g2 = g.copy()
g2.shifts = numpy.zeros((6,3)) g2.shifts = numpy.zeros((6,3))
g2.grids = [numpy.zeros(g.shape) for _ in range(6)] g2data = g2.allocate(0)
epsilon = [g.grids[0],] * 3 epsilon = [gdata[0],] * 3
reciprocal_lattice = numpy.diag(1000/numpy.array([x_period, y_period, z_period])) #cols are vectors reciprocal_lattice = numpy.diag(1000/numpy.array([x_period, y_period, z_period])) #cols are vectors
pyfftw_load_wisdom(WISDOM_FILEPATH) pyfftw_load_wisdom(WISDOM_FILEPATH)
@ -93,8 +93,8 @@ for k0x in [.25]:
z = 0 z = 0
e = v2e(v[0]) e = v2e(v[0])
for i in range(3): for i in range(3):
g2.grids[i] += numpy.real(e[i]) g2data[i] += numpy.real(e[i])
g2.grids[i+3] += numpy.imag(e[i]) g2data[i+3] += numpy.imag(e[i])
f = numpy.sqrt(numpy.real(numpy.abs(n))) # TODO f = numpy.sqrt(numpy.real(numpy.abs(n))) # TODO
print('k0x = {:3g}\n eigval = {}\n f = {}\n'.format(k0x, n, f)) print('k0x = {:3g}\n eigval = {}\n f = {}\n'.format(k0x, n, f))

42
examples/fdfd.py
View File

@ -44,14 +44,17 @@ def test0(solver=generic_solver):
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords] edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
# #### Create the grid, mask, and draw the device #### # #### Create the grid, mask, and draw the device ####
grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3) grid = gridlock.Grid(edge_coords)
grid.draw_cylinder(surface_normal=gridlock.Direction.z, epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
grid.draw_cylinder(epsilon,
surface_normal=2,
center=center, center=center,
radius=max(radii), radius=max(radii),
thickness=th, thickness=th,
eps=n_ring**2, eps=n_ring**2,
num_points=24) num_points=24)
grid.draw_cylinder(surface_normal=gridlock.Direction.z, grid.draw_cylinder(epsilon,
surface_normal=2,
center=center, center=center,
radius=min(radii), radius=min(radii),
thickness=th*1.1, thickness=th*1.1,
@ -64,7 +67,7 @@ def test0(solver=generic_solver):
dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega, dxes = meanas.fdfd.scpml.stretch_with_scpml(dxes, axis=a, polarity=p, omega=omega,
thickness=pml_thickness) thickness=pml_thickness)
J = [numpy.zeros_like(grid.grids[0], dtype=complex) for _ in range(3)] J = [numpy.zeros_like(epsilon[0], dtype=complex) for _ in range(3)]
J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1 J[1][15, grid.shape[1]//2, grid.shape[2]//2] = 1
@ -74,11 +77,11 @@ def test0(solver=generic_solver):
sim_args = { sim_args = {
'omega': omega, 'omega': omega,
'dxes': dxes, 'dxes': dxes,
'epsilon': vec(grid.grids), 'epsilon': vec(epsilon),
} }
x = solver(J=vec(J), **sim_args) x = solver(J=vec(J), **sim_args)
A = operators.e_full(omega, dxes, vec(grid.grids)).tocsr() A = operators.e_full(omega, dxes, vec(epsilon)).tocsr()
b = -1j * omega * vec(J) b = -1j * omega * vec(J)
print('Norm of the residual is ', norm(A @ x - b)) print('Norm of the residual is ', norm(A @ x - b))
@ -117,8 +120,9 @@ def test1(solver=generic_solver):
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords] edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
# #### Create the grid and draw the device #### # #### Create the grid and draw the device ####
grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3) grid = gridlock.Grid(edge_coords)
grid.draw_cuboid(center=center, dimensions=[8e3, w, th], eps=n_wg**2) epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], eps=n_wg**2)
dxes = [grid.dxyz, grid.autoshifted_dxyz()] dxes = [grid.dxyz, grid.autoshifted_dxyz()]
for a in (0, 1, 2): for a in (0, 1, 2):
@ -139,18 +143,18 @@ def test1(solver=generic_solver):
'polarity': +1, 'polarity': +1,
} }
wg_results = waveguide_3d.solve_mode(mode_number=0, omega=omega, epsilon=grid.grids, **wg_args) wg_results = waveguide_3d.solve_mode(mode_number=0, omega=omega, epsilon=epsilon, **wg_args)
J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'], J = waveguide_3d.compute_source(E=wg_results['E'], wavenumber=wg_results['wavenumber'],
omega=omega, epsilon=grid.grids, **wg_args) omega=omega, epsilon=epsilon, **wg_args)
e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args) e_overlap = waveguide_3d.compute_overlap_e(E=wg_results['E'], wavenumber=wg_results['wavenumber'], **wg_args)
pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) pecg = numpy.zeros_like(epsilon)
# pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) # pecg.draw_cuboid(pecg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
# pecg.visualize_isosurface() # pecg.visualize_isosurface(pecg)
pmcg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) pmcg = numpy.zeros_like(epsilon)
# pmcg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) # grid.draw_cuboid(pmcg, center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
# pmcg.visualize_isosurface() # grid.visualize_isosurface(pmcg)
def pcolor(v): def pcolor(v):
vmax = numpy.max(numpy.abs(v)) vmax = numpy.max(numpy.abs(v))
@ -171,9 +175,9 @@ def test1(solver=generic_solver):
sim_args = { sim_args = {
'omega': omega, 'omega': omega,
'dxes': dxes, 'dxes': dxes,
'epsilon': vec(grid.grids), 'epsilon': vec(epsilon),
'pec': vec(pecg.grids), 'pec': vec(pecg),
'pmc': vec(pmcg.grids), 'pmc': vec(pmcg),
} }
x = solver(J=vec(J), **sim_args) x = solver(J=vec(J), **sim_args)

11
examples/fdtd.py
View File

@ -105,22 +105,23 @@ def main():
edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords] edge_coords = [numpy.hstack((-h[::-1], h)) for h in half_edge_coords]
# #### Create the grid, mask, and draw the device #### # #### Create the grid, mask, and draw the device ####
grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3) grid = gridlock.Grid(edge_coords)
grid.draw_slab(surface_normal=gridlock.Direction.z, epsilon = grid.allocate(n_air**2, dtype=dtype)
grid.draw_slab(epsilon,
surface_normal=2,
center=[0, 0, 0], center=[0, 0, 0],
thickness=th, thickness=th,
eps=n_slab**2) eps=n_slab**2)
mask = perturbed_l3(a, r) mask = perturbed_l3(a, r)
grid.draw_polygons(surface_normal=gridlock.Direction.z, grid.draw_polygons(epsilon,
surface_normal=2,
center=[0, 0, 0], center=[0, 0, 0],
thickness=2 * th, thickness=2 * th,
eps=n_air**2, eps=n_air**2,
polygons=mask.as_polygons()) polygons=mask.as_polygons())
print(grid.shape) print(grid.shape)
# #### Create the simulation grid ####
epsilon = [eps.astype(dtype) for eps in grid.grids]
dt = .99/numpy.sqrt(3) dt = .99/numpy.sqrt(3)
e = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)] e = [numpy.zeros_like(epsilon[0], dtype=dtype) for _ in range(3)]

7
examples/tcyl.py
View File

@ -42,8 +42,9 @@ def test1(solver=generic_solver):
edge_coords[0] = numpy.array([-dx, dx]) edge_coords[0] = numpy.array([-dx, dx])
# #### Create the grid and draw the device #### # #### Create the grid and draw the device ####
grid = gridlock.Grid(edge_coords, initial=n_air**2, num_grids=3) grid = gridlock.Grid(edge_coords)
grid.draw_cuboid(center=center, dimensions=[8e3, w, th], eps=n_wg**2) epsilon = grid.allocate(n_air**2, dtype=numpy.float32)
grid.draw_cuboid(epsilon, center=center, dimensions=[8e3, w, th], eps=n_wg**2)
dxes = [grid.dxyz, grid.autoshifted_dxyz()] dxes = [grid.dxyz, grid.autoshifted_dxyz()]
for a in (1, 2): for a in (1, 2):
@ -54,7 +55,7 @@ def test1(solver=generic_solver):
wg_args = { wg_args = {
'omega': omega, 'omega': omega,
'dxes': [(d[1], d[2]) for d in dxes], 'dxes': [(d[1], d[2]) for d in dxes],
'epsilon': vec(g.transpose([1, 2, 0]) for g in grid.grids), 'epsilon': vec(g.transpose([1, 2, 0]) for g in epsilon),
'r0': r0, 'r0': r0,
} }