Cross section and source post-processing (#167)

* Automatically write cross sections to a file when writing a model
* Added script to visualize cross sections
* Added neutron source extraction post-processing function

CMakeLists.txt

@@ -335,6 +335,7 @@ set(UM2_SOURCES
     "src/physics/cmfd.cpp"
     "src/mpact/model.cpp"
     "src/mpact/powers.cpp"
+    "src/mpact/source.cpp"
     "src/gmsh/base_gmsh_api.cpp"
     "src/gmsh/io.cpp"
     "src/gmsh/model.cpp"

include/um2/mpact/source.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include <um2/config.hpp>
+#include <um2/mesh/polytope_soup.hpp>
+#include <um2/physics/material.hpp>
+#include <um2/stdlib/vector.hpp>
+
+namespace um2::mpact
+{
+
+//==============================================================================
+// getSource
+//==============================================================================
+// Given an MPACT FSR output and a list of materials, return a vector
+// containing the fission and scattering sources for each cell.
+
+PURE [[nodiscard]] auto
+getSource(PolytopeSoup const & soup, Vector<Material> const & materials) -> Vector<Vec2F>;
+
+} // namespace um2::mpact

models/c5g7/2d/c5g7_extract_source.cpp

@@ -0,0 +1,54 @@
+// NOLINTBEGIN(misc-include-cleaner)
+
+#include <um2.hpp>
+#include <um2/mpact/source.hpp>
+#include <um2/stdlib/numeric/iota.hpp>
+
+auto
+main(int argc, char ** argv) -> int
+{
+  um2::initialize();
+
+  // Check the number of arguments
+  if (argc != 2) {
+    um2::String const exec_name(argv[0]);
+    um2::logger::error("Usage: ", exec_name, " <FSRmesh>");
+    return 1;
+  }
+
+  // Get the FSR file name
+  um2::String const filename(argv[1]);
+
+  // Get the model's materials
+  um2::Vector<um2::Material> const materials = um2::getC5G7Materials();
+
+  // Import the FSR mesh
+  um2::PolytopeSoup soup(filename);
+
+  // Extract the source
+  um2::Vector<um2::Vec2F> const source = um2::mpact::getSource(soup, materials);
+
+  // Compute the scattering and fission source given k_eff
+  Float constexpr k_eff = 1.1864063;
+  Float constexpr lambda = 1 / k_eff;
+  um2::Vector<Float> scattering_source(source.size());
+  um2::Vector<Float> fission_source(source.size());
+  um2::Vector<Float> total_source(source.size());
+  for (Int i = 0; i < source.size(); ++i) {
+    scattering_source[i] = source[i][0];
+    fission_source[i] = lambda * source[i][1];
+    total_source[i] = scattering_source[i] + fission_source[i];
+  }
+
+  // Add each of the sources as an elset to the soup
+  um2::Vector<Int> ids(source.size());
+  um2::iota(ids.begin(), ids.end(), 0);
+  soup.addElset("scattering_source", ids, scattering_source);
+  soup.addElset("fission_source", ids, fission_source);
+  soup.addElset("total_source", ids, total_source);
+
+  soup.write("c5g7_source.xdmf");
+  um2::finalize();
+  return 0;
+}
+// NOLINTEND(misc-include-cleaner)

models/c5g7/CMakeLists.txt

@@ -1,3 +1,4 @@
 um2_add_executable(./2d/c5g7_2d_variable_grid.cpp)
 um2_add_executable(./2d/c5g7_2d_advanced_mesh.cpp)
+um2_add_executable(./2d/c5g7_extract_source.cpp)
 um2_add_executable(./3d/c5g7_3d_standard_moc.cpp)

models/crocus/CMakeLists.txt

@@ -1,2 +1,3 @@
 um2_add_executable(./crocus_2d_variable_grid.cpp)
 um2_add_executable(./crocus_2d_advanced_mesh.cpp)
+um2_add_executable(./crocus_2d_extract_source.cpp)

models/crocus/crocus_2d_extract_source.cpp

@@ -0,0 +1,108 @@
+// NOLINTBEGIN(misc-include-cleaner)
+
+#include <um2.hpp>
+#include <um2/mpact/source.hpp>
+#include <um2/stdlib/numeric/iota.hpp>
+
+auto
+main(int argc, char ** argv) -> int
+{
+  um2::initialize();
+
+  // Check the number of arguments
+  if (argc != 2) {
+    um2::String const exec_name(argv[0]);
+    um2::logger::error("Usage: ", exec_name, " <FSRmesh>");
+    return 1;
+  }
+
+  // Get the FSR file name
+  um2::String const filename(argv[1]);
+
+  //============================================================================
+  // Materials
+  //============================================================================
+  um2::XSLibrary const xslib(um2::settings::xs::library_path + "/" +
+                             um2::mpact::XSLIB_51G);
+
+  // NOTE: number densities should be computed from the source, but I have simply
+  // ripped them from another CROCUS model for now.
+
+  Float constexpr temp = 293.15; // K. pg. 744 Sec. 3.1
+
+  // UO2
+  //---------------------------------------------------------------------------
+  um2::Material uo2;
+  uo2.setName("UO2");
+  uo2.setDensity(10.556); // pg. 742 Sec. 2.3
+  uo2.setTemperature(temp);
+  uo2.setColor(um2::orange); // Match Fig. 4
+  uo2.addNuclide(92235, 4.30565e-04);
+  uo2.addNuclide(92238, 2.31145e-02);
+  uo2.addNuclide(8016, 4.70902e-02);
+  uo2.populateXSec(xslib);
+
+  // Clad
+  //---------------------------------------------------------------------------
+  um2::Material clad;
+  clad.setName("Clad");
+  clad.setDensity(2.70); // pg. 743 Table 1
+  clad.setTemperature(temp);
+  clad.setColor(um2::slategray);
+  clad.addNuclide(13027, 6.02611e-02);
+  clad.populateXSec(xslib);
+
+  // Umetal
+  //---------------------------------------------------------------------------
+  um2::Material umetal;
+  umetal.setName("Umetal");
+  umetal.setDensity(18.677); // pg. 742 Sec. 2.3
+  umetal.setTemperature(temp);
+  umetal.setColor(um2::red);
+  umetal.addNuclide(92235, 4.53160e-04);
+  umetal.addNuclide(92238, 4.68003e-02);
+  umetal.populateXSec(xslib);
+
+  // Water
+  //---------------------------------------------------------------------------
+  um2::Material water;
+  water.setName("Water");
+  water.setDensity(0.9983); // pg. 743 Table 1
+  water.setTemperature(temp);
+  water.setColor(um2::blue);
+  water.addNuclide(1001, 6.67578e-02);
+  water.addNuclide(8016, 3.33789e-02);
+  water.populateXSec(xslib);
+
+  um2::Vector<um2::Material> const materials = {uo2, clad, umetal, water};
+
+  // Import the FSR mesh
+  um2::PolytopeSoup soup(filename);
+
+  // Extract the source
+  um2::Vector<um2::Vec2F> const source = um2::mpact::getSource(soup, materials);
+
+  // Compute the scattering and fission source given k_eff
+  Float constexpr k_eff = 1.0305077;
+  Float constexpr lambda = 1 / k_eff;
+  um2::Vector<Float> scattering_source(source.size());
+  um2::Vector<Float> fission_source(source.size());
+  um2::Vector<Float> total_source(source.size());
+  for (Int i = 0; i < source.size(); ++i) {
+    scattering_source[i] = source[i][0];
+    fission_source[i] = lambda * source[i][1];
+    total_source[i] = scattering_source[i] + fission_source[i];
+  }
+
+  // Add each of the sources as an elset to the soup
+  um2::Vector<Int> ids(source.size());
+  um2::iota(ids.begin(), ids.end(), 0);
+  soup.addElset("scattering_source", ids, scattering_source);
+  soup.addElset("fission_source", ids, fission_source);
+  soup.addElset("total_source", ids, total_source);
+
+  soup.write("source.xdmf");
+  um2::finalize();
+  return 0;
+}
+// NOLINTEND(misc-include-cleaner)

models/fnr/CMakeLists.txt

@@ -1,3 +1,4 @@
 um2_add_executable(./fnr_heu_2d_single_asy.cpp)
 um2_add_executable(./fnr_leu_2d_variable_grid.cpp) 
 um2_add_executable(./fnr_leu_2d_advanced_mesh.cpp)
+um2_add_executable(./fnr_leu_2d_extract_source.cpp)

models/fnr/fnr_leu_2d_advanced_mesh.cpp

@@ -692,7 +692,7 @@ main(int argc, char ** argv) -> int
 
   um2::gmsh::model::mesh::setMeshFieldFromKnudsenNumber(2, model.materials(), target_kn,
                                                         mfp_threshold, mfp_scale);
-  um2::gmsh::model::mesh::generateMesh(um2::MeshType::Tri);
+  um2::gmsh::model::mesh::generateMesh(um2::MeshType::QuadraticTri);
   um2::gmsh::write("fnr_2d.inp");
 
   //===========================================================================

models/fnr/fnr_leu_2d_extract_source.cpp

@@ -0,0 +1,158 @@
+// NOLINTBEGIN(misc-include-cleaner)
+
+#include <um2.hpp>
+#include <um2/mpact/source.hpp>
+#include <um2/stdlib/numeric/iota.hpp>
+
+auto
+main(int argc, char ** argv) -> int
+{
+  um2::initialize();
+
+  // Check the number of arguments
+  if (argc != 2) {
+    um2::String const exec_name(argv[0]);
+    um2::logger::error("Usage: ", exec_name, " <FSRmesh>");
+    return 1;
+  }
+
+  // Get the FSR file name
+  um2::String const filename(argv[1]);
+
+  //============================================================================
+  // Materials
+  //============================================================================
+  um2::XSLibrary const xslib(um2::settings::xs::library_path + "/" +
+                             um2::mpact::XSLIB_51G);
+
+  // Aluminum
+  //---------------------------------------------------------------------------
+  um2::Material aluminum;
+  aluminum.setName("Aluminum");
+  aluminum.setDensity(2.7); // g/cm^3
+  aluminum.setTemperature(300.0);
+  aluminum.setColor(um2::gray);
+  aluminum.addNuclideWt("Al27", 0.9725);
+  aluminum.addNuclideWt("Mg00", 0.01);
+  aluminum.addNuclideWt("Si00", 0.006);
+  aluminum.addNuclideWt("Fe00", 0.0035);
+  aluminum.addNuclideWt("Cu63", 0.00205437585615717);
+  aluminum.addNuclideWt("Cu65", 0.00094562414384283);
+  aluminum.addNuclideWt("Cr00", 0.003);
+  // No Zinc
+  aluminum.addNuclideWt("Ti00", 0.0005);
+  aluminum.addNuclideWt("Mn55", 0.0005);
+  aluminum.populateXSec(xslib);
+
+  // Water
+  //---------------------------------------------------------------------------
+  um2::Material h2o;
+  h2o.setName("Water");
+  h2o.setDensity(0.99821); // g/cm^3
+  h2o.setTemperature(300.0);
+  h2o.setColor(um2::blue);
+  h2o.addNuclidesAtomPercent({"H1", "O16"}, {2.0 / 3.0, 1.0 / 3.0});
+  h2o.populateXSec(xslib);
+
+  // Heavy water
+  //---------------------------------------------------------------------------
+  um2::Material d2o;
+  d2o.setName("HeavyWater");
+  d2o.setDensity(1.11); // g/cm^3
+  d2o.setTemperature(300.0);
+  d2o.setColor(um2::darkblue);
+  d2o.addNuclidesAtomPercent({"H1", "H2", "O16"}, {0.005 / 3.0, 1.995 / 3.0, 1.0 / 3.0});
+  d2o.populateXSec(xslib);
+
+  // LEU fuel
+  //---------------------------------------------------------------------------
+  um2::Material fuel;
+  fuel.setName("Fuel");
+  fuel.setDensity(2.89924);
+  fuel.setTemperature(300.0);
+  fuel.setColor(um2::red);
+  // Weights from (1), pg. 347
+  Float const u235_wt = 140.61;
+  Float const u234_wt = 1.51;
+  Float const u236_wt = 0.75;
+  Float const u238_wt = 8.04;
+  Float const aluminum_wt = 908;
+  Float const iron_wt = 3.7;
+  Float const heu_wt = u235_wt + u234_wt + u236_wt + u238_wt + aluminum_wt + iron_wt;
+  fuel.addNuclideWt("U235", u235_wt / heu_wt);
+  fuel.addNuclideWt("U234", u234_wt / heu_wt);
+  fuel.addNuclideWt("U236", u236_wt / heu_wt);
+  fuel.addNuclideWt("U238", u238_wt / heu_wt);
+  fuel.addNuclideWt("Al27", 0.9725 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Mg00", 0.01 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Si00", 0.006 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Cu63", 0.00205437585615717 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Cu65", 0.00094562414384283 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Cr00", 0.003 * aluminum_wt / heu_wt);
+  // No Zinc
+  fuel.addNuclideWt("Ti00", 0.0005 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Mn55", 0.0005 * aluminum_wt / heu_wt);
+  fuel.addNuclideWt("Fe00", iron_wt / heu_wt + 0.0035 * aluminum_wt / heu_wt);
+  fuel.populateXSec(xslib);
+
+  // Borated steel
+  //---------------------------------------------------------------------------
+  um2::Material borated_steel;
+  borated_steel.setName("BoratedSteel");
+  borated_steel.setDensity(8.0369);
+  borated_steel.setTemperature(300.0);
+  borated_steel.setColor(um2::black);
+  // From (1), pg. 346
+  borated_steel.addNuclide("B10", 0.001108);
+  borated_steel.addNuclide("B11", 0.005184);
+  borated_steel.addNuclide("Fe00", 0.05644);
+  borated_steel.addNuclide("Ni00", 0.0113);
+  borated_steel.addNuclide("Cr00", 0.0164);
+  borated_steel.populateXSec(xslib);
+
+  // Steel
+  //---------------------------------------------------------------------------
+  um2::Material steel;
+  steel.setName("Steel");
+  steel.setDensity(7.85);
+  steel.setTemperature(300.0);
+  steel.setColor(um2::darkgray);
+  // Same as above without boron
+  steel.addNuclide("Fe00", 0.05644);
+  steel.addNuclide("Ni00", 0.0113);
+  steel.addNuclide("Cr00", 0.0164);
+  steel.populateXSec(xslib);
+
+  um2::Vector<um2::Material> const materials = {aluminum, h2o,           d2o,
+                                                fuel,     borated_steel, steel};
+
+  // Import the FSR mesh
+  um2::PolytopeSoup soup(filename);
+
+  // Extract the source
+  um2::Vector<um2::Vec2F> const source = um2::mpact::getSource(soup, materials);
+
+  // Compute the scattering and fission source given k_eff
+  Float constexpr k_eff = 1.1165101;
+  Float constexpr lambda = 1 / k_eff;
+  um2::Vector<Float> scattering_source(source.size());
+  um2::Vector<Float> fission_source(source.size());
+  um2::Vector<Float> total_source(source.size());
+  for (Int i = 0; i < source.size(); ++i) {
+    scattering_source[i] = source[i][0];
+    fission_source[i] = lambda * source[i][1];
+    total_source[i] = scattering_source[i] + fission_source[i];
+  }
+
+  // Add each of the sources as an elset to the soup
+  um2::Vector<Int> ids(source.size());
+  um2::iota(ids.begin(), ids.end(), 0);
+  soup.addElset("scattering_source", ids, scattering_source);
+  soup.addElset("fission_source", ids, fission_source);
+  soup.addElset("total_source", ids, total_source);
+
+  soup.write("source.xdmf");
+  um2::finalize();
+  return 0;
+}
+// NOLINTEND(misc-include-cleaner)

models/vera/CMakeLists.txt

@@ -3,3 +3,4 @@ um2_add_executable(./problem_1/vera_1a_nogap.cpp)
 um2_add_executable(./problem_4/2d/vera_4_unrodded_2d_variable_grid.cpp)
 um2_add_executable(./problem_4/2d/vera_4_rodded_2d_variable_grid.cpp)
 um2_add_executable(./problem_4/2d/vera_4_rodded_2d_advanced_mesh.cpp)
+um2_add_executable(./problem_4/2d/vera_4_rodded_extract_source.cpp)