deprecate

scipy CG doesn't seem to converge as well as L-BFGS... worth looking
into later

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.

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 ))
+

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',
+                                     method='L-BFGS-B',
-                                     tol=1e-5,
+                                     tol=1e-20,
-                                     options={'maxiter': 30, 'disp':True})
+                                     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',
+                                     method='L-BFGS-B',
-                                     tol=1e-13,
+                                     tol=1e-20,
-                                     options={'maxiter': 100, 'disp':True})
+                                     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]
 
 

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