Model updates + material volumes and energy spectrum post-processing (#…

…166)

*Updated FNR and Nuscale models

*A material volumes file is now written when an MPACT model is written

*Added post-processing utilities for computing the energy spectrum

models/fnr/CMakeLists.txt

@@ -1,2 +1,3 @@
 um2_add_executable(./fnr_heu_2d_single_asy.cpp)
-#um2_add_executable(./fnr_heu_2d.cpp)
+um2_add_executable(./fnr_leu_2d_variable_grid.cpp) 
+um2_add_executable(./fnr_leu_2d_advanced_mesh.cpp)

models/nuscale/nuscale_2d.cpp

@@ -36,17 +36,29 @@ main(int argc, char ** argv) -> int
   //===========================================================================
 
   // Check the number of arguments
-  if (argc != 2) {
+  if (argc != 4) {
     um2::String const exec_name(argv[0]);
-    um2::logger::error("Usage: ", exec_name, " num_coarse_cells");
+    um2::logger::error("Usage: ", exec_name, " target_kn mfp_threshold mfp_scale");
     return 1;
   }
 
   char * end = nullptr;
-  Int const num_coarse_cells = um2::strto<Int>(argv[1], &end);
+  Float const target_kn = um2::strto<Float>(argv[1], &end);
   ASSERT(end != nullptr);
-  ASSERT(num_coarse_cells > 0);
+  ASSERT(target_kn > 0);
 
+  Float const mfp_threshold = um2::strto<Float>(argv[2], &end);
+  ASSERT(end != nullptr);
+  end = nullptr;
+
+  Float const mfp_scale = um2::strto<Float>(argv[3], &end);
+  ASSERT(end != nullptr);
+
+  um2::logger::info("Target Knudsen number: ", target_kn);
+  um2::logger::info("MFP threshold: ", mfp_threshold);
+  um2::logger::info("MFP scale: ", mfp_scale);
+
+  Int constexpr num_coarse_cells = 119;
   um2::String const model_name = "nuscale_2d_" + um2::String(num_coarse_cells) + ".brep";
 
   //============================================================================
@@ -682,7 +694,8 @@ main(int argc, char ** argv) -> int
   // Generate the mesh
   //===========================================================================
 
-  um2::gmsh::model::mesh::setGlobalMeshSize(pin_pitch / 6);
+  um2::gmsh::model::mesh::setMeshFieldFromKnudsenNumber(2, model.materials(), target_kn,
+                                                        mfp_threshold, mfp_scale);
   um2::gmsh::model::mesh::generateMesh(um2::MeshType::QuadraticTri);
   um2::gmsh::write("nuscale_2d.inp");
 
@@ -691,8 +704,7 @@ main(int argc, char ** argv) -> int
   //===========================================================================
 
   model.importCoarseCellMeshes("nuscale_2d.inp");
-  // model.writeCMFDInfo("nuscale_2d_cmfd_info.xdmf");
-  model.write("nuscale_2d.xdmf", /*write_knudsen_data=*/true);
+  model.write("nuscale_2d.xdmf");
   um2::finalize();
   return 0;
 }

post-process/CMakeLists.txt

@@ -4,3 +4,4 @@ um2_add_executable(./extract_meanchordlength.cpp)
 um2_add_executable(./extract_edge_length.cpp)
 um2_add_executable(./extract_cmfd_info.cpp)
 um2_add_executable(./extract_powers.cpp)
+um2_add_executable(./extract_spectrum.cpp)

post-process/extract_spectrum.cpp

@@ -0,0 +1,124 @@
+#include <um2.hpp>
+
+// NOLINTBEGIN(misc-include-cleaner)
+#include <um2/stdlib/algorithm/is_sorted.hpp>
+
+PURE [[nodiscard]] auto
+getSpectrum(um2::String const & filename) -> um2::Vector<Float>;
+
+template <Int P, Int N>
+PURE [[nodiscard]] auto
+getSpectrum(um2::PolytopeSoup const & soup) -> um2::Vector<Float>;
+
+auto
+main(int argc, char ** argv) -> int
+{
+  um2::initialize();
+
+  // Get the file name from the command line
+  if (argc != 2) {
+    um2::String const exec_name(argv[0]);
+    LOG_ERROR("Usage: ", exec_name, " <filename>");
+    return 1;
+  }
+
+  // Get the spectrum
+  um2::String const filename(argv[1]);
+  auto const spectrum = getSpectrum(filename);
+
+  // Write the spectrum to a file
+  FILE * file = fopen("spectrum.txt", "w");
+  if (file == nullptr) {
+    LOG_ERROR("Could not open file spectrum.txt");
+    return 1;
+  }
+
+  for (Int g = 0; g < spectrum.size(); ++g) {
+    // Write group number and flux to file
+    int const ret = fprintf(file, "%d, %.16f\n", g, spectrum[g]);
+    if (ret < 0) {
+      LOG_ERROR("Could not write to file spectrum.txt");
+      return 1;
+    }
+  }
+  int const ret = fclose(file);
+  if (ret != 0) {
+    LOG_ERROR("Could not close file spectrum.txt");
+    return 1;
+  }
+  return 0;
+}
+
+PURE [[nodiscard]] auto
+getSpectrum(um2::String const & filename) -> um2::Vector<Float>
+{
+  um2::PolytopeSoup const soup(filename);
+  // For now, we assume that all the elements are the same type.
+  auto const elem_types = soup.getElemTypes();
+  if (elem_types.size() != 1) {
+    LOG_ERROR("Expected only one element type, but found ", elem_types.size());
+    return {};
+  }
+  switch (elem_types[0]) {
+  case um2::VTKElemType::Triangle:
+    LOG_INFO("FSR mesh type: Triangle");
+    return getSpectrum<1, 3>(soup);
+  case um2::VTKElemType::Quad:
+    LOG_INFO("FSR mesh type: Quad");
+    return getSpectrum<1, 4>(soup);
+  case um2::VTKElemType::QuadraticTriangle:
+    LOG_INFO("FSR mesh type: QuadraticTriangle");
+    return getSpectrum<2, 6>(soup);
+  case um2::VTKElemType::QuadraticQuad:
+    LOG_INFO("FSR mesh type: QuadraticQuad");
+    return getSpectrum<2, 8>(soup);
+  default:
+    LOG_ERROR("Unsupported element type");
+    return {};
+  }
+}
+
+template <Int P, Int N>
+PURE [[nodiscard]] auto
+getSpectrum(um2::PolytopeSoup const & soup) -> um2::Vector<Float>
+{
+  um2::FaceVertexMesh<P, N> const fvm(soup, /*validate=*/false);
+  Int const num_faces = fvm.numFaces();
+
+  // Get the area of each face
+  um2::Vector<Float> face_areas(num_faces);
+  for (Int i = 0; i < num_faces; ++i) {
+    face_areas[i] = fvm.getFace(i).area();
+  }
+
+  um2::Vector<Float> spectrum;
+  um2::Vector<Int> ids;
+  um2::Vector<Float> flux;
+  um2::String elset_name("flux_001");
+  while (true) {
+    flux.clear();
+    soup.getElset(elset_name, ids, flux);
+    if (flux.empty()) {
+      break;
+    }
+    ASSERT(ids.size() == flux.size());
+    ASSERT(ids.size() == num_faces);
+    ASSERT(um2::is_sorted(ids.cbegin(), ids.cend()));
+    Float group_flux = 0;
+    Float local_accum = 0;
+    for (Int i = 0; i < ids.size(); ++i) {
+      local_accum += flux[i] * face_areas[i];
+      // Increment group flux every 128 faces
+      if (i % 128 == 127) {
+        group_flux += local_accum;
+        local_accum = 0;
+      }
+    }
+    group_flux += local_accum;
+    spectrum.emplace_back(group_flux);
+    um2::mpact::incrementASCIINumber(elset_name);
+  }
+  return spectrum;
+}
+
+// NOLINTEND(misc-include-cleaner)

post-process/plot_spectrum.py

@@ -0,0 +1,21 @@
+import matplotlib.pyplot as plt
+import numpy as np
+import sys
+
+# Get the file name from the command line
+if len(sys.argv) != 2:
+    print('Usage: python plot_spectrum.py <file_name>')
+    sys.exit(1)
+
+file_name = sys.argv[1]
+
+# Read file with data of the form: group, neutrons/sec
+data = np.loadtxt(file_name, delimiter=',')
+neutrons = data[:, 1]
+
+total_neutrons = np.sum(neutrons)
+plt.stairs(neutrons/total_neutrons)
+plt.xlabel('Group')
+plt.ylabel('Normalized Neutrons/sec')
+plt.title('Normalized Energy Spectrum')
+plt.show()

src/mpact/model.cpp

@@ -40,6 +40,7 @@
 #include <algorithm>
 #include <cstring>
 #include <fstream>
+#include <iomanip>
 #include <numeric>
 #include <string>
 
@@ -2044,6 +2045,126 @@ writeInputFile(String const & filename, Model const & model)
 
 } // writeInputFile
 
+void
+// NOLINTNEXTLINE(readability-function-cognitive-complexity)
+writeMaterialVolumes(String const & filename, Model const & model)
+{
+
+  // Strip the ending from filename
+  Int const last_period = filename.find_last_of('.');
+  String const base_filename = filename.substr(0, last_period);
+  String const volumes_filename = base_filename + ".volumes";
+  LOG_INFO("Writing MPACT material volumes file: ", volumes_filename);
+
+  // Open the file
+  std::ofstream file(std::string(volumes_filename.data()));
+  if (!file.is_open()) {
+    logger::error("Failed to open file: ", volumes_filename);
+    return;
+  }
+
+  auto const & core = model.core();
+  auto const & lattices = model.lattices();
+  auto const & rtms = model.rtms();
+  auto const & coarse_cells = model.coarseCells();
+  auto const & materials = model.materials();
+  auto const & tris = model.triMeshes();
+  auto const & quads = model.quadMeshes();
+  auto const & tri6s = model.tri6Meshes();
+  auto const & quad8s = model.quad8Meshes();
+
+  Int const nmats = materials.size();
+  Vector<Vector<Float>> volumes(nmats);
+  // For each assembly in the core
+  for (auto const & asy_id : core.children()) {
+    // Get the assembly and determine the number of lattices
+    auto const & assembly = model.getAssembly(asy_id);
+    auto const num_lats = assembly.children().size();
+    // For each of the lattices in the assembly
+    for (Int ilat = 0; ilat < num_lats; ++ilat) {
+      auto const lat_id = assembly.getChild(ilat);
+      // Get the dz of the lattice
+      Float const low_z = assembly.grid().divs(0)[ilat];
+      Float const high_z = assembly.grid().divs(0)[ilat + 1];
+      Float const dz = high_z - low_z;
+      // For each RTM in the lattice
+      for (auto const & rtm_id : lattices[lat_id].children()) {
+        // For each coarse cell in the RTM
+        for (auto const & cc_id : rtms[rtm_id].children()) {
+          // Get the coarse cell
+          auto const & cc = coarse_cells[cc_id];
+          // Based on the mesh type, find the area of each face
+          // NOLINTBEGIN(bugprone-signed-char-misuse,cert-str34-c)
+          switch (cc.mesh_type) {
+          case MeshType::Tri: {
+            auto const & mesh = tris[cc.mesh_id];
+            Int const num_faces = mesh.numFaces();
+            for (Int i = 0; i < num_faces; ++i) {
+              auto const mat_id = static_cast<Int>(cc.material_ids[i]);
+              volumes[mat_id].emplace_back(mesh.getFace(i).area() * dz);
+            }
+            break;
+          }
+          case MeshType::Quad: {
+            auto const & mesh = quads[cc.mesh_id];
+            Int const num_faces = mesh.numFaces();
+            for (Int i = 0; i < num_faces; ++i) {
+              auto const mat_id = static_cast<Int>(cc.material_ids[i]);
+              volumes[mat_id].emplace_back(mesh.getFace(i).area() * dz);
+            }
+            break;
+          }
+          case MeshType::QuadraticTri: {
+            auto const & mesh = tri6s[cc.mesh_id];
+            Int const num_faces = mesh.numFaces();
+            for (Int i = 0; i < num_faces; ++i) {
+              auto const mat_id = static_cast<Int>(cc.material_ids[i]);
+              volumes[mat_id].emplace_back(mesh.getFace(i).area() * dz);
+            }
+            break;
+          }
+          case MeshType::QuadraticQuad: {
+            auto const & mesh = quad8s[cc.mesh_id];
+            Int const num_faces = mesh.numFaces();
+            for (Int i = 0; i < num_faces; ++i) {
+              auto const mat_id = static_cast<Int>(cc.material_ids[i]);
+              volumes[mat_id].emplace_back(mesh.getFace(i).area() * dz);
+            }
+            break;
+          }
+          default:
+            logger::error("Unsupported mesh type");
+          }
+          // NOLINTEND(bugprone-signed-char-misuse,cert-str34-c)
+        }
+      }
+    }
+  }
+
+  // Reduce the volumes to the total volume of each material
+  Vector<Float> total_volumes(nmats, 0);
+  for (Int imat = 0; imat < nmats; ++imat) {
+    // Sort and use the sum function for minimal floating point error
+    // (sum uses pairwise summation)
+    std::sort(volumes[imat].begin(), volumes[imat].end());
+    total_volumes[imat] = um2::sum(volumes[imat].cbegin(), volumes[imat].cend());
+  }
+
+  Int p = 6;
+  if constexpr (std::same_as<Float, double>) {
+    p = 15;
+  }
+
+  for (Int imat = 0; imat < nmats; ++imat) {
+    file << materials[imat].getName().data() << " " << std::setprecision(p)
+         << total_volumes[imat] << '\n';
+  }
+
+  // Close the file
+  file.close();
+
+} // writeMaterialVolumes
+
 } // namespace
 
 //==============================================================================
@@ -2061,6 +2182,8 @@ Model::write(String const & filename, bool const write_knudsen_data,
   }
   // Write the MPACT input file
   writeInputFile(filename, *this);
+  // Write the volumes of each material
+  writeMaterialVolumes(filename, *this);
 }
 
 //==============================================================================