jan/fdfd_tools
2
2
Fork 1
Code Issues Pull Requests Releases Wiki Activity

in progress

Browse Source
This commit is contained in:
Jan Petykiewicz 2019-08-30 22:03:54 -07:00
parent a4c2239ad9
commit 8eac9df76e
4 changed files with 46 additions and 81 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

25
examples/fdfd.py
View File

@ -137,7 +137,9 @@ def test1(solver=generic_solver):
wg_results = waveguide_mode.solve_waveguide_mode(mode_number=0, **wg_args) wg_results = waveguide_mode.solve_waveguide_mode(mode_number=0, **wg_args)
J = waveguide_mode.compute_source(**wg_args, E=wg_results['E'], wavenumber=wg_results['wavenumber']) J = waveguide_mode.compute_source(**wg_args, E=wg_results['E'], wavenumber=wg_results['wavenumber'])
H_overlap = waveguide_mode.compute_overlap_e(**wg_args, **wg_results) H_overlap, slices = waveguide_mode.compute_overlap_ce(E=wg_results['E'], wavenumber=wg_results['wavenumber'],
dxes=dxes, axis=src_axis, polarity=wg_args['polarity'],
slices=wg_args['slices'])
pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3) pecg = gridlock.Grid(edge_coords, initial=0.0, num_grids=3)
# pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1) # pecg.draw_cuboid(center=[700, 0, 0], dimensions=[80, 1e8, 1e8], eps=1)
@ -153,6 +155,13 @@ def test1(solver=generic_solver):
pyplot.axis('equal') pyplot.axis('equal')
pyplot.colorbar() 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)
pyplot.show(block=True)
''' '''
Solve! Solve!
''' '''
@ -186,12 +195,14 @@ def test1(solver=generic_solver):
pyplot.subplot(2, 2, 4) pyplot.subplot(2, 2, 4)
def poyntings(E): def poyntings(E):
e = vec(E) H = functional.e2h(omega, dxes)(E)
h = operators.e2h(omega, dxes) @ e poynting = 0.5 * fdtd.poynting(e=E, h=H.conj()) * dx * dx
cross1 = operators.poynting_e_cross(e, dxes) @ h.conj() cross1 = operators.poynting_e_cross(vec(E), dxes) @ vec(H).conj()
cross2 = operators.poynting_h_cross(h.conj(), dxes) @ e # cross2 = operators.poynting_h_cross(h.conj(), dxes) @ e
s1 = unvec(0.5 * numpy.real(cross1), grid.shape) s1 = unvec(0.5 * numpy.real(cross1), grid.shape)
s2 = unvec(0.5 * numpy.real(-cross2), grid.shape) # s2 = unvec(0.5 * numpy.real(-cross2), grid.shape)
s2 = poynting.real
# s2 = poynting.imag
return s1, s2 return s1, s2
s1x, s2x = poyntings(E) s1x, s2x = poyntings(E)
@ -202,7 +213,7 @@ def test1(solver=generic_solver):
q = [] q = []
for i in range(-5, 30): for i in range(-5, 30):
H_rolled = [numpy.roll(h, i, axis=0) for h in H_overlap] H_rolled = [numpy.roll(h, i, axis=0) for h in H_overlap]
q += [numpy.abs(vec(E) @ vec(H_rolled))] q += [numpy.abs(vec(E) @ vec(H_rolled).conj())]
pyplot.figure() pyplot.figure()
pyplot.plot(q) pyplot.plot(q)
pyplot.title('Overlap with mode') pyplot.title('Overlap with mode')

20
meanas/fdfd/operators.py
View File

@ -453,8 +453,7 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
""" """
shape = [len(dx) for dx in dxes[0]] shape = [len(dx) for dx in dxes[0]]
fx, fy, fz = [avgf(i, shape) for i in range(3)] bx, by, bz = [rotation(i, shape, -1) for i in range(3)]
bx, by, bz = [avgb(i, shape) for i in range(3)]
dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')] dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')]
dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C')) dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C'))
@ -463,12 +462,11 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
Ex, Ey, Ez = [sparse.diags(ei * da) for ei, da in zip(numpy.split(e, 3), dxag)] Ex, Ey, Ez = [sparse.diags(ei * da) for ei, da in zip(numpy.split(e, 3), dxag)]
n = numpy.prod(shape) n = numpy.prod(shape)
zero = sparse.csr_matrix((n, n))
P = sparse.bmat( P = sparse.bmat(
[[ zero, -fx @ Ez @ bz @ dbgy, fx @ Ey @ by @ dbgz], [[ None, -bx @ Ez @ dbgy, bx @ Ey @ dbgz],
[ fy @ Ez @ bz @ dbgx, zero, -fy @ Ex @ bx @ dbgz], [ by @ Ez @ dbgx, None, -by @ Ex @ dbgz],
[-fz @ Ey @ by @ dbgx, fz @ Ex @ bx @ dbgy, zero]]) [-bz @ Ey @ dbgx, bz @ Ex @ dbgy, None]])
return P return P
@ -482,7 +480,7 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix:
""" """
shape = [len(dx) for dx in dxes[0]] shape = [len(dx) for dx in dxes[0]]
fx, fy, fz = [avgf(i, shape) for i in range(3)] fx, fy, fz = [avgf(i, shape) for i in range(3)] #TODO
bx, by, bz = [avgb(i, shape) for i in range(3)] bx, by, bz = [avgb(i, shape) for i in range(3)]
dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')] dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')]
@ -545,4 +543,12 @@ def e_boundary_source(mask: vfield_t,
r3 = sparse.block_diag((r, r, r)) r3 = sparse.block_diag((r, r, r))
jmask = numpy.logical_or(jmask, numpy.abs(r3 @ mask)) jmask = numpy.logical_or(jmask, numpy.abs(r3 @ mask))
# jmask = ((numpy.roll(mask, -1, axis=0) != mask) |
# (numpy.roll(mask, +1, axis=0) != mask) |
# (numpy.roll(mask, -1, axis=1) != mask) |
# (numpy.roll(mask, +1, axis=1) != mask) |
# (numpy.roll(mask, -1, axis=2) != mask) |
# (numpy.roll(mask, +1, axis=2) != mask))
return sparse.diags(jmask.astype(int)) @ full return sparse.diags(jmask.astype(int)) @ full

2
meanas/fdfd/waveguide.py
View File

@ -186,7 +186,7 @@ def _normalized_fields(e: numpy.ndarray,
norm_amplitude = 1 / numpy.sqrt(P) norm_amplitude = 1 / numpy.sqrt(P)
norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric norm_angle = -numpy.angle(e[energy.argmax()]) # Will randomly add a negative sign when mode is symmetric
# Try to break symmetry to assign a consistent sign [experimental] # Try to break symmetry to assign a consistent sign [experimental TODO]
E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape) E_weighted = unvec(e * energy * numpy.exp(1j * norm_angle), shape)
sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum())

80
meanas/fdfd/waveguide_mode.py
View File

@ -112,6 +112,8 @@ def solve_waveguide_mode(mode_number: int,
Apply corrections and expand to 3D Apply corrections and expand to 3D
''' '''
# Correct wavenumber to account for numerical dispersion. # Correct wavenumber to account for numerical dispersion.
print(fields_2d['wavenumber'] / (2/dx_prop * numpy.arcsin(fields_2d['wavenumber'] * dx_prop/2)))
print(fields_2d['wavenumber'].real / (2/dx_prop * numpy.arcsin(fields_2d['wavenumber'].real * dx_prop/2)))
fields_2d['wavenumber'] = 2/dx_prop * numpy.arcsin(fields_2d['wavenumber'] * dx_prop/2) fields_2d['wavenumber'] = 2/dx_prop * numpy.arcsin(fields_2d['wavenumber'] * dx_prop/2)
# Adjust for propagation direction # Adjust for propagation direction
@ -179,20 +181,16 @@ def compute_source(E: field_t,
def compute_overlap_e(E: field_t, def compute_overlap_e(E: field_t,
H: field_t,
wavenumber: complex, wavenumber: complex,
omega: complex,
dxes: dx_lists_t, dxes: dx_lists_t,
axis: int, axis: int,
polarity: int, polarity: int,
slices: List[slice], slices: List[slice],
epsilon: field_t, # TODO unused?? ) -> field_t: # TODO DOCS
mu: field_t = None,
) -> field_t:
""" """
Given an eigenmode obtained by solve_waveguide_mode, calculates overlap_e for the Given an eigenmode obtained by solve_waveguide_mode, calculates overlap_e for the
mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn) mode orthogonality relation Integrate(((E x H_mode) + (E_mode x H)) dot dn)
[assumes reflection symmetry]. [assumes reflection symmetry].i
overlap_e makes use of the e2h operator to collapse the above expression into overlap_e makes use of the e2h operator to collapse the above expression into
(vec(E) @ vec(overlap_e)), allowing for simple calculation of the mode overlap. (vec(E) @ vec(overlap_e)), allowing for simple calculation of the mode overlap.
@ -211,45 +209,20 @@ def compute_overlap_e(E: field_t,
""" """
slices = tuple(slices) slices = tuple(slices)
cross_plane = [slice(None)] * 4 Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes,
cross_plane[axis + 1] = slices[axis] axis=axis, polarity=polarity, slices=slices)
cross_plane = tuple(cross_plane)
# Determine phase factors for parallel slices start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
a_shape = numpy.roll([-1, 1, 1], axis)
a_E = numpy.real(dxes[0][axis]).cumsum()
a_H = numpy.real(dxes[1][axis]).cumsum()
iphi = -polarity * 1j * wavenumber
phase_E = numpy.exp(iphi * (a_E - a_E[slices[axis]])).reshape(a_shape)
phase_H = numpy.exp(iphi * (a_H - a_H[slices[axis]])).reshape(a_shape)
# Expand our slice to the entire grid using the calculated phase factors slices2 = list(slices)
Ee = phase_E * E[cross_plane] slices2[axis] = slice(start, stop)
He = phase_H * H[cross_plane] slices2 = (slice(None), *slices2)
Etgt = numpy.zeros_like(Ee)
Etgt[slices2] = Ee[slices2]
# Write out the operator product for the mode orthogonality integral Etgt /= (Etgt.conj() * Etgt).sum()
domain = numpy.zeros_like(E[0], dtype=int) return Etgt.conj()
domain[slices] = 1
npts = E[0].size
dn = numpy.zeros(npts * 3, dtype=int)
dn[0:npts] = 1
dn = numpy.roll(dn, npts * axis)
e2h = operators.e2h(omega, dxes, mu)
ds = sparse.diags(vec([domain]*3))
h_cross_ = operators.poynting_h_cross(vec(He), dxes)
e_cross_ = operators.poynting_e_cross(vec(Ee), dxes)
overlap_e = dn @ ds @ (-h_cross_ + e_cross_ @ e2h)
# Normalize
dx_forward = dxes[0][axis][slices[axis]]
norm_factor = numpy.abs(overlap_e @ vec(Ee))
overlap_e /= norm_factor * dx_forward
return unvec(overlap_e, E[0].shape)
def solve_waveguide_mode_cylindrical(mode_number: int, def solve_waveguide_mode_cylindrical(mode_number: int,
@ -307,31 +280,6 @@ def solve_waveguide_mode_cylindrical(mode_number: int,
return fields return fields
def compute_overlap_ce(E: field_t,
wavenumber: complex,
dxes: dx_lists_t,
axis: int,
polarity: int,
slices: List[slice],
) -> field_t:
slices = tuple(slices)
Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes,
axis=axis, polarity=polarity, slices=slices)
start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity))
slices2 = list(slices)
slices2[axis] = slice(start, stop)
slices2 = (slice(None), *slices2)
Etgt = numpy.zeros_like(Ee)
Etgt[slices2] = Ee[slices2]
Etgt /= (Etgt.conj() * Etgt).sum()
return Etgt, slices2
def expand_wgmode_e(E: field_t, def expand_wgmode_e(E: field_t,
wavenumber: complex, wavenumber: complex,
dxes: dx_lists_t, dxes: dx_lists_t,