From 89c8ad03188a1edfc8babe9c66106f9335027fed Mon Sep 17 00:00:00 2001 From: Forgejo Actions Date: Wed, 28 Oct 2026 18:29:21 -0700 Subject: [PATCH] [bloch] solve transverse magnetic flux with a consistent energy metric --- docs/api/fdfd.md | 19 ++ meanas/fdfd/bloch.py | 356 ++++++++++++------------- meanas/test/test_bloch_foundations.py | 51 ++++ meanas/test/test_bloch_interactions.py | 38 ++- 4 files changed, 284 insertions(+), 180 deletions(-) diff --git a/docs/api/fdfd.md b/docs/api/fdfd.md index 20b126f..9b71a4e 100644 --- a/docs/api/fdfd.md +++ b/docs/api/fdfd.md @@ -14,6 +14,25 @@ ::: meanas.fdfd.farfield +## Bloch magnetic materials + +Bloch state vectors store the transverse Fourier coefficients of magnetic flux +`B = mu * H`. They equal `H` when permeability is omitted. Pass the same `mu` +to `bloch.hmn_2_hxyz(...)` and `bloch.hmn_2_exyz(...)` when converting solved +vectors to fields. The eigensolver uses the projected `mu**-1` energy metric +for finite positive real permeability, including anisotropic and spatially +varying profiles. Its preconditioner remains finite at the Gamma point. + +## Far-field normalization + +The near-to-far transform includes the cell area `dx * dy`, so normalized +far-field amplitudes describe a Fourier integral and converge under spatial +refinement. This changes amplitudes previously returned as unweighted FFT sums. +The inverse applies the matching Fourier quadrature and returns spacings for +the requested output shape. Both transforms center their origins at index +`size // 2`; padding and cropping preserve that sample. The inverse does not +modify input spectra. + ## SCPML conventions `stretch_with_scpml` now uses `-1` for the low-index face and `+1` for the diff --git a/meanas/fdfd/bloch.py b/meanas/fdfd/bloch.py index 1e15c3a..338f1c7 100644 --- a/meanas/fdfd/bloch.py +++ b/meanas/fdfd/bloch.py @@ -1,97 +1,44 @@ """ -Bloch eigenmode solver/operators +Bloch eigenmode solvers and operators in a spatial Fourier basis. -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 +The two stored coefficients at each wavevector represent transverse magnetic +flux B = mu H. They coincide with H when mu is omitted. Using B ensures +k dot B = 0 even in anisotropic or spatially varying magnetic materials; H +itself need not be transverse. The historical h_mn function names are retained. +Pass the same mu to hmn_2_hxyz and hmn_2_exyz when reconstructing these modes. - 1/mu * curl(1/eps * curl(H_eigenmode)) = (w/c)^2 H_eigenmode +For each Fourier wavevector, generate_kmn constructs unit vectors m and n with +k cross m = |k| n. The Cartesian Fourier field is B_k = b_m m + b_n n. +The Maxwell operator applies, in order: - into + H_k = fftn(ifftn(B_k) / mu) + E_k = fftn(ifftn(k cross H_k) / epsilon) + A(B_k) = k cross E_k - conv(1/mu_k, ik x conv(1/eps_k, ik x H_k)) = (w/c)^2 H_k +Material division is performed in real space and curls in Fourier space. +The explicit factors of i are omitted, giving the historical negative squared +frequency eigenvalues. Wavevectors and frequencies use cycles per length: +frequency = sqrt(-real(eigenvalue)). Reconstructed electric fields retain the +historical overall frequency scaling and are not power-normalized. - where: +For finite positive real permeability, the flux operator is self-adjoint in +the projected magnetic energy metric M = P mu^-1 P, where P projects onto the +transverse Fourier subspace. The block variational solve uses that metric in +its Rayleigh quotients, gradients, line searches, and final Ritz extraction. +The inverse-curl preconditioner regularizes zero wavevectors so static modes +remain accessible at Gamma. Nonzero-frequency modes have no static Fourier +flux component at Gamma; the two stored zero-wavevector basis vectors choose +two representatives of the static sector. - - 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 +Typical usage: - 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 as follows: - - - `h` is shorthand for `H_k` - - `(...)_xyz` denotes the `(x, y, z)` basis - - `(...)_kmn` denotes the `(k, m, n)` basis - - `hm` is the component of `h` in the `m` direction, etc. - - We know - - k @ h = kx hx + ky hy + kz hz = 0 = hk - h = hk + hm + hn = hm + hn - k = kk + km + kn = kk = |k| - - We can write - - k x h = (ky hz - kz hy, - kz hx - kx hz, - kx hy - ky hx)_xyz - = ((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 - - which gives us a straightforward way to perform the cross product - while simultaneously transforming into the `_kmn` basis. - We can also write - - 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)_xyz + (kk hm)(nx, ny, nz)_xyz - = |k| (hm * (nx, ny, nz)_xyz - - hn * (mx, my, mz)_xyz) - - which gives us a way to perform the cross product while simultaneously - trasnforming back into the `_xyz` basis. - - 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 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. - - === - - 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), - direction=[1, 0, 0], - G_matrix=recip_lattice, - epsilon=epsilon, - band=0) + values, vectors = eigsolve(5, k0, reciprocal_lattice, epsilon, mu=mu) + frequency = numpy.sqrt(-numpy.real(values[0])) + electric = hmn_2_exyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0]) + magnetic = hmn_2_hxyz(k0, reciprocal_lattice, epsilon, mu=mu)(vectors[0]) +The block iteration follows the approach of S. G. Johnson and J. D. +Joannopoulos, Optics Express 8, 173-190 (2001). """ from typing import Any, cast @@ -174,7 +121,7 @@ def generate_kmn( 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) + xy_non0 = numpy.any(k_xyz[:, :, :, :2] != 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] @@ -195,13 +142,20 @@ def maxwell_operator( mu: fdfield | None = None ) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]: """ - Generate the Maxwell operator + Generate the Maxwell operator on transverse magnetic-flux coefficients. - conv(1/mu_k, ik x conv(1/eps_k, ik x ___)) + k x conv(1/eps_k, k x conv(1/mu_k, ___)) - which is the spatial-frequency-space representation of + The corresponding real-space curl-curl problem is - 1/mu * curl(1/eps * curl(___)) + curl(1/eps * curl(1/mu * ___)) + + with the overall factor ``-(2*pi)**2`` omitted for the module's + cycles-per-length frequency convention. + + The state represents B = mu H, whose Fourier components are transverse. + It equals H when mu is omitted. Eigenvalues retain the historical negative + squared-frequency convention (the two explicit factors of i are omitted). The operator is a function that acts on a vector h_mn of size `2 * epsilon[0].size` @@ -230,23 +184,28 @@ def maxwell_operator( """ Maxwell operator for Bloch eigenmode simulation. - h is complex 2-field in (m, n) basis, vectorized + h contains the two transverse magnetic-flux coefficients. Args: - h: Raveled h_mn; size `2 * epsilon[0].size`. - Altered in-place. + h: Raveled b_mn; size `2 * epsilon[0].size`. Not modified. Returns: - Raveled conv(1/mu_k, ik x conv(1/eps_k, ik x h_mn)), returned - and overwritten in-place of `h`. + Raveled transverse flux image under the Maxwell operator. """ 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 + # Convert transverse B into H before taking either curl. H itself need + # not be transverse for anisotropic or spatially varying permeability. + if mu is None: + h_m, h_n = hin_m, hin_n + else: + b_xyz = m * hin_m + n * hin_n + h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3)) + h_m = numpy.sum(h_xyz * m, axis=3, keepdims=True) + h_n = numpy.sum(h_xyz * n, axis=3, keepdims=True) + d_xyz = (n * h_m - m * h_n) * k_mag # divide by epsilon temp = ifftn(d_xyz, axes=range(3)) # reuses d_xyz if using pyfftw @@ -257,23 +216,7 @@ def maxwell_operator( b_m = numpy.sum(e_xyz * n, axis=3, keepdims=True) * -k_mag b_n = numpy.sum(e_xyz * m, axis=3, keepdims=True) * +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 - + n * b_n) # noqa - - # divide by mu - temp = ifftn(b_xyz, axes=range(3)) - temp /= mu - h_xyz = fftn(temp, 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 _assemble_hmn_vector(h_m, h_n) + return _assemble_hmn_vector(b_m, b_n) return operator @@ -282,12 +225,12 @@ def hmn_2_exyz( k0: ArrayLike, G_matrix: ArrayLike, epsilon: fdfield, + mu: fdfield | None = None, ) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]: """ - Generate an operator which converts a vectorized spatial-frequency-space - `h_mn` into an E-field distribution, i.e. + Convert transverse Fourier flux coefficients into an unnormalized E field. - ifft(conv(1/eps_k, ik x h_mn)) + ifft(conv(1/eps_k, k x conv(1/mu_k, b_mn))) The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`. @@ -298,17 +241,26 @@ def hmn_2_exyz( G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. epsilon: Dielectric constant distribution for the simulation. All fields are sampled at cell centers (i.e., NOT Yee-gridded) + mu: Permeability used by the solve. The input coefficients describe B; + divide by mu to recover H before computing the electric field. Returns: Function for converting `h_mn` into `E_xyz` """ shape = epsilon[0].shape + (1,) epsilon = numpy.moveaxis(epsilon, 0, -1) + if mu is not None: + mu = numpy.moveaxis(mu, 0, -1) k_mag, m, n = generate_kmn(k0, G_matrix, shape) def operator(h: NDArray[numpy.complex128]) -> cfdfield_t: hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2)) + if mu is not None: + b_xyz = m * hin_m + n * hin_n + h_xyz = fftn(ifftn(b_xyz, axes=range(3)) / mu, axes=range(3)) + hin_m = numpy.sum(h_xyz * m, axis=3, keepdims=True) + hin_n = numpy.sum(h_xyz * n, axis=3, keepdims=True) d_xyz = (n * hin_m - m * hin_n) * k_mag @@ -323,12 +275,12 @@ def hmn_2_hxyz( k0: ArrayLike, G_matrix: ArrayLike, epsilon: fdfield, + mu: fdfield | None = None, ) -> Callable[[NDArray[numpy.complex128]], cfdfield_t]: """ - Generate an operator which converts a vectorized spatial-frequency-space - `h_mn` into an H-field distribution, i.e. + Convert transverse Fourier flux coefficients into a physical H field. - ifft(h_mn) + ifft(b_mn) / mu The operator is a function that acts on a vector `h_mn` of size `2 * epsilon[0].size`. @@ -339,6 +291,8 @@ def hmn_2_hxyz( G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. epsilon: Dielectric constant distribution for the simulation. Only `epsilon[0].shape` is used. + mu: Permeability used by the solve. Input coefficients describe B; + the returned physical H is B / mu. Default is unit permeability. Returns: Function for converting `h_mn` into `H_xyz` @@ -350,7 +304,10 @@ def hmn_2_hxyz( hin_m, hin_n = (hi.reshape(shape) for hi in numpy.split(h, 2)) h_xyz = (m * hin_m + n * hin_n) - return cfdfield_t(numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)])) + field = numpy.array([ifftn(hi) for hi in numpy.moveaxis(h_xyz, 3, 0)]) + if mu is not None: + field /= mu + return cfdfield_t(field) return operator @@ -362,11 +319,16 @@ def inverse_maxwell_operator_approx( mu: fdfield | None = None, ) -> Callable[[NDArray[numpy.complex128]], NDArray[numpy.complex128]]: """ - Generate an approximate inverse of the Maxwell operator, + Generate a positive approximate inverse of minus the flux Maxwell operator, - ik x conv(eps_k, ik x conv(mu_k, ___)) + project_mn(conv(mu_k, cross_inv(conv(eps_k, cross_inv(___))))) - which can be used to improve the speed of ARPACK in shift-invert mode. + At k = 0, inverse crosses use the smallest nonzero grid wavevector scale + (or one on a single-point grid). This finite regularization keeps the + preconditioner positive on static directions, which must remain accessible + to the eigensolver at Gamma. + + It is used for the preconditioned block variational solve. See the `meanas.fdfd.bloch` docstring for more information. @@ -387,6 +349,9 @@ def inverse_maxwell_operator_approx( mu = numpy.moveaxis(mu, 0, -1) k_mag, m, n = generate_kmn(k0, G_matrix, shape) + nonzero = k_mag[k_mag != 0] + reference = float(nonzero.min()) if nonzero.size else 1.0 + k_inverse = numpy.divide(1, k_mag, out=numpy.full_like(k_mag, 1 / reference), where=k_mag != 0) def operator(h: NDArray[numpy.complex128]) -> NDArray[numpy.complex128]: """ @@ -398,31 +363,14 @@ def inverse_maxwell_operator_approx( h: Raveled h_mn; size `2 * epsilon[0].size`. Returns: - Raveled ik x conv(eps_k, ik x conv(mu_k, h_mn)) + Raveled, regularized approximate inverse flux image. """ 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 - + n * hin_n) # noqa - - # multiply by mu - temp = ifftn(h_xyz, axes=range(3)) - temp *= mu - b_xyz = fftn(temp, axes=range(3)) - - # transform back to mn - b_m = numpy.sum(b_xyz * m, axis=3, keepdims=True) - b_n = numpy.sum(b_xyz * n, axis=3, keepdims=True) - # cross product and transform into xyz basis - e_xyz = (n * b_m - - m * b_n) / k_mag + e_xyz = (n * hin_m - m * hin_n) * k_inverse # multiply by epsilon temp = ifftn(e_xyz, axes=range(3)) @@ -430,14 +378,41 @@ def inverse_maxwell_operator_approx( d_xyz = fftn(temp, axes=range(3)) # cross product and transform into mn basis crossinv_t2c - h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) / +k_mag - h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) / -k_mag + h_m = numpy.sum(d_xyz * n, axis=3, keepdims=True) * k_inverse + h_n = numpy.sum(d_xyz * m, axis=3, keepdims=True) * -k_inverse + if mu is not None: + h_xyz = m * h_m + n * h_n + b_xyz = fftn(ifftn(h_xyz, axes=range(3)) * mu, axes=range(3)) + h_m = numpy.sum(b_xyz * m, axis=3, keepdims=True) + h_n = numpy.sum(b_xyz * n, axis=3, keepdims=True) return _assemble_hmn_vector(h_m, h_n) return operator +def _transverse_material_operator( + k0: ArrayLike, + G_matrix: ArrayLike, + material: fdfield, + *, + inverse: bool, + ) -> Callable[[NDArray], NDArray]: + """Apply P material^(+/-1) P on the transverse Fourier subspace.""" + shape = material[0].shape + (1,) + _k_mag, m, n = generate_kmn(k0, G_matrix, shape) + values = numpy.moveaxis(material, 0, -1) + + def apply(h: NDArray) -> NDArray: + hm, hn = (part.reshape(shape) for part in numpy.split(h, 2)) + field = ifftn(m * hm + n * hn, axes=range(3)) + field = field / values if inverse else field * values + transformed = fftn(field, axes=range(3)) + return _assemble_hmn_vector(numpy.sum(transformed * m, axis=3), numpy.sum(transformed * n, axis=3)) + + return apply + + def find_k( frequency: float, tolerance: float, @@ -528,18 +503,20 @@ def eigsolve( G_matrix: 3x3 matrix, with reciprocal lattice vectors as columns. epsilon: Dielectric constant distribution for the simulation. All fields are sampled at cell centers (i.e., NOT Yee-gridded) - mu: Magnetic permability distribution for the simulation. - Default `None` (1 everywhere). - tolerance: Solver stops when fractional change in the objective - `trace(Z.H @ A @ Z @ inv(Z Z.H))` is smaller than the tolerance - max_iters: TODO - reset_iters: TODO - callback: TODO - y0: TODO, initial guess + mu: Finite positive real diagonal permeability. Default is one. + tolerance: Relative projected-gradient tolerance in the magnetic energy + metric. Flat objective traces trigger a conjugate-direction restart. + max_iters: Maximum iteration count; returns Ritz pairs from the final + subspace if the limit is reached. + reset_iters: Periodic conjugate-direction restart interval. + callback: Optional no-argument callback after each subspace update. + y0: Initial flux coefficients with shape (num_modes, 2 * grid_size). + Copied before use. Returns: `(eigenvalues, eigenvectors)` where `eigenvalues[i]` corresponds to the - vector `eigenvectors[i, :]` + transverse flux vector `eigenvectors[i, :]`. Pass mu to the field + converters when it was supplied to this solve. """ k0 = numpy.asarray(k0) @@ -556,6 +533,24 @@ def eigsolve( 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) + metric = inverse_metric = None + if mu is not None: + if numpy.iscomplexobj(mu) or not numpy.isfinite(mu).all() or numpy.any(mu <= 0): + raise ValueError('Bloch variational solving requires finite positive real permeability') + metric = spalg.LinearOperator( + shape=(h_size, h_size), dtype=complex, + matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=True), + ) + inverse_metric = spalg.LinearOperator( + shape=(h_size, h_size), dtype=complex, + matvec=_transverse_material_operator(k0, G_matrix, mu, inverse=False), + ) + + def mass(vectors: NDArray) -> NDArray: + # P mu^-1 P makes the flux operator self-adjoint in the magnetic + # energy inner product, even though it is not Euclidean-Hermitian. + return vectors if metric is None else metric @ vectors + y_shape = (h_size, num_modes) prev_E = 0.0 @@ -569,11 +564,11 @@ def eigsolve( if y0 is None: Z = rng.random(y_shape) + 1j * rng.random(y_shape) else: - Z = numpy.asarray(y0).T + Z = numpy.array(y0, dtype=complex, copy=True).T while True: Z *= num_modes / norm(Z) - ZtZ = Z.conj().T @ Z + ZtZ = Z.conj().T @ mass(Z) try: U = numpy.linalg.inv(ZtZ) except numpy.linalg.LinAlgError: @@ -592,27 +587,28 @@ def eigsolve( for i in range(max_iters): Zt = numpy.conj(Z.T, out=Zt) - ZtZ = Zt @ Z + ZtZ = Zt @ mass(Z) U = numpy.linalg.inv(ZtZ) AZ = scipy_op @ Z.copy() - ZtAZ = Zt @ AZ + ZtAZ = Zt @ mass(AZ) ZtAZU = ZtAZ @ U E_signed = real(trace(ZtAZU)) sgn = numpy.sign(E_signed) E = numpy.abs(E_signed) # G = AZU projected onto the space orthonormal to Z via (1 - ZUZt) - G = (AZ @ U - Z @ U @ ZtAZU) * sgn + G = mass(AZ @ U - Z @ U @ ZtAZU) * sgn if i > 0 and abs(E - prev_E) < tolerance * 0.5 * (E + prev_E + 1e-7): - logger.info( - 'Optimization succeded: ' - f'[change in trace] {abs(E - prev_E)} - 5e-8 ' - f'< {tolerance} [tolerance] * {(E + prev_E) / 2} [value of trace]' - ) - break + # A flat line search alone does not imply a stationary subspace. + # Restart with steepest descent if conjugacy has been lost. + prev_traceGtKG = 0 - KG = scipy_iop @ G # Preconditioned steepest descent direction + if norm(G) <= tolerance * max(norm(mass(AZ) @ U), 1e-7): + break + # Approximate the inverse energy Hessian by a symmetric sandwich of + # P mu P around the inverse-curl / epsilon preconditioner. + KG = scipy_iop @ (G if inverse_metric is None else inverse_metric @ G) traceGtKG = _rtrace_AtB(G, KG) if prev_traceGtKG == 0 or i % reset_iters == 0: @@ -642,11 +638,11 @@ def eigsolve( # Qi = inv(Q) = U' AD = scipy_op @ D.copy() - DtD = D.conj().T @ D - DtAD = D.conj().T @ AD + DtD = D.conj().T @ mass(D) + DtAD = D.conj().T @ mass(AD) - symZtD = _symmetrize(Zt @ D) - symZtAD = _symmetrize(Zt @ AD) + symZtD = _symmetrize(Zt @ mass(D)) + symZtAD = _symmetrize(Zt @ mass(AD)) Qi_memo: list[float | None] = [None, None] @@ -737,8 +733,10 @@ def eigsolve( result = scipy.optimize.minimize_scalar( trace_func, method='bounded', - bounds=(0, pi), - options={'xatol': tolerance}, + # Center the period on zero so the line search can resolve small + # steps without the relative angular floor near theta = pi. + bounds=(-pi / 2, pi / 2), + options={'xatol': min(tolerance, 1e-12)}, ) new_E = result.fun theta = result.x @@ -757,9 +755,9 @@ def eigsolve( # # Recover eigenvectors from Z # - U = numpy.linalg.inv(ZtZ) + U = numpy.linalg.inv(Z.conj().T @ mass(Z)) Y = Z @ scipy.linalg.sqrtm(U).astype(numpy.complex128) - W = Y.conj().T @ (scipy_op @ Y.copy()) + W = Y.conj().T @ mass(scipy_op @ Y.copy()) eigvals, W_eigvecs = numpy.linalg.eig(W) eigvecs = Y @ W_eigvecs @@ -769,10 +767,10 @@ def eigsolve( n = eigvals[i] v /= norm(v) Av = numpy.asarray(scipy_op @ v.copy()).reshape(-1) - eigness = norm(Av - (v.conj() @ Av) * v) - f = numpy.sqrt(-numpy.real(n)) - df = numpy.sqrt(-numpy.real(n) + eigness) - neff_err = kmag * (1 / df - 1 / f) + eigness = norm(Av - n * v) + f = numpy.sqrt(max(0.0, -numpy.real(n))) + df = numpy.sqrt(max(0.0, -numpy.real(n) + eigness)) + neff_err = kmag * (1 / df - 1 / f) if f > 0 and df > 0 else numpy.inf logger.info(f'eigness {i}: {eigness}\n neff_err: {neff_err}') order = numpy.argsort(numpy.abs(eigvals)) diff --git a/meanas/test/test_bloch_foundations.py b/meanas/test/test_bloch_foundations.py index 02a547c..0c680ff 100644 --- a/meanas/test/test_bloch_foundations.py +++ b/meanas/test/test_bloch_foundations.py @@ -1,10 +1,61 @@ import numpy +import pytest from numpy.linalg import norm from ..fdfd import bloch from ._bloch_case import EPSILON, G_MATRIX, H_MN, K0_AXIAL, K0_GENERAL, MU, SHAPE, ZERO_H_MN from .utils import assert_close + +@pytest.mark.parametrize('k0', [[0.2, 0, 0], [0, 0.2, 0], [0, 0, 0.2], [0, 0.2, 0.1], [0, 0, 0]]) +def test_kmn_basis_is_transverse_on_axes_and_at_gamma(k0) -> None: + k0 = numpy.asarray(k0) + shape = (2, 3, 2) + _magnitude, m, n = bloch.generate_kmn(k0, numpy.eye(3), shape) + reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid( + *(numpy.fft.fftfreq(size, 1 / size) for size in shape), indexing='ij', + )), 0, -1) + k = k0 - reciprocal_indices + assert_close(numpy.sum(k * m, axis=-1), 0, atol=1e-14) + assert_close(numpy.sum(k * n, axis=-1), 0, atol=1e-14) + assert_close(numpy.sum(m * n, axis=-1), 0, atol=1e-14) + assert_close(norm(m, axis=-1), 1) + assert_close(norm(n, axis=-1), 1) + + +@pytest.mark.parametrize('mu', [None, MU]) +def test_inverse_preconditioner_is_finite_at_gamma(mu) -> None: + operator = bloch.inverse_maxwell_operator_approx(numpy.zeros(3), G_MATRIX, EPSILON, mu) + with numpy.errstate(divide='raise', invalid='raise', over='raise'): + result = operator(H_MN) + zero = operator(ZERO_H_MN) + assert numpy.isfinite(result).all() + assert_close(zero, 0) + + +def test_anisotropic_magnetic_bloch_modes_obey_full_fourier_maxwell_equations() -> None: + shape = EPSILON.shape[1:] + size = 2 * numpy.prod(shape) + operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU) + matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(size, dtype=complex)]) + values, vectors = numpy.linalg.eig(matrix) + reciprocal_indices = numpy.moveaxis(numpy.array(numpy.meshgrid( + *(numpy.fft.fftfreq(n, 1 / n) for n in shape), indexing='ij', + )), 0, -1) + k = K0_GENERAL - reciprocal_indices + epsilon, mu = numpy.moveaxis(EPSILON, 0, -1), numpy.moveaxis(MU, 0, -1) + h_converter = bloch.hmn_2_hxyz(K0_GENERAL, G_MATRIX, EPSILON, MU) + e_converter = bloch.hmn_2_exyz(K0_GENERAL, G_MATRIX, EPSILON, MU) + for value, vector in zip(values, vectors.T, strict=True): + h = numpy.moveaxis(h_converter(vector), 0, -1) + b = mu * h + h_fourier = numpy.fft.fftn(h, axes=(0, 1, 2)) + electric = numpy.fft.ifftn(numpy.cross(k, h_fourier), axes=(0, 1, 2)) / epsilon + curl_electric = numpy.fft.ifftn(numpy.cross(k, numpy.fft.fftn(electric, axes=(0, 1, 2))), axes=(0, 1, 2)) + assert_close(curl_electric, value * b, atol=1e-13, rtol=1e-11) + assert_close(numpy.sum(k * numpy.fft.fftn(b, axes=(0, 1, 2)), axis=-1), 0, atol=1e-13) + assert_close(numpy.moveaxis(e_converter(vector), 0, -1), electric, atol=1e-13) + def test_generate_kmn_general_case_returns_orthonormal_basis() -> None: k_mag, m_vecs, n_vecs = bloch.generate_kmn(K0_GENERAL, G_MATRIX, SHAPE) diff --git a/meanas/test/test_bloch_interactions.py b/meanas/test/test_bloch_interactions.py index 0628a55..822e11c 100644 --- a/meanas/test/test_bloch_interactions.py +++ b/meanas/test/test_bloch_interactions.py @@ -4,7 +4,7 @@ from numpy.testing import assert_allclose from types import SimpleNamespace from ..fdfd import bloch -from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, Y0, Y0_TWO_MODE, build_overlap_fixture +from ._bloch_case import EPSILON, G_MATRIX, H_SIZE, K0_X, K0_GENERAL, MU, Y0, Y0_TWO_MODE, build_overlap_fixture from .utils import assert_close @@ -81,6 +81,42 @@ def test_eigsolve_returns_finite_modes_with_small_residual() -> None: assert callback_count > 0 +@pytest.mark.parametrize('count', [1, 2]) +def test_iteration_limit_returns_consistent_ritz_pairs_without_mutating_guess(count) -> None: + rng = numpy.random.default_rng(11) + guess = rng.normal(size=(count, H_SIZE)) + 1j * rng.normal(size=(count, H_SIZE)) + original = guess.copy() + values, vectors = bloch.eigsolve(count, K0_X, G_MATRIX, EPSILON, y0=guess, max_iters=1) + operator = bloch.maxwell_operator(K0_X, G_MATRIX, EPSILON) + for value, vector in zip(values, vectors, strict=True): + rayleigh = numpy.vdot(vector, operator(vector).ravel()) / numpy.vdot(vector, vector) + assert_close(value, rayleigh, atol=1e-13) + assert_allclose(vectors.conj() @ vectors.T, numpy.eye(count), atol=1e-12) + numpy.testing.assert_array_equal(guess, original) + + +@pytest.mark.parametrize('mu', [None, MU]) +def test_eigsolve_at_gamma_returns_a_static_mode_with_small_residual(mu) -> None: + with numpy.errstate(divide='raise', invalid='raise', over='raise'): + values, vectors = bloch.eigsolve(1, numpy.zeros(3), G_MATRIX, EPSILON, mu, y0=Y0, max_iters=200, tolerance=1e-9) + residual = bloch.maxwell_operator(numpy.zeros(3), G_MATRIX, EPSILON, mu)(vectors[0]).ravel() - values[0] * vectors[0] + assert abs(values[0]) < 1e-12 + assert numpy.linalg.norm(residual) < 1e-9 + + +def test_anisotropic_magnetic_eigensolve_matches_dense_reference() -> None: + operator = bloch.maxwell_operator(K0_GENERAL, G_MATRIX, EPSILON, MU) + matrix = numpy.column_stack([operator(column).ravel() for column in numpy.eye(H_SIZE, dtype=complex)]) + reference = numpy.linalg.eigvals(matrix) + reference = reference[numpy.argsort(numpy.abs(reference))[:3]] + rng = numpy.random.default_rng(9) + guess = rng.normal(size=(3, H_SIZE)) + 1j * rng.normal(size=(3, H_SIZE)) + values, vectors = bloch.eigsolve(3, K0_GENERAL, G_MATRIX, EPSILON, MU, y0=guess, max_iters=300, tolerance=1e-10) + assert_allclose(values, reference, atol=1e-11, rtol=1e-9) + for value, vector in zip(values, vectors, strict=True): + assert numpy.linalg.norm(matrix @ vector - value * vector) < 1e-9 + + def test_eigsolve_without_initial_guess_returns_finite_modes() -> None: eigvals, eigvecs = bloch.eigsolve( 1,