Implement Material

CMakeLists.txt

@@ -1,15 +1,16 @@
 cmake_minimum_required(VERSION 3.4...3.31) # specifying a minimum CMake version
 project(demeter) # set the project name
 
-set(CMAKE_CXX_STANDARD 14) # specifying the C++ Standard
+set(CMAKE_CXX_STANDARD 17) # specifying the C++ Standard
 set(GCC_COVERAGE_COMPILE_FLAGS "-Wall -Wextra -pedantic")
 set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${GCC_COVERAGE_COMPILE_FLAGS}")
 # add_compile_options(-Wall -Wextra -pedantic)
 
 set(CMAKE_CXX_FLAGS_RELEASE "-O3 -march=native -funroll-loops -funroll-all-loops -DNDEBUG")
-set(CMAKE_CXX_FLAGS_DEBUG "-g -fsanitize=address,undefined")
+set(CMAKE_CXX_FLAGS_DEBUG "-g -fsanitize=address,undefined") # -fanalizer if gcc
 
 # include_directories(include)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/src)
 
 # OpenMP
 find_package(OpenMP)
@@ -24,7 +25,7 @@ set(LINK_LIBS ${LINK_LIBS} Eigen3::Eigen)
 
 # Verbose
 include(CMakePrintHelpers)
-cmake_print_variables(PROJECT_SOURCE_DIR CMAKE_BUILD_TYPE)
+cmake_print_variables(PROJECT_SOURCE_DIR CMAKE_BUILD_TYPE CMAKE_CURRENT_SOURCE_DIR)
 cmake_print_variables(CMAKE_CXX_COMPILER CMAKE_CXX_COMPILER_ID CMAKE_CXX_COMPILER_VERSION)
 cmake_print_variables(CMAKE_CXX_FLAGS COMPILE_DEFINITIONS)
 cmake_print_variables(LINK_LIBS)
@@ -39,8 +40,6 @@ include(cmake/clang-cxx-dev-tools.cmake) # make clang-format and make clang-tidy
 add_executable(main examples/main.cpp) # create an executable using the specified file
 target_link_libraries(main PUBLIC ${LINK_LIBS})
 
-add_executable(simple examples/simple.cpp) # create an executable using the specified file
-target_link_libraries(simple PUBLIC ${LINK_LIBS})
-
-add_executable(test_omp examples/test_omp.cpp) # create an executable using the specified file
-target_link_libraries(test_omp PUBLIC ${LINK_LIBS})
+# Tests
+# enable_testing()
+# add_subdirectory(tests) # add the tests subdirectory

examples/main.cpp

@@ -1,6 +1,39 @@
 #include <iostream>
+#include <vector>
+#include <omp.h>
+#include <Eigen/Dense>
+#include "demeter/model/material.hpp"
 
 int main() {
-  std::cout << "Bonjour !" << std::endl;
+#if defined(_OPENMP)
+  std::cerr << "OpenMP version " << _OPENMP << '\n';
+#endif
+
+#pragma omp parallel for
+  for (int i = 0; i < 9; ++i) {
+    std::cerr << "rank = " << omp_get_thread_num()
+              << " size = " << omp_get_num_threads() << '\n';
+  }
+
+  using Eigen::MatrixXd;
+
+  MatrixXd m(2, 2);
+  m(0, 0) = 3;
+  m(1, 0) = 2.5;
+  m(0, 1) = -1;
+  m(1, 1) = m(1, 0) + m(0, 1);
+  std::cout << m << std::endl;
+
+  // using vec = Eigen::ArrayXd;
+  // using mat = Eigen::ArrayXXd;
+
+  // vec sigma_t = {0.222222, 0.833333};
+  // vec D = 1. / (3. * sigma_t);
+  // vec sigma_a{0.010120, 0.080032};
+  // vec nusigma_f{0., 0.135};
+  // vec sigma_f{0., 0.135};
+  // vec chi{1., 0.};
+  // mat sigma_s{{0.00, 0.02}, {0.00, 0.00}};
+
   return 0;
-}
+}

src/demeter/model/material.cpp

@@ -0,0 +1,110 @@
+#include "material.hpp"
+
+namespace Demeter {
+
+/**
+ * @brief Construct a Material from a set of cross-sections.
+ *
+ * @param[in] sigma_t total cross-section
+ * @param[in] sigma_s scattering cross-section matrix
+ * @param[in] sigma_a absorption cross-section
+ * @param[in] sigma_f fission cross-section
+ * @param[in] nu_sigma_f nu*fission
+ * @param[in] chi fission spectrum
+ * @param[in] name material name (optional)
+ *
+ *
+ * @note
+ * To use it from Python, you can do something like
+ *
+ * >>> sigma_t = np.array([0.1, 0.2, 0.3])
+ * >>> sigma_s = np.array([[0.1, 0.2, 0.3], [0.3, 0.2, 0.1], [0.1, 0.3, 0.2]])
+ * >>> sigma_a = np.array([0.01, 0.02, 0.03])
+ * >>> sigma_f = np.array([0.01, 0.02, 0.03])
+ * >>> nu_sigma_f = np.array([0.01, 0.02, 0.03])
+ * >>> chi = np.array([0.1, 0.2, 0.3])
+ * >>> uo2 = Material(sigma_t, sigma_s, sigma_a, sigma_f, nu_sigma_f, chi,
+ * 'UO2')
+ */
+Material::Material(ArrayXd& sigma_t, ArrayXXd& sigma_s, ArrayXd& sigma_a,
+                   ArrayXd& sigma_f, ArrayXd& nu_sigma_f, ArrayXd& chi,
+                   std::string_view name)
+    : sigma_t_(sigma_t),
+      sigma_s_(sigma_s),
+      sigma_a_(sigma_a),
+      sigma_f_(sigma_f),
+      nu_sigma_f_(nu_sigma_f),
+      chi_(chi),
+      name_(name) {
+  check();
+  fissile_ = sigma_f.maxCoeff() > 0.0;
+}
+
+/**
+ * @brief Construct a Material from a set of cross-sections.
+ * Constructor from rvalue references, it takes ownership of the input arrays.
+ *
+ * @param[in] sigma_t total cross-section
+ * @param[in] sigma_s scattering cross-section matrix
+ * @param[in] sigma_a absorption cross-section
+ * @param[in] sigma_f fission cross-section
+ * @param[in] nu_sigma_f nu*fission
+ * @param[in] chi fission spectrum
+ * @param[in] name material name (optional)
+ *
+ *
+ * @note
+ * To use it from Python, you can do something like
+ *
+ * >>> sigma_t = np.array([0.1, 0.2, 0.3])
+ * >>> sigma_s = np.array([[0.1, 0.2, 0.3], [0.3, 0.2, 0.1], [0.1, 0.3, 0.2]])
+ * >>> sigma_a = np.array([0.01, 0.02, 0.03])
+ * >>> sigma_f = np.array([0.01, 0.02, 0.03])
+ * >>> nu_sigma_f = np.array([0.01, 0.02, 0.03])
+ * >>> chi = np.array([0.1, 0.2, 0.3])
+ * >>> uo2 = Material(sigma_t, sigma_s, sigma_a, sigma_f, nu_sigma_f, chi,
+ * 'UO2')
+ */
+Material::Material(ArrayXd&& sigma_t, ArrayXXd&& sigma_s, ArrayXd&& sigma_a,
+                   ArrayXd&& sigma_f, ArrayXd&& nu_sigma_f, ArrayXd&& chi,
+                   std::string_view name)
+    : sigma_t_(std::move(sigma_t)),
+      sigma_s_(std::move(sigma_s)),
+      sigma_a_(std::move(sigma_a)),
+      sigma_f_(std::move(sigma_f)),
+      nu_sigma_f_(std::move(nu_sigma_f)),
+      chi_(std::move(chi)),
+      num_groups_(sigma_t.size()),
+      name_(name) {
+  check();
+  fissile_ = sigma_f_.maxCoeff() > 0.0;
+}
+
+/**
+ * @brief Move constructor
+ *
+ * @param[in] other other Material object
+ */
+Material::Material(Material&& other)
+    : sigma_t_(std::move(other.sigma_t_)),
+      sigma_s_(std::move(other.sigma_s_)),
+      sigma_a_(std::move(other.sigma_a_)),
+      sigma_f_(std::move(other.sigma_f_)),
+      nu_sigma_f_(std::move(other.nu_sigma_f_)),
+      chi_(std::move(other.chi_)),
+      num_groups_(other.num_groups_),
+      name_(other.name_),
+      fissile_(other.fissile_) {}
+
+void Material::check() const {
+  assert(num_groups_ > 0);
+  assert(sigma_t_.rows() == num_groups_);
+  assert(sigma_s_.rows() == num_groups_);
+  assert(sigma_s_.cols() == num_groups_);
+  assert(sigma_a_.size() == num_groups_);
+  assert(sigma_f_.size() == num_groups_);
+  assert(nu_sigma_f_.size() == num_groups_);
+  assert(chi_.size() == num_groups_);
+}
+
+}  // namespace Demeter

src/demeter/model/material.hpp

@@ -0,0 +1,90 @@
+#pragma once
+
+#include <vector>
+#include <string>
+#include <string_view>
+#include <Eigen/Core>
+
+// TODO add optional xs and name
+// TODO some material maybe not need all xs, make them optional
+namespace Demeter {
+
+/**
+ * @brief Describes the material properties. The cross-sections (xs) are:
+ * total, scatering, absorption and fission.
+ */
+class Material {
+  using ArrayXd = Eigen::ArrayXd;  // TODO move to a global definition header
+  using ArrayXXd = Eigen::ArrayXXd;
+
+ public:
+  Material(ArrayXd& sigma_t, ArrayXXd& sigma_s, ArrayXd& sigma_a,
+           ArrayXd& sigma_f, ArrayXd& nu_sigma_f, ArrayXd& chi,
+           std::string_view name = "");
+
+  // rvalue reference constructor
+  Material(ArrayXd&& sigma_t, ArrayXXd&& sigma_s, ArrayXd&& sigma_a,
+           ArrayXd&& sigma_f, ArrayXd&& nu_sigma_f, ArrayXd&& chi,
+           std::string_view name = "");
+
+  // move constructor
+  Material(Material&& other);
+
+  // copy constructor
+  Material(const Material& other);
+
+  auto name() const { return name_; }
+  auto NumEnergyGroups() const { return num_groups_; }
+  const auto& SigmaT() const { return sigma_t_; }
+  auto MaxSigmaT() const { return sigma_t_.maxCoeff(); };
+  const auto& SigmaS() const { return sigma_s_; }
+  const auto& SigmaA() const { return sigma_a_; }
+  const auto& SigmaF() const { return sigma_f_; }
+  const auto& NuSigmaF() const { return nu_sigma_f_; }
+  const auto& Chi() const { return chi_; }
+  auto SigmaT(size_t group) const { return sigma_t_(check(group)); }
+  auto SigmaS(size_t from, size_t to) const {
+    return sigma_s_(check(from), check(to));
+  };
+  auto SigmaA(size_t group) const { return sigma_a_(check(group)); }
+  auto SigmaF(size_t group) const { return sigma_f_(check(group)); };
+  auto NuSigmaF(size_t group) const { return nu_sigma_f_(check(group)); };
+  auto Chi(size_t group) const { return chi_(check(group)); };
+  bool fissile() const { return fissile_; }
+
+ private:
+  void check() const;
+  size_t check(size_t group) const {
+    assert((group >= 0) and (group < num_groups_));
+    return group;
+  }
+
+ private:
+  /* A name for the Material */
+  std::string name_;
+
+  /* The number of energy groups */
+  size_t num_groups_;
+
+  /* The total xs for each energy group */
+  ArrayXd sigma_t_;
+
+  /* A 2D array of the scattering cross-section matrix from/into each group */
+  ArrayXXd sigma_s_;
+
+  /* The absorption xs for each energy group */
+  ArrayXd sigma_a_;
+
+  /* The fission xs for each energy group */
+  ArrayXd sigma_f_;
+
+  /* The fission xs multiplied by nu for each energy group */
+  ArrayXd nu_sigma_f_;
+
+  /* The chi values for each energy group */
+  ArrayXd chi_;
+
+  /* The Material is fissile if it contains a non-zero fission xs */
+  bool fissile_;
+};
+}  // namespace Demeter