diff --git a/examples/fdfd.py b/examples/fdfd.py index f5130be..cf3d206 100644 --- a/examples/fdfd.py +++ b/examples/fdfd.py @@ -137,7 +137,9 @@ def test1(solver=generic_solver): 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']) - 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.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.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! ''' @@ -186,12 +195,14 @@ def test1(solver=generic_solver): pyplot.subplot(2, 2, 4) def poyntings(E): - e = vec(E) - h = operators.e2h(omega, dxes) @ e - cross1 = operators.poynting_e_cross(e, dxes) @ h.conj() - cross2 = operators.poynting_h_cross(h.conj(), dxes) @ e + H = functional.e2h(omega, dxes)(E) + poynting = 0.5 * fdtd.poynting(e=E, h=H.conj()) * dx * dx + cross1 = operators.poynting_e_cross(vec(E), dxes) @ vec(H).conj() +# cross2 = operators.poynting_h_cross(h.conj(), dxes) @ e 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 s1x, s2x = poyntings(E) @@ -202,7 +213,7 @@ def test1(solver=generic_solver): q = [] for i in range(-5, 30): 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.plot(q) pyplot.title('Overlap with mode') diff --git a/meanas/fdfd/operators.py b/meanas/fdfd/operators.py index c2df8ab..47cdf14 100644 --- a/meanas/fdfd/operators.py +++ b/meanas/fdfd/operators.py @@ -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]] - fx, fy, fz = [avgf(i, shape) for i in range(3)] - bx, by, bz = [avgb(i, shape) for i in range(3)] + bx, by, bz = [rotation(i, shape, -1) for i in range(3)] dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')] 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)] n = numpy.prod(shape) - zero = sparse.csr_matrix((n, n)) P = sparse.bmat( - [[ zero, -fx @ Ez @ bz @ dbgy, fx @ Ey @ by @ dbgz], - [ fy @ Ez @ bz @ dbgx, zero, -fy @ Ex @ bx @ dbgz], - [-fz @ Ey @ by @ dbgx, fz @ Ex @ bx @ dbgy, zero]]) + [[ None, -bx @ Ez @ dbgy, bx @ Ey @ dbgz], + [ by @ Ez @ dbgx, None, -by @ Ex @ dbgz], + [-bz @ Ey @ dbgx, bz @ Ex @ dbgy, None]]) 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]] - 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)] 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)) 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 + diff --git a/meanas/fdfd/waveguide.py b/meanas/fdfd/waveguide.py index 63ecb98..f84bbfe 100644 --- a/meanas/fdfd/waveguide.py +++ b/meanas/fdfd/waveguide.py @@ -186,7 +186,7 @@ def _normalized_fields(e: numpy.ndarray, norm_amplitude = 1 / numpy.sqrt(P) 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) sign = numpy.sign(E_weighted[:, :max(shape[0]//2, 1), :max(shape[1]//2, 1)].real.sum()) diff --git a/meanas/fdfd/waveguide_mode.py b/meanas/fdfd/waveguide_mode.py index 5c02655..cfa888a 100644 --- a/meanas/fdfd/waveguide_mode.py +++ b/meanas/fdfd/waveguide_mode.py @@ -112,6 +112,8 @@ def solve_waveguide_mode(mode_number: int, Apply corrections and expand to 3D ''' # 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) # Adjust for propagation direction @@ -179,20 +181,16 @@ def compute_source(E: field_t, def compute_overlap_e(E: field_t, - H: field_t, wavenumber: complex, - omega: complex, dxes: dx_lists_t, axis: int, polarity: int, slices: List[slice], - epsilon: field_t, # TODO unused?? - mu: field_t = None, - ) -> field_t: + ) -> field_t: # TODO DOCS """ 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) - [assumes reflection symmetry]. + [assumes reflection symmetry].i 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. @@ -211,45 +209,20 @@ def compute_overlap_e(E: field_t, """ slices = tuple(slices) - cross_plane = [slice(None)] * 4 - cross_plane[axis + 1] = slices[axis] - cross_plane = tuple(cross_plane) + Ee = expand_wgmode_e(E=E, wavenumber=wavenumber, dxes=dxes, + axis=axis, polarity=polarity, slices=slices) - # Determine phase factors for parallel slices - 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) + start, stop = sorted((slices[axis].start, slices[axis].start - 2 * polarity)) - # Expand our slice to the entire grid using the calculated phase factors - Ee = phase_E * E[cross_plane] - He = phase_H * H[cross_plane] + slices2 = list(slices) + slices2[axis] = slice(start, stop) + slices2 = (slice(None), *slices2) + Etgt = numpy.zeros_like(Ee) + Etgt[slices2] = Ee[slices2] - # Write out the operator product for the mode orthogonality integral - domain = numpy.zeros_like(E[0], dtype=int) - 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) + Etgt /= (Etgt.conj() * Etgt).sum() + return Etgt.conj() def solve_waveguide_mode_cylindrical(mode_number: int, @@ -307,31 +280,6 @@ def solve_waveguide_mode_cylindrical(mode_number: int, 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, wavenumber: complex, dxes: dx_lists_t,