subpixel averaging using adaptive cubature

[oskooi]

The routines fallback_chi1inv_row in src/meepgeom.cpp and normal_vector in src/anisotropic_averaging.cpp provide fallback methods using adaptive cubature for computing the effective permittivity tensor and interface normal for subpixel averaging in cases not involving GeometricObjects such as a material-function or epsilon-input-file. Because these numerical integration routines tend to be rather slow, these functions have been disabled and no smoothing is performed for material-function or epsilon-input-file.

It would be useful to enable the fallback adaptive cubature routines as an option mainly for e.g. testing purposes. Currently, they seem to be inaccessible from the interface. As an example, the following script which is adapted from the tutorial example for the ring resonator creates the ring geometry using a material function rather than two Cylinders.

import meep as mp

n = 3.4                 # index of waveguide                                                                                                              
w = 1                   # width of waveguide                                                                                                              
r = 1                   # inner radius of ring                                                                                                            
pad = 4                 # padding between waveguide and edge of PML                                                                                       
dpml = 2                # thickness of PML                                                                                                                
sxy = 2*(r+w+pad+dpml)  # cell size                                                                                                                       

def ring_resonator(p):
    rr = (p.x**2+p.y**2)**0.5
    if (rr > r) and (rr < r+w):
        return mp.Medium(index=n)
    return mp.air

geometry = [mp.Block(center=mp.Vector3(),
                     size=mp.Vector3(sxy,sxy),
                     material=ring_resonator)]

fcen = 0.175             # pulse center frequency                                                                                                         
df = 0.05                # pulse frequency width                                                                                                          

src = [mp.Source(mp.GaussianSource(fcen, fwidth=df),
                 component=mp.Ez,
                 center=mp.Vector3(r+0.1))]

sim = mp.Simulation(cell_size=mp.Vector3(sxy,sxy),
                    geometry=geometry,
                    sources=src,
                    resolution=20,
                    eps_averaging=True,
                    symmetries=[mp.Mirror(mp.Y)],
                    boundary_layers=[mp.PML(dpml)])

sim.run(mp.after_sources(mp.Harminv(mp.Ez, mp.Vector3(r+0.1), fcen, df)),
        until_after_sources=300)

The output from Harminv is the same regardless of whether eps_averaging is True or False.

Using MPI version 3.1, 1 processes
-----------
Initializing structure...
Halving computational cell along direction y
Working in 2D dimensions.
Computational cell is 16 x 16 x 0 with resolution 20
     block, center = (0,0,0)
          size (16,16,0)
          axes (1,0,0), (0,1,0), (0,0,1)
time for set_epsilon = 0.992031 s
-----------
harminv0:, frequency, imag. freq., Q, |amp|, amplitude, error
harminv0:, 0.17584282569060494, -6.472394727883501e-05, 1358.4062243072299, 0.03933965687736313, 0.03682037000791641-0.013851676999872385i, 3.0340847418159567e-09+0.0i
run 0 finished at t = 500.0 (20000 timesteps)

Elapsed run time = 4.6469 s

[stevengj]

To enable this, basically I think you need two things:

• At this line, check if it is a material function. If so, you can optionally call fallback_chi1inv_row instead of goto trivial

• There will have to be a way to optionally turn cubature of material functions on (off by default), without affecting anything else. One possibility would be an extra average=true/false flag in the material-function type itself, that could be checked on the line mentioned above.

In principle, we could also do subpixel averaging from a high-resolution epsilon-input-file, but I suspect that will be insanely slow because the corresponding bilinear interpolation ε(x) is only piecewise linear, which will destroy the convergence of the adaptive cubature scheme we are using.

In the long run, there is the potential for a more efficient scheme:

1. Approximate the surface locally by a plane. This can be done analytically for a geometric object like a block or a sphere. For a material function (or ε input file), it can be done by 3 bisection searches to find the intersections of the surface with three edges of the voxel.

2. Compute the intersection of the half-space (one side of the plane) with the voxel using analytical formulas. The formulas are rather messy, unfortunately. @wsshin already implemented this in Julia for subpixel averaging in his FDFD code: https://github.com/stevengj/GeometryPrimitives.jl/blob/master/src/vxlcut.jl

The reason this works is that, for smooth surfaces, in the limit of small pixels the intersection with the pixels is always approximately planar (plus higher-order corrections that don't affect 2nd-order accuracy).

[oskooi]

Beyond just enabling a do_averaging property in #771, it would be useful to specify the tol and maxeval parameters for material_function::eff_chi1inv_row. Currently, these parameters have fixed default values of 1e-4 and 100000.

[stevengj]

Instead of a do_averaging property, maybe we should have an averaging_maxeval property that defaults to 1, which you can increase if you want to do averaging.

[ChristopherHogan]

tol and maxeval can already be set in the Simulation constructor via the keyword arguments subpixel_tol and subpixel_maxeval. These get passed from the structure constructor to eff_chi1inv_row.

[stevengj]

An optimization: in normal_vector, if the first 2ᵈ (1 << d in C, where d = number_of_directions(v.dim)) points in the for (i ...) loop have the same value of chi1p1(ft, pt), then it should give up and return zero_vec(v.dim), which will cause the fallback to go to the trivial case.

This will eliminate the averaging from regions where ε is uniform. (In principle, normal_vector should return 0 in this case already, but I think it may not due to roundoff errors.)

[stevengj]

It looks like the

    for (int i = 0; i < 3; ++i)
      chi1inv_row[i] = n[rownum] * n[i] * (minveps - 1 / meps);
    chi1inv_row[rownum] += 1 / meps;

code in meepgeom.cpp is wrong. In particular, the correct formula for the effective ε⁻¹ in the scalar case is minveps*P + (1/meps)*(I-P) where P is the matrix whose entries are Pᵢⱼ = n[i]*n[j]. This corresponds to:

    for (int i = 0; i < 3; ++i)
      chi1inv_row[i] = n[rownum] * n[i] * (minveps - 1 / meps);
    chi1inv_row[rownum] = 1 / meps + n[rownum] * n[rownum] * (minveps - 1 / meps);

Update: I take that back. The previous code was correct. The += on the diagonal term gives the missing term.

[stevengj]

The ε integrand (ceps_func etcetera) looks a little odd to me for cylindrical coordinates, so we should eventually double check that. (It might be okay. It is including a radius factor via the scale factor r, but there is a similar factor in vol so it cancels.)

[stevengj]

I just noticed a bug: these lines should use chi1p1_inv on the right-hand-side, not chi1p1. Similarly for these lines, which are currently never called.

[stevengj]

Also, we can change the meep::abs(gradient) == 0 to meep::abs(gradient) < 1e-8