diff --git a/src/config.cpp b/src/config.cpp
index 7f534ed6..dbd50466 100644
--- a/src/config.cpp
+++ b/src/config.cpp
@@ -38,7 +38,7 @@ void XDGConfig::initialize() {
   // libMesh respects the OpenMP settings of the host application
   if (n_threads() == -1) {
   #ifdef XDG_HAVE_OPENMP
-    set_n_threads(omp_get_num_threads());
+    set_n_threads(omp_get_max_threads());
   #else
     set_n_threads(1);
   #endif
diff --git a/tools/CMakeLists.txt b/tools/CMakeLists.txt
index 08004ab3..50a0c650 100644
--- a/tools/CMakeLists.txt
+++ b/tools/CMakeLists.txt
@@ -1,6 +1,7 @@
 set(TOOL_NAMES
 particle_sim
 ray_fire
+ray_benchmark
 find_volume
 point_in_volume
 overlap_check
diff --git a/tools/ray_benchmark.cpp b/tools/ray_benchmark.cpp
new file mode 100644
index 00000000..dea28351
--- /dev/null
+++ b/tools/ray_benchmark.cpp
@@ -0,0 +1,322 @@
+#include <cstdint>
+#include <cstdlib>
+#include <filesystem>
+#include <iostream>
+#include <memory>
+#include <stdexcept>
+#include <string>
+#include <vector>
+
+#include "argparse/argparse.hpp"
+#include <fmt/ranges.h>
+
+#include "xdg/config.h"
+#include "xdg/constants.h"
+#include "xdg/error.h"
+#include "xdg/timer.h"
+#include "xdg/vec3da.h"
+#include "xdg/xdg.h"
+
+#include "ray_benchmark.h"
+
+
+using namespace xdg;
+
+int main(int argc, char** argv)
+{
+  argparse::ArgumentParser args("XDG Raytracing throughput benchmarking tool",
+                                "1.0",
+                                argparse::default_arguments::help);
+
+  args.add_argument("filename")
+    .help("Path to the input file");
+
+  args.add_argument("volume")
+    .help("Volume ID to query")
+    .scan<'i', int>();
+
+  args.add_argument("-n", "--num-rays")
+    .default_value<std::uint32_t>(10'000'000)
+    .help("Number of rays to cast for the benchmark")
+    .scan<'u', std::uint32_t>();
+
+  args.add_argument("-s", "--seed")
+    .default_value<std::uint32_t>(12345)
+    .help("Seed for random ray generation")
+    .scan<'u', std::uint32_t>();
+
+  args.add_argument("-o", "-p", "--origin", "--position")
+    .help("Ray origin/position. Defaults to the center of the model bounding box")
+    .scan<'g', double>()
+    .nargs(3);
+
+  args.add_argument("--volume-center")
+    .default_value(false)
+    .implicit_value(true)
+    .help("Use the queried volume bounding box center as the ray origin");
+
+  args.add_argument("-m", "--mesh-library")
+    .help("Mesh library to use. One of (MOAB, LIBMESH)")
+    .default_value("MOAB");
+
+  args.add_argument("-rt", "--rt-library")
+    .help("Ray tracing library to use. Currently implemented: EMBREE")
+    .default_value("EMBREE");
+
+  args.add_argument("-l", "--list")
+    .default_value(false)
+    .implicit_value(true)
+    .help("List all volumes in the file and exit");
+
+  args.add_argument("-sr", "--source-radius")
+    .default_value(0.0)
+    .help("Radius of a scattered source around the origin")
+    .scan<'g', double>();
+
+  args.add_argument("--format")
+    .default_value("human")
+    .choices("human", "csv")
+    .help("stdout format. Human readable (default) or csv");
+
+  args.add_description(
+    "Benchmarks ray-fire throughput for a selected mesh volume. A source "
+    "position is provided and ray directions are randomly generated from it.");
+
+  try {
+    args.parse_args(argc, argv);
+  }
+  catch (const std::runtime_error& err) {
+    std::cout << err.what() << std::endl;
+    std::cout << args;
+    exit(0);
+  }
+
+  std::string mesh_str = args.get<std::string>("--mesh-library");
+  std::string rt_str = args.get<std::string>("--rt-library");
+  std::string rt_label = rt_str;
+
+  RTLibrary rt_lib;
+  if (rt_str == "EMBREE") {
+    rt_lib = RTLibrary::EMBREE;
+  } else {
+    fatal_error("Ray tracing library '{}' is not implemented in this benchmark tool yet", rt_str);
+  }
+
+  MeshLibrary mesh_lib;
+  if (mesh_str == "MOAB") {
+    mesh_lib = MeshLibrary::MOAB;
+  } else if (mesh_str == "LIBMESH") {
+    mesh_lib = MeshLibrary::LIBMESH;
+  } else {
+    fatal_error("Invalid mesh library '{}' specified", mesh_str);
+  }
+
+  const MeshID volume = args.get<int>("volume");
+  const std::string model_filename = args.get<std::string>("filename");
+  const std::string model_name = std::filesystem::path(model_filename).filename().string();
+  const std::size_t num_rays = args.get<std::uint32_t>("--num-rays");
+  const std::uint32_t seed = args.get<std::uint32_t>("--seed");
+  const double source_radius = args.get<double>("--source-radius");
+  const std::string output_format = args.get<std::string>("--format");
+
+  Timer wall_timer;
+  Timer setup_timer;
+  Timer generation_timer;
+  Timer trace_timer;
+
+  wall_timer.start();
+
+  // XDG setup and ray tracer initialisation
+  setup_timer.start();
+  std::shared_ptr<XDG> xdg = XDG::create(mesh_lib, rt_lib);
+  const auto& mesh_manager = xdg->mesh_manager();
+  mesh_manager->load_file(model_filename);
+  mesh_manager->init();
+
+  if (args.get<bool>("--list")) {
+    std::cout << "[" << fmt::format("{}", fmt::join(mesh_manager->volumes(), ", ")) << "]\n";
+    return 0;
+  }
+
+  const auto origin_arg = args.present<std::vector<double>>("--origin");
+  bool use_volume_center = args.get<bool>("--volume-center");
+  if (origin_arg && use_volume_center) {
+    warning("--volume-center enabled but an explicit origin was also provided. The explicit origin will be used and the volume center will be ignored.");
+    use_volume_center = false;
+  }
+
+  Position origin = mesh_manager->global_bounding_box().center();
+  if (origin_arg) {
+    origin = Position(origin_arg.value());
+  } else if (use_volume_center) {
+    origin = mesh_manager->volume_bounding_box(volume).center();
+  }
+
+  xdg->prepare_volume_for_raytracing(volume);
+  xdg->ray_tracing_interface()->init();
+
+  const bool origin_in_volume = xdg->point_in_volume(volume, origin);
+  if (!origin_in_volume) {
+    const std::string origin_label = source_radius > 0.0 ? "source center" : "ray origin";
+    warning(fmt::format("The {} ({}, {}, {}) is not inside volume {}. Results may be skewed by fewer intersections with the BVH.",
+                        origin_label, origin.x, origin.y, origin.z, volume));
+  } else {
+    const double nearest_surface_distance = xdg->closest_distance(volume, origin);
+    if (source_radius > nearest_surface_distance) {
+      warning(fmt::format("The source radius ({}) is larger than the nearest surface distance ({}) from the source center to volume {}. "
+                          "Some sampled source points may be outside the volume.",
+                          source_radius, nearest_surface_distance, volume));
+    }
+  }
+
+  setup_timer.stop();
+
+  if (rt_lib == RTLibrary::EMBREE) {
+    rt_label += " (" + std::to_string(XDGConfig::config().n_threads())
+             + " CPU threads)";
+  }
+
+  const auto num_faces = mesh_manager->num_volume_faces(volume);
+
+
+  // Generate random rays from source
+  generation_timer.start();
+  std::vector<Position> origins(num_rays);
+  std::vector<Direction> directions(num_rays);
+
+  #pragma omp parallel for schedule(runtime)
+  for (std::size_t i = 0; i < num_rays; ++i) {
+    std::uint32_t state = seed ^ static_cast<std::uint32_t>(i);
+    auto sample = tools::benchmark::random_spherical_source(origin.x,
+                                                            origin.y,
+                                                            origin.z,
+                                                            state,
+                                                            source_radius);
+    origins[i] = Position(sample.position[0],
+                          sample.position[1],
+                          sample.position[2]);
+    directions[i] = Direction(sample.direction[0],
+                              sample.direction[1],
+                              sample.direction[2]);
+  }
+  generation_timer.stop();
+
+  // Trace rays
+  trace_timer.start();
+
+  std::size_t num_hits = 0;
+
+  #pragma omp parallel for schedule(runtime) reduction(+:num_hits)
+  for (std::size_t i = 0; i < num_rays; ++i) {
+    const auto hit = xdg->ray_fire(volume, origins[i], directions[i]);
+    if (hit.second != ID_NONE) num_hits++;
+  }
+
+  trace_timer.stop();
+
+  const std::size_t num_misses = num_rays - num_hits;
+  const double hit_fraction = num_rays > 0
+    ? static_cast<double>(num_hits) / static_cast<double>(num_rays)
+    : 0.0;
+  const double generation_time = generation_timer.elapsed();
+  const double trace_time = trace_timer.elapsed();
+  const double end_to_end_time = generation_time + trace_time;
+  const double setup_time = setup_timer.elapsed();
+  const double trace_only_rps = trace_time > 0.0
+    ? static_cast<double>(num_rays) / trace_time
+    : 0.0;
+  const double end_to_end_rps = end_to_end_time > 0.0
+    ? static_cast<double>(num_rays) / end_to_end_time
+    : 0.0;
+
+  wall_timer.stop();
+  const double wall_time = wall_timer.elapsed();
+
+  const std::vector<std::string> csv_columns {
+    "model",
+    "mesh_library",
+    "rt_library",
+    "volume",
+    "num_faces",
+    "num_rays",
+    "num_hits",
+    "num_misses",
+    "hit_fraction",
+    "seed",
+    "source_radius",
+    "origin_x",
+    "origin_y",
+    "origin_z",
+    "n_threads",
+    "initialisation_time_s",
+    "generation_time_s",
+    "trace_time_s",
+    "generation_trace_time_s",
+    "end_to_end_throughput_rays_per_s",
+    "trace_only_throughput_rays_per_s",
+    "wall_time_s"
+  };
+
+  const std::vector<std::string> csv_values {
+    model_name,
+    mesh_str,
+    rt_str,
+    fmt::format("{}", volume),
+    fmt::format("{}", num_faces),
+    fmt::format("{}", num_rays),
+    fmt::format("{}", num_hits),
+    fmt::format("{}", num_misses),
+    fmt::format("{}", hit_fraction),
+    fmt::format("{}", seed),
+    fmt::format("{}", source_radius),
+    fmt::format("{}", origin.x),
+    fmt::format("{}", origin.y),
+    fmt::format("{}", origin.z),
+    fmt::format("{}", XDGConfig::config().n_threads()),
+    fmt::format("{}", setup_time),
+    fmt::format("{}", generation_time),
+    fmt::format("{}", trace_time),
+    fmt::format("{}", end_to_end_time),
+    fmt::format("{}", end_to_end_rps),
+    fmt::format("{}", trace_only_rps),
+    fmt::format("{}", wall_time)
+  };
+
+  if (output_format == "csv") {
+    std::cout << fmt::format("{}\n", fmt::join(csv_columns, ","));
+    std::cout << fmt::format("{}\n", fmt::join(csv_values, ","));
+  } else {
+    std::cout << "\nXDG ray benchmark results\n";
+    std::cout << "----------------------------------------\n";
+    std::cout << "Model                 : " << model_name << "\n";
+    std::cout << "Mesh library          : " << mesh_str << "\n";
+    std::cout << "Ray tracing library   : " << rt_label << "\n";
+    std::cout << "Volume                : " << volume << "\n";
+    std::cout << "Volume faces          : " << num_faces << "\n";
+    std::cout << "Seed                  : " << seed << "\n";
+    std::cout << "Rays                  : " << num_rays << "\n";
+    if (source_radius != 0.0) {
+      std::cout << "Source center         : "
+                << origin.x << ", " << origin.y << ", " << origin.z << "\n";
+      std::cout << "Source radius         : " << source_radius << "\n";
+    } else {
+      std::cout << "Origin (fixed)        : "
+                << origin.x << ", " << origin.y << ", " << origin.z << "\n";
+    }
+    std::cout << "----------------------------------------\n";
+    std::cout << "Hits                  : " << num_hits << "\n";
+    std::cout << "Misses                : " << num_misses << "\n";
+    std::cout << "Hit fraction          : " << hit_fraction << "\n";
+    std::cout << "----------------------------------------\n";
+    std::cout << "Initialisation time   : " << setup_time << " s\n";
+    std::cout << "Ray generation time   : " << generation_time << " s\n";
+    std::cout << "Ray tracing time      : " << trace_time << " s\n";
+    std::cout << "Generation + tracing  : " << end_to_end_time << " s\n";
+    std::cout << "Full wall-clock time  : " << wall_time << " s\n";
+    std::cout << "----------------------------------------\n";
+    std::cout << "End-to-end throughput : " << end_to_end_rps << " rays/s\n";
+    std::cout << "Trace-only throughput : " << trace_only_rps << " rays/s\n";
+  }
+
+  return 0;
+}
diff --git a/tools/ray_benchmark.h b/tools/ray_benchmark.h
new file mode 100644
index 00000000..65a48339
--- /dev/null
+++ b/tools/ray_benchmark.h
@@ -0,0 +1,71 @@
+#ifndef _XDG_RAY_BENCHMARK_H
+#define _XDG_RAY_BENCHMARK_H
+
+#include <cmath>
+#include <cstdint>
+
+namespace xdg::tools::benchmark {
+
+struct SourceSample {
+  double position[3];
+  double direction[3];
+};
+
+#ifdef _OPENMP
+#pragma omp declare target
+#endif
+
+inline double rand01(std::uint32_t& state)
+{
+  state = state * 1664525u + 1013904223u;
+  return static_cast<double>(state) * (1.0 / 4294967296.0);
+}
+
+inline void random_unit_dir_lcg(std::uint32_t& state, double direction[3])
+{
+  double x1;
+  double x2;
+  double s;
+
+  do {
+    x1 = rand01(state) * 2.0 - 1.0;
+    x2 = rand01(state) * 2.0 - 1.0;
+    s = x1 * x1 + x2 * x2;
+  } while (s <= 0.0 || s >= 1.0);
+
+  const double t = 2.0 * std::sqrt(1.0 - s);
+  direction[0] = x1 * t;
+  direction[1] = x2 * t;
+  direction[2] = 1.0 - 2.0 * s;
+}
+
+inline SourceSample random_spherical_source(double origin_x,
+                                            double origin_y,
+                                            double origin_z,
+                                            std::uint32_t state,
+                                            double source_radius)
+{
+  SourceSample sample;
+  random_unit_dir_lcg(state, sample.direction);
+
+  sample.position[0] = origin_x;
+  sample.position[1] = origin_y;
+  sample.position[2] = origin_z;
+
+  if (source_radius > 0.0) {
+    const double radius = source_radius * std::cbrt(rand01(state));
+    sample.position[0] += sample.direction[0] * radius;
+    sample.position[1] += sample.direction[1] * radius;
+    sample.position[2] += sample.direction[2] * radius;
+  }
+
+  return sample;
+}
+
+#ifdef _OPENMP
+#pragma omp end declare target
+#endif
+
+} // namespace xdg::tools::benchmark
+
+#endif // _XDG_RAY_BENCHMARK_H