Getting solution checking to work with Cuda.

packages/MC/CMakeLists.txt

@@ -123,6 +123,7 @@ LIST(APPEND SOURCES ${MC_SOURCES})
 
 FILE(GLOB DRIVER_HEADERS mc_driver/*.hh)
 SET(DRIVER_SOURCES
+  mc_driver/Manager_Base.cc
   mc_driver/Manager.pt.cc
   mc_driver/Manager_Builder.cc
   mc_driver/Problem_Builder.pt.cc

packages/MC/cuda_mc/Manager_Cuda.hh

@@ -104,6 +104,15 @@ class Manager_Cuda : public mc::Manager_Base
     // Output.
     void output();
 
+  protected:
+
+    // Get keff
+    void get_keff(double& keff_mean, double& keff_var);
+
+    // Get flux tally result
+    void get_flux(std::vector<double>& flux_val,
+                  std::vector<double>& flux_std);
+
   private:
     // >>> IMPLEMENTATION
 

packages/MC/cuda_mc/Manager_Cuda.t.cuh

 //----------------------------------*-C++-*----------------------------------//
-/*!
- * \file   cuda_mc/Manager_Cuda.t.cuh
+/*!  * \file   cuda_mc/Manager_Cuda.t.cuh
  * \author Steven Hamilton
  * \date   Wed Jun 18 11:21:16 2014
  * \brief  Manager_Cuda member definitions.
@@ -443,6 +442,36 @@ void Manager_Cuda<Geometry_DMM>::build_physics(RCP_ParameterList master_db)
     d_physics = SDP_Physics( physics_host );
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * \brief Get keff
+ */
+template <class Geometry_DMM>
+void Manager_Cuda<Geometry_DMM>::get_keff(double& keff_mean, double& keff_var)
+{
+    REQUIRE(d_keff_solver);
+    auto keff_tally = d_keff_solver->keff_tally();
+    CHECK(keff_tally);
+    keff_mean = keff_tally->mean();
+    keff_var = keff_tally->variance();
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * \brief Extract flux from tally
+ */
+template <class Geometry_DMM>
+void Manager_Cuda<Geometry_DMM>::get_flux(
+        std::vector<double>& flux_val,
+        std::vector<double>& flux_std)
+{
+    REQUIRE(d_tallier_dmm);
+
+    // Get mean and std dev from cell tally
+    auto cell_tally = d_tallier_dmm->cell_tally();
+    flux_val = cell_tally->results();
+    flux_std = cell_tally->std_dev();
+}
 
 } // end namespace cuda_mc
 

packages/MC/cuda_mc/Tallier.cuh

@@ -156,6 +156,9 @@ class Tallier_DMM : public cuda::Device_Memory_Manager<Tallier<Geometry>>
     // Add keff tally.
     void add_keff_tally(SP_Keff_Tally_DMM tally);
 
+    // Access cell tally
+    SP_Cell_Tally_DMM cell_tally() const {return d_cell_tally_dmm;}
+
     // >>> HOST-SIDE TALLY OPERATIONS
 
     // Tell the tallies to begin active kcode cycles

packages/MC/mc_driver/Manager.hh

@@ -95,6 +95,15 @@ class Manager : public Manager_Base
     // Output.
     void output();
 
+  protected:
+
+    // Get keff
+    void get_keff(double& keff_mean, double& keff_var);
+
+    // Get computed flux solution
+    void get_flux(std::vector<double>& flux_val,
+                  std::vector<double>& flux_std);
+
   private:
     // >>> IMPLEMENTATION
 

packages/MC/mc_driver/Manager.t.hh

@@ -206,63 +206,14 @@ void Manager<Geometry>::solve()
 
 //---------------------------------------------------------------------------//
 /*!
- * \brief Check solution for correctness
+ * \brief Get flux from cell tally
  */
 template <class Geometry>
-void Manager<Geometry>::check_solution(double      k_ref,
-                                       double      k_std_ref,
-                                       std::string flux_ref_file)
+void Manager<Geometry>::get_flux(std::vector<double>& flux_val,
+                                 std::vector<double>& flux_std)
 {
     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();
@@ -276,78 +227,32 @@ void Manager<Geometry>::check_solution(double      k_ref,
     }
     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;
+    int num_cells = results.size();
+
+    // Copy results into vectors
+    flux_val.resize(num_cells);
+    flux_std.resize(num_cells);
+    int cell = 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++;
+        flux_val[cell] = val.second.first;
+        flux_std[cell] = val.second.second;
+        cell++;
     }
-    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 Get keff from solver
+ */
+template <class Geometry>
+void Manager<Geometry>::get_keff(double& keff_mean, double& keff_var)
+{
+    REQUIRE(d_keff_solver);
+    auto keff_tally = d_keff_solver->keff_tally();
+    CHECK(keff_tally);
+    keff_mean = keff_tally->mean();
+    keff_var  = keff_tally->variance();
 }
 
 //---------------------------------------------------------------------------//

packages/MC/mc_driver/Manager_Base.cc

+//---------------------------------*-C++-*-----------------------------------//
+/*!
+ * \file   MC/mc_driver/Manager_Base.cc
+ * \author Steven Hamilton
+ * \date   Thu Aug 30 13:55:07 2018
+ * \brief  Manager_Base class definitions.
+ * \note   Copyright (c) 2018 Oak Ridge National Laboratory, UT-Battelle, LLC.
+ */
+//---------------------------------------------------------------------------//
+
+#include "Manager_Base.hh"
+
+#include "Utils/comm/global.hh"
+#include "Utils/harness/DBC.hh"
+
+namespace mc
+{
+//---------------------------------------------------------------------------//
+/*!
+ * \brief Check solution for correctness
+ */
+void Manager_Base::check_solution(double      k_ref,
+                                  double      k_std_ref,
+                                  std::string flux_ref_file)
+{
+    if (profugus::node() > 0)
+        return;
+
+    // Store existing iostream state
+    std::ios io_state(nullptr);
+    io_state.copyfmt(std::cout);
+
+    // Compare computed eigenvalue to reference
+    double keff_mean, keff_var;
+    this->get_keff(keff_mean, keff_var);
+    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
+    std::vector<double> computed_flux;
+    std::vector<double> computed_std;
+    this->get_flux(computed_flux, computed_std);
+    CHECK(computed_flux.size() == ref_flux.size());
+    CHECK(computed_std.size()  == ref_flux.size());
+    int num_cells = ref_flux.size();
+
+    std::cout << "Flux comparison" << std::endl;
+    for (int cell = 0; cell < num_cells; ++cell)
+    {
+        std::cout << computed_flux[cell] << " " << ref_flux[cell] << " " << computed_std[cell] << " " << ref_flux_std[cell] << std::endl;
+    }
+
+    // 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 (int cell = 0; cell < num_cells; ++cell)
+    {
+        double flux_val = computed_flux[cell];
+        double flux_std = computed_std[cell];
+
+        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;
+}
+
+//---------------------------------------------------------------------------//
+} // end namespace mc
+
+//---------------------------------------------------------------------------//
+// end of MC/mc_driver/Manager_Base.cc
+//---------------------------------------------------------------------------//

packages/MC/mc_driver/Manager_Base.hh

@@ -43,13 +43,22 @@ class Manager_Base
       virtual void solve() = 0;
 
       // Check solution for correctness
-      virtual void check_solution(
+      void check_solution(
           double      k_ref,
           double      k_std_ref,
-          std::string flux_file) = 0;
+          std::string flux_file);
 
       // Output data.
       virtual void output() = 0;
+
+  protected:
+
+      // Get keff
+      virtual void get_keff(double& keff_mean, double& keff_var) = 0;
+
+      // Get computed flux solution
+      virtual void get_flux(std::vector<double>& flux_val,
+                            std::vector<double>& flux_std) = 0;
 };
 
 } // end namespace mc