Correct the weights for dipole sources in cylindrical coordinates

src/meep.hpp

@@ -1524,6 +1524,10 @@ class fields_chunk {
   }
 
   double last_source_time();
+
+  // sources.cpp
+  std::complex<double> sum_sources(field_type ft);
+
   // monitor.cpp
   std::complex<double> get_field(component, const ivec &) const;
 
@@ -1761,6 +1765,9 @@ class fields {
   void change_m(double new_m);
   bool has_nonlinearities(bool parallel = true) const;
 
+  // sources.coo
+  std::complex<double> sum_sources(field_type ft);
+
   // time.cpp
   std::vector<double> time_spent_on(time_sink sink);
   double mean_time_spent_on(time_sink);

src/sources.cpp

@@ -270,7 +270,13 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is
     loc += shift * (0.5 * inva);
 
     vec rel_loc = loc - data->center;
-    amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc);
+
+    // in cylindrical coordinates, we need to account for the r term in the integral.
+    if (fc->gv.dim == Dcyl)
+      amps_array[idx_vol] =
+          IVEC_LOOP_WEIGHT(s0, s1, e0, e1, dV0 + dV1 * loop_i2) * amp * data->A(rel_loc);
+    else
+      amps_array[idx_vol] = IVEC_LOOP_WEIGHT(s0, s1, e0, e1, 1) * amp * data->A(rel_loc);
 
     // check for invalid sources at r=0 in cylindrical coordinates
     if (fc->gv.dim == Dcyl && loc.r() == 0 && amps_array[idx_vol] != 0.0) {
@@ -478,9 +484,16 @@ void fields::add_volume_source(component c, const src_time &src, const volume &w
   src_vol_chunkloop_data data;
   data.A = A ? A : one;
   data.amp = amp;
+  // correct units for J delta-function amplitude
   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
+    if (where.in_direction(d) == 0.0 && !nosize_direction(d)) { // delta-fun
+      if (d == R && where.center().in_direction(d) > 0)
+        // in cylindrical coordinates, we need to scale by $r$ for a dipole.
+        data.amp *= 1.0 / where.center().in_direction(d);
+      else if (gv.dim != Dcyl)
+        // we handle the resolution scaling at the chunk level for cylindrical coordinates.
+        data.amp *= gv.a;
+    }
   }
   sources = src.add_to(sources, &data.src);
   data.center = (where.get_min_corner() + where.get_max_corner()) * 0.5;
@@ -755,4 +768,29 @@ diffractedplanewave::diffractedplanewave(int g_[3], double axis_[3], std::comple
   p = p_;
 };
 
+/*
+Occasionally, it's useful to check the sum of all the actual source amplitudes
+(e.g. to ensure that the boundary weights are being computed properly).
+`sum_sources()` provides a convenient interface to do just that.
+*/
+std::complex<double> fields::sum_sources(field_type ft) {
+  std::complex<double> total_sources = 0.0;
+
+  for (int i = 0; i < num_chunks; i++)
+    if (chunks[i]->is_mine()) total_sources += chunks[i]->sum_sources(ft);
+
+  return sum_to_all(total_sources);
+}
+
+std::complex<double> fields_chunk::sum_sources(field_type ft) {
+  std::complex<double> total_sources = 0.0;
+  for (const src_vol &sv : sources[ft]) {
+    for (size_t j = 0; j < sv.num_points(); j++) {
+      meep::master_printf("src %f\n", sv.amplitude_at(j));
+      total_sources += sv.amplitude_at(j);
+    }
+  }
+  return total_sources;
+}
+
 } // namespace meep