From 7cbbaedcdb6ceea0352bf5be3a0267119810b899 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Wed, 19 Apr 2017 23:18:09 -0700 Subject: [PATCH 01/29] bump version number --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 3da94ed..4a3441a 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ from setuptools import setup, find_packages setup(name='fdfd_tools', - version='0.2', + version='0.3', description='FDFD Electromagnetic simulation tools', author='Jan Petykiewicz', author_email='anewusername@gmail.com', From 43d14642582151e726665d4f8d36dba276e4add9 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sat, 20 May 2017 21:21:50 -0700 Subject: [PATCH 02/29] Remote pyplot.hold It's deprecated now --- examples/test.py | 1 - 1 file changed, 1 deletion(-) diff --git a/examples/test.py b/examples/test.py index a40309f..a7e1746 100644 --- a/examples/test.py +++ b/examples/test.py @@ -190,7 +190,6 @@ def test1(solver=generic_solver): s1x, s2x = poyntings(E) pyplot.plot(s1x[0].sum(axis=2).sum(axis=1)) - pyplot.hold(True) pyplot.plot(s2x[0].sum(axis=2).sum(axis=1)) pyplot.show() From 50334723429e94660a7835884aa97d923f430a5d Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sat, 20 May 2017 21:22:28 -0700 Subject: [PATCH 03/29] Use ravel instead of flatten where possible --- fdfd_tools/operators.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/fdfd_tools/operators.py b/fdfd_tools/operators.py index 1691f36..c7fc578 100644 --- a/fdfd_tools/operators.py +++ b/fdfd_tools/operators.py @@ -304,7 +304,7 @@ def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmat i_ind = numpy.arange(n) j_ind = numpy.ravel_multi_index(ijk, shape, order='C') - vij = (numpy.ones(n), (i_ind, j_ind.flatten(order='C'))) + vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) d = sparse.csr_matrix(vij, shape=(n, n)) @@ -348,7 +348,7 @@ def shift_with_mirror(axis: int, shape: List[int], shift_distance: int=1) -> spa if len(shape) == 3: j_ind += ijk[2] * shape[0] * shape[1] - vij = (numpy.ones(n), (i_ind, j_ind.flatten(order='C'))) + vij = (numpy.ones(n), (i_ind, j_ind.ravel(order='C'))) d = sparse.csr_matrix(vij, shape=(n, n)) return d @@ -369,7 +369,7 @@ def deriv_forward(dx_e: List[numpy.ndarray]) -> List[sparse.spmatrix]: def deriv(axis): return rotation(axis, shape, 1) - sparse.eye(n) - Ds = [sparse.diags(+1 / dx.flatten(order='C')) @ deriv(a) + Ds = [sparse.diags(+1 / dx.ravel(order='C')) @ deriv(a) for a, dx in enumerate(dx_e_expanded)] return Ds @@ -390,7 +390,7 @@ def deriv_back(dx_h: List[numpy.ndarray]) -> List[sparse.spmatrix]: def deriv(axis): return rotation(axis, shape, -1) - sparse.eye(n) - Ds = [sparse.diags(-1 / dx.flatten(order='C')) @ deriv(a) + Ds = [sparse.diags(-1 / dx.ravel(order='C')) @ deriv(a) for a, dx in enumerate(dx_h_expanded)] return Ds @@ -461,8 +461,8 @@ def poynting_e_cross(e: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: fx, fy, fz = [avgf(i, shape) for i in range(3)] bx, by, bz = [avgb(i, shape) for i in range(3)] - dxag = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')] - dbgx, dbgy, dbgz = [sparse.diags(dx.flatten(order='C')) + dxag = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[0], indexing='ij')] + dbgx, dbgy, dbgz = [sparse.diags(dx.ravel(order='C')) for dx in numpy.meshgrid(*dxes[1], indexing='ij')] Ex, Ey, Ez = [sparse.diags(ei * da) for ei, da in zip(numpy.split(e, 3), dxag)] @@ -490,8 +490,8 @@ def poynting_h_cross(h: vfield_t, dxes: dx_lists_t) -> sparse.spmatrix: fx, fy, fz = [avgf(i, shape) for i in range(3)] bx, by, bz = [avgb(i, shape) for i in range(3)] - dxbg = [dx.flatten(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')] - dagx, dagy, dagz = [sparse.diags(dx.flatten(order='C')) + dxbg = [dx.ravel(order='C') for dx in numpy.meshgrid(*dxes[1], indexing='ij')] + dagx, dagy, dagz = [sparse.diags(dx.ravel(order='C')) for dx in numpy.meshgrid(*dxes[0], indexing='ij')] Hx, Hy, Hz = [sparse.diags(hi * db) for hi, db in zip(numpy.split(h, 3), dxbg)] From 9d337444276cd7b2367fc313b39e572bfbdd5b28 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sat, 20 May 2017 21:22:43 -0700 Subject: [PATCH 04/29] Fix docstring for rotation --- fdfd_tools/operators.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/fdfd_tools/operators.py b/fdfd_tools/operators.py index c7fc578..bfcfcdc 100644 --- a/fdfd_tools/operators.py +++ b/fdfd_tools/operators.py @@ -284,7 +284,8 @@ def m2j(omega: complex, def rotation(axis: int, shape: List[int], shift_distance: int=1) -> sparse.spmatrix: """ - Utility operator for performing a circular shift along a specified axis by 1 element. + Utility operator for performing a circular shift along a specified axis by a + specified number of elements. :param axis: Axis to shift along. x=0, y=1, z=2 :param shape: Shape of the grid being shifted From 6748181f8ff8d43e9792b33287b670bee07786cb Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sat, 20 May 2017 21:23:18 -0700 Subject: [PATCH 05/29] use logging module for progress reports --- fdfd_tools/solvers.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/fdfd_tools/solvers.py b/fdfd_tools/solvers.py index bb230f2..066725c 100644 --- a/fdfd_tools/solvers.py +++ b/fdfd_tools/solvers.py @@ -3,6 +3,7 @@ Solvers for FDFD problems. """ from typing import List, Callable, Dict, Any +import logging import numpy from numpy.linalg import norm @@ -11,6 +12,9 @@ import scipy.sparse.linalg from . import operators +logger = logging.getLogger(__name__) + + def _scipy_qmr(A: scipy.sparse.csr_matrix, b: numpy.ndarray, **kwargs @@ -29,20 +33,20 @@ def _scipy_qmr(A: scipy.sparse.csr_matrix, ''' iter = 0 - def print_residual(xk): + def log_residual(xk): nonlocal iter iter += 1 if iter % 100 == 0: - print('Solver residual at iteration', iter, ':', norm(A @ xk - b)) + logger.info('Solver residual at iteration {} : {}'.format(iter, norm(A @ xk - b))) if 'callback' in kwargs: def augmented_callback(xk): - print_residual(xk) + log_residual(xk) kwargs['callback'](xk) kwargs['callback'] = augmented_callback else: - kwargs['callback'] = print_residual + kwargs['callback'] = log_residual ''' Run the actual solve @@ -83,7 +87,7 @@ def generic(omega: complex, b: numpy.ndarray x: numpy.ndarray Default is a wrapped version of scipy.sparse.linalg.qmr() - which doesn't return convergence info and prints the residual + which doesn't return convergence info and logs the residual every 100 iterations. :param matrix_solver_opts: Passed as kwargs to matrix_solver(...) :return: E-field which solves the system. From c14298484db30b20a02e9aaca50b2f0767ce9c17 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 24 Sep 2017 19:11:56 -0700 Subject: [PATCH 06/29] Fix eigenvalue solver for complex matrices --- fdfd_tools/waveguide_mode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index df62143..51a6266 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -64,7 +64,7 @@ def solve_waveguide_mode_2d(mode_number: int, eigval = None for _ in range(40): - eigval = v @ A @ v + eigval = v.conj() @ A @ v if numpy.linalg.norm(A @ v - eigval * v) < 1e-13: break w = spalg.spsolve(A - eigval * sparse.eye(A.shape[0]), v) From d3c22006bdd4f7bb51a06d0782ea1ac4513c8833 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 24 Sep 2017 19:12:48 -0700 Subject: [PATCH 07/29] ie -> i.e. (docs) --- fdfd_tools/operators.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fdfd_tools/operators.py b/fdfd_tools/operators.py index bfcfcdc..02c1197 100644 --- a/fdfd_tools/operators.py +++ b/fdfd_tools/operators.py @@ -182,10 +182,10 @@ def eh_full(omega: complex, :param mu: Vectorized magnetic permeability (default 1 everywhere) :param pec: Vectorized mask specifying PEC cells. Any cells where pec != 0 are interpreted as containing a perfect electrical conductor (PEC). - The PEC is applied per-field-component (ie, pec.size == epsilon.size) + The PEC is applied per-field-component (i.e., pec.size == epsilon.size) :param pmc: Vectorized mask specifying PMC cells. Any cells where pmc != 0 are interpreted as containing a perfect magnetic conductor (PMC). - The PMC is applied per-field-component (ie, pmc.size == epsilon.size) + The PMC is applied per-field-component (i.e., pmc.size == epsilon.size) :return: Sparse matrix containing the wave operator """ if numpy.any(numpy.equal(pec, None)): From 7342c8efd7b246ae6c3182b5970dcddb0e86417f Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 24 Sep 2017 19:13:10 -0700 Subject: [PATCH 08/29] Use ravel instead of flatten for vec() --- fdfd_tools/vectorization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/vectorization.py b/fdfd_tools/vectorization.py index 2377d39..1c2134a 100644 --- a/fdfd_tools/vectorization.py +++ b/fdfd_tools/vectorization.py @@ -27,7 +27,7 @@ def vec(f: field_t) -> vfield_t: """ if numpy.any(numpy.equal(f, None)): return None - return numpy.hstack(tuple((fi.flatten(order='C') for fi in f))) + return numpy.hstack(tuple((fi.ravel(order='C') for fi in f))) def unvec(v: vfield_t, shape: numpy.ndarray) -> field_t: From 17fa2aa3d370635a1bb8b9c713b7b93979c46ba3 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 24 Sep 2017 19:13:37 -0700 Subject: [PATCH 09/29] In-place normalization during eigensolve --- fdfd_tools/waveguide_mode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index 51a6266..fac43e2 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -67,8 +67,8 @@ def solve_waveguide_mode_2d(mode_number: int, eigval = v.conj() @ A @ v if numpy.linalg.norm(A @ v - eigval * v) < 1e-13: break - w = spalg.spsolve(A - eigval * sparse.eye(A.shape[0]), v) - v = w / numpy.linalg.norm(w) + v = spalg.spsolve(A - eigval * sparse.eye(A.shape[0]), v) + v /= numpy.linalg.norm(v) # Calculate the wave-vector (force the real part to be positive) wavenumber = numpy.sqrt(eigval) From 001bf1e2eff28702084aa439b9b9a9fd648ef524 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 24 Sep 2017 19:14:30 -0700 Subject: [PATCH 10/29] Clarify eigensolver documentation --- fdfd_tools/waveguide_mode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index fac43e2..c6b9900 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -45,7 +45,7 @@ def solve_waveguide_mode_2d(mode_number: int, ''' Shift by the absolute value of the largest eigenvalue, then find a few of the - largest (shifted) eigenvalues. The shift ensures that we find the largest + largest-magnitude (shifted) eigenvalues. The shift ensures that we find the largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval) ''' From bacc6fea3f570f87046de80e3003908c52fd4183 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 24 Sep 2017 22:28:08 -0700 Subject: [PATCH 11/29] Move eigensolver code out to separate module --- fdfd_tools/eigensolvers.py | 95 ++++++++++++++++++++++++++++++++++++ fdfd_tools/waveguide.py | 2 +- fdfd_tools/waveguide_mode.py | 39 +++------------ 3 files changed, 102 insertions(+), 34 deletions(-) create mode 100644 fdfd_tools/eigensolvers.py diff --git a/fdfd_tools/eigensolvers.py b/fdfd_tools/eigensolvers.py new file mode 100644 index 0000000..6c06296 --- /dev/null +++ b/fdfd_tools/eigensolvers.py @@ -0,0 +1,95 @@ +""" +Solvers for eigenvalue / eigenvector problems +""" +from typing import Tuple, List +import numpy +from numpy.linalg import norm +from scipy import sparse +import scipy.sparse.linalg as spalg + + +def power_iteration(operator: sparse.spmatrix, + guess_vector: numpy.ndarray = None, + iterations: int = 20, + ) -> Tuple[complex, numpy.ndarray]: + """ + Use power iteration to estimate the dominant eigenvector of a matrix. + + :param operator: Matrix to analyze. + :param guess_vector: Starting point for the eigenvector. Default is a randomly chosen vector. + :param iterations: Number of iterations to perform. Default 20. + :return: (Largest-magnitude eigenvalue, Corresponding eigenvector estimate) + """ + if numpy.any(numpy.equal(guess_vector, None)): + v = numpy.random.rand(operator.shape[0]) + else: + v = guess_vector + + for _ in range(iterations): + v = operator @ v + v /= norm(v) + + lm_eigval = v.conj() @ operator @ v + return lm_eigval, v + + +def rayleigh_quotient_iteration(operator: sparse.spmatrix, + guess_vector: numpy.ndarray, + iterations: int = 40, + tolerance: float = 1e-13, + ) -> Tuple[complex, numpy.ndarray]: + """ + Use Rayleigh quotient iteration to refine an eigenvector guess. + + :param operator: Matrix to analyze. + :param guess_vector: Eigenvector to refine. + :param iterations: Maximum number of iterations to perform. Default 40. + :param tolerance: Stop iteration if (A - I*eigenvalue) @ v < tolerance. + Default 1e-13. + :return: (eigenvalue, eigenvector) + """ + v = guess_vector + for _ in range(iterations): + eigval = v.conj() @ operator @ v + if norm(operator @ v - eigval * v) < tolerance: + break + v = spalg.spsolve(operator - eigval * sparse.eye(operator.shape[0]), v) + v /= norm(v) + + return eigval, v + + +def signed_eigensolve(operator: sparse.spmatrix, + how_many: int, + negative: bool = False, + ) -> Tuple[numpy.ndarray, numpy.ndarray]: + """ + Find the largest-magnitude positive-only (or negative-only) eigenvalues and + eigenvectors of the provided matrix. + + :param operator: Matrix to analyze. + :param how_many: How many eigenvalues to find. + :param negative: Whether to find negative-only eigenvalues. + Default False (positive only). + :return: (sorted list of eigenvalues, 2D ndarray of corresponding eigenvectors) + eigenvectors[:, k] corresponds to the k-th eigenvalue + """ + # Use power iteration to estimate the dominant eigenvector + lm_eigval, _ = power_iteration(operator) + + ''' + Shift by the absolute value of the largest eigenvalue, then find a few of the + largest-magnitude (shifted) eigenvalues. A positive shift ensures that we find the + largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the + range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval) + ''' + if negative: + shifted_operator = operator - abs(lm_eigval) * sparse.eye(operator.shape[0]) + else: + shifted_operator = operator + abs(lm_eigval) * sparse.eye(operator.shape[0]) + + eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50) + + k = eigenvalues.argsort() + return eigenvalues[k], eigenvectors[:, k] + diff --git a/fdfd_tools/waveguide.py b/fdfd_tools/waveguide.py index a8ae1f2..0725bac 100644 --- a/fdfd_tools/waveguide.py +++ b/fdfd_tools/waveguide.py @@ -23,7 +23,7 @@ import numpy from numpy.linalg import norm import scipy.sparse as sparse -from . import unvec, dx_lists_t, field_t, vfield_t +from . import vec, unvec, dx_lists_t, field_t, vfield_t from . import operators diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index c6b9900..1670991 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -1,10 +1,10 @@ from typing import Dict, List import numpy import scipy.sparse as sparse -import scipy.sparse.linalg as spalg from . import vec, unvec, dx_lists_t, vfield_t, field_t from . import operators, waveguide, functional +from .eigensolvers import signed_eigensolve, rayleigh_quotient_iteration def solve_waveguide_mode_2d(mode_number: int, @@ -12,12 +12,12 @@ def solve_waveguide_mode_2d(mode_number: int, dxes: dx_lists_t, epsilon: vfield_t, mu: vfield_t = None, - wavenumber_correction: bool = True + wavenumber_correction: bool = True, ) -> Dict[str, complex or field_t]: """ Given a 2d region, attempts to solve for the eigenmode with the specified mode number. - :param mode_number: Number of the mode, 0-indexed + :param mode_number: Number of the mode, 0-indexed. :param omega: Angular frequency of the simulation :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header :param epsilon: Dielectric constant @@ -29,46 +29,19 @@ def solve_waveguide_mode_2d(mode_number: int, ''' Solve for the largest-magnitude eigenvalue of the real operator - by using power iteration. ''' dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] - A_r = waveguide.operator(numpy.real(omega), dxes_real, numpy.real(epsilon), numpy.real(mu)) - # Use power iteration for 20 steps to estimate the dominant eigenvector - v = numpy.random.rand(A_r.shape[0]) - for _ in range(20): - v = A_r @ v - v /= numpy.linalg.norm(v) - - lm_eigval = v @ A_r @ v - - ''' - Shift by the absolute value of the largest eigenvalue, then find a few of the - largest-magnitude (shifted) eigenvalues. The shift ensures that we find the largest - _positive_ eigenvalues, since any negative eigenvalues will be shifted to the range - 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval) - ''' - shifted_A_r = A_r + abs(lm_eigval) * sparse.eye(A_r.shape[0]) - eigvals, eigvecs = spalg.eigs(shifted_A_r, which='LM', k=mode_number + 3, ncv=50) - - # Pick the eigenvalue we want from the few we found - k = eigvals.argsort()[-(mode_number+1)] - v = eigvecs[:, k] + eigvals, eigvecs = signed_eigensolve(A_r, mode_number+3) + v = eigvecs[:, -(mode_number + 1)] ''' Now solve for the eigenvector of the full operator, using the real operator's eigenvector as an initial guess for Rayleigh quotient iteration. ''' A = waveguide.operator(omega, dxes, epsilon, mu) - - eigval = None - for _ in range(40): - eigval = v.conj() @ A @ v - if numpy.linalg.norm(A @ v - eigval * v) < 1e-13: - break - v = spalg.spsolve(A - eigval * sparse.eye(A.shape[0]), v) - v /= numpy.linalg.norm(v) + v, eigval = rayleigh_quotient_iteration(A, v) # Calculate the wave-vector (force the real part to be positive) wavenumber = numpy.sqrt(eigval) From a4616982ca5d986a8a12e00a86f7e87e686a20ee Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Sun, 24 Sep 2017 22:28:39 -0700 Subject: [PATCH 12/29] Add cylindrical coordinate 2D modesolver code --- fdfd_tools/waveguide.py | 59 ++++++++++++++++++++++++++++++++ fdfd_tools/waveguide_mode.py | 66 ++++++++++++++++++++++++++++++++++++ 2 files changed, 125 insertions(+) diff --git a/fdfd_tools/waveguide.py b/fdfd_tools/waveguide.py index 0725bac..85e96bf 100644 --- a/fdfd_tools/waveguide.py +++ b/fdfd_tools/waveguide.py @@ -307,3 +307,62 @@ def e_err(e: vfield_t, op = ch @ mu_inv @ ce @ e - omega ** 2 * (epsilon * e) return norm(op) / norm(e) + + +def cylindrical_operator(omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + r0: float, + ) -> sparse.spmatrix: + """ + Cylindrical coordinate waveguide operator of the form + + TODO + + for use with a field vector of the form [E_r, E_y]. + + This operator can be used to form an eigenvalue problem of the form + A @ [E_r, E_y] = wavenumber**2 * [E_r, E_y] + + which can then be solved for the eigenmodes of the system (an exp(-i * wavenumber * theta) + theta-dependence is assumed for the fields). + + :param omega: The angular frequency of the system + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header (2D) + :param epsilon: Vectorized dielectric constant grid + :param r0: Radius of curvature for the simulation. This should be the minimum value of + r within the simulation domain. + :return: Sparse matrix representation of the operator + """ + + Dfx, Dfy = operators.deriv_forward(dxes[0]) + Dbx, Dby = operators.deriv_back(dxes[1]) + + rx = r0 + numpy.cumsum(dxes[0][0]) + ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0]) + tx = 1 + rx/r0 + ty = 1 + ry/r0 + + Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1))) + Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1))) + + eps_parts = numpy.split(epsilon, 3) + eps_x = sparse.diags(eps_parts[0]) + eps_y = sparse.diags(eps_parts[1]) + eps_z_inv = sparse.diags(1 / eps_parts[2]) + + pa = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dbx, Dby)) + pb = sparse.vstack((Dfx, Dfy)) @ Tx @ eps_z_inv @ sparse.hstack((Dby, Dbx)) + a0 = Ty @ eps_x + omega**-2 * Dby @ Ty @ Dfy + a1 = Tx @ eps_y + omega**-2 * Dbx @ Ty @ Dfx + b0 = Dbx @ Ty @ Dfy + b1 = Dby @ Ty @ Dfx + + diag = sparse.block_diag + op = (omega**2 * diag((Tx, Ty)) + pa) @ diag((a0, a1)) + \ + - (sparse.bmat(((None, Ty), (Tx, None))) + omega**-2 * pb) @ diag((b0, b1)) + + return op + + + diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index 1670991..7c7bc5c 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -272,3 +272,69 @@ def compute_overlap_e(E: field_t, overlap_e /= norm_factor * dx_forward return unvec(overlap_e, E[0].shape) + + +def solve_waveguide_mode_cylindrical(mode_number: int, + omega: complex, + dxes: dx_lists_t, + epsilon: vfield_t, + r0: float, + wavenumber_correction: bool = True, + ) -> Dict[str, complex or field_t]: + """ + Given a 2d (r, y) slice of epsilon, attempts to solve for the eigenmode + of the bent waveguide with the specified mode number. + + :param mode_number: Number of the mode, 0-indexed + :param omega: Angular frequency of the simulation + :param dxes: Grid parameters [dx_e, dx_h] as described in fdfd_tools.operators header. + The first coordinate is assumed to be r, the second is y. + :param epsilon: Dielectric constant + :param r0: Radius of curvature for the simulation. This should be the minimum value of + r within the simulation domain. + :param wavenumber_correction: Whether to correct the wavenumber to + account for numerical dispersion (default True) + :return: {'E': List[numpy.ndarray], 'H': List[numpy.ndarray], 'wavenumber': complex} + """ + + ''' + Solve for the largest-magnitude eigenvalue of the real operator + ''' + dxes_real = [[numpy.real(dx) for dx in dxi] for dxi in dxes] + + A_r = waveguide.cylindrical_operator(numpy.real(omega), dxes_real, numpy.real(epsilon), r0) + eigvals, eigvecs = signed_eigensolve(A_r, mode_number + 3) + v = eigvecs[:, -(mode_number+1)] + + ''' + Now solve for the eigenvector of the full operator, using the real operator's + eigenvector as an initial guess for Rayleigh quotient iteration. + ''' + A = waveguide.cylindrical_operator(omega, dxes, epsilon, r0) + eigval, v = rayleigh_quotient_iteration(A, v) + + # Calculate the wave-vector (force the real part to be positive) + wavenumber = numpy.sqrt(eigval) + wavenumber *= numpy.sign(numpy.real(wavenumber)) + + ''' + Perform correction on wavenumber to account for numerical dispersion. + + See Numerical Dispersion in Taflove's FDTD book. + This correction term reduces the error in emitted power, but additional + error is introduced into the E_err and H_err terms. This effect becomes + more pronounced as beta increases. + ''' + if wavenumber_correction: + wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - numpy.real(wavenumber) + + shape = [d.size for d in dxes[0]] + v = numpy.hstack((v, numpy.zeros(shape[0] * shape[1]))) + fields = { + 'wavenumber': wavenumber, + 'E': unvec(v, shape), +# 'E': unvec(e, shape), +# 'H': unvec(h, shape), + } + + return fields From ea04fc42bec11a30430b7ade52bbbba90103fdd7 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Mon, 25 Sep 2017 00:10:01 -0700 Subject: [PATCH 13/29] Fix switched args --- fdfd_tools/waveguide_mode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index 7c7bc5c..cffc916 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -41,7 +41,7 @@ def solve_waveguide_mode_2d(mode_number: int, eigenvector as an initial guess for Rayleigh quotient iteration. ''' A = waveguide.operator(omega, dxes, epsilon, mu) - v, eigval = rayleigh_quotient_iteration(A, v) + eigval, v = rayleigh_quotient_iteration(A, v) # Calculate the wave-vector (force the real part to be positive) wavenumber = numpy.sqrt(eigval) From a4cc96395397748b5ff047227693efc505d55e1e Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Tue, 17 Oct 2017 12:58:15 -0700 Subject: [PATCH 14/29] bump version number --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 4a3441a..ef1df08 100644 --- a/setup.py +++ b/setup.py @@ -3,7 +3,7 @@ from setuptools import setup, find_packages setup(name='fdfd_tools', - version='0.3', + version='0.4', description='FDFD Electromagnetic simulation tools', author='Jan Petykiewicz', author_email='anewusername@gmail.com', From 73e3fa18b11c962dffcbed80a03d522bd38e3fb2 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Tue, 17 Oct 2017 13:00:46 -0700 Subject: [PATCH 15/29] fix cylindrical operator --- fdfd_tools/waveguide.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/fdfd_tools/waveguide.py b/fdfd_tools/waveguide.py index 85e96bf..1566a74 100644 --- a/fdfd_tools/waveguide.py +++ b/fdfd_tools/waveguide.py @@ -340,8 +340,8 @@ def cylindrical_operator(omega: complex, rx = r0 + numpy.cumsum(dxes[0][0]) ry = r0 + dxes[0][0]/2.0 + numpy.cumsum(dxes[1][0]) - tx = 1 + rx/r0 - ty = 1 + ry/r0 + tx = rx/r0 + ty = ry/r0 Tx = sparse.diags(vec(tx[:, None].repeat(dxes[0][1].size, axis=1))) Ty = sparse.diags(vec(ty[:, None].repeat(dxes[1][1].size, axis=1))) From 4bf862761100bdb103a7e23034b02a0b5b47249d Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 5 Nov 2017 14:50:30 -0800 Subject: [PATCH 16/29] clarify beta -> wavenumber --- fdfd_tools/waveguide_mode.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/waveguide_mode.py b/fdfd_tools/waveguide_mode.py index cffc916..e50cc6a 100644 --- a/fdfd_tools/waveguide_mode.py +++ b/fdfd_tools/waveguide_mode.py @@ -323,7 +323,7 @@ def solve_waveguide_mode_cylindrical(mode_number: int, See Numerical Dispersion in Taflove's FDTD book. This correction term reduces the error in emitted power, but additional error is introduced into the E_err and H_err terms. This effect becomes - more pronounced as beta increases. + more pronounced as the wavenumber increases. ''' if wavenumber_correction: wavenumber -= 2 * numpy.sin(numpy.real(wavenumber / 2)) - numpy.real(wavenumber) From 6503b488ced14a023c2b5a8345efcf73458ce36f Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 5 Nov 2017 16:33:19 -0800 Subject: [PATCH 17/29] add farfield.py --- fdfd_tools/farfield.py | 220 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 220 insertions(+) create mode 100644 fdfd_tools/farfield.py diff --git a/fdfd_tools/farfield.py b/fdfd_tools/farfield.py new file mode 100644 index 0000000..84a04ba --- /dev/null +++ b/fdfd_tools/farfield.py @@ -0,0 +1,220 @@ +""" +Functions for performing near-to-farfield transformation (and the reverse). +""" +from typing import Dict, List +import numpy +from numpy.fft import fft2, fftshift, fftfreq, ifft2, ifftshift +from numpy import pi + + +def near_to_farfield(E_near: List[numpy.ndarray], + H_near: List[numpy.ndarray], + dx: float, + dy: float, + padded_size: List[int] = None + ) -> Dict[str]: + """ + Compute the farfield, i.e. the distribution of the fields after propagation + through several wavelengths of uniform medium. + + The input fields should be complex phasors. + + :param E_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse + E fields (e.g. [Ex, Ey] for calculating the farfield toward the z-direction). + :param H_near: List of 2 ndarrays containing the 2D phasor field slices for the transverse + H fields (e.g. [Hx, hy] for calculating the farfield towrad the z-direction). + :param dx: Cell size along x-dimension, in units of wavelength. + :param dy: Cell size along y-dimension, in units of wavelength. + :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`. + Powers of 2 are most efficient for FFT computation. + Default is the smallest power of 2 larger than the input, for each axis. + :returns: Dict with keys + 'E_far': Normalized E-field farfield; multiply by + (i k exp(-i k r) / (4 pi r)) to get the actual field value. + 'H_far': Normalized H-field farfield; multiply by + (i k exp(-i k r) / (4 pi r)) to get the actual field value. + 'kx', 'ky': Wavevector values corresponding to the x- and y- axes in E_far and H_far, + normalized to wavelength (dimensionless). + 'dkx', 'dky': step size for kx and ky, normalized to wavelength. + 'theta': arctan2(ky, kx) corresponding to each (kx, ky). + This is the angle in the x-y plane, counterclockwise from above, starting from +x. + 'phi': arccos(kz / k) corresponding to each (kx, ky). + This is the angle away from +z. + """ + + if not len(E_near) == 2: + raise Exception('E_near must be a length-2 list of ndarrays') + if not len(H_near) == 2: + raise Exception('H_near must be a length-2 list of ndarrays') + + s = E_near[0].shape + if not all(s == f.shape for f in E_near + H_near): + raise Exception('All fields must be the same shape!') + + if padded_size is None: + padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int) + if not hasattr(padded_size, '__len__'): + padded_size = (padded_size, padded_size) + + En_fft = [fftshift(fft2(fftshift(Eni), s=padded_size)) for Eni in E_near] + Hn_fft = [fftshift(fft2(fftshift(Hni), s=padded_size)) for Hni in H_near] + + # Propagation vectors kx, ky + k = 2 * pi + kxx = 2 * pi * fftshift(fftfreq(padded_size[0], dx)) + kyy = 2 * pi * fftshift(fftfreq(padded_size[1], dy)) + + kx, ky = numpy.meshgrid(kxx, kyy, indexing='ij') + kxy2 = kx * kx + ky * ky + kxy = numpy.sqrt(kxy2) + kz = numpy.sqrt(k * k - kxy2) + + sin_th = ky / kxy + cos_th = kx / kxy + cos_phi = kz / k + + sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0 + cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1 + + # Normalized vector potentials N, L + N = [-Hn_fft[1] * cos_phi * cos_th + Hn_fft[0] * cos_phi * sin_th, + Hn_fft[1] * sin_th + Hn_fft[0] * cos_th] + L = [ En_fft[1] * cos_phi * cos_th - En_fft[0] * cos_phi * sin_th, + -En_fft[1] * sin_th - En_fft[0] * cos_th] + + E_far = [-L[1] - N[0], + L[0] - N[1]] + H_far = [-E_far[1], + E_far[0]] + + theta = numpy.arctan2(ky, kx) + phi = numpy.arccos(cos_phi) + theta[numpy.logical_and(kx == 0, ky == 0)] = 0 + phi[numpy.logical_and(kx == 0, ky == 0)] = 0 + + # Zero fields beyond valid (phi, theta) + invalid_ind = kxy2 >= k * k + theta[invalid_ind] = 0 + phi[invalid_ind] = 0 + for i in range(2): + E_far[i][invalid_ind] = 0 + H_far[i][invalid_ind] = 0 + + outputs = { + 'E': E_far, + 'H': H_far, + 'dkx': kx[1]-kx[0], + 'dky': ky[1]-ky[0], + 'kx': kx, + 'ky': ky, + 'theta': theta, + 'phi': phi, + } + + return outputs + + + +def far_to_nearfield(E_far: List[numpy.ndarray], + H_far: List[numpy.ndarray], + dkx: float, + dky: float, + padded_size: List[int] = None + ) -> Dict[str]: + """ + Compute the farfield, i.e. the distribution of the fields after propagation + through several wavelengths of uniform medium. + + The input fields should be complex phasors. + + :param E_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse + E fields (e.g. [Ex, Ey] for calculating the nearfield toward the z-direction). + Fields should be normalized so that + E_far = E_far_actual / (i k exp(-i k r) / (4 pi r)) + :param H_far: List of 2 ndarrays containing the 2D phasor field slices for the transverse + H fields (e.g. [Hx, hy] for calculating the nearfield toward the z-direction). + Fields should be normalized so that + H_far = H_far_actual / (i k exp(-i k r) / (4 pi r)) + :param dkx: kx discretization, in units of wavelength. + :param dky: ky discretization, in units of wavelength. + :param padded_size: Shape of the output. A single integer `n` will be expanded to `(n, n)`. + Powers of 2 are most efficient for FFT computation. + Default is the smallest power of 2 larger than the input, for each axis. + :returns: Dict with keys + 'E': E-field nearfield + 'H': H-field nearfield + 'dx', 'dy': spatial discretization, normalized to wavelength (dimensionless) + """ + + if not len(E_far) == 2: + raise Exception('E_far must be a length-2 list of ndarrays') + if not len(H_far) == 2: + raise Exception('H_far must be a length-2 list of ndarrays') + + s = E_far[0].shape + if not all(s == f.shape for f in E_far + H_far): + raise Exception('All fields must be the same shape!') + + if padded_size is None: + padded_size = (2**numpy.ceil(numpy.log2(s))).astype(int) + if not hasattr(padded_size, '__len__'): + padded_size = (padded_size, padded_size) + + + k = 2 * pi + kxs = fftshift(fftfreq(s[0], 1/(s[0] * dkx))) + kys = fftshift(fftfreq(s[0], 1/(s[1] * dky))) + + kx, ky = numpy.meshgrid(kxs, kys, indexing='ij') + kxy2 = kx * kx + ky * ky + kxy = numpy.sqrt(kxy2) + + kz = numpy.sqrt(k * k - kxy2) + + sin_th = ky / kxy + cos_th = kx / kxy + cos_phi = kz / k + + sin_th[numpy.logical_and(kx == 0, ky == 0)] = 0 + cos_th[numpy.logical_and(kx == 0, ky == 0)] = 1 + + # Zero fields beyond valid (phi, theta) + invalid_ind = kxy2 >= k * k + theta[invalid_ind] = 0 + phi[invalid_ind] = 0 + for i in range(2): + E_far[i][invalid_ind] = 0 + H_far[i][invalid_ind] = 0 + + + # Normalized vector potentials N, L + L = [0.5 * E_far[1], + -0.5 * E_far[0]] + N = [L[1], + -L[0]] + + En_fft = [-( L[0] * sin_th + L[1] * cos_phi * cos_th)/cos_phi, + -(-L[0] * cos_th + L[1] * cos_phi * sin_th)/cos_phi] + + Hn_fft = [( N[0] * sin_th + N[1] * cos_phi * cos_th)/cos_phi, + (-N[0] * cos_th + N[1] * cos_phi * sin_th)/cos_phi] + + for i in range(2): + En_fft[i][cos_phi == 0] = 0 + Hn_fft[i][cos_phi == 0] = 0 + + E_near = [ifftshift(ifft2(ifftshift(Ei), s=padded_size)) for Ei in En_fft] + H_near = [ifftshift(ifft2(ifftshift(Hi), s=padded_size)) for Hi in Hn_fft] + + dx = 2 * pi / (s[0] * dkx) + dy = 2 * pi / (s[0] * dky) + + outputs = { + 'E': E_near, + 'H': H_near, + 'dx': dx, + 'dy': dy, + } + + return outputs + From 4aa2d07cef0e8793bbe48bfcf4ac246b193d467f Mon Sep 17 00:00:00 2001 From: jan Date: Sat, 9 Dec 2017 18:21:37 -0800 Subject: [PATCH 18/29] Add Bloch eigenproblem --- fdfd_tools/bloch.py | 406 ++++++++++++++++++++++++++++++++++++ fdfd_tools/eigensolvers.py | 19 +- fdfd_tools/vectorization.py | 2 +- 3 files changed, 420 insertions(+), 7 deletions(-) create mode 100644 fdfd_tools/bloch.py diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py new file mode 100644 index 0000000..d7a2ffa --- /dev/null +++ b/fdfd_tools/bloch.py @@ -0,0 +1,406 @@ +''' +Bloch eigenmode solver/operators + +This module contains functions for generating and solving the + 3D Bloch eigenproblem. The approach is to transform the problem + into the (spatial) fourier domain, transforming the equation + 1/mu * curl(1/eps * curl(H)) = (w/c)^2 H + into + conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k + where: + - the _k subscript denotes a 3D fourier transformed field + - each component of H_k corresponds to a plane wave with wavevector k + - x is the cross product + - conv denotes convolution + + Since k and H are orthogonal for each plane wave, we can use each + k to create an orthogonal basis (k, m, n), with k x m = n, and + |m| = |n| = 1. The cross products are then simplified with + + k @ h = kx hx + ky hy + kz hz = 0 = hk + h = hk + hm + hn = hm + hn + k = kk + km + kn = kk = |k| + + k x h = (ky hz - kz hy, + kz hx - kx hz, + kx hy - ky hx) + = ((k x h) @ k, (k x h) @ m, (k x h) @ n)_kmn + = (0, (m x k) @ h, (n x k) @ h)_kmn # triple product ordering + = (0, kk (-n @ h), kk (m @ h))_kmn # (m x k) = -|k| n, etc. + = |k| (0, -h @ n, h @ m)_kmn + + k x h = (km hn - kn hm, + kn hk - kk hn, + kk hm - km hk)_kmn + = (0, -kk hn, kk hm)_kmn + = (-kk hn)(mx, my, mz) + (kk hm)(nx, ny, nz) + = |k| (hm * (nx, ny, nz) - hn * (mx, my, mz)) + + where h is shorthand for H_k, (...)_kmn deontes the (k, m, n) basis, + and e.g. hm is the component of h in the m direction. + + We can also simplify conv(X_k, Y_k) as fftn(X * ifftn(Y_k)). + + Using these results and storing H_k as h = (hm, hn), we have + e_xyz = fftn(1/eps * ifftn(|k| (hm * n - hn * m))) + b_mn = |k| (-e_xyz @ n, e_xyz @ m) + h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n)) + which forms the operator from the left side of the equation. + + We can then use ARPACK in shift-invert mode (via scipy.linalg.eigs) + to find the eigenvectors for this operator. + + This approach is similar to the one used in MPB and derived at the start of + SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods + for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001) + + === + + Typically you will want to do something like + + recip_lattice = numpy.diag(1/numpy.array(epsilon[0].shape * dx)) + n, v = bloch.eigsolve(5, k0, recip_lattice, epsilon) + f = numpy.sqrt(-numpy.real(n[0])) + n_eff = norm(recip_lattice @ k0) / f + + v2e = bloch.hmn_2_exyz(k0, recip_lattice, epsilon) + e_field = v2e(v[0]) + + k, f = find_k(frequency=1/1550, + tolerance=(1/1550 - 1/1551), + search_direction=[1, 0, 0], + G_matrix=recip_lattice, + epsilon=epsilon, + band=0) + +''' + +from typing import List, Tuple, Callable, Dict +import numpy +from numpy.fft import fftn, ifftn, fftfreq +import scipy +from scipy.linalg import norm +import scipy.sparse.linalg as spalg + +from . import field_t + + +def generate_kmn(k0: numpy.ndarray, + G_matrix: numpy.ndarray, + shape: numpy.ndarray + ) -> Tuple[numpy.ndarray, numpy.ndarray, numpy.ndarray]: + """ + Generate a (k, m, n) orthogonal basis for each k-vector in the simulation grid. + + :param k0: [k0x, k0y, k0z], Bloch wavevector, in G basis. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param shape: [nx, ny, nz] shape of the simulation grid. + :return: (|k|, m, n) where |k| has shape tuple(shape) + (1,) + and m, n have shape tuple(shape) + (3,). + All are given in the xyz basis (e.g. |k|[0,0,0] = norm(G_matrix @ k0)). + """ + k0 = numpy.array(k0) + + Gi_grids = numpy.meshgrid(*(fftfreq(n, 1/n) for n in shape[:3]), indexing='ij') + Gi = numpy.stack(Gi_grids, axis=3) + + k_G = k0[None, None, None, :] - Gi + k_xyz = numpy.rollaxis(G_matrix @ numpy.rollaxis(k_G, 3, 2), 3, 2) + + m = numpy.broadcast_to([0, 1, 0], tuple(shape[:3]) + (3,)).astype(float) + n = numpy.broadcast_to([0, 0, 1], tuple(shape[:3]) + (3,)).astype(float) + + xy_non0 = numpy.any(k_xyz[:, :, :, 0:1] != 0, axis=3) + if numpy.any(xy_non0): + u = numpy.cross(k_xyz[xy_non0], [0, 0, 1]) + m[xy_non0, :] = u / norm(u, axis=1)[:, None] + + z_non0 = numpy.any(k_xyz != 0, axis=3) + if numpy.any(z_non0): + v = numpy.cross(k_xyz[z_non0], m[z_non0]) + n[z_non0, :] = v / norm(v, axis=1)[:, None] + + k_mag = norm(k_xyz, axis=3)[:, :, :, None] + return k_mag, m, n + + +def maxwell_operator(k0: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t, + mu: field_t = None + ) -> Callable[[numpy.ndarray], numpy.ndarray]: + """ + Generate the Maxwell operator + conv(1/mu_k, ik x conv(1/eps_k, ik x ___)) + which is the spatial-frequency-space representation of + 1/mu * curl(1/eps * curl(___)) + + The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size) + + See the module-level docstring for more information. + + :param k0: Bloch wavevector, [k0x, k0y, k0z]. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + All fields are sampled at cell centers (i.e., NOT Yee-gridded) + :param mu: Magnetic permability distribution for the simulation. + Default None (1 everywhere). + :return: Function which applies the maxwell operator to h_mn. + """ + + shape = epsilon[0].shape + (1,) + k_mag, m, n = generate_kmn(k0, G_matrix, shape) + + epsilon = numpy.stack(epsilon, 3) + if mu is not None: + mu = numpy.stack(mu, 3) + + def operator(h: numpy.ndarray): + """ + Maxwell operator for Bloch eigenmode simulation. + + h is complex 2-field in (m, n) basis, vectorized + + :param h: Raveled h_mn; size (2 * epsilon[0].size). + :return: Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)). + """ + hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] + + #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis + + # cross product and transform into xyz basis + d_xyz = (n * hin_m - + m * hin_n) * k_mag + + # divide by epsilon + e_xyz = ifftn(fftn(d_xyz, axes=range(3)) / epsilon, axes=range(3)) + + # cross product and transform into mn basis + b_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] * -k_mag + b_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] * +k_mag + + if mu is None: + h_m, h_n = b_m, b_n + else: + # transform from mn to xyz + b_xyz = (m * b_m[:, :, :, None] + + n * b_n[:, :, :, None]) + + # divide by mu + h_xyz = ifftn(fftn(b_xyz, axes=range(3)) / mu, axes=range(3)) + + # transform back to mn + h_m = numpy.sum(h_xyz * m, axis=3) + h_n = numpy.sum(h_xyz * n, axis=3) + return numpy.hstack((h_m.ravel(), h_n.ravel())) + + return operator + + +def hmn_2_exyz(k0: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t, + ) -> Callable[[numpy.ndarray], field_t]: + """ + Generate an operator which converts a vectorized spatial-frequency-space + h_mn into an E-field distribution, i.e. + ifft(conv(1/eps_k, ik x h_mn)) + + The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size) + + See the module-level docstring for more information. + + :param k0: Bloch wavevector, [k0x, k0y, k0z]. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + All fields are sampled at cell centers (i.e., NOT Yee-gridded) + :return: Function for converting h_mn into E_xyz + """ + shape = epsilon[0].shape + (1,) + epsilon = numpy.stack(epsilon, 3) + + k_mag, m, n = generate_kmn(k0, G_matrix, shape) + + def operator(h: numpy.ndarray) -> field_t: + hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] + d_xyz = (n * hin_m - + m * hin_n) * k_mag + + # divide by epsilon + return [ei for ei in numpy.rollaxis(fftn(d_xyz, axes=range(3)) / epsilon, 3)] + + return operator + + +def hmn_2_hxyz(k0: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t + ) -> Callable[[numpy.ndarray], field_t]: + """ + Generate an operator which converts a vectorized spatial-frequency-space + h_mn into an H-field distribution, i.e. + ifft(h_mn) + + The operator is a function that acts on a vector h_mn of size (2 * epsilon[0].size) + + See the module-level docstring for more information. + + :param k0: Bloch wavevector, [k0x, k0y, k0z]. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + Only epsilon[0].shape is used. + :return: Function for converting h_mn into H_xyz + """ + shape = epsilon[0].shape + (1,) + k_mag, m, n = generate_kmn(k0, G_matrix, shape) + + def operator(h: numpy.ndarray): + hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] + h_xyz = (m * hin_m + + n * hin_n) + return [fftn(hi) for hi in numpy.rollaxis(h_xyz, 3)] + + return operator + + +def inverse_maxwell_operator_approx(k0: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t, + mu: field_t = None + ) -> Callable[[numpy.ndarray], numpy.ndarray]: + """ + Generate an approximate inverse of the Maxwell operator, + ik x conv(eps_k, ik x conv(mu_k, ___)) + which can be used to improve the speed of ARPACK in shift-invert mode. + + See the module-level docstring for more information. + + :param k0: Bloch wavevector, [k0x, k0y, k0z]. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + All fields are sampled at cell centers (i.e., NOT Yee-gridded) + :param mu: Magnetic permability distribution for the simulation. + Default None (1 everywhere). + :return: Function which applies the approximate inverse of the maxwell operator to h_mn. + """ + shape = epsilon[0].shape + (1,) + epsilon = numpy.stack(epsilon, 3) + + k_mag, m, n = generate_kmn(k0, G_matrix, shape) + + if mu is not None: + mu = numpy.stack(mu, 3) + + def operator(h: numpy.ndarray): + """ + Approximate inverse Maxwell operator for Bloch eigenmode simulation. + + h is complex 2-field in (m, n) basis, vectorized + + :param h: Raveled h_mn; size (2 * epsilon[0].size). + :return: Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn)) + """ + hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] + + #{d,e,h}_xyz fields are complex 3-fields in (1/x, 1/y, 1/z) basis + + if mu is None: + b_m, b_n = hin_m, hin_n + else: + # transform from mn to xyz + h_xyz = (m * hin_m[:, :, :, None] + + n * hin_n[:, :, :, None]) + + # multiply by mu + b_xyz = ifftn(fftn(h_xyz, axes=range(3)) * mu, axes=range(3)) + + # transform back to mn + b_m = numpy.sum(b_xyz * m, axis=3) + b_n = numpy.sum(b_xyz * n, axis=3) + + # cross product and transform into xyz basis + e_xyz = (n * b_m - + m * b_n) / k_mag + + # multiply by epsilon + d_xyz = ifftn(fftn(e_xyz, axes=range(3)) * epsilon, axes=range(3)) + + # cross product and transform into mn basis crossinv_t2c + h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag + h_n = numpy.sum(e_xyz * m, axis=3)[:, :, :, None] / -k_mag + + return numpy.hstack((h_m.ravel(), h_n.ravel())) + + return operator + + +def eigsolve(num_modes: int, + k0: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t, + mu: field_t = None + ) -> Tuple[numpy.ndarray, numpy.ndarray]: + """ + Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector + k0 of the specified structure. + + :param k0: Bloch wavevector, [k0x, k0y, k0z]. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + All fields are sampled at cell centers (i.e., NOT Yee-gridded) + :param mu: Magnetic permability distribution for the simulation. + Default None (1 everywhere). + :return: (eigenvalues, eigenvectors) where eigenvalues[i] corresponds to the + vector eigenvectors[i, :] + """ + h_size = 2 * epsilon[0].size + + mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu) + imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu) + + scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop) + scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop) + + _eigvals, eigvecs = spalg.eigs(scipy_op, num_modes, sigma=0, OPinv=scipy_iop, which='LM') + eigvals = numpy.sum(eigvecs * (scipy_op @ eigvecs), axis=0) / numpy.sum(eigvecs * eigvecs, axis=0) + order = numpy.argsort(-eigvals) + return eigvals[order], eigvecs.T[order] + + +def find_k(frequency: float, + tolerance: float, + search_direction: numpy.ndarray, + G_matrix: numpy.ndarray, + epsilon: field_t, + mu: field_t = None, + band: int = 0 + ) -> Tuple[numpy.ndarray, float]: + """ + Search for a bloch vector that has a given frequency. + + :param frequency: Target frequency. + :param tolerance: Target frequency tolerance. + :param search_direction: k-vector direction to search along. + :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. + :param epsilon: Dielectric constant distribution for the simulation. + All fields are sampled at cell centers (i.e., NOT Yee-gridded) + :param mu: Magnetic permability distribution for the simulation. + Default None (1 everywhere). + :param band: Which band to search in. Default 0 (lowest frequency). + return: (k, actual_frequency) The found k-vector and its frequency + """ + + search_direction = numpy.array(search_direction) / norm(search_direction) + + def get_f(k0_mag: float, band: int = 0): + k0 = search_direction * k0_mag + n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon) + f = numpy.sqrt(numpy.abs(numpy.real(n[band]))) + return f + + res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency), 0.25, + method='Bounded', bounds=(0, 0.5), + options={'xatol': abs(tolerance)}) + return res.x * search_direction, res.fun + frequency + + diff --git a/fdfd_tools/eigensolvers.py b/fdfd_tools/eigensolvers.py index 6c06296..24d7339 100644 --- a/fdfd_tools/eigensolvers.py +++ b/fdfd_tools/eigensolvers.py @@ -29,7 +29,7 @@ def power_iteration(operator: sparse.spmatrix, v = operator @ v v /= norm(v) - lm_eigval = v.conj() @ operator @ v + lm_eigval = v.conj() @ (operator @ v) return lm_eigval, v @@ -59,7 +59,7 @@ def rayleigh_quotient_iteration(operator: sparse.spmatrix, return eigval, v -def signed_eigensolve(operator: sparse.spmatrix, +def signed_eigensolve(operator: sparse.spmatrix or spalg.LinearOperator, how_many: int, negative: bool = False, ) -> Tuple[numpy.ndarray, numpy.ndarray]: @@ -83,12 +83,19 @@ def signed_eigensolve(operator: sparse.spmatrix, largest _positive_ eigenvalues, since any negative eigenvalues will be shifted to the range 0 >= neg_eigval + abs(lm_eigval) > abs(lm_eigval) ''' + shift = numpy.abs(lm_eigval) if negative: - shifted_operator = operator - abs(lm_eigval) * sparse.eye(operator.shape[0]) - else: - shifted_operator = operator + abs(lm_eigval) * sparse.eye(operator.shape[0]) + shift *= -1 - eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50) + # Try to combine, use general LinearOperator if we fail + try: + shifted_operator = operator + shift * sparse.eye(operator.shape[0]) + except TypeError: + shifted_operator = operator + spalg.LinearOperator(shape=operator.shape, + matvec=lambda v: shift * v) + + shifted_eigenvalues, eigenvectors = spalg.eigs(shifted_operator, which='LM', k=how_many, ncv=50) + eigenvalues = shifted_eigenvalues - shift k = eigenvalues.argsort() return eigenvalues[k], eigenvectors[:, k] diff --git a/fdfd_tools/vectorization.py b/fdfd_tools/vectorization.py index 1c2134a..e7b9645 100644 --- a/fdfd_tools/vectorization.py +++ b/fdfd_tools/vectorization.py @@ -1,7 +1,7 @@ """ Functions for moving between a vector field (list of 3 ndarrays, [f_x, f_y, f_z]) and a 1D array representation of that field [f_x0, f_x1, f_x2,... f_y0,... f_z0,...]. -Vectorized versions of the field use column-major (ie., Fortran, Matlab) ordering. +Vectorized versions of the field use row-major (ie., C-style) ordering. """ From d09eff990f9160e70bc837174990305a7be791db Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 17 Dec 2017 20:51:34 -0800 Subject: [PATCH 19/29] Update Rayleigh quotient iteration to allow arbitrary linear operators --- fdfd_tools/eigensolvers.py | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/fdfd_tools/eigensolvers.py b/fdfd_tools/eigensolvers.py index 24d7339..e348cd7 100644 --- a/fdfd_tools/eigensolvers.py +++ b/fdfd_tools/eigensolvers.py @@ -33,10 +33,11 @@ def power_iteration(operator: sparse.spmatrix, return lm_eigval, v -def rayleigh_quotient_iteration(operator: sparse.spmatrix, +def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperator, guess_vector: numpy.ndarray, iterations: int = 40, tolerance: float = 1e-13, + solver=None, ) -> Tuple[complex, numpy.ndarray]: """ Use Rayleigh quotient iteration to refine an eigenvector guess. @@ -46,16 +47,33 @@ def rayleigh_quotient_iteration(operator: sparse.spmatrix, :param iterations: Maximum number of iterations to perform. Default 40. :param tolerance: Stop iteration if (A - I*eigenvalue) @ v < tolerance. Default 1e-13. + :param solver: Solver function of the form x = solver(A, b). + By default, use scipy.sparse.spsolve for sparse matrices and + scipy.sparse.bicgstab for general LinearOperator instances. :return: (eigenvalue, eigenvector) """ + try: + _test = operator - sparse.eye(operator.shape) + shift = lambda eigval: eigval * sparse.eye(operator.shape[0]) + if solver is None: + solver = spalg.spsolve + except TypeError: + shift = lambda eigval: spalg.LinearOperator(shape=operator.shape, + dtype=operator.dtype, + matvec=lambda v: eigval * v) + if solver is None: + solver = lambda A, b: spalg.bicgstab(A, b)[0] + v = guess_vector + v /= norm(v) for _ in range(iterations): - eigval = v.conj() @ operator @ v + eigval = v.conj() @ (operator @ v) if norm(operator @ v - eigval * v) < tolerance: break - v = spalg.spsolve(operator - eigval * sparse.eye(operator.shape[0]), v) - v /= norm(v) + shifted_operator = operator - shift(eigval) + v = solver(shifted_operator, v) + v /= norm(v) return eigval, v From 000cfabd788275be7de7ead2a4d85c8a95e734c6 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 17 Dec 2017 21:32:29 -0800 Subject: [PATCH 20/29] switch fft, ifft --- fdfd_tools/bloch.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index d7a2ffa..df63e6e 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -173,7 +173,7 @@ def maxwell_operator(k0: numpy.ndarray, m * hin_n) * k_mag # divide by epsilon - e_xyz = ifftn(fftn(d_xyz, axes=range(3)) / epsilon, axes=range(3)) + e_xyz = fftn(ifftn(d_xyz, axes=range(3)) / epsilon, axes=range(3)) # cross product and transform into mn basis b_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] * -k_mag @@ -187,7 +187,7 @@ def maxwell_operator(k0: numpy.ndarray, n * b_n[:, :, :, None]) # divide by mu - h_xyz = ifftn(fftn(b_xyz, axes=range(3)) / mu, axes=range(3)) + h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3)) # transform back to mn h_m = numpy.sum(h_xyz * m, axis=3) @@ -227,7 +227,7 @@ def hmn_2_exyz(k0: numpy.ndarray, m * hin_n) * k_mag # divide by epsilon - return [ei for ei in numpy.rollaxis(fftn(d_xyz, axes=range(3)) / epsilon, 3)] + return [ei for ei in numpy.rollaxis(ifftn(d_xyz, axes=range(3)) / epsilon, 3)] return operator @@ -258,7 +258,7 @@ def hmn_2_hxyz(k0: numpy.ndarray, hin_m, hin_n = [hi.reshape(shape) for hi in numpy.split(h, 2)] h_xyz = (m * hin_m + n * hin_n) - return [fftn(hi) for hi in numpy.rollaxis(h_xyz, 3)] + return [ifftn(hi) for hi in numpy.rollaxis(h_xyz, 3)] return operator @@ -312,7 +312,7 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray, n * hin_n[:, :, :, None]) # multiply by mu - b_xyz = ifftn(fftn(h_xyz, axes=range(3)) * mu, axes=range(3)) + b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3)) # transform back to mn b_m = numpy.sum(b_xyz * m, axis=3) @@ -323,7 +323,7 @@ def inverse_maxwell_operator_approx(k0: numpy.ndarray, m * b_n) / k_mag # multiply by epsilon - d_xyz = ifftn(fftn(e_xyz, axes=range(3)) * epsilon, axes=range(3)) + d_xyz = fftn(ifftn(e_xyz, axes=range(3)) * epsilon, axes=range(3)) # cross product and transform into mn basis crossinv_t2c h_m = numpy.sum(e_xyz * n, axis=3)[:, :, :, None] / +k_mag From 4a9596921f507e79a3ea59ac2d0d6129a49b3c08 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 17 Dec 2017 21:32:59 -0800 Subject: [PATCH 21/29] rename search_direction to direction --- fdfd_tools/bloch.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index df63e6e..33e9255 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -68,7 +68,7 @@ This module contains functions for generating and solving the k, f = find_k(frequency=1/1550, tolerance=(1/1550 - 1/1551), - search_direction=[1, 0, 0], + direction=[1, 0, 0], G_matrix=recip_lattice, epsilon=epsilon, band=0) @@ -369,7 +369,7 @@ def eigsolve(num_modes: int, def find_k(frequency: float, tolerance: float, - search_direction: numpy.ndarray, + direction: numpy.ndarray, G_matrix: numpy.ndarray, epsilon: field_t, mu: field_t = None, @@ -380,7 +380,7 @@ def find_k(frequency: float, :param frequency: Target frequency. :param tolerance: Target frequency tolerance. - :param search_direction: k-vector direction to search along. + :param direction: k-vector direction to search along. :param G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. :param epsilon: Dielectric constant distribution for the simulation. All fields are sampled at cell centers (i.e., NOT Yee-gridded) @@ -390,10 +390,10 @@ def find_k(frequency: float, return: (k, actual_frequency) The found k-vector and its frequency """ - search_direction = numpy.array(search_direction) / norm(search_direction) + direction = numpy.array(direction) / norm(direction) def get_f(k0_mag: float, band: int = 0): - k0 = search_direction * k0_mag + k0 = direction * k0_mag n, _v = eigsolve(band + 1, k0, G_matrix=G_matrix, epsilon=epsilon) f = numpy.sqrt(numpy.abs(numpy.real(n[band]))) return f @@ -401,6 +401,6 @@ def find_k(frequency: float, res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency), 0.25, method='Bounded', bounds=(0, 0.5), options={'xatol': abs(tolerance)}) - return res.x * search_direction, res.fun + frequency + return res.x * direction, res.fun + frequency From 39979edc440a7e3e547f81921b0bfba8bf8ba895 Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 17 Dec 2017 21:33:53 -0800 Subject: [PATCH 22/29] implement eigenvalue algorithm from Johnson paper. Could also use arpack + refinement, but that's also slow. --- fdfd_tools/bloch.py | 109 ++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 100 insertions(+), 9 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index 33e9255..8d85b19 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -47,12 +47,10 @@ This module contains functions for generating and solving the h_mn = fftn(1/mu * ifftn(b_m * m + b_n * n)) which forms the operator from the left side of the equation. - We can then use ARPACK in shift-invert mode (via scipy.linalg.eigs) - to find the eigenvectors for this operator. - - This approach is similar to the one used in MPB and derived at the start of - SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods + We can then use a preconditioned block Rayleigh iteration algorithm, as in + SG Johnson and JD Joannopoulos, Block-iterative frequency-domain methods for Maxwell's equations in a planewave basis, Optics Express 8, 3, 173-190 (2001) + (similar to that used in MPB) to find the eigenvectors for this operator. === @@ -76,14 +74,19 @@ This module contains functions for generating and solving the ''' from typing import List, Tuple, Callable, Dict +import logging import numpy from numpy.fft import fftn, ifftn, fftfreq import scipy +import scipy.optimize from scipy.linalg import norm import scipy.sparse.linalg as spalg +from .eigensolvers import rayleigh_quotient_iteration from . import field_t +logger = logging.getLogger(__name__) + def generate_kmn(k0: numpy.ndarray, G_matrix: numpy.ndarray, @@ -338,7 +341,8 @@ def eigsolve(num_modes: int, k0: numpy.ndarray, G_matrix: numpy.ndarray, epsilon: field_t, - mu: field_t = None + mu: field_t = None, + tolerance = 1e-8, ) -> Tuple[numpy.ndarray, numpy.ndarray]: """ Find the first (lowest-frequency) num_modes eigenmodes with Bloch wavevector @@ -355,15 +359,102 @@ def eigsolve(num_modes: int, """ h_size = 2 * epsilon[0].size + ''' + Generate the operators + ''' mop = maxwell_operator(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu) imop = inverse_maxwell_operator_approx(k0=k0, G_matrix=G_matrix, epsilon=epsilon, mu=mu) scipy_op = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=mop) scipy_iop = spalg.LinearOperator(dtype=complex, shape=(h_size, h_size), matvec=imop) - _eigvals, eigvecs = spalg.eigs(scipy_op, num_modes, sigma=0, OPinv=scipy_iop, which='LM') - eigvals = numpy.sum(eigvecs * (scipy_op @ eigvecs), axis=0) / numpy.sum(eigvecs * eigvecs, axis=0) - order = numpy.argsort(-eigvals) + y_shape = (h_size, num_modes) + + def rayleigh_quotient(Z: numpy.ndarray, approx_grad: bool = True): + """ + Absolute value of the block Rayleigh quotient, and the associated gradient. + + See Johnson and Joannopoulos, Opt. Expr. 8, 3 (2001) for details (full + citation in module docstring). + + === + + Notes on my understanding of the procedure: + + Minimize f(Y) = |trace((Y.H @ A @ Y)|, making use of Y = Z @ inv(Z.H @ Z)^(1/2) + (a polar orthogonalization of Y). This gives f(Z) = |trace(Z.H @ A @ Z @ U)|, + where U = inv(Z.H @ Z). We minimize the absolute value to find the eigenvalues + with smallest magnitude. + + The gradient is P @ (A @ Z @ U), where P = (1 - Z @ U @ Z.H) is a projection + onto the space orthonormal to Z. If approx_grad is True, the approximate + inverse of the maxwell operator is used to precondition the gradient. + """ + z = Z.reshape(y_shape) + U = numpy.linalg.inv(z.conj().T @ z) + zU = z @ U + AzU = scipy_op @ zU + zTAzU = z.conj().T @ AzU + f = numpy.real(numpy.trace(zTAzU)) + if approx_grad: + df_dy = scipy_iop @ (AzU - zU @ zTAzU) + else: + df_dy = (AzU - zU @ zTAzU) + return numpy.abs(f), numpy.sign(f) * df_dy.ravel() + + ''' + Use the conjugate gradient method and the approximate gradient calculation to + quickly find approximate eigenvectors. + ''' + result = scipy.optimize.minimize(rayleigh_quotient, + numpy.random.rand(*y_shape), + jac=True, + method='CG', + tol=1e-5, + options={'maxiter': 30, 'disp':True}) + + result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False), + result.x, + jac=True, + method='CG', + tol=1e-13, + options={'maxiter': 100, 'disp':True}) + + z = result.x.reshape(y_shape) + + ''' + Recover eigenvectors from Z + ''' + U = numpy.linalg.inv(z.conj().T @ z) + y = z @ scipy.linalg.sqrtm(U) + w = y.conj().T @ (scipy_op @ y) + + eigvals, w_eigvecs = numpy.linalg.eig(w) + eigvecs = y @ w_eigvecs + + for i in range(len(eigvals)): + v = eigvecs[:, i] + n = eigvals[i] + v /= norm(v) + logger.info('eigness {}: {}'.format(i, norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v ))) + + ev2 = eigvecs.copy() + for i in range(len(eigvals)): + logger.info('Refining eigenvector {}'.format(i)) + eigvals[i], ev2[:, i] = rayleigh_quotient_iteration(scipy_op, + guess_vector=eigvecs[:, i], + iterations=40, + tolerance=tolerance * numpy.real(numpy.sqrt(eigvals[i])) * 2, + solver = lambda A, b: spalg.bicgstab(A, b, maxiter=200)[0]) + eigvecs = ev2 + order = numpy.argsort(numpy.abs(eigvals)) + + for i in range(len(eigvals)): + v = eigvecs[:, i] + n = eigvals[i] + v /= norm(v) + logger.info('eigness {}: {}'.format(i, norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v ))) + return eigvals[order], eigvecs.T[order] From f312d735039668527a55f628cb548c342511b39a Mon Sep 17 00:00:00 2001 From: jan Date: Sun, 17 Dec 2017 22:55:55 -0800 Subject: [PATCH 23/29] Return real part of the gradient --- fdfd_tools/bloch.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index 8d85b19..2d343ce 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -400,7 +400,7 @@ def eigsolve(num_modes: int, df_dy = scipy_iop @ (AzU - zU @ zTAzU) else: df_dy = (AzU - zU @ zTAzU) - return numpy.abs(f), numpy.sign(f) * df_dy.ravel() + return numpy.abs(f), numpy.sign(f) * numpy.real(df_dy).ravel() ''' Use the conjugate gradient method and the approximate gradient calculation to From 16f97d7f6be0ccef2bf440bfa58e258857fe22c2 Mon Sep 17 00:00:00 2001 From: jan Date: Mon, 18 Dec 2017 00:13:29 -0800 Subject: [PATCH 24/29] Add ability to set bounds for find_k --- fdfd_tools/bloch.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index 2d343ce..ad7c555 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -464,7 +464,9 @@ def find_k(frequency: float, G_matrix: numpy.ndarray, epsilon: field_t, mu: field_t = None, - band: int = 0 + band: int = 0, + k_min: float = 0, + k_max: float = 0.5, ) -> Tuple[numpy.ndarray, float]: """ Search for a bloch vector that has a given frequency. @@ -489,8 +491,10 @@ def find_k(frequency: float, f = numpy.sqrt(numpy.abs(numpy.real(n[band]))) return f - res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency), 0.25, - method='Bounded', bounds=(0, 0.5), + res = scipy.optimize.minimize_scalar(lambda x: abs(get_f(x, band) - frequency), + (k_min + k_max) / 2, + method='Bounded', + bounds=(k_min, k_max), options={'xatol': abs(tolerance)}) return res.x * direction, res.fun + frequency From 85030448c3917ebeadc9b2c4a690ffc5d6b6f46b Mon Sep 17 00:00:00 2001 From: jan Date: Thu, 21 Dec 2017 20:11:30 -0800 Subject: [PATCH 25/29] Use L-BFGS instead of CG, and remove rayleigh iteration refinement scipy CG doesn't seem to converge as well as L-BFGS... worth looking into later --- fdfd_tools/bloch.py | 46 ++++++++++++++++++++++----------------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index ad7c555..739cc61 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -359,6 +359,8 @@ def eigsolve(num_modes: int, """ h_size = 2 * epsilon[0].size + kmag = norm(G_matrix @ k0) + ''' Generate the operators ''' @@ -409,16 +411,26 @@ def eigsolve(num_modes: int, result = scipy.optimize.minimize(rayleigh_quotient, numpy.random.rand(*y_shape), jac=True, - method='CG', - tol=1e-5, - options={'maxiter': 30, 'disp':True}) + method='L-BFGS-B', + tol=1e-20, + options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30}) + result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False), result.x, jac=True, - method='CG', - tol=1e-13, - options={'maxiter': 100, 'disp':True}) + method='L-BFGS-B', + tol=1e-20, + options={'maxiter': 2000, 'ptol':1e-18, 'disp':True}) + + for i in range(20): + result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False), + result.x, + jac=True, + method='L-BFGS-B', + tol=1e-20, + options={'maxiter': 70, 'gtol':1e-18, 'disp':True}) + z = result.x.reshape(y_shape) @@ -436,25 +448,13 @@ def eigsolve(num_modes: int, v = eigvecs[:, i] n = eigvals[i] v /= norm(v) - logger.info('eigness {}: {}'.format(i, norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v ))) + eigness = norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v ) + f = numpy.sqrt(-numpy.real(n)) + df = numpy.sqrt(-numpy.real(n + eigness)) + neff_err = kmag * (1/df - 1/f) + logger.info('eigness {}: {}\n neff_err: {}'.format(i, eigness, neff_err)) - ev2 = eigvecs.copy() - for i in range(len(eigvals)): - logger.info('Refining eigenvector {}'.format(i)) - eigvals[i], ev2[:, i] = rayleigh_quotient_iteration(scipy_op, - guess_vector=eigvecs[:, i], - iterations=40, - tolerance=tolerance * numpy.real(numpy.sqrt(eigvals[i])) * 2, - solver = lambda A, b: spalg.bicgstab(A, b, maxiter=200)[0]) - eigvecs = ev2 order = numpy.argsort(numpy.abs(eigvals)) - - for i in range(len(eigvals)): - v = eigvecs[:, i] - n = eigvals[i] - v /= norm(v) - logger.info('eigness {}: {}'.format(i, norm(scipy_op @ v - (v.conj() @ (scipy_op @ v)) * v ))) - return eigvals[order], eigvecs.T[order] From a70687f5e39b27bf9a579529f51f566fead577e1 Mon Sep 17 00:00:00 2001 From: jan Date: Thu, 21 Dec 2017 20:11:42 -0800 Subject: [PATCH 26/29] add bloch example --- examples/bloch.py | 68 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) create mode 100644 examples/bloch.py diff --git a/examples/bloch.py b/examples/bloch.py new file mode 100644 index 0000000..8bbd30e --- /dev/null +++ b/examples/bloch.py @@ -0,0 +1,68 @@ +import numpy, scipy, gridlock, fdfd_tools +from fdfd_tools import bloch +from numpy.linalg import norm +import logging + +logging.basicConfig(level=logging.DEBUG) +logger = logging.getLogger(__name__) + + +dx = 40 +x_period = 400 +y_period = z_period = 2000 +g = gridlock.Grid([numpy.arange(-x_period/2, x_period/2, dx), + numpy.arange(-1000, 1000, dx), + numpy.arange(-1000, 1000, dx)], + shifts=numpy.array([[0,0,0]]), + initial=1.445**2, + periodic=True) + +g.draw_cuboid([0,0,0], [200e8, 220, 220], eps=3.47**2) + +#x_period = y_period = z_period = 13000 +#g = gridlock.Grid([numpy.arange(3), ]*3, +# shifts=numpy.array([[0, 0, 0]]), +# initial=2.0**2, +# periodic=True) + +g2 = g.copy() +g2.shifts = numpy.zeros((6,3)) +g2.grids = [numpy.zeros(g.shape) for _ in range(6)] + +epsilon = [g.grids[0],] * 3 +reciprocal_lattice = numpy.diag(1e6/numpy.array([x_period, y_period, z_period])) #cols are vectors + +#print('Finding k at 1550nm') +#k, f = bloch.find_k(frequency=1/1550, +# tolerance=(1/1550 - 1/1551), +# direction=[1, 0, 0], +# G_matrix=reciprocal_lattice, +# epsilon=epsilon, +# band=0) +# +#print("k={}, f={}, 1/f={}, k/f={}".format(k, f, 1/f, norm(reciprocal_lattice @ k) / f )) + +print('Finding f at [0.25, 0, 0]') +for k0x in [.25]: + k0 = numpy.array([k0x, 0, 0]) + + kmag = norm(reciprocal_lattice @ k0) + tolerance = (1e6/1550) * 1e-4/1.5 # df = f * dn_eff / n + logger.info('tolerance {}'.format(tolerance)) + + n, v = bloch.eigsolve(4, k0, G_matrix=reciprocal_lattice, epsilon=epsilon, tolerance=tolerance) + v2e = bloch.hmn_2_exyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon) + v2h = bloch.hmn_2_hxyz(k0, G_matrix=reciprocal_lattice, epsilon=epsilon) + ki = bloch.generate_kmn(k0, reciprocal_lattice, g.shape) + + z = 0 + e = v2e(v[0]) + for i in range(3): + g2.grids[i] += numpy.real(e[i]) + g2.grids[i+3] += numpy.imag(e[i]) + + f = numpy.sqrt(numpy.real(numpy.abs(n))) # TODO + print('k0x = {:3g}\n eigval = {}\n f = {}\n'.format(k0x, n, f)) + n_eff = norm(reciprocal_lattice @ k0) / f + print('kmag/f = n_eff = {} \n wl = {}\n'.format(n_eff, 1/f )) + From 66712efd49d99fd547196933b4846e022461fbea Mon Sep 17 00:00:00 2001 From: jan Date: Wed, 27 Dec 2017 01:44:45 -0800 Subject: [PATCH 27/29] scipy L-BFGS silently converts to float, so view as floats when dealing with it.' --- fdfd_tools/bloch.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index 739cc61..b172e21 100644 --- a/fdfd_tools/bloch.py +++ b/fdfd_tools/bloch.py @@ -392,7 +392,7 @@ def eigsolve(num_modes: int, onto the space orthonormal to Z. If approx_grad is True, the approximate inverse of the maxwell operator is used to precondition the gradient. """ - z = Z.reshape(y_shape) + z = Z.view(dtype=complex).reshape(y_shape) U = numpy.linalg.inv(z.conj().T @ z) zU = z @ U AzU = scipy_op @ zU @@ -402,26 +402,35 @@ def eigsolve(num_modes: int, df_dy = scipy_iop @ (AzU - zU @ zTAzU) else: df_dy = (AzU - zU @ zTAzU) - return numpy.abs(f), numpy.sign(f) * numpy.real(df_dy).ravel() + + df_dy_flat = df_dy.view(dtype=float).ravel() + return numpy.abs(f), numpy.sign(f) * df_dy_flat ''' Use the conjugate gradient method and the approximate gradient calculation to quickly find approximate eigenvectors. ''' result = scipy.optimize.minimize(rayleigh_quotient, - numpy.random.rand(*y_shape), + numpy.random.rand(*y_shape, 2), jac=True, method='L-BFGS-B', tol=1e-20, options={'maxiter': 2000, 'gtol':0, 'ftol':1e-20 , 'disp':True})#, 'maxls':80, 'm':30}) + result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, True), + result.x, + jac=True, + method='L-BFGS-B', + tol=1e-20, + options={'maxiter': 2000, 'gtol':0, 'disp':True}) + result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False), result.x, jac=True, method='L-BFGS-B', tol=1e-20, - options={'maxiter': 2000, 'ptol':1e-18, 'disp':True}) + options={'maxiter': 2000, 'gtol':0, 'disp':True}) for i in range(20): result = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False), @@ -429,10 +438,13 @@ def eigsolve(num_modes: int, jac=True, method='L-BFGS-B', tol=1e-20, - options={'maxiter': 70, 'gtol':1e-18, 'disp':True}) + options={'maxiter': 70, 'gtol':0, 'disp':True}) + if result.nit == 0: + # We took 0 steps, so re-running won't help + break - z = result.x.reshape(y_shape) + z = result.x.view(dtype=complex).reshape(y_shape) ''' Recover eigenvectors from Z From 47dd0df8bc919c9289e9a6434eff5f2e2b4b4452 Mon Sep 17 00:00:00 2001 From: jan Date: Sat, 6 Jan 2018 13:51:42 -0800 Subject: [PATCH 28/29] fix operator test --- fdfd_tools/eigensolvers.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fdfd_tools/eigensolvers.py b/fdfd_tools/eigensolvers.py index e348cd7..d0ee541 100644 --- a/fdfd_tools/eigensolvers.py +++ b/fdfd_tools/eigensolvers.py @@ -53,7 +53,7 @@ def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperato :return: (eigenvalue, eigenvector) """ try: - _test = operator - sparse.eye(operator.shape) + _test = operator - sparse.eye(operator.shape[0]) shift = lambda eigval: eigval * sparse.eye(operator.shape[0]) if solver is None: solver = spalg.spsolve From 03fc9e6d700c7d7840ae9ebd66fcdc9d980c8e91 Mon Sep 17 00:00:00 2001 From: Jan Petykiewicz Date: Wed, 19 Feb 2020 19:07:57 -0800 Subject: [PATCH 29/29] deprecate --- README.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/README.md b/README.md index 5a3f49c..35a190c 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,11 @@ # fdfd_tools +** DEPRECATED ** + +The functionality in this module is now provided by [meanas](https://mpxd.net/code/jan/meanas). + +----------------------- + **fdfd_tools** is a python package containing utilities for creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD) electromagnetic simulations.