fix deadlock in add_srcdata via new require_source_components() funct…

…ion (#1521)

* fix deadlock in add_srcdata via new require_source_components() function

* typo

* rm redundant reduction

* timing

* don't call self.fields.require_source_components if there were no sources

* actually we should call init_sim in add_sources, in case only some processes have sources

python/simulation.py

@@ -1930,8 +1930,7 @@ def use_real(self):
             v = Vector3(self.k_point.x, self.k_point.y) if self.special_kz else self.k_point
             self.fields.use_bloch(py_v3_to_vec(self.dimensions, v, self.is_cylindrical))
 
-        for s in self.sources:
-            self.add_source(s)
+        self.add_sources()
 
         for hook in self.init_sim_hooks:
             hook()
@@ -2368,6 +2367,13 @@ def add_source(self, src):
             else:
                 add_vol_src(src.amplitude * 1.0)
 
+    def add_sources(self):
+        if self.fields is None:
+            self.init_sim() # in case only some processes have IndexedSources
+        for s in self.sources:
+            self.add_source(s)
+        self.fields.require_source_components() # needed by IndexedSource objects
+
     def _evaluate_dft_objects(self):
         for dft in self.dft_objects:
             if dft.swigobj is None:
@@ -3567,8 +3573,7 @@ def change_sources(self, new_sources):
         self.sources = new_sources
         if self.fields:
             self.fields.remove_sources()
-            for s in self.sources:
-                self.add_source(s)
+            self.add_sources()
 
     def reset_meep(self):
         """

src/fields.cpp

@@ -472,15 +472,30 @@ bool fields_chunk::alloc_f(component c) {
   return changed;
 }
 
-void fields::require_component(component c) {
-  if (!gv.has_field(c))
-    abort("cannot require a %s component in a %s grid", component_name(c), dimension_name(gv.dim));
-
-  if (beta != 0 && gv.dim != D2) abort("Nonzero beta unsupported in dimensions other than 2.");
+// allocate fields for components required by any source on any process
+// ... this is needed after calling the low-level fields::add_srcdata
+void fields::require_source_components() {
+  int needed[NUM_FIELD_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;
+    }
+  }
+  int allneeded[NUM_FIELD_COMPONENTS];
+  am_now_working_on(MpiAllTime);
+  or_to_all(needed, allneeded, NUM_FIELD_COMPONENTS);
+  finished_working();
 
-  components_allocated = true;
+  bool aniso2d = is_aniso2d();
+  for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
+    if (allneeded[c])
+      _require_component(component(c), aniso2d);
+}
 
-  // check if we are in 2d but anisotropy couples xy with z
+// check if we are in 2d but anisotropy couples xy with z
+bool fields::is_aniso2d() {
   bool aniso2d = false;
   if (gv.dim == D2) {
     int i;
@@ -494,9 +509,18 @@ void fields::require_component(component c) {
     aniso2d = or_to_all(i < num_chunks);
     finished_working();
   }
+  else if (beta != 0)
+    abort("Nonzero beta unsupported in dimensions other than 2.");
   if (aniso2d && beta != 0 && is_real)
     abort("Nonzero beta need complex fields when mu/epsilon couple TE and TM");
-  aniso2d = aniso2d || (beta != 0); // beta couples TE/TM
+  return aniso2d || (beta != 0); // beta couples TE/TM
+}
+
+void fields::_require_component(component c, bool aniso2d) {
+  if (!gv.has_field(c))
+    abort("cannot require a %s component in a %s grid", component_name(c), dimension_name(gv.dim));
+
+  components_allocated = true;
 
   // allocate fields if they haven't been allocated yet for this component
   int need_to_reconnect = 0;
@@ -506,10 +530,10 @@ void fields::require_component(component c) {
         if (chunks[i]->alloc_f(c_alloc)) need_to_reconnect++;
   }
 
-  if (need_to_reconnect) figure_out_step_plan();
-  am_now_working_on(MpiAllTime);
-  if (sum_to_all(need_to_reconnect)) chunk_connections_valid = false;
-  finished_working();
+  if (need_to_reconnect) {
+    figure_out_step_plan();
+    chunk_connections_valid = false; // connect_chunks() will synchronize this for us
+  }
 }
 
 void fields_chunk::remove_sources() {

src/meep.hpp

@@ -1728,7 +1728,10 @@ class fields {
   void add_volume_source(component c, const src_time &src, const volume &,
                          std::complex<double> amp = 1.0);
   void add_volume_source(const src_time &src, const volume &, gaussianbeam beam);
-  void require_component(component c);
+  bool is_aniso2d();
+  void require_source_components();
+  void _require_component(component c, bool aniso2d);
+  void require_component(component c) { _require_component(c, is_aniso2d()); }
   void add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex<double>* amp_arr);
 
   // mpb.cpp

src/sources.cpp

@@ -319,7 +319,10 @@ void fields::add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, st
   fields_chunk *fc = chunks[cur_data.fc_idx];
   if (!fc->is_mine()) abort("wrong fields chunk");
   fc->sources[ft] = tmp->add_to(fc->sources[ft]);
-  require_component(c);
+  // 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()
+  //     after all add_srcdata calls are complete.
 }
 
 static realnum *amp_func_data_re = NULL;