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.

Sign up for GitHub

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

Jump to bottom

Update integration weights to account for cylindrical coordinates #2470

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Update integration weights to account for cylindrical coordinates #2470

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
145 changes: 145 additions & 0 deletions python/tests/test_cyl_weights.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
"""Tests that the integration weights are computed properly in cylindrical coordinates.

"""

import unittest
from typing import Tuple
import meep as mp
import numpy as np


resolution = 80 # pixels/μm
n = 2.4 # refractive index of dielectric layer
wvl = 1.0 # wavelength (in vacuum)
fcen = 1 / wvl # center frequency of source/monitor


def led_flux(dmat: float, h: float, rpos: float, m: int) -> Tuple[float, float]:
"""Computes the radiated and total flux of a point source embedded
within a dielectric layer above a lossless ground plane.

Args:
dmat: thickness of dielectric layer.
h: height of dipole above ground plane as a fraction of dmat.
rpos: position of source in radial direction.
m: angular φ dependence of the fields exp(imφ).

Returns:
The radiated and total flux as a 2-Tuple.
"""
L = 20 # length of non-PML region in radial direction
dair = 1.0 # thickness of air padding
dpml = 1.0 # PML thickness
sr = L + dpml
sz = dmat + dair + dpml
cell_size = mp.Vector3(sr, 0, sz)

boundary_layers = [
mp.PML(dpml, direction=mp.R),
mp.PML(dpml, direction=mp.Z, side=mp.High),
]

src_pt = mp.Vector3(rpos, 0, -0.5 * sz + h * dmat)
sources = [
mp.Source(
src=mp.GaussianSource(fcen, fwidth=0.1 * fcen),
component=mp.Er,
center=src_pt,
),
]

geometry = [
mp.Block(
material=mp.Medium(index=n),
center=mp.Vector3(0, 0, -0.5 * sz + 0.5 * dmat),
size=mp.Vector3(mp.inf, mp.inf, dmat),
)
]

sim = mp.Simulation(
resolution=resolution,
cell_size=cell_size,
dimensions=mp.CYLINDRICAL,
m=m,
boundary_layers=boundary_layers,
sources=sources,
geometry=geometry,
)

flux_air_mon = sim.add_flux(
fcen,
0,
1,
mp.FluxRegion(
center=mp.Vector3(0.5 * L, 0, 0.5 * sz - dpml),
size=mp.Vector3(L, 0, 0),
),
mp.FluxRegion(
center=mp.Vector3(L, 0, 0.5 * sz - dpml - 0.5 * dair),
size=mp.Vector3(0, 0, dair),
),
)

sim.run(
mp.dft_ldos(fcen, 0, 1),
until_after_sources=mp.stop_when_fields_decayed(
50.0,
mp.Er,
src_pt,
1e-8,
),
)

flux_air = mp.get_fluxes(flux_air_mon)[0]

if rpos == 0:
dV = np.pi / (resolution**3)
else:
dV = 2 * np.pi * rpos / (resolution**2)

# total flux from point source via LDOS
flux_src = -np.real(sim.ldos_Fdata[0] * np.conj(sim.ldos_Jdata[0])) * dV

print(f"flux-cyl:, {rpos:.2f}, {m:3d}, {flux_src:.6f}, {flux_air:.6f}")

return flux_air, flux_src

class TestCylWeights(unittest.TestCase):
def test_weights(self):
layer_thickness = 0.7 * wvl / n
dipole_height = 0.5

rpos = [0.5, 1.5]
flux_tol = 1e-5 # threshold flux to determine when to truncate expansion
normalized_flux = []
for rp in rpos:
# analytic upper bound on m based on coupling to free-space modes
# in light cone of source medium
cutoff_M = int(rp * 2 * np.pi * fcen * n)
ms = range(cutoff_M + 1)
flux_src_tot = 0
flux_air_tot = 0
flux_air_max = 0
for m in ms:
flux_air, flux_src = led_flux(
layer_thickness,
dipole_height,
rp,
m,
)
flux_air_tot += flux_air if m == 0 else 2 * flux_air
flux_src_tot += flux_src if m == 0 else 2 * flux_src
if flux_air > flux_air_max:
flux_air_max = flux_air
if m > 0 and (flux_air / flux_air_max) < flux_tol:
break

normalized_flux.append(flux_air_tot / rp**2)

print("--------------------------")
print(normalized_flux)
# TODO
# check for equivalence

if __name__ == "__main__":
unittest.main()
128 changes: 102 additions & 26 deletions src/loop_in_chunks.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -257,39 +257,100 @@ static inline int iabs(int i) { return (i < 0 ? -i : i); }
void compute_boundary_weights(grid_volume gv, const volume &where, ivec &is, ivec &ie,
bool snap_empty_dimensions, vec &s0, vec &e0, vec &s1, vec &e1) {
LOOP_OVER_DIRECTIONS(gv.dim, d) {
double w0, w1;
w0 = 1. - where.in_direction_min(d) * gv.a + 0.5 * is.in_direction(d);
w1 = 1. + where.in_direction_max(d) * gv.a - 0.5 * ie.in_direction(d);
// w0 = (i+1) - a (from above)
double w0 = 1. - where.in_direction_min(d) * gv.a + 0.5 * is.in_direction(d);
// w1 = b - (i+3) (from above)
double w1 = 1. + where.in_direction_max(d) * gv.a - 0.5 * ie.in_direction(d);

// only needed for d==R cases
double idx_i = is.in_direction(d);
double idx_in = ie.in_direction(d);

// Case 1 (a and b separated by at least 2 grid points)
if (ie.in_direction(d) >= is.in_direction(d) + 3 * 2) {
s0.set_direction(d, w0 * w0 / 2);
s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2);
e0.set_direction(d, w1 * w1 / 2);
e1.set_direction(d, 1 - (1 - w1) * (1 - w1) / 2);
// cylindrical coordinates in R direction require extra care
if (d == R) {
// include dV in the weight to prevent division by zero
if (is.in_direction(R) == 0)
s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2);
else
s0.set_direction(d, (-w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2) / idx_i );

s1.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0 + (idx_i+1)/2 + 1.0 / 6.0) / (idx_i+2));
e0.set_direction(d, (w1 * w1 * w1 / 3 + (idx_in-2) * w1 * w1 / 2) / (idx_in));
e1.set_direction(d, (-w1 * w1 * w1 / 3 - (idx_in-4) * w1 * w1 / 2 + (idx_in-2) * w1 + (idx_in-2)/2 - 1/6) / (idx_in-2));
} else {
s0.set_direction(d, w0 * w0 / 2);
s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2);
e0.set_direction(d, w1 * w1 / 2);
e1.set_direction(d, 1 - (1 - w1) * (1 - w1) / 2);
}

}
// Case 2 (one grid point between a and b)
else if (ie.in_direction(d) == is.in_direction(d) + 2 * 2) {
s0.set_direction(d, w0 * w0 / 2);
s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2 - (1 - w1) * (1 - w1) / 2);
e0.set_direction(d, w1 * w1 / 2);
e1.set_direction(d, s1.in_direction(d));
// cylindrical coordinates in R direction require extra care
if (d == R) {
// include dV in the weight to prevent division by zero
if (is.in_direction(R) == 0)
s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2);
else
s0.set_direction(d, -(w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2) / idx_i );

s1.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0
- w1 * w1 * w1 / 3 - (idx_i) * w1 * w1 / 2 + (idx_i+2) * w1) / (idx_i+2));
e0.set_direction(d, (w1 * w1 * w1 / 3 + (idx_i+2) * w1 * w1 / 2) / (idx_i+2));
e1.set_direction(d, s1.in_direction(d));
} else {
s0.set_direction(d, w0 * w0 / 2);
s1.set_direction(d, 1 - (1 - w0) * (1 - w0) / 2 - (1 - w1) * (1 - w1) / 2);
e0.set_direction(d, w1 * w1 / 2);
e1.set_direction(d, s1.in_direction(d));
}
}
// Case 4 (as (3), but a = b: interpolation, not integration)
// to make the logic easier, we check for 4 before 3...
else if (where.in_direction_min(d) == where.in_direction_max(d)) {
if (snap_empty_dimensions) {
if (w0 > w1)
ie.set_direction(d, is.in_direction(d));
// cylindrical coordinates should be the same as cartesian, with the
// exception at r=0 (due to our convention when defining dv1).
if (snap_empty_dimensions) {
if (w0 > w1)
ie.set_direction(d, is.in_direction(d));
else
is.set_direction(d, ie.in_direction(d));
w0 = w1 = 1.0;
}
// include dV in the weight to prevent division by zero
if ((d == R) && (is.in_direction(R) == 0))
s0.set_direction(d, w0);
else
is.set_direction(d, ie.in_direction(d));
w0 = w1 = 1.0;
s0.set_direction(d, w0);
s1.set_direction(d, w1);
e0.set_direction(d, w1);
e1.set_direction(d, w0);
}
s0.set_direction(d, w0);
s1.set_direction(d, w1);
e0.set_direction(d, w1);
e1.set_direction(d, w0);
}
// Case 3 (no grid points between a and b)
else if (ie.in_direction(d) == is.in_direction(d) + 1 * 2) {
s0.set_direction(d, w0 * w0 / 2 - (1 - w1) * (1 - w1) / 2);
e0.set_direction(d, w1 * w1 / 2 - (1 - w0) * (1 - w0) / 2);
s1.set_direction(d, e0.in_direction(d));
e1.set_direction(d, s0.in_direction(d));
// cylindrical coordinates in R direction require extra care
if (d == R) {
// include dV in the weight to prevent division by zero
if (is.in_direction(R) == 0)
s0.set_direction(d, -w0 * w0 * w0 / 3 + w0 * w0 / 2 - w1 * w1 * w1 / 3 + w1 *w1 /2 - 1/6);
else
s0.set_direction(d, (-w0 * w0 * w0 / 3 + (idx_i+2) * w0 * w0 / 2
- w1 * w1 * w1 / 3 - (idx_i-2) * w1 * w1 / 2 + idx_i * w1 - idx_i/2 - 1/6) / idx_i );

e0.set_direction(d, (w0 * w0 * w0 / 3 - (idx_i+4) * w0 * w0 / 2 + (idx_i+2) * w0
+ w1 * w1 * w1 /3 + (idx_i) * w1 * w1 / 2 - (idx_i+2) / 2 + 1/6) / (idx_i+2));
s1.set_direction(d, e0.in_direction(d));
e1.set_direction(d, s0.in_direction(d));

} else {
s0.set_direction(d, w0 * w0 / 2 - (1 - w1) * (1 - w1) / 2);
e0.set_direction(d, w1 * w1 / 2 - (1 - w0) * (1 - w0) / 2);
s1.set_direction(d, e0.in_direction(d));
e1.set_direction(d, s0.in_direction(d));
}
}
else
meep::abort("bug: impossible(?) looping boundaries");
Expand Down Expand Up @@ -352,9 +413,19 @@ void fields::loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, con
vec yee_c(gv.yee_shift(Centered) - gv.yee_shift(cgrid));
ivec iyee_c(gv.iyee_shift(Centered) - gv.iyee_shift(cgrid));
volume wherec(where + yee_c);

/* Integrating zero crossings in cylindrical coordinates is really hairy...
let's just punt for now...*/
if ((gv.dim == Dcyl) && (wherec.get_min_corner().in_direction(R) < 0))
meep::abort("Monitors and sources must not extend below r=0.");

ivec is(vec2diel_floor(wherec.get_min_corner(), gv.a, zero_ivec(gv.dim)) - iyee_c);
ivec ie(vec2diel_ceil(wherec.get_max_corner(), gv.a, zero_ivec(gv.dim)) - iyee_c);

//In cylindrical coordinates, we don't need the extra pad point beyond r=0
if ((gv.dim == Dcyl) && (is.in_direction(R) < 0))
is.set_direction(R,is.in_direction(R) + 2); // add a yee voxel width

vec s0(gv.dim), e0(gv.dim), s1(gv.dim), e1(gv.dim);
compute_boundary_weights(gv, where, is, ie, snap_empty_dimensions, s0, e0, s1, e1);

Expand Down Expand Up @@ -509,7 +580,12 @@ void fields::loop_in_chunks(field_chunkloop chunkloop, void *chunkloop_data, con
}
if (gv.dim == Dcyl) {
dV1 = dV0 * 2 * pi * gv.inva;
dV0 *= 2 * pi * fabs((S.transform(chunks[i]->gv[isc], sn) + shift).in_direction(R));
// we include the dV term in the computation of the s0 weight when
// r=0 so that we don't have to worry about dividing by zero.
if (is.in_direction(R) == 0)
dV0 = 1;
else
dV0 *= 2 * pi * fabs((S.transform(chunks[i]->gv[isc], sn) + shift).in_direction(R));
}

chunkloop(chunks[i], i, cS, isc, iec, s0c, s1c, e0c, e1c, dV0, dV1, shifti, ph, S, sn,
Expand Down