Merge branch 'issue-53'

Signed-off-by: Rodrigo Tobar <[email protected]>

src/CMakeLists.txt

@@ -14,6 +14,7 @@ set(VR_SOURCES
     bgfield.cxx
     buildandsortarrays.cxx
     endianutils.cxx
+    exceptions.cxx
     fofalgo.cxx
     gadgetio.cxx
     "${git_revision_cxx}"

src/exceptions.cxx

@@ -0,0 +1,45 @@
+#include <sstream>
+
+#ifdef USEOPENMP
+#include <omp.h>
+#endif // USEOPENMP
+
+#include "exceptions.h"
+
+#include "allvars.h"
+
+namespace vr
+{
+
+static std::string build_error_message(const Particle &part, const std::string &function)
+{
+    std::ostringstream os;
+    os << "Particle density not positive, cannot continue.\n\n"
+       << "Particle information: "
+       << "id=" << part.GetID() << ", pid=" << part.GetPID()
+       << ", type=" << part.GetType()
+       << ", pos=(" << part.X() << ", " << part.Y() << ", " << part.Z()
+       << "), vel=(" << part.Vx() << ", " << part.Vy() << ", " << part.Vz()
+       << "), mass=" << part.GetMass()
+       << ", density=" << part.GetDensity() << '\n';
+    os << "Execution context: MPI enabled=";
+#ifdef USEMPI
+    os << "yes, rank=" << ThisTask << '/' << NProcs;
+#else
+    os << "no";
+#endif // USEMPI
+    os << ", OpenMP enabled=";
+#ifdef USEOPENMP
+    os << "yes, thread=" << omp_get_thread_num() << '/' << omp_get_num_threads();
+#else
+    os << "no";
+#endif // USEOPENMP
+    return os.str();
+}
+
+non_positive_density::non_positive_density(const Particle &part, const std::string &function)
+  : exception(build_error_message(part, function))
+{
+}
+
+}  // namespace vr

src/exceptions.h

@@ -0,0 +1,37 @@
+#include <exception>
+#include <string>
+
+// Forward declaration
+namespace NBody
+{
+class Particle;
+}
+
+namespace vr
+{
+
+/// Base class for all exceptions in velociraptor
+class exception : public std::exception
+{
+public:
+    exception(std::string what)
+      : m_what(std::move(what))
+    {}
+
+    const char *what() const noexcept override
+    {
+        return m_what.c_str();
+    }
+
+private:
+    std::string m_what;
+};
+
+/// Exception thrown when a particle with non-positive density is found
+class non_positive_density : public exception
+{
+public:
+    non_positive_density(const NBody::Particle &part, const std::string &function);
+};
+
+} // namespace vr

src/localbgcomp.cxx

@@ -3,6 +3,7 @@
     \todo once ratio is calculated, must figure out best way to do global mpi fit. Probably best to combine the binned distribution and fit that
  */
 
+#include "exceptions.h"
 #include "stf.h"
 
 /*! This calculates the logarithmic ratio of the measured velocity density and the expected velocity density assuming a bg muiltivariate gaussian distribution
@@ -87,29 +88,7 @@ if (nbodies > ompsubsearchnum)
     for (i=0;i<nbodies;i++)
     {
         if (!(Part[i].GetDensity() > 0)) {
-            const auto &part = Part[i];
-            std::ostringstream os;
-            os << "Particle density not positive, cannot continue.\n\n"
-               << " Information for particle " << i << '/' << nbodies
-               << ": id=" << part.GetID() << ", pid=" << part.GetPID()
-               << ", type=" << part.GetType()
-               << ", pos=(" << part.X() << ", " << part.Y() << ", " << part.Z()
-               << "), vel=(" << part.Vx() << ", " << part.Vy() << ", " << part.Vz()
-               << "), mass=" << part.GetMass()
-               << ", density=" << part.GetDensity() << '\n';
-            os << "Information for execution context: MPI enabled=";
-#ifdef USEMPI
-            os << "yes, rank=" << ThisTask << '/' << NProcs;
-#else
-            os << "no";
-#endif // USEMPI
-            os << ", OpenMP enabled=";
-#ifdef USEOPENMP
-            os << "yes, thread=" << omp_get_thread_num() << '/' << nthreads;
-#else
-            os << "no";
-#endif // USEOPENMP
-            throw std::runtime_error(os.str());
+            throw vr::non_positive_density(Part[i], __PRETTY_FUNCTION__);
         }
         tempdenv=Part[i].GetDensity()/opt.Nsearch;
 #ifdef USEOPENMP

src/localfield.cxx

@@ -4,6 +4,7 @@
 
 //--  Local Velocity density routines
 
+#include "exceptions.h"
 #include "logging.h"
 #include "stf.h"
 #include "swiftinterface.h"
@@ -952,7 +953,6 @@ reduction(+:nprocessed,ntot)
 
     MEMORY_USAGE_REPORT(debug, opt);
 
-    if (nimport>0) {
 #ifdef USEOPENMP
 #pragma omp parallel default(shared) \
 private(id,v2,nnids,nnr2,nnidsneighbours,nnr2neighbours,weight,pqx,pqv,Pval,pid2)
@@ -992,12 +992,14 @@ reduction(+:nprocessed)
                 pqx->Push(nnids[j], nnr2[j]);
             }
         }
-        //search neighbouring domain and update priority queue
-        treeneighbours->FindNearestPos(leafnodes[i].cm,nnidsneighbours,nnr2neighbours,nimportsearch);
-        for (auto j = 0; j < nimportsearch; j++) {
-            if (nnr2neighbours[j] < pqx->TopPriority()){
-                pqx->Pop();
-                pqx->Push(nnidsneighbours[j]+nbodies, nnr2neighbours[j]);
+        //search neighbouring domain and update priority queue if any neighbours imported
+        if (nimport > 0) {
+            treeneighbours->FindNearestPos(leafnodes[i].cm,nnidsneighbours,nnr2neighbours,nimportsearch);
+            for (auto j = 0; j < nimportsearch; j++) {
+                if (nnr2neighbours[j] < pqx->TopPriority()){
+                    pqx->Pop();
+                    pqx->Push(nnidsneighbours[j]+nbodies, nnr2neighbours[j]);
+                }
             }
         }
         for (auto j = 0; j < opt.Nsearch; j++) {
@@ -1044,7 +1046,6 @@ reduction(+:nprocessed)
 }
 #endif
         delete treeneighbours;
-    }
     delete[] PartDataIn;
     delete[] PartDataGet;
     delete[] NNDataIn;
@@ -1057,4 +1058,16 @@ reduction(+:nprocessed)
     //free memory
     if (itreeflag) delete tree;
     if (period!=NULL) delete[] period;
+
+    // Double-check that valid densities have been set in all particles
+    for (int i = 0; i < nbodies; i++)
+    {
+        const auto &part = Part[i];
+#ifdef STRUCDEN
+        if (part.GetType()<=0) continue;
+#endif
+        if (!(part.GetDensity() > 0)) {
+            throw vr::non_positive_density(Part[i], __PRETTY_FUNCTION__);
+        }
+    }
 }