Update integration weights to account for cylindrical coordinates

[smartalecH]

closes #2467.

This PR attempts to compute the proper integration weights for $r$ in cylindrical coordinates. There are a few things to be aware of:

1. @oskooi I don't think this fixes the issue with dipoles near $r=0$. Primarily because the interpolation formula for a dipole in cylindrical coordinates appears to be the same as Cartesian (at least according to the in-repo notes). The integration is obviously different, so that affects other cases involving lines/planes, but not points in $r$. So while this PR should help with some of the edge elements of monitors, it's not the fix we're looking for.
2. I added an abort if the user attempts to place a source or monitor that spans $r<0$. Physically, I don't even know what this means (a double source/monitor?) and it's a pain to deal with in the code (as described in the notes). So we punt.
3. Dealing with the weights right at $r=0$ is tricky. We don't want to divide by zero, so I include the dV term in the weight itself here. This required some additional modification to the loop_in_chunks() function, but doesn't seem to affect the IVEC_LOOP_WEIGHT macro (contrary to what the notes suggest) so I left that alone
4. Getting the indexing and resolution just right is a bookkeeping nightmare... so we should double check this.
5. The current test isn't a suitable test, since this doesn't fix that issue (plus it takes far too long to run with all of the m simulations sequentially, even for just two points relatively close to $r=0$.

I think this PR is more in the noise floor when it comes to comparing accuracy of a cylindrical simulation and a true 3D simulation with dipole sources.

@stevengj, since the issue isn't with source grid weights themselves, can you think of any other factors that would impact a dipole source in cylindrical coordinates?

[stevengj]

Primarily because the interpolation formula for a dipole in cylindrical coordinates appears to be the same as Cartesian (at least according to the in-repo notes).

I think that might be wrong — for a dipole source, what you want to do is not quite interpolation, I think. Really you should think of it as a certain integral you need to match (e.g. a finite-width delta approximation you are trying to sample, such that the integral remains the same).

[stevengj]

Would be good to check that it computes ∫ₐᵇf(r)rdr exactly for any affine f(r). It should suffice to check that it integrates 1 and r exactly for various intervals [a,b]

[smartalecH]

Hmm I'm still not convinced that there's a problem with the current implementation of a dipole in cylindrical coordinates.

First, consider the problem in cartesian coordinates:

| -------o-----|
a        x     b 

where $a$ and $b$ are grid points, and $x$ is the location of our dipole in continuous space. The weights we use at grid points $a$ and $b$ we'll call $w_a$ and $w_b$ respectively.

Meep currently defines the weights using simple linear interpolation, such that:
$$\Delta x = b-a$$
$$w_a + w_b = 1$$
$$w_a=\frac{1}{\Delta x}\left(b-x\right)$$
$$w_b=\frac{1}{\Delta x}\left(x-a\right)$$

(The $\Delta x$ term isn't applied until the IVEC_LOOP_WEIGHT multiplies by dV).

Now, consider the problem in cylindrical coordinates:

| -------o-----|
a        r     b 

where the only difference here is $r$. But now, meep does a slightly different operation to determine the weights because of the dV0 and dV1 terms used to compute the $r$ in the integral:
$$w_a=\frac{1}{\Delta r}a\left(b-r\right)$$
$$w_a=\frac{1}{\Delta r}b\left(r-a\right)$$

Let's see what happens when we "integrate" over the weights by summing them together:
$$w_a + w_b$$
$$\frac{1}{\Delta r}a\left(b-r\right)+\frac{1}{\Delta r}b\left(r-a\right)$$
$$=\frac{b-a}{\Delta r} r$$
$$=r$$

...which is exactly what we want! Whereas in cartesian coordinates, we expect the integral of our delta function to sum to 1, we expect the integral in cylindrical coordinates for our delta function to sum to $r$, because we don't incorporate that factor by default in meep (and it's part of the definition of a delta function in cylindrical coordinates).

So again, I'm not sure the issue is with the weights of the dipole source...

[smartalecH]

So again, I'm not sure the issue is with the weights of the dipole source...

... and yet, experimentally, the code doesn't seem to be behaving as intended. I ran a simple test where I have a 1D cylindrical simulation, and I pull the weights of a dipole that I "slide" along $R$. Ideally, the sum of the weights should be $r$ itself... but there's a lot of noise...

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

sim = mp.Simulation(
    resolution=10,
    cell_size=mp.Vector3(100,0,0),
    dimensions=mp.CYLINDRICAL,
    m=1,
)
sim.init_sim()
r_bank = np.linspace(0.5,100,50)
w_bank = []
for r in r_bank:
    (x,y,z,w) = sim.get_array_metadata(mp.Volume(center=mp.Vector3(r,0,0),size=mp.Vector3(0,0,0)),)
    w_bank.append(w)


plt.figure()
plt.plot(r_bank,r_bank,color='k',label='ideal weights')
plt.plot(r_bank,w_bank, color='r',label='measured weights')
plt.xlabel("$r$ postion")
plt.ylabel("sum of weights")
plt.legend()
plt.savefig("weights.png")
plt.show()

(Now I'm using get_array_metadata() to pull the weights, which automatically "collapses" the weights of a dipole by summing it together. There could be a bug here rather than the actual weight computation routine?)

[oskooi]

As a useful reference, here is a schematic of ΔV in cylindrical coordinates from slide 2 of these lecture notes:

[oskooi]

As I previously commented in #2108, I think there could be a bug in fields::add_volume_source related to the normalization of the $\delta$-function source in cylindrical coordinates:

meep/src/sources.cpp

Lines 481 to 484 in aedc794

	LOOP_OVER_DIRECTIONS(gv.dim, d) {
	if (where.in_direction(d) == 0.0 && !nosize_direction(d)) // delta-fun
	data.amp *= gv.a; // correct units for J delta-function amplitude
	}

The normalization is $\Delta r * \Delta z$ which is the same as 2D: $\Delta x * \Delta y$. A $\delta(r)$ source in cylindrical coordinates would seem to require an additional $r$ term in the normalization which is missing.

As a test, I modified the point source in python/examples/point_dipole_cyl.py to include amplitude=1/rpos for rpos > 0. However, this did not affect the plot of flux_air_tot / rpos**2 flux vs. rpos. The normalized radiated flux is still not a constant independent of rp.

[smartalecH]

The normalization is $\Delta r * \Delta z$ which is the same as 2D: $\Delta x * \Delta y$. A $\delta(r)$ source in cylindrical coordinates would seem to require an additional $r$ term in the normalization which is missing.

Wait, I thought this was already known (per our previous discussion). In fact, I thought this was why the test divides by $r^2$ in the first place (to account for the missing $r$ in the sources)?

This isn't a "bug", just a matter of convention. My above comment discusses this.

(By the way, we do have an extra $r$ term when computing dr. See here:

meep/src/loop_in_chunks.cpp

Lines 308 to 316 in aedc794

	We also pass CHUNKLOOP() dV0 and dV1, such that the integration
	"grid_volume" dV is dV0 + dV1 * iloopR, where iloopR is the loop
	variable (starting from 0 at the starting integer coord and
	incrementing by 1) corresponding to the direction R. Note that, in
	the LOOP_OVER_IVECS macro, iloopR corresponds to the loop variable
	loop_i2 in Dcyl (cylindrical coordinates). In other coordinates,
	dV1 is 0. Note also that by "grid_volume" dV we mean the integration
	unit corresponding to the dimensionality of WHERE (e.g. an area if
	WHERE is 2d, etc.)

and it is applied to the source... just within the weights themselves:

meep/src/sources.cpp

Line 273 in aedc794

amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc);

meep/src/meep/vec.hpp

Lines 371 to 383 in aedc794

	// integration weight for using LOOP_OVER_IVECS with field::integrate
	#define IVEC_LOOP_WEIGHT1x(s0, s1, e0, e1, i, n, dir) \
	((i > 1 && i < n - 2) \
	? 1.0 \
	: (i == 0 ? (s0).in_direction(meep::direction(dir)) \
	: (i == 1 ? (s1).in_direction(meep::direction(dir)) \
	: i == n - 1 ? (e0).in_direction(meep::direction(dir)) \
	: (i == n - 2 ? (e1).in_direction(meep::direction(dir)) : 1.0))))
	#define IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, k) \
	IVEC_LOOP_WEIGHT1x(s0, s1, e0, e1, loop_i##k, loop_n##k, loop_d##k)
	#define IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV) \
	(IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 3) * \
	(IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 2) * ((dV)*IVEC_LOOP_WEIGHT1(s0, s1, e0, e1, 1))))

)

EDIT: wait, @oskooi, I think you're on to something.

Notice that the weight update passes in 1 for the integration weight, (presumably because of the factor you are referencing). Instead, it should look like this:

IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2);

and the "global" scaling you are describing shouldn't exist at all. It should be done at the chunk level.

[smartalecH]

@oskooi you can try this patch to see if this resolves the behavior you're seeing: https://github.com/smartalecH/meep/tree/cyl_dipoles

[oskooi]

you can try this patch to see if this resolves the behavior you're seeing:

There are still rpos-dependent values for flux_air_tot / rpos**2 in the modified test based on point_dipole_cyl.py using your cyl_dipoles branch. The extraction-efficiency values are slightly different compared to the master branch which is unexpected. Finally, three of the unit tests in python/tests are failing using the cyl_dipoles branch: harmonics.py, known_results.py, and near2far.py.

[smartalecH]

There are still rpos-dependent values for flux_air_tot / rpos**2 in the modified test based on point_dipole_cyl.py using your cyl_dipoles

Bummer.

The extraction-efficiency values are slightly different compared to the master branch which is unexpected

Why is this unexpected? I would expect them to be different if the dipole is changing.

Finally, three of the unit tests in python/tests are failing using the cyl_dipoles branch

Yes this is expected. I didn't write any logic to check for dipoles in the chunk code (which is much harder). So tests involving source volumes (not dipoles) or going to have an additional scale factor with this patch. I wanted to test if this fixed the issue before addressing all corner cases.

[oskooi]

Why is this unexpected? I would expect them to be different if the dipole is changing.

Since the extraction efficiency is a dimensionless quantity, scaling the amplitude of a point source by a constant factor should not affect its result. The scale factor would only affect dimensionful quantities such as the radiated flux (numerator) and source power (denominator) which indeed it does.

[smartalecH]

scaling the amplitude of a point source by a constant factor should not affect its result.

But the point of this is that the point source is changing. We aren't just changing the relative amplitude, we're changing the profile. LEE is highly dependent on the type of source (position, size, polarization, etc). So the LEE should in fact change (otherwise, we wouldn't bother fixing the bug).

(Recall that nothing in meep is a true point source. It's a projection onto the closest grid points)

[stevengj]

To be clear, there are two separate issues here:

• Do we define a cylindrical point source as δ(r) or δ(r)/r (or δ(r)/2πr)? This is mostly a matter of convention. As long as we document what we are doing, the user can easily multiply or divide by r as needed (in the source amplitude or in post-processing). Currently, we define it as δ(r); we may want to change this in the future, but I suggest that we stick with the current convention for now.
• How do we discretize δ(r) onto the grid? Currently, we are using just linear interpolation, but I think that introduces a first-order error in cylindrical coordinates. That is what we would like to fix, and is independent of whether there is a 1/r in the amplitude convention.

(A symptom of the fact that the latter is just about discretization is that you should find that it goes away as you crank up the resolution. i.e. try your test from #2108 at a higher resolution, and you should find that it approaches a constant at a smaller radius.)

[smartalecH]

How do we discretize δ(r) onto the grid? Currently, we are using just linear interpolation, but I think that introduces a first-order error in cylindrical coordinates. That is what we would like to fix, and is independent of whether there is a 1/r in the amplitude convention.

But from my comment above, I think we already interpolating in cylindrical coordinates correctly. So the issue doesn't appear to be the method (or math) but something else (perhaps the implementation?)

[oskooi]

A symptom of the fact that the latter is just about discretization is that you should find that it goes away as you crank up the resolution. i.e. try your test from #2108 at a higher resolution, and you should find that it approaches a constant at a smaller radius.

Looks like the discretization error for a point source $\delta(r-r_0)$ does not go away with increasing resolution.

As a demonstration, here are results for dipole emission in vacuum (same test in #2108) plotting radiated flux / $r_0^2$ vs. $r_0$ (dipole position) for three different resolutions. This is using the master branch (not #2483). The results in the plot seem to be converged for the lowest resolution of 50 since there is no significant change at the higher resolutions.

Simulation layout with dipole at $r_0 = 0$.