Add documentation for bfast_scaled_k parameter of Simulation class constructor

[oskooi]

Adds documentation for the API of the recently added BFAST feature from #2609.

There is a formula in the BFAST paper for an upper bound for the Courant parameter:

In testing, I found that it is often necessary to specify a value much less than this upper bound (i.e., adding an additional factor of $1/n^2$ where $n$ is the refractive index of the lossless medium of the source as in the updated unit test in this PR). As a result, I think it would be better not to specify a default value as suggested in item 3 of #2609 (comment).

Note: we will need to regenerate the Markdown pages for the user manual from the docstrings after this is merged in order for these changes to appear in ReadTheDocs.

cc @Dan2357

[stevengj]

This bound come (I think) from a Von Neumann stability analysis, which should be a tight bound for homogeneous medium, and is usually pretty close to tight even for inhomogeneous media.

It should have a factor of 1/n for a medium with index n, but I don't see why it would be 1/n^2. Can you double check this for a homogeneous medium?

(If there is any doubt here it would be good to dig back through the derivation.)

[oskooi]

Reproduce the results in this paper first for vacuum and then for a homogeneous medium whereby $c$ is replaced by $c/n$:

[oskooi]

It should have a factor of 1/n for a medium with index n, but I don't see why it would be 1/n^2. Can you double check this for a homogeneous medium?

For the case of homogeneous media of index $n$, the necessary Courant factor to ensure stability is $S_n = (1 - \sin(\theta)) / (n\sqrt{D})$. This is consistent with equation (31) of the BFAST paper for which $S = c\Delta t / \Delta x$ and $c$ is replaced by $c/n$.

However, for the case heterogeneous media involving two media of index $n_1$ and $n_2$ separated by a flat interface with the planewave source in $n_1$, my experiments have shown that for $n_1=1.4$ and $n_2=3.5$ a Courant factor of $(1 - \sin(\theta)) / (n_1\sqrt{D})$ is only stable up to $\theta \approx 48.0^\circ$. For $\theta> 48.0^\circ$, it is necessary to use $(1 - \sin(\theta)) / (n_1^2\sqrt{D})$.

[oskooi]

Addendum to my previous comment. The tests involving homogeneous media that I referenced all involved $n = 1.4$. These simulations were stable for large $\theta$ approaching $90^\circ$ using the formula for $S_n$.

However, if I change the media to just vacuum ($n = 1$) the simulation is only stable up to about $\theta = 20^\circ$ using $S = (1 - \sin(\theta))/\sqrt{D}$. This means we are not able to reproduce the results from Fig. 1 of the BFAST paper using equation 31.

[stevengj]

Maybe try a 2d calculation with a finite period and the s polarization, to more closely match what @Dan2357 tried?

[oskooi]

Maybe try a 2d calculation with a finite period and the s polarization,

Launching an oblique planewave in 2d is complicated by the fact that the amplitude function of the source applies to a single frequency component. Nevertheless, I tried setting this up using the amplitude function $e^{i k_{BFAST} \cdot r}$ where $k_{BFAST} = n(\cos(\theta), \sin(\theta), 0)$ for a source medium of index $n$ (= 1.5) and incident angle $\theta$ (= 40.0°). The center wavelength of the pulse is 1 μm. A snapshot of $E_z$ of the pulsed source shows various spectral components but is not symmetric about the line source at the origin as would be expected. The fields eventually blow up even when using the Courant factor from the BFAST paper of $(1 - \sin(\theta)) / (n * \sqrt{2})$.

Note that the unit test for the BFAST feature in tests/test_refl_angular.py uses a 1D cell with a non-zero k_point which results in a 3D simulation. The source is a point and thus the amplitude function is not used.

The main issue here is: how to generate a pulsed oblique planewave for a broad bandwidth in 2d (and 3d). It seems this is not currently possible using the conventional approach of an amplitude function.

"""Launch a planewave at oblique incidence in 2d using the BFAST feature."""

import cmath
import math

import matplotlib.pyplot as plt
import meep as mp
import numpy as np


mp.verbosity(2)

resolution_um = 50
pml_um = 3.0
size_um = 10.0
cell_size = mp.Vector3(size_um + 2 * pml_um, size_um, 0)
pml_layers = [mp.PML(thickness=pml_um, direction=mp.X)]

# Incident angle of planewave. 0 is +x with rotation in                                                           
# counter clockwise (CCW) direction around z axis.                                                                
incident_angle = np.radians(40.0)

wavelength_um = 1.0
frequency = 1 / wavelength_um

n_mat = 1.5  # refractive index of homogeneous material                                                           
default_material = mp.Medium(index=n_mat)

bfast_scaled_k = (
    n_mat * np.cos(incident_angle),
    n_mat * np.sin(incident_angle),
    0
)
print(f"bfast_scaled_k = {bfast_scaled_k}")

if incident_angle == 0:
    eig_parity = mp.EVEN_Y + mp.ODD_Z
    symmetries = [mp.Mirror(mp.Y)]
else:
    eig_parity = mp.ODD_Z
    symmetries = []

def pw_amp(k, x0):
    def _pw_amp(x):
	return cmath.exp(1j * 2 * math.pi * k.dot(x + x0))
    return _pw_amp

src_pt = mp.Vector3()

sources = [
    mp.Source(
	src=mp.GaussianSource(frequency, fwidth=0.2 * frequency),
	center=src_pt,
	size=mp.Vector3(0, size_um, 0),
	component=mp.Ez,
	amp_func=pw_amp(
            mp.Vector3(bfast_scaled_k[0], bfast_scaled_k[1], 0),
            src_pt
	)
    )
]

Courant = (1 - math.sin(incident_angle)) / (math.sqrt(2) * n_mat)

sim = mp.Simulation(
    cell_size=cell_size,
    resolution=resolution_um,
    Courant=Courant,
    boundary_layers=pml_layers,
    sources=sources,
    k_point=mp.Vector3(),
    bfast_scaled_k=bfast_scaled_k,
    default_material=default_material,
    symmetries=symmetries,
)

def ez_snapshot(sim):
    fig, ax = plt.subplots()
    sim.plot2D(
        ax=ax,
        fields=mp.Ez,
        field_parameters={"colorbar":True}
    )
    ax.set_title(f"t = {sim.meep_time():.2f}")
    fig.savefig(
        f"planewave_source_t{sim.meep_time():.2f}.png",
        dpi=150,
        bbox_inches="tight"
    )
    plt.close()

fig, ax = plt.subplots()
sim.plot2D()
fig.savefig(
    f"planewave_source_layout.png",
    dpi=150,
    bbox_inches="tight"
)

sim.run(
    mp.at_every(3.56, ez_snapshot),
    until_after_sources=mp.stop_when_fields_decayed(10.0, mp.Ez, src_pt, 1e-6)
)

[stevengj]

When you are testing stability, it doesn't matter what source you use. Just a point source should be fine.

If you want test response to a planewave in 2d, remember that the fields in the BFAST approach are not the physical fields — they have the exp(ikx) factored out. The same applies to the currents. So, you don't need a pw_amp source, and can just use a constant amplitude.

[stevengj]

(You might try replacing the PML with a scalar absorber, to check whether there is a bug in the PML terms causing an instability.)

[oskooi]

Fields from a point source are not symmetric:

[stevengj]

You don't expect the fields from a point source to be symmetric, because the implicit exp(ikx) breaks the mirror symmetry. You effectively have a phased array of point sources.

[stevengj]

Note that when we use this feature, k_point should be set to zero — the exp(ikx) dependence has been factored out of the fields, and what is left is periodic.