From 85030448c3917ebeadc9b2c4a690ffc5d6b6f46b Mon Sep 17 00:00:00 2001 From: jan Date: Thu, 21 Dec 2017 20:11:30 -0800 Subject: [PATCH 1/5] 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 2/5] 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 3/5] 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 4/5] 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 5/5] 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.