make greencyl_tol a parameter (#3064)

* make greencyl_tol a parameter

* add greencyl_tol to wrappers

* work around SWIG bug: make greencyl_tol a required argument

Author: stevengj

python/adjoint/objective.py

@@ -443,6 +443,7 @@ def __init__(
         nperiods: Optional[int] = 1,
         decimation_factor: Optional[int] = 0,
         norm_near_fields: Optional[NearToFarData] = None,
+        greencyl_tol: float = 1e-3,
     ):
         """Initialize an instance of differentiable Fourier fields instance.
 
@@ -470,6 +471,7 @@ def __init__(
         self.decimation_factor = decimation_factor
         self.norm_near_fields = norm_near_fields
         self.nperiods = nperiods
+        self.greencyl_tol = greencyl_tol
 
     def register_monitors(self, frequencies):
         self._frequencies = np.asarray(frequencies)
@@ -501,7 +503,7 @@ def place_adjoint_source(self, dJ):
         )
 
         all_nearsrcdata = self._monitor.swigobj.near_sourcedata(
-            far_pt_vec, farpt_list, self._nfar_pts, dJ
+            far_pt_vec, farpt_list, self._nfar_pts, dJ, self.greencyl_tol
         )
         for near_data in all_nearsrcdata:
             cur_comp = near_data.near_fd_comp

python/meep.i

@@ -329,12 +329,12 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i
 }
 
 // Wrapper around meep::dft_near2far::farfield
-PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) {
+PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v, double greencyl_tol) {
     // Return value: New reference
     Py_ssize_t len = f->freq.size() * 6;
     PyObject *res = PyList_New(len);
 
-    std::complex<double> *ff_arr = f->farfield(v);
+    std::complex<double> *ff_arr = f->farfield(v, greencyl_tol);
 
     for (Py_ssize_t i = 0; i < len; i++) {
         PyList_SetItem(res, i, PyComplex_FromDoubles(ff_arr[i].real(), ff_arr[i].imag()));
@@ -347,13 +347,13 @@ PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v) {
 
 // Wrapper around meep::dft_near2far::get_farfields_array
 PyObject *_get_farfields_array(meep::dft_near2far *n2f, const meep::volume &where,
-                               double resolution) {
+                               double resolution, double greencyl_tol) {
     // Return value: New reference
     size_t dims[4] = {1, 1, 1, 1};
     int rank = 0;
     size_t N = 1;
 
-    double *EH = n2f->get_farfields_array(where, rank, dims, N, resolution);
+    double *EH = n2f->get_farfields_array(where, rank, dims, N, resolution, greencyl_tol);
 
     if (!EH) return PyArray_SimpleNew(0, 0, NPY_CDOUBLE);
 
@@ -662,8 +662,8 @@ PyObject *py_do_harminv(PyObject *vals, double dt, double f_min, double f_max, i
                      double spectral_density, double Q_thresh, double rel_err_thresh,
                      double err_thresh, double rel_amp_thresh, double amp_thresh);
 
-PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v);
-PyObject *_get_farfields_array(meep::dft_near2far *n2f, const meep::volume &where, double resolution);
+PyObject *_get_farfield(meep::dft_near2far *f, const meep::vec & v, double greencyl_tol);
+PyObject *_get_farfields_array(meep::dft_near2far *n2f, const meep::volume &where, double resolution, double greencyl_tol);
 PyObject *_dft_ldos_ldos(meep::dft_ldos *f);
 PyObject *_dft_ldos_F(meep::dft_ldos *f);
 PyObject *_dft_ldos_J(meep::dft_ldos *f);

python/simulation.py

@@ -3215,7 +3215,7 @@ def load_minus_energy(self, fname: str, energy: DftEnergy):
         self.load_energy(fname, energy)
         energy.scale_dfts(-1.0)
 
-    def get_farfield(self, near2far, x):
+    def get_farfield(self, near2far, x, greencyl_tol: float = 1e-3):
         """
         Given a `Vector3` point `x` which can lie anywhere outside the near-field surface,
         including outside the cell and a `near2far` object, returns the computed
@@ -3228,6 +3228,7 @@ def get_farfield(self, near2far, x):
         return mp._get_farfield(
             near2far.swigobj,
             py_v3_to_vec(self.dimensions, x, is_cylindrical=self.is_cylindrical),
+            greencyl_tol,
         )
 
     def get_farfields(
@@ -3237,6 +3238,7 @@ def get_farfields(
         where: Volume = None,
         center: Vector3Type = None,
         size: Vector3Type = None,
+        greencyl_tol: float = 1e-3,
     ):
         """
         Like `output_farfields` but returns a dictionary of NumPy arrays instead of
@@ -3253,7 +3255,9 @@ def get_farfields(
             self.init_sim()
         vol = self._volume_from_kwargs(where, center, size)
         self.fields.am_now_working_on(mp.GetFarfieldsTime)
-        result = mp._get_farfields_array(near2far.swigobj, vol, resolution)
+        result = mp._get_farfields_array(
+            near2far.swigobj, vol, resolution, greencyl_tol
+        )
         self.fields.finished_working()
         res_ex = complexarray(result[0], result[1])
         res_ey = complexarray(result[2], result[3])

src/meep.hpp

@@ -1360,17 +1360,17 @@ class dft_near2far {
   dft_near2far(const dft_near2far &f);
 
   /* return an array (Ex,Ey,Ez,Hx,Hy,Hz) x Nfreq of the far fields at x */
-  std::complex<double> *farfield(const vec &x);
+  std::complex<double> *farfield(const vec &x, double greencyl_tol = 1e-3);
 
   /* like farfield, but requires F to be Nfreq*6 preallocated array, and
      does *not* perform the reduction over processes...an MPI allreduce
      summation by the caller is required to get the final result ... used
      by other output routine to efficiently get far field on a grid of pts */
-  void farfield_lowlevel(std::complex<double> *F, const vec &x);
+  void farfield_lowlevel(std::complex<double> *F, const vec &x, double greencyl_tol = 1e-3);
 
   /* Return a newly allocated array with all far fields */
   double *get_farfields_array(const volume &where, int &rank, size_t *dims, size_t &N,
-                              double resolution);
+                              double resolution, double greencyl_tol = 1e-3);
 
   /* output far fields on a grid to an HDF5 file */
   void save_farfields(const char *fname, const char *prefix, const volume &where,
@@ -1400,7 +1400,7 @@ class dft_near2far {
   double periodic_k[2], period[2];
 
   std::vector<sourcedata> near_sourcedata(const vec &x_0, double *farpt_list, size_t nfar_pts,
-                                          const std::complex<double> *dJ);
+                                          const std::complex<double> *dJ, double greencyl_tol);
 };
 
 /* Class to compute local-density-of-states spectra: the power spectrum

src/near2far.cpp

@@ -348,7 +348,7 @@ void greencyl(std::complex<double> *EH, const vec &x, double freq, double eps, d
   }
 }
 
-void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
+void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x, double greencyl_tol) {
   if (x.dim != D3 && x.dim != D2 && x.dim != Dcyl)
     meep::abort("only 2d or 3d or cylindrical far-field computation is supported");
   greenfunc green = x.dim == D2 ? green2d : green3d;
@@ -384,7 +384,7 @@ void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
             std::complex<double> cphase = std::polar(1.0, phase);
             if (x.dim == Dcyl)
               greencyl(EH6, x, freq[i], eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i], f->fc->m,
-                       1e-3);
+                       greencyl_tol);
             else
               green(EH6, x, freq[i], eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]);
             for (int j = 0; j < 6; ++j)
@@ -397,19 +397,19 @@ void dft_near2far::farfield_lowlevel(std::complex<double> *EH, const vec &x) {
   }
 }
 
-std::complex<double> *dft_near2far::farfield(const vec &x) {
+std::complex<double> *dft_near2far::farfield(const vec &x, double greencyl_tol) {
   std::complex<double> *EH, *EH_local;
   const size_t Nfreq = freq.size();
   EH_local = new std::complex<double>[6 * Nfreq];
-  farfield_lowlevel(EH_local, x);
+  farfield_lowlevel(EH_local, x, greencyl_tol);
   EH = new std::complex<double>[6 * Nfreq];
   sum_to_all(EH_local, EH, 6 * Nfreq);
   delete[] EH_local;
   return EH;
 }
 
 double *dft_near2far::get_farfields_array(const volume &where, int &rank, size_t *dims, size_t &N,
-                                          double resolution) {
+                                          double resolution, double greencyl_tol) {
   /* compute output grid size etc. */
   double dx[3] = {0, 0, 0};
   direction dirs[3] = {X, Y, Z};
@@ -457,7 +457,7 @@ double *dft_near2far::get_farfields_array(const volume &where, int &rank, size_t
           start = t;
           last_point = this_point;
         }
-        farfield_lowlevel(EH1, x);
+        farfield_lowlevel(EH1, x, greencyl_tol);
         if (verbosity > 1) all_wait(); // Allow consistent progress updates from master
         ptrdiff_t idx = (i0 * dims[1] + i1) * dims[2] + i2;
         for (size_t i = 0; i < Nfreq; ++i)
@@ -655,7 +655,8 @@ dft_near2far fields::add_dft_near2far(const volume_list *where, const double *fr
 // Modified from farfield_lowlevel
 std::vector<struct sourcedata> dft_near2far::near_sourcedata(const vec &x_0, double *farpt_list,
                                                              size_t nfar_pts,
-                                                             const std::complex<double> *dJ) {
+                                                             const std::complex<double> *dJ,
+                                                             double greencyl_tol) {
   if (x_0.dim != D3 && x_0.dim != D2 && x_0.dim != Dcyl)
     meep::abort("only 2d or 3d or cylindrical far-field computation is supported");
   greenfunc green = x_0.dim == D2 ? green2d : green3d;
@@ -697,7 +698,7 @@ std::vector<struct sourcedata> dft_near2far::near_sourcedata(const vec &x_0, dou
             for (size_t ipt = 0; ipt < nfar_pts; ++ipt) {
               vec x = vec(farpt_list[3 * ipt], farpt_list[3 * ipt + 1], farpt_list[3 * ipt + 2]);
               if (x_0.dim == Dcyl)
-                greencyl(EH6, x, freq[i], eps, mu, xs, c0, w, f->fc->m, 1e-3);
+                greencyl(EH6, x, freq[i], eps, mu, xs, c0, w, f->fc->m, greencyl_tol);
               else
                 green(EH6, x, freq[i], eps, mu, xs, c0, w);
               for (int j = 0; j < 6; ++j)