Radiation pattern of a dipole emitter does not converge with resolution in cylindrical coordinates

[oskooi]

As part of a convergence analysis of the radiation pattern of a dipole emitter in cylindrical coordinates (based on the calculation procedure described in #2472), it seems the results do not converge with resolution. Not sure whether this is an indication of a bug somewhere or the way the calculation has been set up.

The test structure involves a point source at $r = 0$ in a finite-size "disc" (a modification of the structure in the LED tutorial which is infinitely extended). There are two different configurations of the closed box of near-field monitors (see schematics below). By Poynting's theorem, the two monitor configurations should produce identical results in the limit of infinite resolution.

While the extraction efficiency over all angles computed using the Poynting flux is in fact identical (= 1.0), there are small but noticeable differences in the radiation pattern computed using a near-to-far field transformation. These discrepancies do not go away with increasing resolution. The script used to generate the results is here.

Here is a side-by-side comparison of the extraction efficiency (really the fractional flux $P_\theta$ in the far field) within four angular cones at a resolution of 200. The final column is the absolute difference in the results (i.e., a measure of the error).

half angle of collection cone	Config. 1 ext. eff.	Config. 2 ext. eff.	\	ee1-ee2\
10°	0.4539	0.4421	0.0118
20°	0.7451	0.7480	0.0029
30°	0.8683	0.8642	0.0041
40°	0.9143	0.9154	0.0011

Same results at a resolution of 400. The difference in the results (final column) is not decreasing as would be expected.

half angle of collection cone	Config. 1 ext. eff.	Config. 2 ext. eff.	\	ee1-ee2\
10°	0.4524	0.4399	0.0125
20°	0.7420	0.7446	0.0026
30°	0.8644	0.8601	0.0043
40°	0.9110	0.9118	0.0008

Configuration 1

resolution = 200

resolution = 400

Configuration 2

resolution = 200

resolution = 400

[oskooi]

There is an even larger disparity in the radiation pattern for Configurations 1 and 2 for a dipole source at $r = 0.5$ μm. These results are for a resolution of 200 pixels/μm. For λ = 1 μm and a disc radius of 0.625 μm, this is equivalent to 200 pixels and 125 pixels, respectively, which should be sufficient for convergence.

These results seem to indicate that there is likely a bug somewhere, possibly in the near-to-far field transformation in cylindrical coordinates?

Configuration 1

Configuration 2

[smartalecH]

Does the total LEE itself converge?

If so, then you can narrow the bug down to the near2far.

[oskooi]

Does the total LEE itself converge?

The total EE is converged as I mentioned in my first comment.

The bug is likely in the near to far transformation function greencyl.

[smartalecH]

From your gist:

for i in range(npts):
        ff = sim.get_farfield(n2f_mon,
                              mp.Vector3(r*math.cos(angles[i]),
                                         0,
                                         r*math.sin(angles[i])))
        E[i,:] = [np.conj(ff[j]) for j in range(3)]
        H[i,:] = [ff[j+3] for j in range(3)]

    Pr = np.real(E[:,1]*H[:,2]-E[:,2]*H[:,1])
    Pz = np.real(E[:,0]*H[:,1]-E[:,1]*H[:,0])
    Prz = np.sqrt(np.square(Pr)+np.square(Pz))

Note that in the far field, the solutions are plane waves, and you should simply be able to calculate the Poynting flux using $|E|^2$.

Is 1000 wavelengths far enough away (should be... but still)?

Also, where is the integral that computes the fractional LEE in a cone? I didn't see that in the gist.

[smartalecH]

So I would expect a difference between configuration 1 and configuration 2 because of that vertical line monitor that changes in $r$. Recall that the integration weights have a first order error with $r$. That error increases the closer you get to $r=0$.

But that error should go away with resolution.

[oskooi]

It looks like the radiation pattern obtained using two different configurations of the near-field monitors converges to the same result only when the length $L$ of the cell in the $r$ direction approaches $\infty$. This means that the DFT monitor in the $r$ direction becomes further and further away from $r = 0$.

As a demonstration, here are the radiation patterns generated from a dipole at $r = 0.2$ ($M + 1 = 7$ simulations) inside a "disc" for $L$ of 5, 10, 20, and 40. Also shown are the schematic of the cell for each configuration and two rows of data, one per DFT monitor configuration, listing the total flux and four numbers (in parentheses) for the fractional flux in the far field within an angular cone of 10°, 20°, 30°, and 40° half angles.

Note 1: The size of the disc is $0.2L$ and thus increases with $L$. This is why the radiation pattern changes for each $L$. The position of the dipole source, however, is fixed at $r = 0.2$ for all cases.)

Note 2: I verified that the results are converged with the PML thickness in $r$ and $z$.

1. $L = 5$

exteff-1:, 0.999059, (0.05, 0.22, 0.41, 0.58)
exteff-2:, 0.998968, (0.04, 0.22, 0.41, 0.59)

2. $L = 10$

exteff-1:, 0.2, 0.999256, (0.21, 0.46, 0.67, 0.75)
exteff-2:, 0.2, 0.999203, (0.21, 0.46, 0.68, 0.75)

3. $L = 20$

exteff-1:, 0.2, 0.999384, (0.31, 0.56, 0.70, 0.79)
exteff-2:, 0.2, 0.999346, (0.31, 0.56, 0.70, 0.79)

4. $L = 40$

exteff-1:, 0.2, 0.999691, (0.39, 0.59, 0.71, 0.80)
exteff-2:, 0.2, 0.999606, (0.39, 0.59, 0.71, 0.80)

[smartalecH]

This means that the DFT monitor in the $r$ direction becomes further and further away from $r=0$

Hmm. I'm wondering if this is because the monitors are further away, or if it's because the fields are converging better with the larger cells?

This test is turning too many knobs... it's changing the cell size, the size of the disc, and the locations of the monitors. Can we simplify the test matrix here to exactly narrow down what's going on here?

For example, keep the geometry and inner monitor position at the same position, but increase the cell size (and obviously the outer monitor should move with it). Then, you can safely test two things: (1) Does the inner-monitor response change as the cell size increase? (2) does the difference between inner and outer monitor response decrease with cell size?

[stevengj]

There could be a bug in the near2far calculation for the radial-direction interface?

In particular, the cylindrical near2far is implemented by integrating the Cartesian Green's function over φ. For the ±z surfaces, only the position of the near-field source point changes during this integral. However, for the radial surface, the normal vector r̂ changes direction with φ as well, and maybe we aren't taking this into account? It looks like we do, however?

Maybe there is a bug in these lines? But that would affect the ±z surface as well? It would be good to check that these lines are actually being executed — i.e. that we are passing c0 as cylindrical components and not cartesian?

It could also be as simple as the +r surface missing the r factor? i.e. the weight f0 in greencyl (which comes from the DFT weight) should be multiplied by r.

[oskooi]

Maybe there is a bug in these lines? But that would affect the ±z surface as well? It would be good to check that these lines are actually being executed — i.e. that we are passing c0 as cylindrical components and not cartesian?

Checked. These lines are being executed as expected.

It could also be as simple as the +r surface missing the r factor? i.e. the weight f0 in greencyl (which comes from the DFT weight) should be multiplied by r.

Tried scaling the weight f0 by $r$ as suggested but it does not fix the problem.

[stevengj]

Might be worth doing a simple test of cylindrical near2far, just like we did for Cartesian: put a source somewhere (at the center?), compute the field at some point by explicit FDTD, and then compute the field at the same point using near2far. They should match (in the limit of infinite resolution). This can be a relatively small calculation, for a single m.

[oskooi]

Might be worth doing a simple test of cylindrical near2far, just like we did for Cartesian: put a source somewhere (at the center?), compute the field at some point by explicit FDTD, and then compute the field at the same point using near2far. They should match (in the limit of infinite resolution). This can be a relatively small calculation, for a single m.

Setting up such a test involving a point source within a "disc" at a single m demonstrates that there indeed is a discrepancy in the radiated fields in air above the disc computed using the two different methods of: (1) a DFT monitor and (2) a near-to-far field transformation. The discrepancy does not go away with increasing resolution which would be expected if numerical dispersion was the primary cause. I tested this with different values of m and monitor configurations and the results are similar.

The script used to generate to these results is here.

[oskooi]

After further investigation, it looks like the discrepancy between the DFT fields and the near-to-far field transformation only occurs for m = 0. For m ≠ 0, there is excellent agreement. This can be demonstrated for a point source in vacuum (i.e., no disc) above the perfect-metallic reflector ground plane. The previous setup placed the point source within a dielectric disc which was in fact a weak cavity (the radiated flux in air was nearly 100% of the input power of the dipole). We demonstrated in a separate tutorial that when computing the radiated fields of a cavity using the near-to-far field transformation, the enclosing box must be as close to the cavity as possible in order to fully capture its fields. This is why two different near-field monitor configurations were producing different far-field radiation patterns. The monitor configuration furthest from the disc cavity was less accurate.

With the modified setup, here are the results for the radiated fields for m of 0, 1 and 2. The m = 1, 2 case show nearly exact agreement between the DFT fields and near-to-far field transformation. For m = 0, there is a large discrepancy near r = 0. These results were generated using this script.

[stevengj]

The Er fields for m=0 don't make sense. (They should look similar to the m=2 case.) Can you plot the raw time-domain Er field (e.g. in a movie) for the m=0 case to see what they look like, and whether they are going to zero as r goes to zero?

You could also look at other test cases, e.g. if you put an Ez point source at r=0, with m=0, it should match the field of a 3d dipole (i.e. the 3d Green's function, which we know analytically).

[oskooi]

Can you plot the raw time-domain Er field (e.g. in a movie) for the m=0 case to see what they look like, and whether they are going to zero as r goes to zero?

If we take a snapshot of $|E_r|$ on a logscale over the entire cell at some time $t > 0$ and from this data plot a cross section along $r$ at a $z$ corresponding to the point source, we find that the $E_r$ fields do not go to zero at $r = 0$ as would be expected for this field component. This is a clear indication of a bug.

The presence of a bug is perhaps not surprising since currently there is no unit test for the $m = 0$ case.

(Note: the time-domain results are consistent with the DFT fields above which also show a non-zero $E_r$ at $r = 0$.)

The peak in the second plot is the position of the point source. Note that the fields are decaying to zero in the PML region.

Script used to generate these results is here.

[stevengj]

Can you try a really old version of Meep (e.g. via a C++ simulation) and see if it ever worked?

[oskooi]

Looks like the bug is due to the change to update_eh.cpp in #2383 which affects $E_z$ at $r = 0$. Reverting this change produces the expected result.

[stevengj]

It seems like we really want an r=0 version of step_EDHB, which is specialized to assume that there are no off-diagonal terms (since that would break axisymmetry at r=0, regardless of m). Something like this that you call after step_update_EDHB, which we can then call unconditionally (regardless of m or whether it's a cylindrical calculation), where you pass in little_corner and big_corner for is and ie, and not little_corner0:

void step_update_EDHB0(RPR f, component fc, const grid_volume &gv, const ivec is, const ivec ie,
                      const RPR g, const RPR u,  ptrdiff_t s, const RPR chi2,
                      const RPR chi3, RPR fw, direction dsigw, const RPR sigw, const RPR kapw) {
  (void)fc; // currently unused
  if (!f) return;

  // 1.  check that is is cylindrical and starts at r=0, otherwise return
  // 2. change ie.set_direction(R) = is.set_direction(R)   ... this way it will only loop over r=0

  if (dsigw != NO_DIRECTION) { //////// PML case (with fw) /////////////
    KSTRIDE_DEF(dsigw, kw, is, gv);
    { // diagonal u
      if (chi3) {
        {
          PLOOP_OVER_IVECS(gv, is, ie, i) {
            realnum gs = g[i];
            realnum us = u[i];
            DEF_kw;
            realnum fwprev = fw[i], kapwkw = kapw[kw], sigwkw = sigw[kw];
            fw[i] = (gs * us) * calc_nonlinear_u(gs * gs, gs, us, chi2[i], chi3[i]);
            f[i] += (kapwkw + sigwkw) * fw[i] - (kapwkw - sigwkw) * fwprev;
          }
        }
      }
      else if (u) {
        PLOOP_OVER_IVECS(gv, is, ie, i) {
          realnum gs = g[i];
          realnum us = u[i];
          DEF_kw;
          realnum fwprev = fw[i], kapwkw = kapw[kw], sigwkw = sigw[kw];
          fw[i] = (gs * us);
          f[i] += (kapwkw + sigwkw) * fw[i] - (kapwkw - sigwkw) * fwprev;
        }
      }
      else {
        PLOOP_OVER_IVECS(gv, is, ie, i) {
          DEF_kw;
          realnum fwprev = fw[i], kapwkw = kapw[kw], sigwkw = sigw[kw];
          fw[i] = g[i];
          f[i] += (kapwkw + sigwkw) * fw[i] - (kapwkw - sigwkw) * fwprev;
        }
      }
    }
  }
  else {            /////////////// no PML (no fw) ///////////////////
    { // diagonal u
      if (chi3) {
        {
          PLOOP_OVER_IVECS(gv, is, ie, i) {
            realnum gs = g[i];
            realnum us = u[i];
            f[i] = (gs * us) * calc_nonlinear_u(gs * gs, gs, us, chi2[i], chi3[i]);
          }
        }
      }
      else if (u) {
        PLOOP_OVER_IVECS(gv, is, ie, i) {
          realnum gs = g[i];
          realnum us = u[i];
          f[i] = (gs * us);
        }
      }
      else
        PLOOP_OVER_IVECS(gv, is, ie, i) { f[i] = g[i]; }
    }
  }
}