#include "openmc/mesh.h"

#include <algorithm> // for copy, equal, min, min_element

#include <cassert>

#include <cmath> // for ceil

#include <cstddef> // for size_t

#include <cstdint> // for uint64_t

#include <cstring> // for memcpy

#include <limits>

#include <numeric> // for accumulate

#include <string>

#ifdef _MSC_VER

#include <intrin.h> // for _InterlockedCompareExchange

#endif

#ifdef OPENMC_MPI

#include "mpi.h"

#endif

#include "openmc/tensor.h"

#include <fmt/core.h> // for fmt

#include "openmc/capi.h"

#include "openmc/constants.h"

#include "openmc/container_util.h"

#include "openmc/error.h"

#include "openmc/file_utils.h"

#include "openmc/geometry.h"

#include "openmc/hdf5_interface.h"

#include "openmc/material.h"

#include "openmc/memory.h"

#include "openmc/message_passing.h"

#include "openmc/openmp_interface.h"

#include "openmc/output.h"

#include "openmc/particle_data.h"

#include "openmc/plot.h"

#include "openmc/random_dist.h"

#include "openmc/search.h"

#include "openmc/settings.h"

#include "openmc/string_utils.h"

#include "openmc/tallies/filter.h"

#include "openmc/tallies/tally.h"

#include "openmc/timer.h"

#include "openmc/volume_calc.h"

#include "openmc/xml_interface.h"

#ifdef OPENMC_LIBMESH_ENABLED

#include "libmesh/mesh_modification.h"

#include "libmesh/mesh_tools.h"

#include "libmesh/numeric_vector.h"

#include "libmesh/replicated_mesh.h"

#endif

#ifdef OPENMC_DAGMC_ENABLED

#include "moab/FileOptions.hpp"

#endif

#include "openmc/xdg.h"

namespace openmc {

//==============================================================================

// Global variables

//==============================================================================

// Value used to indicate an empty slot in the hash table. We use -2 because

// the value -1 is used to indicate a void material.

constexpr int32_t EMPTY = -2;

namespace model {

std::unordered_map<int32_t, int32_t> mesh_map;

vector<unique_ptr<Mesh>> meshes;

} // namespace model

#ifdef OPENMC_LIBMESH_ENABLED

namespace settings {

unique_ptr<libMesh::LibMeshInit> libmesh_init;

const libMesh::Parallel::Communicator* libmesh_comm {nullptr};

} // namespace settings

#endif

//==============================================================================

// Helper functions

//==============================================================================

//! Update an intersection point if the given candidate is closer.

//

//! The first 6 arguments are coordinates for the starting point of a particle

//! and its intersection with a mesh surface. If the distance between these

//! two points is shorter than the given `min_distance`, then the `r` argument

//! will be updated to match the intersection point, and `min_distance` will

//! also be updated.

inline bool check_intersection_point(double x1, double x0, double y1, double y0,

double z1, double z0, Position& r, double& min_distance)

{

double dist =

std::pow(x1 - x0, 2) + std::pow(y1 - y0, 2) + std::pow(z1 - z0, 2);

if (dist < min_distance) {

r.x = x1;

r.y = y1;

r.z = z1;

min_distance = dist;

return true;

}

return false;

}

//! Atomic compare-and-swap for signed 32-bit integer

//

//! \param[in,out] ptr Pointer to value to update

//! \param[in,out] expected Value to compare to

//! \param[in] desired If comparison is successful, value to update to

//! \return True if the comparison was successful and the value was updated

inline bool atomic_cas_int32(int32_t* ptr, int32_t& expected, int32_t desired)

{

#if defined(__GNUC__) || defined(__clang__)

// For gcc/clang, use the __atomic_compare_exchange_n intrinsic

return __atomic_compare_exchange_n(

ptr, &expected, desired, false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST);

#elif defined(_MSC_VER)

// For MSVC, use the _InterlockedCompareExchange intrinsic

int32_t old_val =

_InterlockedCompareExchange(reinterpret_cast<volatile long*>(ptr),

static_cast<long>(desired), static_cast<long>(expected));

return (old_val == expected);

#else

#error "No compare-and-swap implementation available for this compiler."

#endif

}

// Helper function equivalent to std::bit_cast in C++20

template<typename To, typename From>

inline To bit_cast_value(const From& value)

{

To out;

std::memcpy(&out, &value, sizeof(To));

return out;

}

inline void atomic_update_double(double* ptr, double value, bool is_min)

{

#if defined(__GNUC__) || defined(__clang__)

using may_alias_uint64_t [[gnu::may_alias]] = uint64_t;

auto* bits_ptr = reinterpret_cast<may_alias_uint64_t*>(ptr);

uint64_t current_bits = __atomic_load_n(bits_ptr, __ATOMIC_SEQ_CST);

double current = bit_cast_value<double>(current_bits);

while (is_min ? (value < current) : (value > current)) {

uint64_t desired_bits = bit_cast_value<uint64_t>(value);

uint64_t expected_bits = current_bits;

if (__atomic_compare_exchange_n(bits_ptr, &expected_bits, desired_bits,

false, __ATOMIC_SEQ_CST, __ATOMIC_SEQ_CST)) {

return;

}

current_bits = expected_bits;

current = bit_cast_value<double>(current_bits);

}

#elif defined(_MSC_VER)

auto* bits_ptr = reinterpret_cast<volatile long long*>(ptr);

long long current_bits = *bits_ptr;

double current = bit_cast_value<double>(current_bits);

while (is_min ? (value < current) : (value > current)) {

long long desired_bits = bit_cast_value<long long>(value);

long long old_bits =

_InterlockedCompareExchange64(bits_ptr, desired_bits, current_bits);

if (old_bits == current_bits) {

return;

}

current_bits = old_bits;

current = bit_cast_value<double>(current_bits);

}

#else

#error "No compare-and-swap implementation available for this compiler."

#endif

}

inline void atomic_max_double(double* ptr, double value)

{

atomic_update_double(ptr, value, false);

}

inline void atomic_min_double(double* ptr, double value)

{

atomic_update_double(ptr, value, true);

}

namespace detail {

//==============================================================================

// MaterialVolumes implementation

//==============================================================================

void MaterialVolumes::add_volume(

int index_elem, int index_material, double volume, const BoundingBox* bbox)

{

// This method handles adding elements to the materials hash table,

// implementing open addressing with linear probing. Consistency across

// multiple threads is handled by with an atomic compare-and-swap operation.

// Ideally, we would use #pragma omp atomic compare, but it was introduced in

// OpenMP 5.1 and is not widely supported yet.

// Loop for linear probing

for (int attempt = 0; attempt < table_size_; ++attempt) {

// Determine slot to check, making sure it is positive

int slot = (index_material + attempt) % table_size_;

if (slot < 0)

slot += table_size_;

int32_t* slot_ptr = &this->materials(index_elem, slot);

// Non-atomic read of current material

int32_t current_val = *slot_ptr;

// Found the desired material; accumulate volume and bbox

if (current_val == index_material) {

#pragma omp atomic

this->volumes(index_elem, slot) += volume;

if (bbox) {

atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);

atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);

atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);

atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);

atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);

atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);

}

return;

}

// Slot appears to be empty; attempt to claim

if (current_val == EMPTY) {

// Attempt compare-and-swap from EMPTY to index_material

int32_t expected_val = EMPTY;

bool claimed_slot =

atomic_cas_int32(slot_ptr, expected_val, index_material);

// If we claimed the slot or another thread claimed it but the same

// material was inserted, proceed to accumulate

if (claimed_slot || (expected_val == index_material)) {

#pragma omp atomic

this->volumes(index_elem, slot) += volume;

if (bbox) {

atomic_min_double(&this->bboxes(index_elem, slot, 0), bbox->min.x);

atomic_min_double(&this->bboxes(index_elem, slot, 1), bbox->min.y);

atomic_min_double(&this->bboxes(index_elem, slot, 2), bbox->min.z);

atomic_max_double(&this->bboxes(index_elem, slot, 3), bbox->max.x);

atomic_max_double(&this->bboxes(index_elem, slot, 4), bbox->max.y);

atomic_max_double(&this->bboxes(index_elem, slot, 5), bbox->max.z);

}

return;

}

}

}

// If table is full, set a flag that can be checked later

table_full_ = true;

}

void MaterialVolumes::add_volume_unsafe(

int index_elem, int index_material, double volume, const BoundingBox* bbox)

{

// Linear probe

for (int attempt = 0; attempt < table_size_; ++attempt) {

// Determine slot to check, making sure it is positive

int slot = (index_material + attempt) % table_size_;

if (slot < 0)

slot += table_size_;

// Read current material

int32_t current_val = this->materials(index_elem, slot);

// Found the desired material; accumulate volume and bbox

if (current_val == index_material) {

this->volumes(index_elem, slot) += volume;

if (bbox) {

this->bboxes(index_elem, slot, 0) =

std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);

this->bboxes(index_elem, slot, 1) =

std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);

this->bboxes(index_elem, slot, 2) =

std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);

this->bboxes(index_elem, slot, 3) =

std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);

this->bboxes(index_elem, slot, 4) =

std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);

this->bboxes(index_elem, slot, 5) =

std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);

}

return;

}

// Claim empty slot

if (current_val == EMPTY) {

this->materials(index_elem, slot) = index_material;

this->volumes(index_elem, slot) += volume;

if (bbox) {

this->bboxes(index_elem, slot, 0) =

std::min(this->bboxes(index_elem, slot, 0), bbox->min.x);

this->bboxes(index_elem, slot, 1) =

std::min(this->bboxes(index_elem, slot, 1), bbox->min.y);

this->bboxes(index_elem, slot, 2) =

std::min(this->bboxes(index_elem, slot, 2), bbox->min.z);

this->bboxes(index_elem, slot, 3) =

std::max(this->bboxes(index_elem, slot, 3), bbox->max.x);

this->bboxes(index_elem, slot, 4) =

std::max(this->bboxes(index_elem, slot, 4), bbox->max.y);

this->bboxes(index_elem, slot, 5) =

std::max(this->bboxes(index_elem, slot, 5), bbox->max.z);

}

return;

}

}

// If table is full, set a flag that can be checked later

table_full_ = true;

}

} // namespace detail

//==============================================================================

// Mesh implementation

//==============================================================================

std::string read_mesh_type(const pugi::xml_node& node)

{

std::string mesh_type;

if (check_for_node(node, "type")) {

mesh_type = get_node_value(node, "type");

} else {

mesh_type = "regular";

}

return std::move(mesh_type);

}

std::string read_mesh_type(hid_t group)

{

std::string mesh_type;

if (object_exists(group, "type")) {

read_dataset(group, "type", mesh_type);

} else {

mesh_type = "regular";

}

return std::move(mesh_type);

}

template<typename T>

const std::unique_ptr<Mesh>& Mesh::create(T dataset)

{

std::string mesh_type = read_mesh_type(dataset);

// Determine mesh type. Add to model vector and map

if (mesh_type == RegularMesh::mesh_type) {

model::meshes.push_back(make_unique<RegularMesh>(dataset));

} else if (mesh_type == RectilinearMesh::mesh_type) {

model::meshes.push_back(make_unique<RectilinearMesh>(dataset));

} else if (mesh_type == CylindricalMesh::mesh_type) {

model::meshes.push_back(make_unique<CylindricalMesh>(dataset));

} else if (mesh_type == SphericalMesh::mesh_type) {

model::meshes.push_back(make_unique<SphericalMesh>(dataset));

} else if (mesh_type == UnstructuredMesh::mesh_type) {

model::meshes.push_back(UnstructuredMesh::create(dataset));

} else {

fatal_error(fmt::format("Invalid mesh type: {}", mesh_type));

}

// Map ID to position in vector

model::mesh_map[model::meshes.back()->id_] = model::meshes.size() - 1;

return model::meshes.back();

}

Mesh::Mesh(pugi::xml_node node)

{

// Read mesh id

id_ = std::stoi(get_node_value(node, "id"));

if (check_for_node(node, "name"))

name_ = get_node_value(node, "name");

}

Mesh::Mesh(hid_t group)

{

// Read mesh ID

read_attribute(group, "id", id_);

// Read mesh name

if (object_exists(group, "name")) {

read_dataset(group, "name", name_);

}

}

void Mesh::set_id(int32_t id)

{

assert(id >= 0 || id == C_NONE);

// Clear entry in mesh map in case one was already assigned

if (id_ != C_NONE) {

model::mesh_map.erase(id_);

id_ = C_NONE;

}

// Ensure no other mesh has the same ID

if (model::mesh_map.find(id) != model::mesh_map.end()) {

throw std::runtime_error {

fmt::format("Two meshes have the same ID: {}", id)};

}

// If no ID is specified, auto-assign the next ID in the sequence

if (id == C_NONE) {

id = 0;

for (const auto& m : model::meshes) {

id = std::max(id, m->id_);

}

++id;

}

// Update ID and entry in the mesh map

id_ = id;

// find the index of this mesh in the model::meshes vector

// (search in reverse because this mesh was likely just added to the vector)

auto it = std::find_if(model::meshes.rbegin(), model::meshes.rend(),

[this](const std::unique_ptr<Mesh>& mesh) { return mesh.get() == this; });

model::mesh_map[id] = std::distance(model::meshes.begin(), it.base()) - 1;

}

vector<double> Mesh::volumes() const

{

vector<double> volumes(n_bins());

for (int i = 0; i < n_bins(); i++) {

volumes[i] = this->volume(i);

}

return volumes;

}

//! Default (Cartesian) axis labels used for surface bin labels.

std::array<const char*, 3> Mesh::axis_labels() const

{

return {"x", "y", "z"};

}

//! Build the surface component of a mesh surface tally bin label.

//! surf_index: 0=out/min, 1=in/min, 2=out/max, 3=in/max

std::string Mesh::surface_bin_label(int surf_index) const

{

auto labels = this->axis_labels();

int dim = surf_index / 4;

int code = surf_index % 4;

bool incoming = (code == 1) || (code == 3);

bool max = (code == 2) || (code == 3);

return fmt::format(" {}, {}-{}", incoming ? "Incoming" : "Outgoing",

labels[dim], max ? "max" : "min");

}

void Mesh::material_volumes(int nx, int ny, int nz, int table_size,

int32_t* materials, double* volumes) const

{

this->material_volumes(nx, ny, nz, table_size, materials, volumes, nullptr);

}

void Mesh::material_volumes(int nx, int ny, int nz, int table_size,

int32_t* materials, double* volumes, double* bboxes) const

{

if (mpi::master) {

header("MESH MATERIAL VOLUMES CALCULATION", 7);

}

write_message(7, "Number of mesh elements = {}", n_bins());

write_message(7, "Number of rays (x) = {}", nx);

write_message(7, "Number of rays (y) = {}", ny);

write_message(7, "Number of rays (z) = {}", nz);

int64_t n_total = static_cast<int64_t>(nx) * ny +

static_cast<int64_t>(ny) * nz +

static_cast<int64_t>(nx) * nz;

write_message(7, "Total number of rays = {}", n_total);

write_message(7, "Table size per mesh element = {}", table_size);

Timer timer;

timer.start();

// Create object for keeping track of materials/volumes

detail::MaterialVolumes result(materials, volumes, bboxes, table_size);

bool compute_bboxes = bboxes != nullptr;

// Determine bounding box

auto bbox = this->bounding_box();

std::array<int, 3> n_rays = {nx, ny, nz};

// Determine effective width of rays

Position width = bbox.max - bbox.min;

width.x = (nx > 0) ? width.x / nx : 0.0;

width.y = (ny > 0) ? width.y / ny : 0.0;

width.z = (nz > 0) ? width.z / nz : 0.0;

#pragma omp parallel

{

// Preallocate vector for mesh indices and length fractions and particle

vector<int> bins;

vector<double> length_fractions;

Particle p;

SourceSite site;

site.E = 1.0;

site.particle = ParticleType::neutron();

bool verbose = settings::verbosity >= 10;

// Save the cells occupied immediately before a boundary crossing.

auto save_cell_state = [&p]() {

for (int j = 0; j < p.n_coord(); ++j) {

p.cell_last(j) = p.coord(j).cell();

}

p.n_coord_last() = p.n_coord();

};

// Initialize cell history after locating a ray inside the model.

auto initialize_cell_state = [&p, &save_cell_state]() {

if (p.cell_born() == C_NONE)

p.cell_born() = p.lowest_coord().cell();

save_cell_state();

};

// Reset a failed coordinate search while preserving position and direction.

auto reset_geometry_state = [&p]() {

Position r = p.r();

Direction u = p.u();

p.init_from_r_u(r, u);

p.coord(0).universe() = model::root_universe;

};

for (int axis = 0; axis < 3; ++axis) {

// Set starting position and direction

site.r = {0.0, 0.0, 0.0};

site.r[axis] = bbox.min[axis];

site.u = {0.0, 0.0, 0.0};

site.u[axis] = 1.0;

// Determine width of rays and number of rays in other directions

int ax1 = (axis + 1) % 3;

int ax2 = (axis + 2) % 3;

double min1 = bbox.min[ax1];

double min2 = bbox.min[ax2];

double d1 = width[ax1];

double d2 = width[ax2];

int n1 = n_rays[ax1];

int n2 = n_rays[ax2];

if (n1 == 0 || n2 == 0) {

continue;

}

// Divide rays in first direction over MPI processes by computing starting

// and ending indices

int min_work = n1 / mpi::n_procs;

int remainder = n1 % mpi::n_procs;

int n1_local = (mpi::rank < remainder) ? min_work + 1 : min_work;

int i1_start = mpi::rank * min_work + std::min(mpi::rank, remainder);

int i1_end = i1_start + n1_local;

// Add the contribution from a ray segment. The positions used here are

// kept separate from the particle position because the latter is moved a

// tiny distance across each surface for robust geometry searches.

auto add_segment = [&](const Position& r0, const Position& r1,

int i_material) {

double distance = r1[axis] - r0[axis];

if (distance <= 0.0)

return;

bins.clear();

length_fractions.clear();

this->bins_crossed(r0, r1, site.u, bins, length_fractions);

double cumulative_frac = 0.0;

for (int i_bin = 0; i_bin < bins.size(); i_bin++) {

int mesh_index = bins[i_bin];

double length = distance * length_fractions[i_bin];

double volume = length * d1 * d2;

if (compute_bboxes) {

double axis_start = r0[axis] + distance * cumulative_frac;

double axis_end = axis_start + length;

cumulative_frac += length_fractions[i_bin];

Position contrib_min = site.r;

Position contrib_max = site.r;

contrib_min[ax1] = site.r[ax1] - 0.5 * d1;

contrib_max[ax1] = site.r[ax1] + 0.5 * d1;

contrib_min[ax2] = site.r[ax2] - 0.5 * d2;

contrib_max[ax2] = site.r[ax2] + 0.5 * d2;

contrib_min[axis] = std::min(axis_start, axis_end);

contrib_max[axis] = std::max(axis_start, axis_end);

BoundingBox contrib_bbox {contrib_min, contrib_max};

contrib_bbox &= bbox;

result.add_volume(mesh_index, i_material, volume, &contrib_bbox);

} else {

result.add_volume(mesh_index, i_material, volume);

}

}

};

// Loop over rays on face of bounding box

#pragma omp for collapse(2)

for (int i1 = i1_start; i1 < i1_end; ++i1) {

for (int i2 = 0; i2 < n2; ++i2) {

site.r[ax1] = min1 + (i1 + 0.5) * d1;

site.r[ax2] = min2 + (i2 + 0.5) * d2;

p.from_source(&site);

// Set the physical endpoint of this ray at the far mesh face.

Position r_mesh_end = site.r;

r_mesh_end[axis] = bbox.max[axis];

// Determine particle's location

bool inside_model = exhaustive_find_cell(p, verbose);

if (inside_model) {

initialize_cell_state();

} else {

// Clear any partial descent into nested universes before searching

// for the first root-universe boundary from undefined space.

reset_geometry_state();

}

// Physical position through which volume has been accumulated. This

// differs by TINY_BIT from p.r() after crossing a surface.

Position r_scored = site.r;

while (r_scored[axis] < r_mesh_end[axis]) {

if (!inside_model) {

// The ray is outside the model. Advance to the next surface of

// any cell in the root universe, as is done for ray-traced

// plots. Undefined space traversed along the way is void.

Position r0 = p.r();

p.advance_to_boundary_from_void();

// If no model surface lies before the mesh edge, score the

// remaining exterior interval as void and finish the ray.

double distance_to_mesh_end = r_mesh_end[axis] - r0[axis];

if (p.boundary().surface() == SURFACE_NONE ||

p.boundary().distance() >= distance_to_mesh_end) {

add_segment(r_scored, r_mesh_end, MATERIAL_VOID);

break;

}

// Determine the physical position of the model boundary.

Position r_boundary = r0 + p.boundary().distance() * p.u();

// Score the exterior interval and record its physical endpoint.

add_segment(r_scored, r_boundary, MATERIAL_VOID);

r_scored = r_boundary;

// Check whether advancing through the surface entered the model.

inside_model = exhaustive_find_cell(p, verbose);

if (inside_model) {

initialize_cell_state();

} else {

// Clear any partial coordinate search before looking for the

// next surface from undefined space.

reset_geometry_state();

}

continue;

}

// Find the distance to the nearest boundary

BoundaryInfo boundary = distance_to_boundary(p);

// Convert the material index to a user-facing ID

int i_material = p.material();

if (i_material != C_NONE) {

i_material = model::materials[i_material]->id();

}

// If no model boundary lies before the mesh edge, score the

// remaining material interval and finish the ray.

double distance_to_mesh_end = r_mesh_end[axis] - p.r()[axis];

if (boundary.distance() >= distance_to_mesh_end) {

add_segment(r_scored, r_mesh_end, i_material);

break;

}

// Determine the physical position of the model boundary.

Position r_boundary = p.r() + boundary.distance() * p.u();

// Score the material interval and record its physical endpoint.

add_segment(r_scored, r_boundary, i_material);

r_scored = r_boundary;

// Cross the next geometric surface. The small forward movement

// and neighbor-list search mirror Ray::trace, allowing a failed

// search to mean that the ray has left the model rather than that

// a transport particle has been lost.

save_cell_state();

// Move just beyond the surface to make the next search robust.

p.move_distance(boundary.distance() + TINY_BIT);

// Set surface that particle is on and adjust coordinate levels

p.surface() = boundary.surface();

p.n_coord() = boundary.coord_level();

// Update the geometry state according to the boundary type.

if (boundary.lattice_translation()[0] != 0 ||

boundary.lattice_translation()[1] != 0 ||

boundary.lattice_translation()[2] != 0) {

// Particle crosses lattice boundary

cross_lattice(p, boundary, verbose);

inside_model = true;

} else {

// Search for the cell on the opposite side of a surface.

inside_model = neighbor_list_find_cell(p, verbose);

}

// Treat a failed cell search as a transition to exterior void.

if (!inside_model) {

// Reset the geometry state so the next iteration can search for

// another disjoint portion of the model.

reset_geometry_state();

}

}

}

}

}

}

// Check for errors

if (result.table_full()) {

throw std::runtime_error("Maximum number of materials for mesh material "

"volume calculation insufficient.");

}

// Compute time for raytracing

double t_raytrace = timer.elapsed();

#ifdef OPENMC_MPI

// Combine results from multiple MPI processes

if (mpi::n_procs > 1) {

int total = this->n_bins() * table_size;

int total_bbox = total * 6;

if (mpi::master) {

// Allocate temporary buffer for receiving data

vector<int32_t> mats(total);

vector<double> vols(total);

vector<double> recv_bboxes;

if (compute_bboxes) {

recv_bboxes.resize(total_bbox);

}

for (int i = 1; i < mpi::n_procs; ++i) {

// Receive material indices and volumes from process i

MPI_Recv(mats.data(), total, MPI_INT32_T, i, i, mpi::intracomm,

MPI_STATUS_IGNORE);

MPI_Recv(vols.data(), total, MPI_DOUBLE, i, i, mpi::intracomm,

MPI_STATUS_IGNORE);

if (compute_bboxes) {

MPI_Recv(recv_bboxes.data(), total_bbox, MPI_DOUBLE, i, i,

mpi::intracomm, MPI_STATUS_IGNORE);

}

// Combine with existing results; we can call thread unsafe version of

// add_volume because each thread is operating on a different element

#pragma omp for

for (int index_elem = 0; index_elem < n_bins(); ++index_elem) {

for (int k = 0; k < table_size; ++k) {

int index = index_elem * table_size + k;

if (mats[index] != EMPTY) {

if (compute_bboxes) {

int bbox_index = index * 6;

BoundingBox slot_bbox {

{recv_bboxes[bbox_index + 0], recv_bboxes[bbox_index + 1],

recv_bboxes[bbox_index + 2]},

{recv_bboxes[bbox_index + 3], recv_bboxes[bbox_index + 4],

recv_bboxes[bbox_index + 5]}};

result.add_volume_unsafe(

index_elem, mats[index], vols[index], &slot_bbox);

} else {

result.add_volume_unsafe(index_elem, mats[index], vols[index]);

}

}

}

}

}

} else {

// Send material indices and volumes to process 0

MPI_Send(materials, total, MPI_INT32_T, 0, mpi::rank, mpi::intracomm);

MPI_Send(volumes, total, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);

if (compute_bboxes) {

MPI_Send(bboxes, total_bbox, MPI_DOUBLE, 0, mpi::rank, mpi::intracomm);

}

}

}

// Report time for MPI communication

double t_mpi = timer.elapsed() - t_raytrace;

#else

double t_mpi = 0.0;

#endif

// Normalize based on known volumes of elements

for (int i = 0; i < this->n_bins(); ++i) {

// Estimated total volume in element i

double volume = 0.0;

for (int j = 0; j < table_size; ++j) {

volume += result.volumes(i, j);

}

// Renormalize volumes based on known volume of element i

double norm = this->volume(i) / volume;

for (int j = 0; j < table_size; ++j) {

result.volumes(i, j) *= norm;

}

}

// Get total time and normalization time

timer.stop();

double t_total = timer.elapsed();

double t_norm = t_total - t_raytrace - t_mpi;

// Show timing statistics

if (settings::verbosity < 7 || !mpi::master)

return;

header("Timing Statistics", 7);

fmt::print(" Total time elapsed = {:.4e} seconds\n", t_total);

fmt::print(" Ray tracing = {:.4e} seconds\n", t_raytrace);

fmt::print(" MPI communication = {:.4e} seconds\n", t_mpi);

fmt::print(" Normalization = {:.4e} seconds\n", t_norm);

fmt::print(" Calculation rate = {:.4e} rays/seconds\n",

n_total / t_raytrace);

fmt::print(" Calculation rate (per thread) = {:.4e} rays/seconds\n",

n_total / (t_raytrace * mpi::n_procs * num_threads()));

std::fflush(stdout);

}

void Mesh::to_hdf5(hid_t group) const

{

// Create group for mesh

std::string group_name = fmt::format("mesh {}", id_);

hid_t mesh_group = create_group(group, group_name.c_str());

// Write mesh type

write_dataset(mesh_group, "type", this->get_mesh_type());

// Write mesh ID

write_attribute(mesh_group, "id", id_);

// Write mesh name

write_dataset(mesh_group, "name", name_);

// Write mesh data

this->to_hdf5_inner(mesh_group);

// Close group

close_group(mesh_group);

}

//==============================================================================

// Structured Mesh implementation

//==============================================================================

std::string StructuredMesh::bin_label(int bin) const

{

MeshIndex ijk = get_indices_from_bin(bin);

if (n_dimension_ > 2) {

return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]);

} else if (n_dimension_ > 1) {

return fmt::format("Mesh Index ({}, {})", ijk[0], ijk[1]);

} else {

return fmt::format("Mesh Index ({})", ijk[0]);

}

}

tensor::Tensor<int> StructuredMesh::get_shape_tensor() const

{

return tensor::Tensor<int>(shape_.data(), static_cast<size_t>(n_dimension_));

}

Position StructuredMesh::sample_element(

const MeshIndex& ijk, uint64_t* seed) const

{

// lookup the lower/upper bounds for the mesh element

double x_min = negative_grid_boundary(ijk, 0);

double x_max = positive_grid_boundary(ijk, 0);

double y_min = (n_dimension_ >= 2) ? negative_grid_boundary(ijk, 1) : 0.0;

double y_max = (n_dimension_ >= 2) ? positive_grid_boundary(ijk, 1) : 0.0;

double z_min = (n_dimension_ == 3) ? negative_grid_boundary(ijk, 2) : 0.0;

double z_max = (n_dimension_ == 3) ? positive_grid_boundary(ijk, 2) : 0.0;

return {x_min + (x_max - x_min) * prn(seed),

y_min + (y_max - y_min) * prn(seed), z_min + (z_max - z_min) * prn(seed)};

}

//==============================================================================

// Unstructured Mesh implementation

//==============================================================================

std::string read_mesh_interface(const pugi::xml_node& node)

{

std::string mesh_interface;

if (check_for_node(node, "interface")) {

mesh_interface = get_node_value(node, "interface");

} else {

mesh_interface = "native";

}

return std::move(mesh_interface);

}

std::string read_mesh_interface(hid_t group)

{

std::string mesh_interface;

if (object_exists(group, "interface")) {

read_dataset(group, "interface", mesh_interface);

} else {

mesh_interface = "native";

}

return std::move(mesh_interface);

}

std::string read_mesh_library(const pugi::xml_node& node)

{

std::string mesh_library;

if (check_for_node(node, "library")) {

mesh_library = get_node_value(node, "library");

} else {

mesh_library = "moab";

}

return std::move(mesh_library);

}

std::string read_mesh_library(hid_t group)

{

std::string mesh_library;

if (object_exists(group, "library")) {

read_dataset(group, "library", mesh_library);

} else {

mesh_library = "moab";

}

return std::move(mesh_library);

}

template<typename T>

std::unique_ptr<UnstructuredMesh> UnstructuredMesh::create(T dataset)

{

std::string mesh_interface = read_mesh_interface(dataset);

std::string mesh_library = read_mesh_library(dataset);

if (mesh_interface == "xdg") {

#ifdef OPENMC_XDG_ENABLED

return make_unique<XDGMesh>(dataset);

#else

fatal_error(fmt::format(

"XDG unstructured mesh interface with library '{}' is not enabled in "

"this build of OpenMC.",

mesh_library));

#endif

}

if (mesh_interface != "native") {

fatal_error(fmt::format(

"Unrecognized unstructured mesh interface '{}'.", mesh_interface));

}

#ifdef OPENMC_DAGMC_ENABLED

if (mesh_library == MOABMesh::mesh_lib_type) {

return make_unique<MOABMesh>(dataset);

}

#endif

#ifdef OPENMC_LIBMESH_ENABLED

if (mesh_library == LibMesh::mesh_lib_type) {

return make_unique<LibMesh>(dataset);

}

#endif

fatal_error(fmt::format("Native unstructured mesh library '{}' is not "

"enabled in this build of OpenMC.",

mesh_library));

}

UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)

{

n_dimension_ = 3;

// check if a length unit multiplier was specified

if (check_for_node(node, "length_multiplier")) {

length_multiplier_ = std::stod(get_node_value(node, "length_multiplier"));

}

// get the filename of the unstructured mesh to load

if (check_for_node(node, "filename")) {

filename_ = get_node_value(node, "filename");

if (!file_exists(filename_)) {

fatal_error("Mesh file '" + filename_ + "' does not exist!");