Merge branch 'develop' into collisionFilter

include/openmc/capi.h

@@ -28,7 +28,8 @@ extern "C" {
   int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end);
   int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end);
   int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end);
-  int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end);
+  int openmc_extend_meshes(int32_t n, const char* type, int32_t* index_start,
+                           int32_t* index_end);
   int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end);
   int openmc_filter_get_id(int32_t index, int32_t* id);
   int openmc_filter_get_type(int32_t index, char* type);
@@ -73,11 +74,7 @@ extern "C" {
   int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
   int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
   int openmc_mesh_get_id(int32_t index, int32_t* id);
-  int openmc_mesh_get_dimension(int32_t index, int** id, int* n);
-  int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
   int openmc_mesh_set_id(int32_t index, int32_t id);
-  int openmc_mesh_set_dimension(int32_t index, int n, const int* dims);
-  int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width);
   int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
   int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
   int openmc_new_filter(const char* type, int32_t* index);
@@ -86,6 +83,15 @@ extern "C" {
   int openmc_plot_geometry();
   int openmc_id_map(const void* slice, int32_t* data_out);
   int openmc_property_map(const void* slice, double* data_out);
+  int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx,
+                       double** grid_y, int* ny, double** grid_z, int* nz);
+  int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x,
+                       const int nx, const double* grid_y, const int ny,
+                       const double* grid_z, const int nz);
+  int openmc_regular_mesh_get_dimension(int32_t index, int** id, int* n);
+  int openmc_regular_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
+  int openmc_regular_mesh_set_dimension(int32_t index, int n, const int* dims);
+  int openmc_regular_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width);
   int openmc_reset();
   int openmc_reset_timers();
   int openmc_run();

include/openmc/error.h

@@ -10,7 +10,7 @@
 #include "openmc/capi.h"
 #include "openmc/settings.h"
 
-#ifdef __GNUC__
+#if defined(__GNUC__) || defined(__clang__)
 #define UNREACHABLE() __builtin_unreachable()
 #else
 #define UNREACHABLE() (void)0

include/openmc/geometry.h

@@ -43,14 +43,11 @@ bool check_cell_overlap(Particle& p, bool error=true);
 //!
 //! \param p A particle to be located.  This function will populate the
 //!   geometry-dependent data fields of the particle.
-//! \param use_neighbor_lists If true, neighbor lists will be used to accelerate
-//!   the geometry search, but this only works if the cell attribute of the
-//!   particle's lowest coordinate level is valid and meaningful.
 //! \return True if the particle's location could be found and ascribed to a
 //!   valid geometry coordinate stack.
 //==============================================================================
-
-bool find_cell(Particle& p, bool use_neighbor_lists);
+bool exhaustive_find_cell(Particle& p);
+bool neighbor_list_find_cell(Particle& p); // Only usable on surface crossings
 
 //==============================================================================
 //! Move a particle into a new lattice tile.

include/openmc/mesh.h

@@ -114,6 +114,14 @@ public:
   void bins_crossed(const Particle& p, std::vector<int>& bins,
                     std::vector<double>& lengths) const override;
 
+  //! Count number of bank sites in each mesh bin / energy bin
+  //
+  //! \param[in] Pointer to bank sites
+  //! \param[in] Number of bank sites
+  //! \param[out] Whether any bank sites are outside the mesh
+  xt::xtensor<double, 1> count_sites(const Particle::Bank* bank,
+                                     int64_t length, bool* outside) const;
+
   //! Get bin given mesh indices
   //
   //! \param[in] Array of mesh indices
@@ -200,16 +208,6 @@ public:
 
   void to_hdf5(hid_t group) const override;
 
-  // New methods
-
-  //! Count number of bank sites in each mesh bin / energy bin
-  //
-  //! \param[in] bank Array of bank sites
-  //! \param[out] Whether any bank sites are outside the mesh
-  //! \return Array indicating number of sites in each mesh/energy bin
-  xt::xtensor<double, 1> count_sites(const Particle::Bank* bank, int64_t length,
-    bool* outside) const;
-
   // Data members
 
   double volume_frac_; //!< Volume fraction of each mesh element
@@ -221,6 +219,7 @@ class RectilinearMesh : public StructuredMesh
 {
 public:
   // Constructors
+  RectilinearMesh() = default;
   RectilinearMesh(pugi::xml_node node);
 
   // Overriden methods
@@ -239,8 +238,9 @@ public:
 
   void to_hdf5(hid_t group) const override;
 
-private:
   std::vector<std::vector<double>> grid_;
+
+  int set_grid();
 };
 
 #ifdef DAGMC
@@ -424,8 +424,6 @@ void read_meshes(pugi::xml_node root);
 //! \param[in] group HDF5 group
 void meshes_to_hdf5(hid_t group);
 
-RegularMesh* get_regular_mesh(int32_t index);
-
 void free_memory_mesh();
 
 } // namespace openmc

include/openmc/plot.h

@@ -182,7 +182,7 @@ T PlotBase::get_map() const {
         p.r()[in_i] = xyz[in_i] + in_pixel * x;
         p.n_coord_ = 1;
         // local variables
-        bool found_cell = find_cell(p, 0);
+        bool found_cell = exhaustive_find_cell(p);
         j = p.n_coord_ - 1;
         if (level >= 0) { j = level; }
         if (found_cell) {

include/openmc/position.h

@@ -51,6 +51,18 @@ struct Position {
     }
   }
 
+  // Access to x, y, or z by compile time known index (specializations below)
+  template<int i>
+  const double& get() const
+  {
+    throw std::out_of_range {"Index in Position must be between 0 and 2."};
+  }
+  template<int i>
+  double& get()
+  {
+    throw std::out_of_range {"Index in Position must be between 0 and 2."};
+  }
+
   // Other member functions
 
   //! Dot product of two vectors
@@ -77,6 +89,38 @@ struct Position {
   double z = 0.;
 };
 
+// Compile-time known member index access functions
+template<>
+inline const double& Position::get<0>() const
+{
+  return x;
+}
+template<>
+inline const double& Position::get<1>() const
+{
+  return y;
+}
+template<>
+inline const double& Position::get<2>() const
+{
+  return z;
+}
+template<>
+inline double& Position::get<0>()
+{
+  return x;
+}
+template<>
+inline double& Position::get<1>()
+{
+  return y;
+}
+template<>
+inline double& Position::get<2>()
+{
+  return z;
+}
+
 // Binary operators
 inline Position operator+(Position a, Position b) { return a += b; }
 inline Position operator+(Position a, double b)   { return a += b; }

include/openmc/settings.h

@@ -88,7 +88,7 @@ extern RunMode run_mode;                 //!< Run mode (eigenvalue, fixed src, e
 extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
 extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
 extern std::unordered_set<int> source_write_surf_id; //!< Surface ids where sources will be written
-extern int64_t max_particles;            //!< maximum number of particles to be banked on surfaces per process
+extern int64_t max_surface_particles;    //!< maximum number of particles to be banked on surfaces per process
 extern TemperatureMethod temperature_method;           //!< method for choosing temperatures
 extern double temperature_tolerance;     //!< Tolerance in [K] on choosing temperatures
 extern double temperature_default;       //!< Default T in [K]

include/openmc/shared_array.h

@@ -73,12 +73,12 @@ public:
   {
     // Atomically capture the index we want to write to
     int64_t idx;
-    #pragma omp atomic capture
+    #pragma omp atomic capture seq_cst
     idx = size_++;
 
     // Check that we haven't written off the end of the array
     if (idx >= capacity_) {
-      #pragma omp atomic write
+      #pragma omp atomic write seq_cst
       size_ = capacity_;
       return -1;
     }