Merge pull request #2 from paulromano/dlopen_source

Fixes to custom source feature

docs/source/io_formats/settings.rst

@@ -142,7 +142,7 @@ materials in the problem and is specified using a :ref:`mesh_element`.
 ----------------------------
 
 Determines whether to use event-based parallelism instead of the default
-history-based parallelism. 
+history-based parallelism.
 
   *Default*: false
 
@@ -467,7 +467,7 @@ attributes/sub-elements:
     more documentation on how to build sources can be found in :ref:`custom_source`
 
     *Deafult*: None
-     
+
   :space:
     An element specifying the spatial distribution of source sites. This element
     has the following attributes:
@@ -605,59 +605,60 @@ attributes/sub-elements:
 Custom Sources
 ++++++++++++++
 
-It is often the case that one may wish to simulate a complex source distribution,
-which may include physics not present within OpenMC or to be phase space complex. It
-is possible to define a complex source with an externally defined source function
-that is loaded at runtime. A simple example source is shown below. 
+It is often the case that one may wish to simulate a complex source
+distribution, which may include physics not present within OpenMC or to be phase
+space complex. It is possible to define a complex source with an externally
+defined source function that is loaded at runtime. A simple example source is
+shown below.
 
 .. code-block:: c++
 
-   #include <iostream>
    #include "openmc/random_lcg.h"
    #include "openmc/source.h"
    #include "openmc/particle.h"
 
    // you must have external C linkage here
-   extern "C" openmc::Particle::Bank sample_source(const int64_t seed) {
+   extern "C" openmc::Particle::Bank sample_source(uint64_t* seed) {
      openmc::Particle::Bank particle;
-      // wgt
-      particle.particle = openmc::Particle::Type::neutron;
-      particle.wgt = 1.0;
-      // position 
-      double angle = 2.*M_PI*openmc::prn();
-      double radius = 3.0;
-      particle.r.x = radius*std::cos(angle);
-      particle.r.y = radius*std::sin(angle);
-      particle.r.z = 0.;
-      // angle
-      particle.u = {1.,0,0};
-      particle.E = 14.08e6;
-      particle.delayed_group = 0;
-      return particle;    
+     // weight
+     particle.particle = openmc::Particle::Type::neutron;
+     particle.wgt = 1.0;
+     // position
+     double angle = 2.0 * M_PI * openmc::prn(seed);
+     double radius = 3.0;
+     particle.r.x = radius * std::cos(angle);
+     particle.r.y = radius * std::sin(angle);
+     particle.r.z = 0.0;
+     // angle
+     particle.u = {1.0, 0.0, 0.0};
+     particle.E = 14.08e6;
+     particle.delayed_group = 0;
+     return particle;
   }
 
 The above source, creates 14.08 MeV neutrons, with an istropic direction
 vector but distributed in a ring with a 3 cm radius. This routine is
-not particular complex, but should serve as an example upon which to build
-more complicated sources. 
+not particularly complex, but should serve as an example upon which to build
+more complicated sources.
+
+  .. note:: The function signature must be declared to be extern "C".
 
-  .. note:: The function signature must be declared to be extern.
-	    
   .. note:: You should only use the openmc::prn() random number generator
 
-In order to build your external source you need the following CMakeLists.txt file
+In order to build your external source you need the following CMakeLists.txt
+file
 
 .. code-block:: cmake
 
    cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
-   project(openmc_sources CXX) 
+   project(openmc_sources CXX)
    add_library(source SHARED source_ring.cpp)
    find_package(OpenMC REQUIRED HINTS <path to openmc>)
    target_link_libraries(source OpenMC::libopenmc)
 
-You will now have a libsouce.so (or .dylib) file in this directory, now point the library
-attribute of the source XML element to this file and you will be able to sample 
-particles. 
+You will now have a libsource.so (or .dylib) file in this directory, now point
+the library attribute of the source XML element to this file and you will be
+able to sample particles.
 
 .. _univariate:
 

examples/xml/custom_source/CMakeLists.txt

@@ -0,0 +1,8 @@
+cmake_minimum_required(VERSION 3.3 FATAL_ERROR)
+project(openmc_sources CXX)
+add_library(source SHARED source_ring.cpp)
+find_package(OpenMC REQUIRED)
+if (OpenMC_FOUND)
+  message(STATUS "Found OpenMC: ${OpenMC_DIR}")
+endif()
+target_link_libraries(source OpenMC::libopenmc)

examples/xml/custom_source/settings.xml

@@ -8,8 +8,7 @@
 
   <!-- Starting source -->
   <source>
-    <!--<library>./source_sample.so</library> -->
-    <library>./source_ring.so</library> 
+    <library>build/libsource.so</library>
   </source>
 
 </settings>

examples/xml/custom_source/source_ring.cpp

@@ -3,17 +3,16 @@
 #include "openmc/source.h"
 #include "openmc/particle.h"
 
-// you must have external C linkage here otherwise 
+// you must have external C linkage here otherwise
 // dlopen will not find the file
-extern "C" openmc::Particle::Bank sample_source() 
+extern "C" openmc::Particle::Bank sample_source(uint64_t* seed)
 {
   openmc::Particle::Bank particle;
   // wgt
   particle.particle = openmc::Particle::Type::neutron;
-  particle.wgt = 1.0; 
+  particle.wgt = 1.0;
   // position
-    
-  double angle = 2. * M_PI * openmc::prn();
+  double angle = 2. * M_PI * openmc::prn(seed);
   double radius = 3.0;
   particle.r.x = radius * std::cos(angle);
   particle.r.y = radius * std::sin(angle);
@@ -22,5 +21,5 @@ extern "C" openmc::Particle::Bank sample_source()
   particle.u = {1.0, 0.0, 0.0};
   particle.E = 14.08e6;
   particle.delayed_group = 0;
-  return particle;    
+  return particle;
 }

include/openmc/source.h

@@ -59,10 +59,6 @@ private:
 //! Initialize source bank from file/distribution
 extern "C" void initialize_source();
 
-// as yet uncreated function to sample a source from a shared object
-// 
-// extern "C" Particle::Bank sample_source();
-
 //! Sample a site from all external source distributions in proportion to their
 //! source strength
 //! \param[inout] seed Pseudorandom seed pointer
@@ -73,7 +69,7 @@ Particle::Bank sample_external_source(uint64_t* seed);
 void fill_source_bank_fixedsource();
 
 //! Fill source bank at the end of a generation for dlopen based source simulation
-void fill_source_bank_dlopen_source();
+void fill_source_bank_custom_source();
 
 void free_memory_source();
 

src/source.cpp

@@ -1,10 +1,13 @@
 #include "openmc/source.h"
 
-#include <algorithm> // for move
-#include <sstream> // for stringstream
-
 #if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__))
-#include <dlfcn.h> // for dlopen
+#define HAS_DYNAMIC_LINKING
+#endif
+
+#include <algorithm> // for move
+
+#ifdef HAS_DYNAMIC_LINKING
+#include <dlfcn.h> // for dlopen, dlsym, dlclose, dlerror
 #endif
 
 #include <fmt/core.h>
@@ -77,13 +80,11 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
         settings::path_source));
     }
   } else if (check_for_node(node, "library")) {
-     settings::path_source_library = get_node_value(node, "library", false, true);
-     // check if it exists
-     if (!file_exists(settings::path_source_library)) {
-      	std::stringstream msg;
-      	msg << "Library file " << settings::path_source_library << "' does not exist.";
-	      fatal_error(msg);
-     }
+    settings::path_source_library = get_node_value(node, "library", false, true);
+    if (!file_exists(settings::path_source_library)) {
+      fatal_error(fmt::format("Source library '{}' does not exist.",
+        settings::path_source_library));
+    }
   } else {
 
     // Spatial distribution for external source
@@ -249,7 +250,7 @@ void initialize_source()
 {
   write_message("Initializing source particles...", 5);
 
-  if (settings::path_source != "") {
+  if (!settings::path_source.empty()) {
     // Read the source from a binary file instead of sampling from some
     // assumed source distribution
 
@@ -273,31 +274,27 @@ void initialize_source()
 
     // Close file
     file_close(file_id);
-  } else if ( settings::path_source_library != "" ) {
-    // Get the source from a library object 
-    std::stringstream msg;
-    msg << "Sampling from library source " << settings::path_source << "...";
-    write_message(msg, 6);
+  } else if (!settings::path_source_library.empty()) {
+
+  #ifdef HAS_DYNAMIC_LINKING
+    // Get the source from a library object
+    write_message(fmt::format("Sampling from library source {}...",
+      settings::path_source), 6);
 
-    #if defined (__unix__) || (defined (__APPLE__) && defined (__MACH__))
     // Open the library
-    void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY);
-    if(!source_library) {
-      std::stringstream msg("Couldn't open source library " + settings::path_source_library);
-      fatal_error(msg);
+    auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
+    if (!source_library) {
+      fatal_error("Couldn't open source library " + settings::path_source_library);
     }
-    #else
-    std::stringstream msg("This feature has not yet been implemented for non POSIX systems");
-    fatal_error(msg);
-    #endif
-
-    // load the symbol
-    typedef Particle::Bank (*sample_t)(uint64_t seed);
 
     // reset errors
     dlerror();
 
-    fill_source_bank_dlopen_source();
+    fill_source_bank_custom_source();
+  #else
+    fatal_error("Custom source libraries have not yet been implemented for "
+      "non-POSIX systems");
+  #endif
 
   } else {
     // Generation source sites from specified distribution in user input
@@ -359,47 +356,44 @@ void free_memory_source()
 }
 
 // fill the source bank from the external source
-void fill_source_bank_dlopen_source() 
+void fill_source_bank_custom_source()
 {
-  std::stringstream msg;
-
+#ifdef HAS_DYNAMIC_LINKING
   // Open the library
-  void* source_library = dlopen(settings::path_source_library.c_str(),RTLD_LAZY);
-  if(!source_library) {
-    std::stringstream msg("Couldn't open source library " + settings::path_source_library);
-    fatal_error(msg);
+  auto source_library = dlopen(settings::path_source_library.c_str(), RTLD_LAZY);
+  if (!source_library) {
+    fatal_error("Couldn't open source library " + settings::path_source_library);
   }
 
-  // load the symbol
-  typedef Particle::Bank (*sample_t)(uint64_t seed);
-
   // reset errors
   dlerror();
 
   // get the function from the library
-  sample_t sample_source = (sample_t) dlsym(source_library, "sample_source");
-  const char *dlsym_error = dlerror();
-    
+  using sample_t = Particle::Bank (*)(uint64_t* seed);
+  auto sample_source = reinterpret_cast<sample_t>(dlsym(source_library, "sample_source"));
+
   // check for any dlsym errors
+  auto dlsym_error = dlerror();
   if (dlsym_error) {
-    std::cout << dlsym_error << std::endl;
     dlclose(source_library);
-    fatal_error("Couldn't open the sample_source symbol");
+    fatal_error(fmt::format("Couldn't open the sample_source symbol: {}", dlsym_error));
   }
 
   // Generation source sites from specified distribution in the
-  // library source 
+  // library source
   for (int64_t i = 0; i < simulation::work_per_rank; ++i) {
     // initialize random number seed
-    int64_t id = simulation::total_gen*settings::n_particles +
-      simulation::work_index[mpi::rank] + i + 1;
+    int64_t id = (simulation::total_gen + overall_generation()) *
+      settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
      uint64_t seed = init_seed(id, STREAM_SOURCE);
+
     // sample external source distribution
-    simulation::source_bank[i] = sample_source(seed);
+    simulation::source_bank[i] = sample_source(&seed);
   }
 
   // release the library
   dlclose(source_library);
+#endif
 }
 
 void fill_source_bank_fixedsource()
@@ -414,8 +408,8 @@ void fill_source_bank_fixedsource()
       // sample external source distribution
       simulation::source_bank[i] = sample_external_source(&seed);
     }
-  } else if (settings::path_source.empty() && !settings::path_source.empty()) {
-     fill_source_bank_dlopen_source();
+  } else if (settings::path_source.empty() && !settings::path_source_library.empty()) {
+     fill_source_bank_custom_source();
   }
 }