Continue modeling

src/model/material.hpp

@@ -0,0 +1,57 @@
+#pragma once
+
+#include <vector>
+#include <string>
+
+// TODO switch to Eigen
+// TODO add optional
+namespace Demeter {
+
+/**
+ * @brief Describes the material properties. The cross-sections (xs) are:
+ * total, scatering, absorption and fission.
+ */
+class Material {
+ public:
+  Material(std::vector<double>& sigma_t, std::vector<double>& sigma_s,
+           std::vector<double>& sigma_a, std::vector<double>& sigma_f,
+           std::vector<double>& nu_sigma_f, std::vector<double>& chi);
+  // rvalue reference constructor
+  Material(std::vector<double>&& sigma_t, std::vector<double>&& sigma_s,
+           std::vector<double>&& sigma_a, std::vector<double>&& sigma_f,
+           std::vector<double>&& nu_sigma_f, std::vector<double>&& chi);
+  // move constructor
+  Material(Material&& other);
+  // copy constructor
+  Material(const Material& other);
+
+ private:
+  /* A name for the Material */
+  std::string name_;
+
+  /* The number of energy groups */
+  int num_groups_;
+
+  /* The total xs for each energy group */
+  std::vector<double> sigma_t_;
+
+  /* A 2D (TODO) array of the scattering cross-section matrix from/into each
+   * group */
+  std::vector<double> sigma_s_;
+
+  /* The absorption xs for each energy group */
+  std::vector<double> sigma_a_;
+
+  /* The fission xs for each energy group */
+  std::vector<double> sigma_f_;
+
+  /* The fission xs multiplied by nu for each energy group */
+  std::vector<double> nu_sigma_f_;
+
+  /* The chi values for each energy group */
+  std::vector<double> chi_;
+
+  /* The Material is fissile if it contains a non-zero fission xs */
+  bool fissile_;
+};
+}  // namespace Demeter

src/core/spec.md → src/model/spec.md

@@ -7,8 +7,8 @@ This also simplifies the specification of the geometry sub-module.
 
 ```python
 # Initialize material cross sections using NumPy data arrays
-sigma_t = np.array([0.1,0.15,0.2,0.25,0.35,0.4,0.45,0.5])
-sigma_f = np.array([0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4])
+sigma_t = np.array([0.45,0.5])
+sigma_f = np.array([0.35,0.4])
 
 # Instantiate an OpenMOC Material class object with an optional string name
 material = Material(sigma_t, sigma_f, name='test material')
@@ -23,12 +23,11 @@ mox43.setMaterial(material)
 Lets the user describe the geometry of the problem.
 
 Cell
-PinCell(dx (sca), dy (sca), radii, mats, pitch, pitch)
-Lattice(dx (vec), dy (vec))
-  symetry 1/2 1/8
+PinCell
+Assembly
+Core
 
-Point
-Edge
+Symetry 1/2 1/8
 
 ## Material and Cross section
 
@@ -40,10 +39,3 @@ Material : store the cross section of a material
   vec nu*fission v*sigma_f
   sca chi
   str name (optional)
-
-## Mesh
-
-Utilisee par le solveur pour les calculs. La geometrie est convertie en maillage.
-
-Point
-Edge

tests/regress/test00.py

@@ -0,0 +1,77 @@
+from demeter import *
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Create materials
+sigma_t = np.array([0.1, 0.2, 0.3])
+sigma_f = np.array([0.01, 0.02, 0.03]),
+sigma_a = np.array([0.001, 0.002, 0.003])
+uO2 = Material(sigma_t, sigma_f, sigma_a, name= 'uO2')
+
+sigma_t = np.array([0.1, 0.2, 0.3])
+sigma_f = np.array([0.01, 0.02, 0.03]),
+sigma_a = np.array([0.001, 0.002, 0.003])
+h20 = Material(sigma_t, sigma_f, sigma_a, name= 'H20')
+
+# Create pin cells
+uO2_cell = PinCell(0.5, uO2)
+h2O_cell = PinCell(0.5, h20)
+
+# Create assembly
+dx = np.array([10., 20.])
+dy = np.array([10., 20.])
+motif = [uO2_cell, uO2_cell,
+        h2O_cell, h2O_cell]
+assembly = Lattice(motif, dx, dy)
+
+# Create core
+core = Lattice([assembly]) # => deduction des tailles
+# core.refineMesh(2)
+# core.getXMin().setBoundaryCondition(BoundaryCondition.Vacuum)
+# core.getXMax().setBoundaryCondition(BoundaryCondition.Reflection)
+
+# Maybe an intermediate geometry class that handles core or assembly
+
+# Create solver
+solver = Solver(core, 4, 10)
+
+# Set boundary conditions
+# solver.setBoundaryCondition(BoundaryCondition.Vacuum, BoundarySide.Left)
+# solver.setBoundaryCondition(BoundaryCondition.Reflection, BoundarySide.Right)
+
+# Solve the transport equation
+solver.solve()
+
+# Access data
+print(f"Core has {len(core.getAssemblies())} assemblies")
+print(f"Assembly has {len(core.getAssemblies()[0].getPinCells())} pin cells")
+print(f"Pin cell radius: {core.getAssemblies()[0].getPinCells()[0].getRadius()}")
+print(f"Material total cross section: {core.getAssemblies()[0].getPinCells()[0].getMaterial().getCrossSection().getTotal()}")
+
+# Get flux data
+flux = solver.getFlux()
+
+# Plot geometry and flux
+fig, ax = plt.subplots()
+
+# Plot pin cells
+for i, assembly in enumerate(core.getAssemblies()):
+  for j, pin_cell in enumerate(assembly.getPinCells()):
+    circle = plt.Circle((i, j), pin_cell.getRadius(), edgecolor='black', facecolor='none')
+    ax.add_artist(circle)
+
+# Plot flux
+flux_array = np.array(flux)
+cax = ax.matshow(flux_array, cmap='viridis')
+
+# Add colorbar
+fig.colorbar(cax)
+
+# Set plot limits and labels
+ax.set_xlim(-1, len(core.getAssemblies()))
+ax.set_ylim(-1, len(core.getAssemblies()[0].getPinCells()))
+ax.set_xlabel('Assembly Index')
+ax.set_ylabel('Pin Cell Index')
+ax.set_title('Flux Distribution')
+
+plt.show()