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. 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 )) + diff --git a/fdfd_tools/bloch.py b/fdfd_tools/bloch.py index ad7c555..b172e21 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 ''' @@ -390,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 @@ -400,27 +402,49 @@ 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='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, 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='CG', - tol=1e-13, - options={'maxiter': 100, 'disp':True}) + method='L-BFGS-B', + tol=1e-20, + options={'maxiter': 2000, 'gtol':0, 'disp':True}) - z = result.x.reshape(y_shape) + 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':0, 'disp':True}) + if result.nit == 0: + # We took 0 steps, so re-running won't help + break + + + z = result.x.view(dtype=complex).reshape(y_shape) ''' Recover eigenvectors from Z @@ -436,25 +460,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] 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