Solver Methods

PhoXonic.jl provides multiple solver methods for different problem sizes:

Memory and Computational Cost

N = number of plane waves = basis.num_pw × ncomponents(wave)

For large-scale calculations with N > 10,000, see Matrix-Free Methods.

Recommendations

Note: For 3D photonic crystals, TransverseEM is strongly recommended over FullVectorEM. See TransverseEM vs FullVectorEM for why.

Note: For 3D phononic crystals, FullElastic is an experimental implementation. See FullElastic (Experimental) for why parameter tuning is challenging.

DenseMethod

The default solver using LAPACK's dense eigenvalue decomposition.

Basic Usage

solver = Solver(TEWave(), geo, (64, 64), DenseMethod(); cutoff=7)

DenseMethod for 3D

For small 3D systems, DenseMethod with shift > 0 can be used:

solver = Solver(FullVectorEM(), geo, (12, 12, 12), DenseMethod(shift=0.01); cutoff=5)  # Use cutoff≥7 for high ε

How it works:

Unlike iterative methods which use shift-and-invert transformation, DenseMethod(shift=σ) computes all eigenvalues using LAPACK, then filters out eigenvalues where ω² < σ. This post-hoc filtering effectively removes the spurious longitudinal modes at ω ≈ 0.

Important difference from iterative methods:

This means Dense and iterative methods may return different bands for the same bands=1:n request if there are intermediate eigenvalues. Both results are valid eigenvalues, just in different order.

When to use DenseMethod for 3D:

• Small systems (N × 3 < 2000)
• Validation and debugging
• Need all eigenvalues in ascending order

KrylovKitMethod

Iterative solver based on KrylovKit.jl using Arnoldi iteration.

Basic Usage

solver = Solver(TEWave(), geo, (128, 128), KrylovKitMethod(); cutoff=15)

Shift-and-Invert for 3D

For 3D calculations, use shift-and-invert to skip longitudinal modes:

solver = Solver(FullVectorEM(), geo, (16, 16, 16), KrylovKitMethod(shift=0.01); cutoff=7)

The shift parameter σ transforms the generalized eigenvalue problem:

Eigenvalues closest to σ become the largest |μ|, making them easiest for iterative solvers to find.

Why is Shift Needed for 3D?

The 3D H-field formulation does not explicitly enforce ∇·H = 0, producing:

• Longitudinal modes: λ ≈ 0 (unphysical, violate transversality)
• Transverse modes: λ > 0 (physical electromagnetic waves)

Without shift, iterative solvers converge to the spurious λ ≈ 0 modes first. With shift σ > 0, eigenvalues λ < σ are effectively filtered out.

Choosing σ

Guidelines for choosing σ:

1. Lower bound: σ must be larger than the longitudinal mode eigenvalues (≈ 0). A small positive value like σ = 0.01 typically works.

2. Upper bound: σ should be smaller than the smallest physical eigenvalue you want to find. For a homogeneous medium: λmin ≈ |k|²/εmax.

   Example: |k| = 0.5, ε = 4 → λ_min ≈ 0.0625. Use σ < 0.06.

3. Safe default: For most 3D photonic calculations, σ = 0.01 works well when |k| ≥ 0.1 (in units of 2π/a).

4. Near k = 0: At the Γ point (k = 0), the lowest transverse modes also have λ → 0. Use a very small shift (σ ~ 0.001) or skip the Γ point.

5. High-contrast materials: For large ε, eigenvalues scale as 1/ε. Reduce σ proportionally: σ ≈ 0.01 / ε_max.

Troubleshooting

Phononic Eigenvalue Scaling

Phononic eigenvalues ω² are typically O(10¹⁰) for common materials (steel, epoxy, etc.), which can cause numerical instability in iterative solvers. PhoXonic.jl automatically scales the eigenvalue problem for SHWave, PSVWave, FullElastic, and Longitudinal1D:

Scaled problem: (A/s) x = (λ/s) B x
where s = c² × |k|² and c² = C₁₁/ρ (longitudinal wave speed squared)

This normalizes eigenvalues to O(1), ensuring stable convergence. The scaling is:

• Automatic: No user configuration needed
• Transparent: Results are unscaled before return
• Photonic-safe: Scale factor = 1.0 for TE/TM/FullVectorEM

Example: Steel/epoxy at |k| = 100 rad/m:

• Without scaling: ω² ~ 10¹⁰ → KrylovKit may fail
• With scaling: ω²/s ~ 1 → stable convergence

KrylovKitMethod Parameters

method = KrylovKitMethod(
    tol = 1e-8,       # Convergence tolerance
    maxiter = 300,    # Maximum iterations
    krylovdim = 30,   # Krylov subspace dimension
    verbosity = 0,    # 0=silent, 1=warnings, 2=info
    shift = 0.0,      # Spectral shift (0 = no shift)
    matrix_free = false  # Use matrix-free operators for shift-and-invert
)

Matrix-Free Shift-and-Invert

For large 3D calculations where memory is limited, the matrix_free=true option avoids building dense matrices:

# Memory-efficient for large 3D systems
solver = Solver(
    FullVectorEM(), geo, (32, 32, 32),
    KrylovKitMethod(shift=0.01, matrix_free=true);
    cutoff=5
)

Comparison:

How it works:

• Default mode: Builds dense matrices LHS, RHS, computes LU factorization of (A - σB)
• Matrix-free mode: Uses FFT-based operators and iterative CG solver for inner system

For most 3D calculations with cutoff ≤ 5, the default dense mode is faster. Use matrix_free=true when memory becomes the limiting factor.

LOBPCGMethod

LOBPCG (Locally Optimal Block Preconditioned Conjugate Gradient) is an iterative eigenvalue solver. A. V. Knyazev, SIAM J. Sci. Comput. 23, 517 (2001). DOI:10.1137/S1064827500366124

When to Use LOBPCG

LOBPCG is effective for:

• Large-scale band structure calculations (cutoff ≥ 12)
• Symmetric generalized eigenvalue problems A x = λ B x
• Computing multiple eigenvalues simultaneously (block method)

Performance comparison (2D PSV wave, Steel/Epoxy):

*Speedup < 1 means Dense is faster

Recommendation:

• cutoff < 10: Use Dense (faster due to BLAS optimization)
• cutoff ≥ 12: Use LOBPCG (speedup increases with problem size)

Basic Usage

# 2D photonic
solver = Solver(TEWave(), geo, (64, 64), LOBPCGMethod(); cutoff=12)

# 2D phononic (LOBPCG effective for large problems)
solver = Solver(PSVWave(), geo, (64, 64), LOBPCGMethod(); cutoff=15)

Shift-and-Invert for 3D

For 3D calculations, LOBPCG requires shift > 0 to skip spurious longitudinal modes:

solver = Solver(FullVectorEM(), geo, (16, 16, 16), LOBPCGMethod(shift=0.01); cutoff=7)

solver = Solver(FullElastic(), geo, (16, 16, 16), LOBPCGMethod(shift=0.01); cutoff=5)

Why is shift needed for 3D?

The 3D H-field formulation does not enforce the transversality constraint ∇·H = 0, leading to spurious longitudinal modes at λ ≈ 0. Without shift:

• LOBPCG tries to find smallest eigenvalues, converging to these spurious modes
• The RHS matrix B becomes singular in this subspace → PosDefException

With shift=σ, the problem transforms to (A - σB)⁻¹ B x = μ x, which only returns eigenvalues λ > σ, effectively filtering out the spurious modes.

LOBPCGMethod Parameters

method = LOBPCGMethod(
    tol = 1e-3,              # Convergence tolerance (default)
    maxiter = 100,           # Maximum iterations (default)
    shift = 0.0,             # Spectral shift (required for 3D)
    warm_start = true,       # Use previous eigenvectors as initial guess
    first_dense = true,      # Solve first k-point with Dense
)

Warm Start for Band Structure Calculations

When computing band structures with compute_bands, LOBPCG uses warm start to speed up calculations. The eigenvectors from the previous k-point serve as the initial guess for the next k-point.

How it works:

1. First k-point is solved with DenseMethod (when first_dense=true)
2. Subsequent k-points use previous eigenvectors as initial guess
3. Since eigenvectors vary smoothly along k-path, LOBPCG converges quickly

Usage:

# Automatic warm start (default)
solver = Solver(PSVWave(), geo, (64, 64), LOBPCGMethod(); cutoff=15)
bands = compute_bands(solver, kpath; bands=1:10)

Limitations

• For small problems (cutoff < 10), Dense is faster due to BLAS optimization
• Phononic problems with band crossings at Γ point may require additional care
• Accuracy is typically 10-100 rad/s compared to Dense (0.1-1% relative error)

LOBPCG vs KrylovKit

For details on matrix-free methods and their resolution constraints, see Matrix-Free Methods.

API Reference

• Solver API - DenseMethod, KrylovKitMethod, LOBPCGMethod