WIP: near2far for cylindrical coordinates. (#1090)

* initial stab at a greencyl function

* whoops

* formatting

* use meep::pi

* fix get_farfields_array in cylindrical coords

* slight code cleanup

* rm redundant 2pi factor

* smaller N0

* add some comments to greencyl

* typo in comment

* whoops, 2/N is integer division in C

* add green functions to meep internals header in case we want to export them some day, and to make testing easier

* whoops, fix evaluation angles

* clarify normalization of greencyl

src/meep_internals.hpp

@@ -151,4 +151,12 @@ void step_beta_stride1(realnum *f, component c, const realnum *g, const grid_vol
       step_beta(f, c, g, gv, betadt, dsig, siginv, fu, dsigu, siginvu, cndinv, fcnd);              \
   } while (0)
 
+// analytical Green's functions from near2far.cpp, which we might want to expose someday
+void green3d(std::complex<double> *EH, const vec &x, double freq, double eps, double mu,
+             const vec &x0, component c0, std::complex<double> f0);
+void green2d(std::complex<double> *EH, const vec &x, double freq, double eps, double mu,
+             const vec &x0, component c0, std::complex<double> f0);
+void greencyl(std::complex<double> *EH, const vec &x, double freq, double eps, double mu,
+              const vec &x0, component c0, std::complex<double> f0, double m, double tol);
+
 } // namespace meep

src/near2far.cpp

@@ -20,7 +20,7 @@
    compute the fields on a "far" surface via the homogeneous-medium Green's
    function in 2d or 3d. */
 
-#include <meep.hpp>
+#include "meep_internals.hpp"
 #include <assert.h>
 #include "config.h"
 #include <math.h>
@@ -245,8 +245,74 @@ void green2d(std::complex<double> *EH, const vec &x, double freq, double eps, do
   }
 }
 
+// cylindrical Green's function constructed by integrating green3d as the source
+// term rotates around the z axis with exp(im*phi) dependence, integrated to a tolerance tol.
+// (note: this is the Green's function divided by 2pi*x0.r(), to compensate for a 2piR factor
+//  in the near2far add_dft weight.)
+void greencyl(std::complex<double> *EH, const vec &x, double freq, double eps, double mu,
+              const vec &x0, component c0, std::complex<double> f0, double m, double tol) {
+  if (x0.dim != Dcyl) abort("wrong dimensionality in greencyl");
+  vec x_3d(x.dim == Dcyl ? x.r() : x.x(), x.y(), x.z());
+  direction d = component_direction(c0);
+  component cx = direction_component(c0, X), cy = direction_component(c0, Y);
+  for (int j = 0; j < 6; ++j)
+    EH[j] = 0;
+
+  /* Perform phi integral.  Since phi integrand is smooth, quadrature with equally spaced points
+     should converge exponentially fast with the number N of quadrature points.  We
+     repeatedly double N until convergence to tol is achieved, re-using previous points. */
+  const int N0 = 4;
+  double dphi = 2.0 / N0; // factor of 2*pi*r is already included in add_dft weight
+  for (int N = N0; N <= 65536; N *= 2) {
+    std::complex<double> EH_sum[6];
+    dphi *= 0.5; // delta phi is halved because N doubles
+    double dphi2pi = dphi * 2*pi;
+    for (int j = 0; j < 6; ++j)
+      EH_sum[j] = EH[j] * 0.5; // re-use previous quadrature points (with halved dphi)
+    /* N-point quadrature points i = 0..N-1.  After the first iteration (N==N0), we
+       only need to sum over odd i, since the even i were summed for the previous N. */
+    for (int i = (N > N0); i < N; i += 1 + (N > N0)) {
+      double phi = i * dphi2pi, c = cos(phi), s = sin(phi);
+      vec x0_phi(x0.r() * c, x0.r() * s, x0.z()); // source point rotated by phi
+      std::complex<double> EH_phi[6], f0_exp_imphi = f0 * std::polar(1.0, m * phi) * dphi;
+      /* if the source direction is in the r or phi directions, then we must rotate
+        the direction of the source current in the xy plane as well */
+      if (d == Z) { // source currents in z direction don't rotate
+        green3d(EH_phi, x_3d, freq, eps, mu, x0_phi, c0, f0_exp_imphi);
+        for (int j = 0; j < 6; ++j)
+          EH_sum[j] += EH_phi[j];
+      }
+      else if (d == R) { // r_hat = c x_hat + s y_hat
+        green3d(EH_phi, x_3d, freq, eps, mu, x0_phi, cx, f0_exp_imphi * c);
+        for (int j = 0; j < 6; ++j)
+          EH_sum[j] += EH_phi[j];
+        green3d(EH_phi, x_3d, freq, eps, mu, x0_phi, cy, f0_exp_imphi * s);
+        for (int j = 0; j < 6; ++j)
+          EH_sum[j] += EH_phi[j];
+      }
+      else { // (d == P):  phi_hat = c y_hat - s x_hat
+        green3d(EH_phi, x_3d, freq, eps, mu, x0_phi, cx, f0_exp_imphi * (-s));
+        for (int j = 0; j < 6; ++j)
+          EH_sum[j] += EH_phi[j];
+        green3d(EH_phi, x_3d, freq, eps, mu, x0_phi, cy, f0_exp_imphi * c);
+        for (int j = 0; j < 6; ++j)
+          EH_sum[j] += EH_phi[j];
+      }
+    }
+    // accumulate the new and old sums and check how much the integral has changed in L1 norm
+    double sumdiff = 0, sumabs = 0;
+    for (int j = 0; j < 6; ++j) {
+      sumdiff += abs(EH[j] - EH_sum[j]);
+      sumabs += abs(EH_sum[j]);
+      EH[j] = EH_sum[j];
+    }
+    if (sumdiff <= sumabs * tol) break; // doubling N changed sum by less than tol
+  }
+}
+
 void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
-  if (x.dim != D3 && x.dim != D2) abort("only 2d or 3d far-field computation is supported");
+  if (x.dim != D3 && x.dim != D2 && x.dim != Dcyl)
+    abort("only 2d or 3d or cylindrical far-field computation is supported");
   greenfunc green = x.dim == D2 ? green2d : green3d;
 
   for (int i = 0; i < 6 * Nfreq; ++i)
@@ -278,7 +344,10 @@ void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
               xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1 * period[1]);
             double phase = phase0 + i1 * periodic_k[1];
             std::complex<double> cphase = std::polar(1.0, phase);
-            green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]);
+            if (x.dim == Dcyl)
+              greencyl(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i], f->fc->m, 1e-3);
+            else
+              green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]);
             for (int j = 0; j < 6; ++j)
               EH[i * 6 + j] += EH6[j] * cphase;
           }
@@ -305,17 +374,19 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_
   double dx[3] = {0, 0, 0};
   direction dirs[3] = {X, Y, Z};
 
+  rank = 0;
+  N = dims[0] = dims[1] = dims[2] = 1;
+
   LOOP_OVER_DIRECTIONS(where.dim, d) {
     dims[rank] = int(floor(where.in_direction(d) * resolution));
-    if (dims[rank] <= 1) {
+    if (dims[rank] <= 1)
       dims[rank] = 1;
-      dx[rank] = 0;
-    }
     else
       dx[rank] = where.in_direction(d) / (dims[rank] - 1);
     N *= dims[rank];
     dirs[rank++] = d;
   }
+  if (where.dim == Dcyl) dirs[2] = P; // otherwise Z is listed twice
 
   if (N * Nfreq < 1) return NULL; /* nothing to output */
 
@@ -327,7 +398,6 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_
   std::complex<double> *EH1 = new std::complex<double>[6 * Nfreq];
 
   double start = wall_time();
-  size_t total_points = dims[0] * dims[1] * dims[2];
   size_t last_point = 0;
 
   vec x(where.dim);
@@ -339,9 +409,9 @@ realnum *dft_near2far::get_farfields_array(const volume &where, int &rank, size_
         x.set_direction(dirs[2], where.in_direction_min(dirs[2]) + i2 * dx[2]);
         double t;
         if (verbosity > 0 && (t = wall_time()) > start + MEEP_MIN_OUTPUT_TIME) {
-          size_t this_point = (dims[1] * dims[2] * i0) + (dims[2] * i1) + i2 + 1;
+          size_t this_point = (dims[1] * i0 + i1) * dims[2] + i2 + 1;
           master_printf("get_farfields_array working on point %zu of %zu (%d%% done), %g s/point\n",
-                        this_point, total_points, (int)((double)this_point / total_points * 100),
+                        this_point, N, (int)((double)this_point / N * 100),
                         (t - start) / (std::max(1, (int)(this_point - last_point))));
           start = t;
           last_point = this_point;