Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

subpixel averaging for boundaries of variable-material objects #1399

Open
wants to merge 4 commits into
base: master
Choose a base branch
from

Conversation

stevengj
Copy link
Collaborator

@stevengj stevengj commented Oct 10, 2020

This PR implements subpixel averaging for the boundaries of objects with variable materials (e.g. user-defined material functions). I realized that we can still handle this case accurately and efficiently, as long as the material function is continuous inside the object.

For example, here is a simulation of transmitted power through a cylinder of linearly-varying ε:
image

def linear_eps(p):
    return 1 + 11 * (p.x + 0.4) / 0.8

def flux2(res):
    sim = mp.Simulation(cell_size=(7,1), resolution=res, boundary_layers=[mp.PML(direction=mp.X, thickness=0.5)],
                    k_point=mp.Vector3(0,0), symmetries=[mp.Mirror(direction=mp.Y)],
                    sources=[mp.Source(mp.GaussianSource(frequency=0.7, fwidth=0.2, is_integrated=True),
                                       component=mp.Ez, center=(-3,0), size=(0,1))],
                    geometry=[mp.Cylinder(radius=0.4, height=mp.inf, epsilon_func=linear_eps)])
    tran = sim.add_flux(0.7, 0.2, 1, mp.FluxRegion(center=(3,0), size=(0,1)))
    sim.run(until=200)
    return tran.flux()[0]

If I compute and plot the convergence with respect to resolution using this code:

import matplotlib.pyplot as plt
import numpy as np
r = [10,20,40,80,160,320,640]
f2 = np.asarray([flux2(r) for r in r])
r = np.asarray(r)
plt.loglog(r[0:-1], np.abs(f2[0:-1] - f2[-1]) / f2[-1], "ro-")
plt.loglog(r[0:-1], 30 / r[0:-1]**2, "k--")
plt.loglog(r[0:-1], 1 / r[0:-1], "b:")
plt.legend(["flux2 error", "$1/r^2$ reference", "$1/r$ reference"])
plt.xlabel("resolution $r$")
plt.ylabel("fractional error")

then before this PR I get first-order convergence due to the lack of subpixel averaging at the boundaries of the cylinder:
image
but after this PR I get second-order convergence:
image

}

/* check for trivial case of only one object/material */
if (material_type_equal(mat, mat_behind)) {
if (is_variable(mat)) eval_material_pt(mat, vec_to_vector3(v.center()));
goto trivial;
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In #1539, this will branch to the new averaging procedure if mat is a variable material.


Requires moderately horrifying logic to figure things out properly,
stolen from MPB. */
static bool get_front_object(const meep::volume &v, geom_box_tree geometry_tree, vector3 &pcenter,
const geometric_object **o_front, vector3 &shiftby_front,
material_type &mat_front, material_type &mat_behind) {
material_type &mat_front, material_type &mat_behind,
vector3 &p_front, vector3 &p_behind) {
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that we added p_front and p_behind here — this is so that we have a point to evaluate the material at if it is a variable materials (grid or user).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

1 participant