Fix memory leaks in src_vol (NanoComp#1644)

The old implementation's weak point was the `add_to` method, which
sometimes transferred memory ownership correctly and sometimes did not.

* Rewrite `src_vol` using STL containers.
* Replace the linked list of source volumes in `fields_chunk` with a
vector. Adjust code that iterates for source volumes accordingly.
* Exclude `fields_chunk` from SWIG wrappers.

python/meep.i

@@ -1470,6 +1470,7 @@ void _get_gradient(PyObject *grad, PyObject *fields_a, PyObject *fields_f, PyObj
 %ignore is_medium;
 %ignore is_metal;
 %ignore meep::choose_chunkdivision;
+%ignore meep::fields_chunk;
 %ignore meep::infinity;
 
 %ignore std::vector<meep::volume>::vector(size_type);

scheme/Makefile.am

@@ -39,7 +39,7 @@ endif
 
 # workaround missing namespace prefix in swig
 meep_renames.i: $(LIBHDRS)
-	(echo "// AUTOMATICALLY GENERATED -- DO NOT EDIT"; perl -pe 's/^ *class +([A-Za-z_0-9:]*)( *| *:[^{]*){.*$$/%rename(meep_\1) meep::\1;/' $(LIBHDRS) | grep "%rename" | sort -u; echo; grep -hv typedef $(LIBHDRS) | perl -pe 's/(inline|const|extern|static) +//g' | perl -pe 's/^[A-Za-z_0-9:<>]+[* ]+([A-Za-z_0-9:]*) *\(.*$$/%rename(meep_\1) meep::\1;/' | grep "%rename" | sed '/choose_chunkdivision/d' | sort -u; ) > $@
+	(echo "// AUTOMATICALLY GENERATED -- DO NOT EDIT"; perl -pe 's/^ *class +([A-Za-z_0-9:]*)( *| *:[^{]*){.*$$/%rename(meep_\1) meep::\1;/' $(LIBHDRS) | grep "%rename" | sed '/meep_fields_chunk/d' | sort -u; echo; grep -hv typedef $(LIBHDRS) | perl -pe 's/(inline|const|extern|static) +//g' | perl -pe 's/^[A-Za-z_0-9:<>]+[* ]+([A-Za-z_0-9:]*) *\(.*$$/%rename(meep_\1) meep::\1;/' | grep "%rename" | sed '/choose_chunkdivision/d' | sort -u; ) > $@
 
 # work around bug in swig, where it doesn't prepend namespace to friend funcs
 meep_swig_bug_workaround.i: $(LIBHDRS)

scheme/meep.i

@@ -281,6 +281,7 @@ static meep::vec my_kpoint_func(double freq, int mode, void *user_data) {
 %warnfilter(302,325,451,503,509);
 
 %ignore meep::choose_chunkdivision;
+%ignore meep::fields_chunk;
 %ignore meep_geom::fragment_stats;
 
 %include "meep_renames.i"

src/array_slice.cpp

@@ -209,18 +209,18 @@ static void get_source_slice_chunkloop(fields_chunk *fc, int ichnk, component cg
   if (has_direction(dim, Y)) NY = ((slice_imax - slice_imin).in_direction(Y) / 2) + 1;
 
   for (int ft = 0; ft < NUM_FIELD_TYPES; ft++)
-    for (src_vol *s = fc->sources[ft]; s; s = s->next) {
+    for (const src_vol &s : fc->get_sources(static_cast<field_type>(ft))) {
       component cS = S.transform(data->source_component, -sn);
-      if (s->c != cS) continue;
+      if (s.c != cS) continue;
 
       // loop over point sources in this src_vol. for each point source,
       // the src_vol stores the amplitude and the global index of the
       // symmetry-parent grid point, from which we need to compute the
       // local index of the symmetry-child grid point within this
       // slice (that is, if it even lies within the slice)
-      for (size_t npt = 0; npt < s->npts; npt++) {
-        complex<double> amp = s->A[npt];
-        ptrdiff_t chunk_index = s->index[npt];
+      for (size_t npt = 0; npt < s.num_points(); npt++) {
+        const complex<double> &amp = s.amplitude_at(npt);
+        ptrdiff_t chunk_index = s.index_at(npt);
         ivec iloc_parent = fc->gv.iloc(Dielectric, chunk_index);
         ivec iloc_child = S.transform(iloc_parent, sn) + shift;
         if (!in_subgrid(slice_imin, iloc_child, slice_imax)) continue; // source point outside slice

src/dft_ldos.cpp

@@ -102,41 +102,41 @@ void dft_ldos::update(fields &f) {
 
   for (int ic = 0; ic < f.num_chunks; ic++)
     if (f.chunks[ic]->is_mine()) {
-      for (src_vol *sv = f.chunks[ic]->sources[D_stuff]; sv; sv = sv->next) {
-        component c = direction_component(Ex, component_direction(sv->c));
+      for (const src_vol &sv : f.chunks[ic]->get_sources(D_stuff)) {
+        component c = direction_component(Ex, component_direction(sv.c));
         realnum *fr = f.chunks[ic]->f[c][0];
         realnum *fi = f.chunks[ic]->f[c][1];
         if (fr && fi) // complex E
-          for (size_t j = 0; j < sv->npts; j++) {
-            const ptrdiff_t idx = sv->index[j];
-            const complex<double> A = sv->A[j];
+          for (size_t j = 0; j < sv.num_points(); j++) {
+            const ptrdiff_t idx = sv.index_at(j);
+            const complex<double>& A = sv.amplitude_at(j);
             EJ += complex<double>(fr[idx], fi[idx]) * conj(A);
             Jsum += abs(A);
           }
         else if (fr) { // E is purely real
-          for (size_t j = 0; j < sv->npts; j++) {
-            const ptrdiff_t idx = sv->index[j];
-            const complex<double> A = sv->A[j];
+          for (size_t j = 0; j < sv.num_points(); j++) {
+            const ptrdiff_t idx = sv.index_at(j);
+            const complex<double>& A = sv.amplitude_at(j);
             EJ += double(fr[idx]) * conj(A);
             Jsum += abs(A);
           }
         }
       }
-      for (src_vol *sv = f.chunks[ic]->sources[B_stuff]; sv; sv = sv->next) {
-        component c = direction_component(Hx, component_direction(sv->c));
+      for (const src_vol& sv : f.chunks[ic]->get_sources(B_stuff)) {
+        component c = direction_component(Hx, component_direction(sv.c));
         realnum *fr = f.chunks[ic]->f[c][0];
         realnum *fi = f.chunks[ic]->f[c][1];
         if (fr && fi) // complex H
-          for (size_t j = 0; j < sv->npts; j++) {
-            const ptrdiff_t idx = sv->index[j];
-            const complex<double> A = sv->A[j];
+          for (size_t j = 0; j < sv.num_points(); j++) {
+            const ptrdiff_t idx = sv.index_at(j);
+            const complex<double>& A = sv.amplitude_at(j);
             HJ += complex<double>(fr[idx], fi[idx]) * conj(A);
             Jsum += abs(A);
           }
         else if (fr) { // H is purely real
-          for (size_t j = 0; j < sv->npts; j++) {
-            const ptrdiff_t idx = sv->index[j];
-            const complex<double> A = sv->A[j];
+          for (size_t j = 0; j < sv.num_points(); j++) {
+            const ptrdiff_t idx = sv.index_at(j);
+            const complex<double>& A = sv.amplitude_at(j);
             HJ += double(fr[idx]) * conj(A);
             Jsum += abs(A);
           }

src/fields.cpp

@@ -15,6 +15,8 @@
 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
 
+#include <algorithm>
+#include <utility>
 #include <stdio.h>
 #include <string.h>
 #include <math.h>
@@ -197,7 +199,6 @@ fields_chunk::~fields_chunk() {
     dft_chunks = nxt;
   }
   FOR_FIELD_TYPES(ft) {
-    delete sources[ft];
     delete[] zeroes[ft];
   }
   FOR_FIELD_TYPES(ft) {
@@ -246,7 +247,6 @@ fields_chunk::fields_chunk(structure_chunk *the_s, const char *od, double m, dou
   }
   doing_solve_cw = false;
   solve_cw_omega = 0.0;
-  FOR_FIELD_TYPES(ft) { sources[ft] = NULL; }
   FOR_COMPONENTS(c) DOCMP2 {
     f[c][cmp] = NULL;
     f_u[c][cmp] = NULL;
@@ -308,7 +308,6 @@ fields_chunk::fields_chunk(const fields_chunk &thef, int chunkidx) : gv(thef.gv)
   }
   doing_solve_cw = thef.doing_solve_cw;
   solve_cw_omega = thef.solve_cw_omega;
-  FOR_FIELD_TYPES(ft) { sources[ft] = NULL; }
   FOR_COMPONENTS(c) DOCMP2 {
     f[c][cmp] = NULL;
     f_u[c][cmp] = NULL;
@@ -484,8 +483,9 @@ void fields::require_source_components() {
   memset(needed, 0, sizeof(needed));
   for (int i = 0; i < num_chunks; i++) {
     FOR_FIELD_TYPES(ft) {
-      for (src_vol *src = chunks[i]->sources[ft]; src; src = src->next)
-        needed[src->c] = 1;
+      for (const auto& src : chunks[i]->get_sources(ft)) {
+        needed[src.c] = 1;
+      }
     }
   }
   int allneeded[NUM_FIELD_COMPONENTS];
@@ -544,10 +544,22 @@ void fields::_require_component(component c, bool aniso2d) {
   }
 }
 
+void fields_chunk::add_source(field_type ft, src_vol &&src) {
+  auto it = std::find_if(sources[ft].begin(), sources[ft].end(), [&src](const src_vol &other) {
+    return src_vol::combinable(src, other);
+  });
+
+  if (it != sources[ft].end()) {
+    it->add_amplitudes_from(src);
+    return;
+  }
+
+  sources[ft].push_back(std::move(src));
+}
+
 void fields_chunk::remove_sources() {
   FOR_FIELD_TYPES(ft) {
-    delete sources[ft];
-    sources[ft] = NULL;
+    sources[ft].clear();
   }
 }
 

src/meep.hpp

@@ -1370,7 +1370,7 @@ class fields_chunk {
   bool zero_fields_near_cylorigin; // fields=0 m pixels near r=0 for stability
   double beta;
   int is_real;
-  src_vol *sources[NUM_FIELD_TYPES];
+  std::vector<src_vol> sources[NUM_FIELD_TYPES];
   structure_chunk *new_s;
   structure_chunk *s;
   const char *outdir;
@@ -1398,6 +1398,10 @@ class fields_chunk {
 
   std::complex<double> get_chi1inv(component, direction, const ivec &iloc,
                                    double frequency = 0) const;
+  // Returns the vector of sources volumes for field type `ft`.
+  const std::vector<src_vol> &get_sources(field_type ft) const { return sources[ft]; }
+  // Adds a source volume of field type `ft` and takes ownership of `src`.
+  void add_source(field_type ft, src_vol &&src);
 
   void backup_component(component c);
   void average_with_backup(component c);

src/meep_internals.hpp

@@ -15,6 +15,7 @@
 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
 
+#include <vector>
 #include <string.h>
 #include "meep.hpp"
 
@@ -53,35 +54,37 @@ inline int rmin_bulk(int m) {
   return r;
 }
 
+// A source volume
+// Moveable and copyable
 class src_vol {
 public:
-  src_vol(component cc, src_time *st, size_t n, ptrdiff_t *ind, std::complex<double> *amps);
-  src_vol(const src_vol &sv);
-  ~src_vol() {
-    delete next;
-    delete[] index;
-    delete[] A;
-  }
-
-  src_time *t;
-  ptrdiff_t *index;        // list of locations of sources in grid (indices)
-  size_t npts;             // number of points in list
-  component c;             // field component the source applies to
-  std::complex<double> *A; // list of amplitudes
-
-  std::complex<double> dipole(size_t j) { return A[j] * t->dipole(); }
-  std::complex<double> current(size_t j) { return A[j] * t->current(); }
-  void update(double time, double dt) { t->update(time, dt); }
-
-  bool operator==(const src_vol &sv) const {
-    // note: don't compare sv.A, since this is used to see if we can just
-    // add one source's amplitudes to another in src_vol::add_to
-    return sv.npts == npts && sv.c == c && sv.t == t &&
-           memcmp(sv.index, index, npts * sizeof(ptrdiff_t)) == 0;
-  }
-
-  src_vol *add_to(src_vol *others);
-  src_vol *next;
+  // Constructs a new source volume. Takes ownership of `ind` and `amps`.
+  // Requirement: ind.size() == amps.size()
+  src_vol(component cc, src_time *st, std::vector<ptrdiff_t> &&ind,
+          std::vector<std::complex<double> > &&amps);
+
+  // Checks whether `a` and `b` are combinable, i.e. have the same indices and point to the same
+  // `src_time` instance, but have potentially different amplitudes.
+  static bool combinable(const src_vol &a, const src_vol &b);
+
+  ptrdiff_t index_at(size_t pos) const { return index[pos]; }
+  const std::complex<double> &amplitude_at(size_t pos) const { return amp[pos]; }
+  size_t num_points() const { return index.size(); };
+  const src_time *t() const { return src_t; };
+
+  std::complex<double> dipole(size_t j) const { return amp[j] * src_t->dipole(); };
+  std::complex<double> current(size_t j) const { return amp[j] * src_t->current(); };
+  void update(double time, double dt) { src_t->update(time, dt); };
+  // Merges the amplitudes from volume source `other` into `this`.
+  // Requirement: other.num_points() == this->num_points().
+  // It is recommended to use `combinable` before calling this method.
+  void add_amplitudes_from(const src_vol &other);
+
+  const component c; // field component the source applies to
+private:
+  src_time *src_t;                        // Not owned by us.
+  std::vector<ptrdiff_t> index;           // locations of sources in grid (indices)
+  std::vector<std::complex<double> > amp; // amplitudes
 };
 
 const int num_bandpts = 32;

src/sources.cpp

@@ -15,6 +15,8 @@
 %  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 */
 
+#include <cassert>
+#include <utility>
 #include <stdio.h>
 #include <stdlib.h>
 #include <math.h>
@@ -152,52 +154,25 @@ bool custom_src_time::is_equal(const src_time &t) const {
 
 /*********************************************************************/
 
-src_vol::src_vol(component cc, src_time *st, size_t n, ptrdiff_t *ind, complex<double> *amps) {
-  c = cc;
-  if (is_D(c)) c = direction_component(Ex, component_direction(c));
-  if (is_B(c)) c = direction_component(Hx, component_direction(c));
-  t = st;
-  next = NULL;
-  npts = n;
-  index = ind;
-  A = amps;
+src_vol::src_vol(component cc, src_time *st, std::vector<ptrdiff_t> &&ind,
+                 std::vector<std::complex<double> > &&amps)
+    : c([](component c) -> component {
+        if (is_D(c)) c = direction_component(Ex, component_direction(c));
+        if (is_B(c)) c = direction_component(Hx, component_direction(c));
+        return c;
+      }(cc)),
+      src_t(st), index(std::move(ind)), amp(std::move(amps)) {
+  assert(index.size() == amp.size());
 }
 
-src_vol::src_vol(const src_vol &sv) {
-  c = sv.c;
-  t = sv.t;
-  npts = sv.npts;
-  index = new ptrdiff_t[npts];
-  A = new complex<double>[npts];
-  for (size_t j = 0; j < npts; j++) {
-    index[j] = sv.index[j];
-    A[j] = sv.A[j];
-  }
-  if (sv.next)
-    next = new src_vol(*sv.next);
-  else
-    next = NULL;
+bool src_vol::combinable(const src_vol &a, const src_vol &b) {
+  return (a.c == b.c) && (a.src_t == b.src_t) && (a.index == b.index);
 }
 
-src_vol *src_vol::add_to(src_vol *others) {
-  if (others) {
-    if (*this == *others) {
-      if (npts != others->npts)
-        meep::abort("Cannot add grid_volume sources with different number of points\n");
-      /* Compare all of the indices...if this ever becomes too slow,
-         we can just compare the first and last indices. */
-      for (size_t j = 0; j < npts; j++) {
-        if (others->index[j] != index[j]) meep::abort("Different indices\n");
-        others->A[j] += A[j];
-      }
-    }
-    else
-      others->next = add_to(others->next);
-    return others;
-  }
-  else {
-    next = others;
-    return this;
+void src_vol::add_amplitudes_from(const src_vol &other) {
+  assert(amp.size() == other.num_points());
+  for (size_t i = 0; i < amp.size(); ++i) {
+    amp[i] += other.amp[i];
   }
 }
 
@@ -267,8 +242,8 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is
 
   size_t npts = 1;
   LOOP_OVER_DIRECTIONS(is.dim, d) { npts *= (ie.in_direction(d) - is.in_direction(d)) / 2 + 1; }
-  ptrdiff_t *index_array = new ptrdiff_t[npts];
-  complex<double> *amps_array = new complex<double>[npts];
+  std::vector<ptrdiff_t> index_array(npts);
+  std::vector<complex<double>> amps_array(npts);
 
   complex<double> amp = data->amp * conj(shift_phase);
 
@@ -299,26 +274,22 @@ static void src_vol_chunkloop(fields_chunk *fc, int ichunk, component c, ivec is
 
   if (idx_vol > npts) meep::abort("add_volume_source: computed wrong npts (%zd vs. %zd)", npts, idx_vol);
 
-  src_vol *tmp = new src_vol(c, data->src, idx_vol, index_array, amps_array);
   field_type ft = is_magnetic(c) ? B_stuff : D_stuff;
-  fc->sources[ft] = tmp->add_to(fc->sources[ft]);
+  fc->add_source(ft, src_vol(c, data->src, std::move(index_array), std::move(amps_array)));
 }
 
 void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex<double>* amp_arr){
   sources = src->add_to(sources, &src);
-  ptrdiff_t* index_arr = new ptrdiff_t[n];
-  std::complex<double>* amp_arr_copy = new std::complex<double>[n];
-  for (size_t i = 0; i < n; i++){
-    index_arr[i] = cur_data.idx_arr[i];
-    amp_arr_copy[i] = amp_arr[i];
-  }
+  std::vector<ptrdiff_t> index_arr(cur_data.idx_arr);
+  std::vector<std::complex<double>> amplitudes(amp_arr, amp_arr+n);
   component c = cur_data.near_fd_comp;
-  src_vol *tmp = new src_vol(c, src, n, index_arr, amp_arr_copy);
+
   field_type ft = is_magnetic(c) ? B_stuff : D_stuff;
   if (0 > cur_data.fc_idx or cur_data.fc_idx >= num_chunks) meep::abort("fields chunk index out of range");
   fields_chunk *fc = chunks[cur_data.fc_idx];
   if (!fc->is_mine()) meep::abort("wrong fields chunk");
-  fc->sources[ft] = tmp->add_to(fc->sources[ft]);
+
+  fc->add_source(ft, src_vol(c, src, std::move(index_arr), std::move(amplitudes)));
   // We can't do require_component(c) since that only works if all processes are adding
   // srcdata for the same components in the same order, which may not be true.
   // ... instead, the caller should call fields::require_source_components()