Fix bug for 3D simulation using 1D cell and single-precision floating point (#3157)

* fix bug for 3D simulation using 1D cell and single-precision floating point

* approximately compare phase to +/- 1

* Apply suggestion from @stevengj

* Update boundaries.cpp: rm redundant thephase

---------

Co-authored-by: Steven G. Johnson <stevenj@alum.mit.edu>
Co-authored-by: Steven G. Johnson <stevenj@mit.edu>

Author: oskooi

src/boundaries.cpp

@@ -350,6 +350,21 @@ bool fields_chunk::needs_W_notowned(component c) {
   return false;
 }
 
+/* approximate comparisons of phases to 1 or -1.   exp(ix) * exp(-ix) is not exactly 1
+   due to roundoff errors, but we want to treat it equal for the purpose of determining
+   the boundary conditions, because we have optimized b.c.'s for phase +1 and -1, and
+   roundoff errors on the order of 1e-13 are not significant relative to the truncation
+   errors in meep anyway.  (the alternative would be to accumulate symmetry phases more
+   accurately by adding exponents rather than multiplying). */
+static bool phase_isclose(std::complex<double> thephase, double realphase) {
+  return fabs(thephase.imag()) < 1e-13 && fabs(thephase.real() - realphase) < 1e-13;
+}
+static connect_phase connect_phase_from_phase(std::complex<double> thephase) {
+  return phase_isclose(thephase, 1.0)    ? CONNECT_COPY
+         : phase_isclose(thephase, -1.0) ? CONNECT_NEGATE
+                                         : CONNECT_PHASE;
+}
+
 void fields::connect_the_chunks() {
   /* For some of the chunks, H==B, and we definitely don't need to
      send B between two such chunks.   We'll still send B when
@@ -405,9 +420,7 @@ void fields::connect_the_chunks() {
               if ((chunks[i]->is_mine() || chunks[j]->is_mine()) && chunks[j]->gv.owns(here) &&
                   !(is_B(corig) && is_B(c) && B_redundant[5 * i + corig - Bx] &&
                     B_redundant[5 * j + c - Bx])) {
-                const connect_phase ip = thephase == 1.0
-                                             ? CONNECT_COPY
-                                             : (thephase == -1.0 ? CONNECT_NEGATE : CONNECT_PHASE);
+                const connect_phase ip = connect_phase_from_phase(thephase);
                 comm_sizes[{type(c), ip, pair_j_to_i}] += num_reals_per_voxel;
 
                 if (needs_W_notowned[corig]) {
@@ -468,9 +481,8 @@ void fields::connect_the_chunks() {
           IVEC_LOOP_ILOC(vi, here);
           component c = corig;
           // We're looking at a border element...
-          std::complex<double> thephase_double;
-          if (locate_component_point(&c, &here, &thephase_double) && !on_metal_boundary(here)) {
-            std::complex<realnum> thephase(thephase_double.real(), thephase_double.imag());
+          std::complex<double> thephase;
+          if (locate_component_point(&c, &here, &thephase) && !on_metal_boundary(here)) {
             for (int j = 0; j < num_chunks; j++) {
               const std::pair<int, int> pair_j_to_i{j, i};
               const bool i_is_mine = chunks[i]->is_mine();
@@ -494,10 +506,7 @@ void fields::connect_the_chunks() {
               if (chunks[j]->gv.owns(here) &&
                   !(is_B(corig) && is_B(c) && B_redundant[5 * i + corig - Bx] &&
                     B_redundant[5 * j + c - Bx])) {
-                const connect_phase ip =
-                    thephase == static_cast<realnum>(1.0)
-                        ? CONNECT_COPY
-                        : (thephase == static_cast<realnum>(-1.0) ? CONNECT_NEGATE : CONNECT_PHASE);
+                const connect_phase ip = connect_phase_from_phase(thephase);
                 const ptrdiff_t m = chunks[j]->gv.index(c, here);
 
                 {

tests/three_d.cpp

@@ -223,6 +223,32 @@ int test_pml_splitting(double eps(const vec &), int splitting) {
   return 1;
 }
 
+// Test that a 3D simulation with a 1D cell (zero size in X and Y)
+// and Bloch-periodic boundaries does not crash. This exercises the
+// boundary connection logic where the connect_phase classification
+// must be consistent between the counting pass and the connection
+// setup pass, regardless of the realnum precision. (See issue #2916.)
+int test_1d_cell_bloch() {
+  double a = 52.0;
+  double pml_um = 1.0;
+  double air_um = 10.0;
+  double size_z = pml_um + air_um + pml_um;
+
+  grid_volume gv = vol3d(0, 0, size_z, a);
+  gv.center_origin();
+  structure s(gv, one, pml(pml_um));
+
+  master_printf("Testing 1D cell with Bloch boundaries...\n");
+  fields f(&s);
+  f.use_bloch(vec(-0.2, -0.2, 0.0));
+  f.add_point_source(Ex, 1.0, 0.1, 0.0, 4.0, vec(0, 0, 0));
+
+  while (f.time() < 5.0)
+    f.step();
+
+  return 1;
+}
+
 int main(int argc, char **argv) {
   initialize mpi(argc, argv);
   verbosity = 0;
@@ -239,5 +265,7 @@ int main(int argc, char **argv) {
   for (int s = 2; s < 4; s++)
     if (!test_pml_splitting(one, s)) meep::abort("error in test_pml_splitting vacuum\n");
 
+  if (!test_1d_cell_bloch()) meep::abort("error in test_1d_cell_bloch\n");
+
   return 0;
 }