fix edge case bug with small negative coefficients in tmat

INSIGNEO · Jun 25, 2019 · 4e30b51 · 4e30b51
1 parent 69d8c7f
commit 4e30b51
Show file tree

Hide file tree

Showing 5 changed files with 39 additions and 16 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -38,6 +38,7 @@ set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -Wall -O3 -march=native
 option(DEBUG_VERBOSE "Verbose debug outputs")
 if(DEBUG_VERBOSE)
   add_definitions(-DDEBUG_VERBOSE)
+  add_definitions(-DTMATRIX_DEBUG)
   message(STATUS "Enabled verbose debugging messages")
 endif(DEBUG_VERBOSE)
 

diff --git a/src/fd_routines.cpp b/src/fd_routines.cpp
@@ -84,21 +84,43 @@ void gradient_existing(const DM& dmda, const Vec& srcvec, Vec& tgtvec, integer d
   j_hi += j_lo;
   k_hi += k_lo;
 
-  intvector ofs = {0, 0, 0};
-  ofs[dim] = 1;
-  for (integer i = i_lo; i < i_hi; i++)
+#ifdef GRADIENT_DEBUG
+  MPI_Comm comm;
+  int rank, commsize;
+  PetscObjectGetComm(reinterpret_cast<PetscObject>(dmda), &comm);
+  MPI_Comm_size(comm, &commsize);
+  MPI_Comm_rank(comm, &rank);
+  for (int irank = 0; irank < commsize; irank++)
   {
-    for (integer j = j_lo; j < j_hi; j++)
+    if (irank == rank)
     {
-      for (integer k = k_lo; k < k_hi; k++)
+#endif // GRADIENT_DEBUG
+      intcoord ofs = {0, 0, 0};
+      ofs[dim] = 1;
+      for (integer i = i_lo; i < i_hi; i++)
       {
-        // Remember c-indexing is backwards because PETSc is odd
-        grad_array[k][j][i] = 0.5
-                              * (img_array[k + ofs[2]][j + ofs[1]][i + ofs[0]]
-                                 - img_array[k - ofs[2]][j - ofs[1]][i - ofs[0]]);
+        for (integer j = j_lo; j < j_hi; j++)
+        {
+          for (integer k = k_lo; k < k_hi; k++)
+          {
+#ifdef GRADIENT_DEBUG
+            // Remember c-indexing is backwards because PETSc is odd
+            std::cout << "(" << i << ", " << j << ", " << k << ") +/- " << ofs << ": "
+                      << img_array[k + ofs[2]][j + ofs[1]][i + ofs[0]] << " - "
+                      << img_array[k - ofs[2]][j - ofs[1]][i - ofs[0]] << std::endl;
+#endif // GRADIENT_DEBUG
+            grad_array[k][j][i] = 0.5
+                                  * (img_array[k + ofs[2]][j + ofs[1]][i + ofs[0]]
+                                     - img_array[k - ofs[2]][j - ofs[1]][i - ofs[0]]);
+          }
+        }
       }
+#ifdef GRADIENT_DEBUG
     }
+    MPI_Barrier(comm);
   }
+#endif // GRADIENT_DEBUG
+
   // Release pointers and allow petsc to cleanup
   perr = DMDAVecRestoreArray(dmda, srcvec, &img_array);
   CHKERRXX(perr);

diff --git a/src/mapbase.cpp b/src/mapbase.cpp
@@ -19,7 +19,7 @@ MapBase::MapBase(const intcoord& shape, const integer& ndim, const intcoord& nod
   calculate_node_locations();
 }
 
-floatcoord MapBase::coord_from_index(intcoord index) const
+floatcoord MapBase::coord_from_index(const intcoord& index) const
 {
   return {_node_locations[0][index[0]], _node_locations[1][index[1]], _node_locations[2][index[2]]};
 }
@@ -137,9 +137,9 @@ std::pair<intcoord, intcoord> MapBase::get_pixel_neighbourhood(const intcoord& m
   intcoord cimhi = {0, 0, 1};
   for (integer idim = 0; idim < this->ndim(); idim++)
   {
-    integer ctr = this->node_locs()[idim][map_node[idim]];
-    cimlo[idim] = ctr - this->spacing()[idim];
-    cimhi[idim] = ctr + this->spacing()[idim];
+    floating ctr = this->node_locs()[idim][map_node[idim]];
+    cimlo[idim] = std::ceil(ctr - this->spacing()[idim]);
+    cimhi[idim] = std::floor(ctr + this->spacing()[idim]);
   }
   clamp_location_lo(cimlo, this->mask().index_min());
   clamp_location_hi(cimhi, this->mask().shape());

diff --git a/src/mapbase.hpp b/src/mapbase.hpp
@@ -47,7 +47,7 @@ class MapBase: public GridVariable {
 
   std::pair<intcoord, intcoord> get_pixel_neighbourhood(const intcoord& map_node) const;
 
-  floatcoord coord_from_index(intcoord index) const;
+  floatcoord coord_from_index(const intcoord& index) const;
 
   //! Return the node at the lower corner of the cell containing loc
   template<typename mathtype>

diff --git a/src/tmatrix.cpp b/src/tmatrix.cpp
@@ -291,15 +291,15 @@ std::pair<floating, floating> calculate_entries(integer row, integer col, intege
         auto& rank_cc = row_coeffs[jrank];
         for (size_t iidx = 0; iidx < rank_reqs.size(); iidx++)
         {
-          std::cout << rank_reqs[iidx] << " [" << idof << "]"
+          std::cout << "post " << rank_reqs[iidx] << " [" << idof << "]"
                     << ", rc = " << rank_rc[iidx] << ", cc = " << rank_cc[iidx]
                     << ", gd = " << rank_vals[iidx].first << ", dd = " << rank_vals[iidx].second
                     << std::endl;
 
         }
       }
       std::cout << std::flush;
-      std::this_thread::sleep_for(std::chrono::milliseconds(10));
+      std::this_thread::sleep_for(std::chrono::milliseconds(2));
     }
     MPI_Barrier(comm);
   }