Align source sampling with parameterised source

Allows the source generator to be run by passing in the parameterised
source.
Also provides an example for how to read the output h5 file.
Requires pugixml to align with OpenMC calls.

Author: DanShort12

.gitignore

@@ -30,6 +30,7 @@
 *.exe
 *.out
 *.app
+source_generator
 
 # Build path
 build/
@@ -41,8 +42,6 @@ __pycache__/
 *.py[cod]
 .venv
 
-source_generator
-
 # Outputs
 *.xml
 *.h5

.gitmodules

@@ -1,3 +1,6 @@
 [submodule "pybind11"]
 	path = pybind11
 	url = https://github.com/pybind/pybind11
+[submodule "pugixml"]
+	path = pugixml
+	url = https://github.com/zeux/pugixml.git

CMakeLists.txt

@@ -66,6 +66,8 @@ add_subdirectory(pybind11)
 
 pybind11_add_module(plasma_source ${plasma_source_pybind_SOURCES})
 
+add_subdirectory(pugixml)
+
 # Build source_generator
 list(APPEND source_generator_SOURCES
   ${SRC_DIR}/source_generator.cpp

README.md

@@ -29,7 +29,7 @@ The package can be built and installed in editable mode by:
 
 ## Usage
 
-The parametric plasma source can be sampled either directly in Python 3 or sampled in an OpenMC simulation.
+The parametric plasma source can be sampled either directly in Python 3, sampled in an OpenMC simulation, or sampled using OpenMC via a standalone executable without simulation.
 
 For a better understanding of the varibles take a look at the [C. Fausser et al](https://www.sciencedirect.com/science/article/pii/S0920379612000853) paper.
 

examples/build_and_run_model.py …ples/openmc_model/build_and_run_model.pyexamples/build_and_run_model.py renamed to examples/openmc_model/build_and_run_model.py examples/build_and_run_model.py renamed to examples/openmc_model/build_and_run_model.py

@@ -0,0 +1,34 @@
+"""
+An example of how to read an initial_source.h5 file, as generated from source_generator
+
+In order to create the initial source, first build the source_generator and then run
+the below, which will create an example source and sample it 1000 times.
+
+../../build/source_generator -l ../../build/source_sampling.so -i "major_radius=9.06, \
+minor_radius=2.92258, elongation=1.557, triangularity=0.27, shafranov_shift=0.44789, \
+pedestal_radius=2.33806, ion_density_pedestal=1.09e+20, ion_density_separatrix=3e+19, \
+ion_density_origin=1.09e+20, ion_density_alpha=1, ion_temperature_pedestal=6.09, \
+ion_temperature_separatrix=0.1, ion_temperature_origin=45.9, \
+ion_temperature_alpha=8.06, ion_temperature_beta=6, plasma_type=plasma, plasma_id=1, \
+number_of_bins=100, minimum_toroidal_angle=0, maximum_toroidal_angle=360" \
+-v 10 -n 1000
+"""
+
+import h5py
+import matplotlib.pyplot as plt
+import numpy as np
+
+path_to_initial_source = "initial_source.h5"
+
+with h5py.File(path_to_initial_source, "r") as f:
+    r = np.sqrt(
+        np.square(f["source_bank"][...]["r"]["x"])
+        + np.square(f["source_bank"][...]["r"]["y"])
+    )
+    z = f["source_bank"][...]["r"]["z"]
+
+plt.scatter(r, z)
+
+plt.ion()
+plt.show()
+plt.gca().set_aspect("equal")

examples/plot_flux.py examples/openmc_model/plot_flux.pyexamples/plot_flux.py renamed to examples/openmc_model/plot_flux.py

@@ -1,6 +1,8 @@
 #include <iostream>
 #include <experimental/filesystem>
 
+#include "pugixml.hpp"
+
 #include "openmc/bank.h"
 #include "openmc/constants.h"
 #include "openmc/message_passing.h"
@@ -10,13 +12,14 @@
 
 namespace source_generator
 {
-  void print_settings()
+  void print_settings(pugi::xml_node &root)
   {
     using namespace openmc;
 
     std::cout << "Settings:" << std::endl;
     std::cout << " Number of particles: " << settings::n_particles << std::endl;
-    std::cout << " Source library: " << settings::path_source_library << std::endl;
+    std::cout << " Source library: " << root.child("source").attribute("library").value() << std::endl;
+    std::cout << " Source parameters: " << root.child("source").attribute("parameters").value() << std::endl;
     std::cout << " Output path: " << settings::path_output << std::endl;
     std::cout << " Verbosity: " << settings::verbosity << std::endl;
     std::cout << std::endl;
@@ -29,6 +32,7 @@ namespace source_generator
     std::cout << std::endl;
     std::cout << "Options:" << std::endl;
     std::cout << " -l,--library      Source library, mandatory" << std::endl;
+    std::cout << " -i,--input        Source definition, mandatory" << std::endl;
     std::cout << " -n,--particles    Number of particles, default 1000" << std::endl;
     std::cout << " -o,--output       Output directory, default {current directory}" << std::endl;
     std::cout << " -v,--verbosity    Verbosity, default 5" << std::endl; 
@@ -51,10 +55,12 @@ namespace source_generator
     settings::verbosity = 5;
   }
 
-  int parse_command_line(int argc, char* argv[])
+  int parse_command_line(int argc, char* argv[], pugi::xml_node &root)
   {
     using namespace openmc;
 
+    root.append_child("source");
+
     for (int i=1; i < argc; ++i)
     {
       std::string arg {argv[i]};
@@ -68,8 +74,14 @@ namespace source_generator
         else if (arg == "-l" || arg == "--library")
         {
           i += 1;
+          root.child("source").append_attribute("library").set_value(argv[i]);
           settings::path_source_library = argv[i];
         }
+        else if (arg == "-i" || arg == "--input")
+        {
+          i += 1;
+          root.child("source").append_attribute("parameters").set_value(argv[i]);
+        }
         else if (arg == "-o" || arg == "--output")
         {
           i += 1;
@@ -88,22 +100,40 @@ namespace source_generator
       }
     }
 
-    if (settings::path_source_library.empty())
+    bool missing_arg = false;
+
+    if (!root.child("source").attribute("library"))
     {
       std::cout << "The --library or -l argument is mandatory and must be set." << std::endl;
+      missing_arg = true;
+    }
+
+    if (!root.child("source").attribute("parameters"))
+    {
+      std::cout << "The --input or -i argument is mandatory and must be set." << std::endl;
+      missing_arg = true;
+    }
+
+    if (missing_arg)
+    {
       source_generator::print_help();
       return -1;
     }
-
-    return 0;
+    else
+    {
+      return 0;
+    }
   }
 }
 
 int main(int argc, char* argv[])
 {
+  pugi::xml_document doc;
+  pugi::xml_node root = doc.append_child("settings");
+
   source_generator::set_defaults();
 
-  int run = source_generator::parse_command_line(argc, argv);
+  int run = source_generator::parse_command_line(argc, argv, root);
 
   if (run < 0)
   {
@@ -112,10 +142,11 @@ int main(int argc, char* argv[])
 
   if (openmc::settings::verbosity >= 5)
   {
-    source_generator::print_settings();
+    source_generator::print_settings(root);
   }
 
   std::cout << "Sampling source:" << std::endl;
+  openmc::SourceDistribution source = openmc::SourceDistribution(root.child("source"));
   openmc::calculate_work();
   openmc::allocate_banks();
   openmc::initialize_source();