Adding inline comparison to reference value.

CMakeLists.txt

@@ -64,7 +64,8 @@ ENDIF()
 # COPY DATA FILES
 SET(DATA_FILES
     data/c5g7_252g.xml
-    data/c5g7_3d.xml)
+    data/c5g7_3d.xml
+    data/c5g7_flux_ref.txt)
 INSTALL(FILES ${DATA_FILES} DESTINATION data)
 
 # Do all of the processing for this Tribits project

packages/MC/mc_driver/Manager.hh

@@ -87,6 +87,11 @@ class Manager : public Manager_Base
     // Solve the problem.
     void solve();
 
+    // Check solution for correctness
+    void check_solution(double      k_ref,
+                        double      k_std_ref,
+                        std::string flux_file);
+
     // Output.
     void output();
 

packages/MC/mc_driver/Manager.t.hh

@@ -11,8 +11,10 @@
 #ifndef mc_driver_Manager_t_hh
 #define mc_driver_Manager_t_hh
 
+#include <sstream>
 #include <fstream>
-#include <iomanip>
+#include <iostream>
+//#include <iomanip>
 
 #include "Teuchos_XMLParameterListHelpers.hpp"
 
@@ -24,6 +26,7 @@
 #include "mc/Fission_Source.hh"
 #include "mc/Uniform_Source.hh"
 #include "mc/KCode_Solver.hh"
+#include "mc/Cell_Tally.hh"
 #include "Manager.hh"
 
 namespace mc
@@ -201,6 +204,152 @@ void Manager<Geometry>::solve()
     }
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * \brief Check solution for correctness
+ */
+template <class Geometry>
+void Manager<Geometry>::check_solution(double      k_ref,
+                                       double      k_std_ref,
+                                       std::string flux_ref_file)
+{
+    INSIST(d_keff_solver, "Can only check solution with keff solver.");
+
+    if (profugus::node() > 0)
+        return;
+
+    // Store existing iostream state
+    std::ios io_state(nullptr);
+    io_state.copyfmt(std::cout);
+
+    // Compare computed eigenvalue to reference
+    auto keff = d_keff_solver->keff_tally();
+    CHECK(keff);
+    auto keff_mean = keff->mean();
+    auto keff_var  = keff->variance();
+    std::cout << "Computed eigenvalue:  "
+        << std::fixed << std::setprecision(6) << keff_mean << " +/- "
+        << std::scientific << std::setprecision(2) << std::sqrt(keff_var)
+        << std::endl;
+    std::cout << "Reference eigenvalue: "
+        << std::fixed << std::setprecision(6) << k_ref << " +/- "
+        << std::scientific << std::setprecision(2) << k_std_ref << std::endl;
+    auto keff_diff = std::abs(keff_mean - k_ref);
+    auto keff_combined_std = std::sqrt(keff_var + k_std_ref * k_std_ref);
+    auto keff_rel_diff = keff_diff / keff_combined_std;
+    std::cout << "Difference in eigenvalue is "
+        << std::scientific << std::setprecision(2) << keff_diff << ", which is "
+        << std::fixed << std::setprecision(2) << keff_rel_diff
+        << " relative combined standard deviations." << std::endl;
+    std::cout << std::endl;
+
+    // Restore iostream state
+    std::cout.copyfmt(io_state);
+
+    // Load reference fluxes from file
+    std::ifstream ref_file(flux_ref_file);
+    CHECK(ref_file.is_open());
+    std::vector<double> ref_flux;
+    std::vector<double> ref_flux_std;
+    std::string buffer;
+    while (std::getline(ref_file, buffer))
+    {
+        std::stringstream line(buffer);
+
+        double flux_val, flux_std;
+        line >> flux_val >> flux_std;
+
+        ref_flux.push_back(flux_val);
+        ref_flux_std.push_back(flux_std);
+    }
+
+    // Get computed flux values
+    using Cell_Tally_t = profugus::Cell_Tally<Geometry>;
+    auto tallier = d_keff_solver->tallier();
+    CHECK(tallier->is_finalized());
+    CHECK(tallier->num_pathlength_tallies() >= 2);
+    std::shared_ptr<Cell_Tally_t> flux_tally;
+    for (auto tally : *tallier)
+    {
+        if (tally->name() == "cell")
+            flux_tally = std::dynamic_pointer_cast<Cell_Tally_t>(tally);
+    }
+    CHECK(flux_tally);
+    auto results = flux_tally->results();
+    CHECK(results.size() == ref_flux.size());
+    int num_cells = ref_flux.size();
+
+    // Compare computed and reference solutions
+    // We print the number that are within 1, 2, and 3 combined standard
+    // deviations of the reference, but we only consider the test to have
+    // failed if any cells are more than 5 sigma from reference.
+    int num_1sigma = 0;
+    int num_2sigma = 0;
+    int num_3sigma = 0;
+    int num_failed = 0;
+    int cell_ind = 0;
+    for (auto& val : results)
+    {
+        double flux_val = val.second.first;
+        double flux_std = val.second.second;
+
+        double diff = std::abs(flux_val - ref_flux[cell_ind]);
+        double combined_std = std::sqrt(
+            flux_std * flux_std + ref_flux_std[cell_ind] * ref_flux_std[cell_ind]);
+        double rel_diff = diff / combined_std;
+
+        if (rel_diff < 1.0)
+            num_1sigma++;
+        if (rel_diff < 2.0)
+            num_2sigma++;
+        if (rel_diff < 3.0)
+            num_3sigma++;
+        if (rel_diff > 5.0)
+            num_failed++;
+
+        cell_ind++;
+    }
+    int expected_1sigma = 0.68  * num_cells;
+    int expected_2sigma = 0.95  * num_cells;
+    int expected_3sigma = 0.997 * num_cells;
+    std::cout << "Comparing computed tally results to reference solution for "
+        << num_cells << " cells" << std::endl;
+    std::cout << num_1sigma << " cells were within 1 combined relative "
+        << "standard deviation of reference solution" << std::endl;
+    std::cout << num_2sigma << " cells were within 2 combined relative "
+        << "standard deviations of reference solution" << std::endl;
+    std::cout << num_3sigma << " cells were within 3 combined relative "
+        << "standard deviations of reference solution" << std::endl;
+    std::cout << "Expected "
+        << expected_1sigma << " within 1 sigma, "
+        << expected_2sigma << " within 2 sigma, and "
+        << expected_3sigma << " within 3 sigma" << std::endl << std::endl;
+
+    // Print finale pass/fail status
+    if (keff_rel_diff < 3.0)
+    {
+        std::cout << "CHECK PASSED: Eigenvalue is within 3 standard "
+            << "deviations of the reference solution." << std::endl;
+    }
+    else
+    {
+        std::cout << "CHECK FAILED: Eigenvalue is more than 3 standard"
+            << " deviations from the reference solution." << std::endl;
+    }
+    if (num_failed == 0)
+    {
+        std::cout << "CHECK PASSED: All cells were within 5 standard "
+            << "deviations of the reference solution." << std::endl;
+    }
+    else
+    {
+        std::cout << "CHECK FAILED: " << num_failed << " cells were more than "
+            << "5 standard deviations from the reference solution"
+            << std::endl;
+    }
+    std::cout << std::endl;
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * \brief Do output.

packages/MC/mc_driver/Manager_Base.hh

@@ -42,6 +42,12 @@ class Manager_Base
       // Solve the problem.
       virtual void solve() = 0;
 
+      // Check solution for correctness
+      virtual void check_solution(
+          double      k_ref,
+          double      k_std_ref,
+          std::string flux_file) = 0;
+
       // Output data.
       virtual void output() = 0;
 };

packages/MC/mc_driver/profugus.cc

@@ -102,16 +102,6 @@ std::tuple<std::string,def::size_type> parse_input_arguments(
     {
         std::stringstream np_stream(*(iter+1));
         np_stream >> Np;
-        if (Np < 0)
-        {
-            if (node == 0)
-            {
-                std::cout << "Invalid particle count " << Np << std::endl
-                          << "Particle count must be positive." << std::endl;
-            }
-            exit(1);
-        }
-
         if (node == 0)
             std::cout << "Executing with Np = " << Np << std::endl;
     }
@@ -135,7 +125,7 @@ int main(int argc, char *argv[])
     nodes = profugus::nodes();
 
     profugus::pcout << "=======================================\n"
-                    << "    Profugus MC Mini-APP               \n"
+                    << "    ProfugusMC Mini-APP                \n"
                     << "    (C) ORNL, Battelle, 2014           \n"
                     << "=======================================\n"
                     << profugus::endl;
@@ -154,9 +144,9 @@ int main(int argc, char *argv[])
     try
     {
         std::string install_dir = CMAKE_INSTALL_PREFIX;
-        std::string data_dir = install_dir + "/data";
-        std::string xml_file = data_dir + "/c5g7_3d.xml";
-        std::string xs_file = data_dir + "/c5g7_252g.xml";
+        std::string data_dir = install_dir + "/data/";
+        std::string xml_file = data_dir + "c5g7_3d.xml";
+        std::string xs_file = data_dir + "c5g7_252g.xml";
 
         if (node ==0)
         {
@@ -189,6 +179,12 @@ int main(int argc, char *argv[])
         // solve the problem
         manager->solve();
 
+        // check solution
+        double k_ref = 1.19383926;
+        double k_std_ref = 6.131761e-06;
+        std::string flux_file = data_dir + "c5g7_flux_ref.txt";
+        manager->check_solution(k_ref, k_std_ref, flux_file);
+
         // output
         manager->output();
     }