Compare commits

...

5 Commits

Author SHA1 Message Date
03fc9e6d70 deprecate 2020-02-19 19:10:18 -08:00
jan
47dd0df8bc fix operator test 2018-01-08 11:59:00 -08:00
jan
66712efd49 scipy L-BFGS silently converts to float, so view as floats when dealing with it.' 2017-12-27 01:44:45 -08:00
jan
a70687f5e3 add bloch example 2017-12-21 20:11:42 -08:00
jan
85030448c3 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
2017-12-21 20:11:30 -08:00
4 changed files with 114 additions and 28 deletions

View File

@ -1,5 +1,11 @@
# fdfd_tools # 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 **fdfd_tools** is a python package containing utilities for
creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD) creating and analyzing 2D and 3D finite-difference frequency-domain (FDFD)
electromagnetic simulations. electromagnetic simulations.

68
examples/bloch.py Normal file
View File

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

View File

@ -359,6 +359,8 @@ def eigsolve(num_modes: int,
""" """
h_size = 2 * epsilon[0].size h_size = 2 * epsilon[0].size
kmag = norm(G_matrix @ k0)
''' '''
Generate the operators 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 onto the space orthonormal to Z. If approx_grad is True, the approximate
inverse of the maxwell operator is used to precondition the gradient. 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) U = numpy.linalg.inv(z.conj().T @ z)
zU = z @ U zU = z @ U
AzU = scipy_op @ zU AzU = scipy_op @ zU
@ -400,27 +402,49 @@ def eigsolve(num_modes: int,
df_dy = scipy_iop @ (AzU - zU @ zTAzU) df_dy = scipy_iop @ (AzU - zU @ zTAzU)
else: else:
df_dy = (AzU - zU @ zTAzU) 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 Use the conjugate gradient method and the approximate gradient calculation to
quickly find approximate eigenvectors. quickly find approximate eigenvectors.
''' '''
result = scipy.optimize.minimize(rayleigh_quotient, result = scipy.optimize.minimize(rayleigh_quotient,
numpy.random.rand(*y_shape), numpy.random.rand(*y_shape, 2),
jac=True, 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 = scipy.optimize.minimize(lambda y: rayleigh_quotient(y, False),
result.x, result.x,
jac=True, 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 Recover eigenvectors from Z
@ -436,25 +460,13 @@ def eigsolve(num_modes: int,
v = eigvecs[:, i] v = eigvecs[:, i]
n = eigvals[i] n = eigvals[i]
v /= norm(v) 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)) 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] return eigvals[order], eigvecs.T[order]

View File

@ -53,7 +53,7 @@ def rayleigh_quotient_iteration(operator: sparse.spmatrix or spalg.LinearOperato
:return: (eigenvalue, eigenvector) :return: (eigenvalue, eigenvector)
""" """
try: try:
_test = operator - sparse.eye(operator.shape) _test = operator - sparse.eye(operator.shape[0])
shift = lambda eigval: eigval * sparse.eye(operator.shape[0]) shift = lambda eigval: eigval * sparse.eye(operator.shape[0])
if solver is None: if solver is None:
solver = spalg.spsolve solver = spalg.spsolve