refactor sync_chunk_connections to remove barrier during fields::step (…

…#1635)

src/boundaries.cpp

@@ -30,6 +30,8 @@ namespace meep {
 void fields::set_boundary(boundary_side b, direction d, boundary_condition cond) {
   if (boundaries[b][d] != cond) {
     boundaries[b][d] = cond;
+    // we don't need to call sync_chunk_connections() since set_boundary()
+    // should always be called on every process
     chunk_connections_valid = false;
   }
 }
@@ -120,11 +122,19 @@ void fields::disconnect_chunks() {
   }
 }
 
-void fields::connect_chunks() {
-  /* make sure all processes agree on chunk_connections_valid to avoid deadlocks */
+// this should be called by any code that might set chunk_connections_valid = false,
+// with the caveat that we need to be careful that we call it on all processes
+void fields::sync_chunk_connections() {
+  /* make sure all processes agree on chunk_connections_valid to avoid deadlocks
+     when we eventually call connect_chunks */
   am_now_working_on(MpiAllTime);
   chunk_connections_valid = and_to_all(chunk_connections_valid);
   finished_working();
+}
+
+void fields::connect_chunks() {
+  // might have invalidated connections in step_db, update_eh, or update_pols:
+  if (changed_materials) sync_chunk_connections();
 
   if (!chunk_connections_valid) {
     am_now_working_on(Connecting);

src/fields.cpp

@@ -74,6 +74,7 @@ fields::fields(structure *s, double m, double beta, bool zero_fields_near_cylori
         boundaries[b][d] = None;
     }
   chunk_connections_valid = false;
+  changed_materials = true;
 
   // unit directions are periodic by default:
   FOR_DIRECTIONS(d) {
@@ -126,6 +127,7 @@ fields::fields(const fields &thef)
   for (int b = 0; b < 2; b++)
     FOR_DIRECTIONS(d) { boundaries[b][d] = thef.boundaries[b][d]; }
   chunk_connections_valid = false;
+  changed_materials = true;
 }
 
 fields::~fields() {
@@ -152,6 +154,9 @@ void fields::use_real_fields() {
   is_real = 1;
   for (int i = 0; i < num_chunks; i++)
     chunks[i]->use_real_fields();
+
+  // don't need to call sync_chunk_connections() since use_real_fields()
+  // should always be called on every process
   chunk_connections_valid = false;
 }
 
@@ -492,6 +497,7 @@ void fields::require_source_components() {
   for (int c = 0; c < NUM_FIELD_COMPONENTS; ++c)
     if (allneeded[c])
       _require_component(component(c), aniso2d);
+  sync_chunk_connections();
 }
 
 // check if we are in 2d but anisotropy couples xy with z
@@ -532,7 +538,9 @@ void fields::_require_component(component c, bool aniso2d) {
 
   if (need_to_reconnect) {
     figure_out_step_plan();
-    chunk_connections_valid = false; // connect_chunks() will synchronize this for us
+    // we will eventually call sync_chunk_connections(), in either require_component(c)
+    // or require_components(), to synchronize this across processes:
+    chunk_connections_valid = false;
   }
 }
 
@@ -566,6 +574,7 @@ void fields_chunk::remove_susceptibilities(bool shared_chunks) {
 }
 
 void fields::remove_susceptibilities() {
+  changed_materials = true;
   for (int i = 0; i < num_chunks; i++)
     chunks[i]->remove_susceptibilities(shared_chunks);
 }
@@ -635,6 +644,7 @@ int fields::phase_in_material(const structure *snew, double time) {
   for (int i = 0; i < num_chunks; i++)
     if (chunks[i]->is_mine()) chunks[i]->phase_in_material(snew->chunks[i]);
   phasein_time = (int)(time / dt);
+  changed_materials = true;
   // FIXME: how to handle changes in susceptibilities?
   return phasein_time;
 }

src/meep.hpp

@@ -1742,7 +1742,7 @@ class fields {
   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 require_component(component c) { _require_component(c, is_aniso2d()); sync_chunk_connections(); }
   void add_srcdata(struct sourcedata cur_data, src_time *src, size_t n, std::complex<double>* amp_arr);
 
   // mpb.cpp
@@ -2039,9 +2039,11 @@ class fields {
   void figure_out_step_plan();
   // boundaries.cpp
   bool chunk_connections_valid;
+  bool changed_materials; // keep track of whether materials have changed (in case field chunk connections need sync'ing)
   void find_metals();
   void disconnect_chunks();
   void connect_chunks();
+  void sync_chunk_connections();
   void connect_the_chunks(); // Intended to be ultra-private...
   bool on_metal_boundary(const ivec &);
   ivec ilattice_vector(direction) const;

src/step.cpp

@@ -94,6 +94,8 @@ void fields::step() {
     synchronized_magnetic_fields = save_synchronized_magnetic_fields;
   }
 
+  changed_materials = false; // any material changes were handled in connect_chunks()
+
   if (!std::isfinite(get_field(D_EnergyDensity, gv.center(), false)))
     meep::abort("simulation fields are NaN or Inf");
 }

src/step_db.cpp

@@ -19,6 +19,7 @@
 #include <stdlib.h>
 #include <math.h>
 #include <string.h>
+#include <assert.h>
 
 #include "meep.hpp"
 #include "meep_internals.hpp"
@@ -32,7 +33,10 @@ namespace meep {
 void fields::step_db(field_type ft) {
   for (int i = 0; i < num_chunks; i++)
     if (chunks[i]->is_mine())
-      if (chunks[i]->step_db(ft)) chunk_connections_valid = false;
+      if (chunks[i]->step_db(ft)) {
+       chunk_connections_valid = false;
+       assert(changed_materials);
+      }
 }
 
 bool fields_chunk::step_db(field_type ft) {

src/update_eh.cpp

@@ -16,6 +16,7 @@
 */
 
 #include <string.h>
+#include <assert.h>
 
 #include "meep.hpp"
 #include "meep_internals.hpp"
@@ -28,8 +29,10 @@ void fields::update_eh(field_type ft, bool skip_w_components) {
   if (ft != E_stuff && ft != H_stuff) meep::abort("update_eh only works with E/H");
   for (int i = 0; i < num_chunks; i++)
     if (chunks[i]->is_mine())
-      if (chunks[i]->update_eh(ft, skip_w_components))
+      if (chunks[i]->update_eh(ft, skip_w_components)) {
         chunk_connections_valid = false; // E/H allocated - reconnect chunks
+        assert(changed_materials);
+      }
 }
 
 bool fields_chunk::needs_W_prev(component c) const {

src/update_pols.cpp

@@ -18,6 +18,7 @@
 #include <stdio.h>
 #include <math.h>
 #include <string.h>
+#include <assert.h>
 
 #include "meep.hpp"
 #include "meep_internals.hpp"
@@ -30,7 +31,10 @@ namespace meep {
 void fields::update_pols(field_type ft) {
   for (int i = 0; i < num_chunks; i++)
     if (chunks[i]->is_mine())
-      if (chunks[i]->update_pols(ft)) chunk_connections_valid = false;
+      if (chunks[i]->update_pols(ft)) {
+        chunk_connections_valid = false;
+        assert(changed_materials);
+      }
 }
 
 bool fields_chunk::update_pols(field_type ft) {