Leaky and Albedo Boundary Conditions Implementation (#2724)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

src/boundary_condition.cpp

@@ -28,6 +28,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
   Direction u = surf.reflect(p.r(), p.u(), &p);
   u /= u.norm();
 
+  // Handle the effects of the surface albedo on the particle's weight.
+  BoundaryCondition::handle_albedo(p, surf);
+
   p.cross_reflective_bc(surf, u);
 }
 
@@ -40,6 +43,9 @@ void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
   Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
   u /= u.norm();
 
+  // Handle the effects of the surface albedo on the particle's weight.
+  BoundaryCondition::handle_albedo(p, surf);
+
   p.cross_reflective_bc(surf, u);
 }
 
@@ -130,6 +136,9 @@ void TranslationalPeriodicBC::handle_particle(
       "hitting a surface, but that surface is not recognized by the BC.");
   }
 
+  // Handle the effects of the surface albedo on the particle's weight.
+  BoundaryCondition::handle_albedo(p, surf);
+
   // Pass the new location and surface to the particle.
   p.cross_periodic_bc(surf, new_r, p.u(), new_surface);
 }
@@ -264,6 +273,9 @@ void RotationalPeriodicBC::handle_particle(
   Direction new_u = {
     cos_theta * u.x - sin_theta * u.y, sin_theta * u.x + cos_theta * u.y, u.z};
 
+  // Handle the effects of the surface albedo on the particle's weight.
+  BoundaryCondition::handle_albedo(p, surf);
+
   // Pass the new location, direction, and surface to the particle.
   p.cross_periodic_bc(surf, new_r, new_u, new_surface);
 }

src/dagmc.cpp

@@ -289,10 +289,10 @@ void DAGUniverse::init_geometry()
         bc_value == "transmission") {
       // set to transmission by default (nullptr)
     } else if (bc_value == "vacuum") {
-      s->bc_ = std::make_shared<VacuumBC>();
+      s->bc_ = make_unique<VacuumBC>();
     } else if (bc_value == "reflective" || bc_value == "reflect" ||
                bc_value == "reflecting") {
-      s->bc_ = std::make_shared<ReflectiveBC>();
+      s->bc_ = make_unique<ReflectiveBC>();
     } else if (bc_value == "periodic") {
       fatal_error("Periodic boundary condition not supported in DAGMC.");
     } else {
@@ -310,7 +310,7 @@ void DAGUniverse::init_geometry()
     // if this surface belongs to the graveyard
     if (graveyard && parent_vols.find(graveyard) != parent_vols.end()) {
       // set graveyard surface BC's to vacuum
-      s->bc_ = std::make_shared<VacuumBC>();
+      s->bc_ = make_unique<VacuumBC>();
     }
 
     // add to global array and map

src/surface.cpp

@@ -80,19 +80,36 @@ Surface::Surface(pugi::xml_node surf_node)
     if (surf_bc == "transmission" || surf_bc == "transmit" || surf_bc.empty()) {
       // Leave the bc_ a nullptr
     } else if (surf_bc == "vacuum") {
-      bc_ = std::make_shared<VacuumBC>();
+      bc_ = make_unique<VacuumBC>();
     } else if (surf_bc == "reflective" || surf_bc == "reflect" ||
                surf_bc == "reflecting") {
-      bc_ = std::make_shared<ReflectiveBC>();
+      bc_ = make_unique<ReflectiveBC>();
     } else if (surf_bc == "white") {
-      bc_ = std::make_shared<WhiteBC>();
+      bc_ = make_unique<WhiteBC>();
     } else if (surf_bc == "periodic") {
-      // periodic BC's are handled separately
+      // Periodic BCs are handled separately
     } else {
       fatal_error(fmt::format("Unknown boundary condition \"{}\" specified "
                               "on surface {}",
         surf_bc, id_));
     }
+
+    if (check_for_node(surf_node, "albedo") && bc_) {
+      double surf_alb = std::stod(get_node_value(surf_node, "albedo"));
+
+      if (surf_alb < 0.0)
+        fatal_error(fmt::format("Surface {} has an albedo of {}. "
+                                "Albedo values must be positive.",
+          id_, surf_alb));
+
+      if (surf_alb > 1.0)
+        warning(fmt::format("Surface {} has an albedo of {}. "
+                            "Albedos greater than 1 may cause "
+                            "unphysical behaviour.",
+          id_, surf_alb));
+
+      bc_->set_albedo(surf_alb);
+    }
   }
 }
 
@@ -154,6 +171,7 @@ void Surface::to_hdf5(hid_t group_id) const
 
     if (bc_) {
       write_string(surf_group, "boundary_type", bc_->type(), false);
+      bc_->to_hdf5(surf_group);
     } else {
       write_string(surf_group, "boundary_type", "transmission", false);
     }
@@ -1156,9 +1174,10 @@ void read_surfaces(pugi::xml_node node)
   }
 
   // Loop over XML surface elements and populate the array.  Keep track of
-  // periodic surfaces.
+  // periodic surfaces and their albedos.
   model::surfaces.reserve(n_surfaces);
   std::set<std::pair<int, int>> periodic_pairs;
+  std::unordered_map<int, double> albedo_map;
   {
     pugi::xml_node surf_node;
     int i_surf;
@@ -1221,6 +1240,12 @@ void read_surfaces(pugi::xml_node node)
       if (check_for_node(surf_node, "boundary")) {
         std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
         if (surf_bc == "periodic") {
+          // Check for surface albedo. Skip sanity check as it is already done
+          // in the Surface class's constructor.
+          if (check_for_node(surf_node, "albedo")) {
+            albedo_map[model::surfaces.back()->id_] =
+              std::stod(get_node_value(surf_node, "albedo"));
+          }
           if (check_for_node(surf_node, "periodic_surface_id")) {
             int i_periodic =
               std::stoi(get_node_value(surf_node, "periodic_surface_id"));
@@ -1284,7 +1309,7 @@ void read_surfaces(pugi::xml_node node)
     periodic_pairs.erase(second_unresolved);
   }
 
-  // Assign the periodic boundary conditions
+  // Assign the periodic boundary conditions with albedos
   for (auto periodic_pair : periodic_pairs) {
     int i_surf = model::surface_map[periodic_pair.first];
     int j_surf = model::surface_map[periodic_pair.second];
@@ -1302,11 +1327,19 @@ void read_surfaces(pugi::xml_node node)
     // planes are parallel which indicates a translational periodic boundary
     // condition.  Otherwise, it is a rotational periodic BC.
     if (std::abs(1.0 - dot_prod) < FP_PRECISION) {
-      surf1.bc_ = std::make_shared<TranslationalPeriodicBC>(i_surf, j_surf);
-      surf2.bc_ = surf1.bc_;
+      surf1.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
+      surf2.bc_ = make_unique<TranslationalPeriodicBC>(i_surf, j_surf);
     } else {
-      surf1.bc_ = std::make_shared<RotationalPeriodicBC>(i_surf, j_surf);
-      surf2.bc_ = surf1.bc_;
+      surf1.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
+      surf2.bc_ = make_unique<RotationalPeriodicBC>(i_surf, j_surf);
+    }
+
+    // If albedo data is present in albedo map, set the boundary albedo.
+    if (albedo_map.count(surf1.id_)) {
+      surf1.bc_->set_albedo(albedo_map[surf1.id_]);
+    }
+    if (albedo_map.count(surf2.id_)) {
+      surf2.bc_->set_albedo(albedo_map[surf2.id_]);
     }
   }
 }