Convert particle track writing to C++

CMakeLists.txt

@@ -327,7 +327,6 @@ add_library(libopenmc SHARED
   src/state_point.F90
   src/stl_vector.F90
   src/string.F90
-  src/track_output.F90
   src/vector_header.F90
   src/xml_interface.F90
   src/tallies/tally.F90
@@ -419,6 +418,7 @@ add_library(libopenmc SHARED
   src/tallies/trigger.cpp
   src/timer.cpp
   src/thermal.cpp
+  src/track_output.cpp
   src/volume_calc.cpp
   src/wmp.cpp
   src/xml_interface.cpp

include/openmc/track_output.h

@@ -0,0 +1,18 @@
+#ifndef OPENMC_TRACK_OUTPUT_H
+#define OPENMC_TRACK_OUTPUT_H
+
+#include "openmc/particle.h"
+
+namespace openmc {
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void add_particle_track();
+void write_particle_track(const Particle& p);
+void finalize_particle_track(const Particle& p);
+
+} // namespace openmc
+
+#endif // OPENMC_TRACK_OUTPUT_H

src/constants.F90

@@ -4,12 +4,6 @@ module constants
 
   implicit none
 
-  ! ============================================================================
-  ! VERSIONING NUMBERS
-
-  ! Version numbers for binary files
-  integer, parameter :: VERSION_TRACK(2) = [2, 0]
-
   ! ============================================================================
   ! ADJUSTABLE PARAMETERS
 

src/particle.cpp

@@ -24,6 +24,7 @@
 #include "openmc/tallies/derivative.h"
 #include "openmc/tallies/tally.h"
 #include "openmc/tallies/tally_scoring.h"
+#include "openmc/track_output.h"
 
 namespace openmc {
 
@@ -132,11 +133,6 @@ Particle::from_source(const Bank* src)
   last_E = E;
 }
 
-extern "C" void initialize_particle_track();
-extern "C" void write_particle_track(const Particle*);
-extern "C" void add_particle_track();
-extern "C" void finalize_particle_track(const Particle*);
-
 void
 Particle::transport()
 {
@@ -160,7 +156,7 @@ Particle::transport()
   }
 
   // Prepare to write out particle track.
-  if (write_track) initialize_particle_track();
+  if (write_track) add_particle_track();
 
   // Every particle starts with no accumulated flux derivative.
   if (!model::active_tallies.empty()) zero_flux_derivs();
@@ -194,7 +190,7 @@ Particle::transport()
     }
 
     // Write particle track.
-    if (write_track) write_particle_track(this);
+    if (write_track) write_particle_track(*this);
 
     if (settings::check_overlaps) check_cell_overlap(this);
 
@@ -396,8 +392,8 @@ Particle::transport()
 
   // Finish particle track output.
   if (write_track) {
-    write_particle_track(this);
-    finalize_particle_track(this);
+    write_particle_track(*this);
+    finalize_particle_track(*this);
   }
 }
 

src/track_output.F90

@@ -1,132 +0,0 @@
-!===============================================================================
-! TRACK_OUTPUT handles output of particle tracks (the paths taken by particles
-! as they are transported through the geometry).
-!===============================================================================
-
-module track_output
-
-  use, intrinsic :: ISO_C_BINDING
-
-  use constants
-  use hdf5_interface
-  use particle_header, only: Particle
-  use settings,        only: path_output
-  use simulation_header
-  use string,          only: to_str
-
-  implicit none
-  private
-
-  type TrackCoordinates
-    real(8), allocatable :: coords(:,:)
-  end type TrackCoordinates
-
-  type(TrackCoordinates), allocatable :: tracks(:)
-!$omp threadprivate(tracks)
-
-  public :: initialize_particle_track
-  public :: write_particle_track
-  public :: add_particle_track
-  public :: finalize_particle_track
-
-contains
-
-!===============================================================================
-! INITIALIZE_PARTICLE_TRACK allocates the array to store particle track
-! information
-!===============================================================================
-
-  subroutine initialize_particle_track() bind(C)
-    allocate(tracks(1))
-  end subroutine initialize_particle_track
-
-!===============================================================================
-! WRITE_PARTICLE_TRACK copies particle position to an array.
-!===============================================================================
-
-  subroutine write_particle_track(p) bind(C)
-    type(Particle), intent(in)  :: p
-    real(8), allocatable :: new_coords(:, :)
-
-    integer :: i
-    integer :: n_tracks
-
-    ! Add another column to coords
-    i = size(tracks)
-    if (allocated(tracks(i)%coords)) then
-      n_tracks = size(tracks(i)%coords, 2)
-      allocate(new_coords(3, n_tracks + 1))
-      new_coords(:, 1:n_tracks) = tracks(i)%coords
-      call move_alloc(FROM=new_coords, TO=tracks(i)%coords)
-    else
-      n_tracks = 0
-      allocate(tracks(i)%coords(3, 1))
-    end if
-
-    ! Write current coordinates into the newest column.
-    n_tracks = n_tracks + 1
-    tracks(i)%coords(:, n_tracks) = p%coord(1)%xyz
-  end subroutine write_particle_track
-
-!===============================================================================
-! ADD_PARTICLE_TRACK creates a new entry in the track coordinates for a
-! secondary particle
-!===============================================================================
-
-  subroutine add_particle_track() bind(C)
-    type(TrackCoordinates), allocatable :: new_tracks(:)
-
-    integer :: i
-
-    ! Determine current number of particle tracks
-    i = size(tracks)
-
-    ! Create array one larger than current
-    allocate(new_tracks(i + 1))
-
-    ! Copy memory and move allocation
-    new_tracks(1:i) = tracks(i)
-    call move_alloc(FROM=new_tracks, TO=tracks)
-
-  end subroutine add_particle_track
-
-!===============================================================================
-! FINALIZE_PARTICLE_TRACK writes the particle track array to disk.
-!===============================================================================
-
-  subroutine finalize_particle_track(p) bind(C)
-    type(Particle), intent(in)  :: p
-
-    integer :: i
-    integer :: n_particle_tracks
-    integer(HID_T) :: file_id
-    character(MAX_FILE_LEN) :: fname
-    integer, allocatable :: n_coords(:)
-
-    fname = trim(path_output) // 'track_' // trim(to_str(current_batch)) &
-         // '_' // trim(to_str(current_gen)) // '_' // trim(to_str(p%id)) &
-         // '.h5'
-
-    ! Determine total number of particles and number of coordinates for each
-    n_particle_tracks = size(tracks)
-    allocate(n_coords(n_particle_tracks))
-    do i = 1, n_particle_tracks
-      n_coords(i) = size(tracks(i)%coords, 2)
-    end do
-
-!$omp critical (FinalizeParticleTrack)
-    file_id = file_open(fname, 'w')
-    call write_attribute(file_id, 'filetype', 'track')
-    call write_attribute(file_id, 'version', VERSION_TRACK)
-    call write_attribute(file_id, 'n_particles', n_particle_tracks)
-    call write_attribute(file_id, 'n_coords', n_coords)
-    do i = 1, n_particle_tracks
-      call write_dataset(file_id, 'coordinates_' // trim(to_str(i)), &
-           tracks(i)%coords)
-    end do
-    call file_close(file_id)
-!$omp end critical (FinalizeParticleTrack)
-    deallocate(tracks)
-  end subroutine finalize_particle_track
-
-end module track_output

src/track_output.cpp

@@ -0,0 +1,80 @@
+#include "openmc/track_output.h"
+
+#include "openmc/constants.h"
+#include "openmc/hdf5_interface.h"
+#include "openmc/position.h"
+#include "openmc/settings.h"
+#include "openmc/simulation.h"
+
+#include "xtensor/xtensor.hpp"
+
+#include <cstddef> // for size_t
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace openmc {
+
+//==============================================================================
+// Global variables
+//==============================================================================
+
+// Forward declaration needed in order to declare tracks as threadprivate
+extern std::vector<std::vector<Position>> tracks;
+#pragma omp threadprivate(tracks)
+
+std::vector<std::vector<Position>> tracks;
+
+//==============================================================================
+// Non-member functions
+//==============================================================================
+
+void add_particle_track()
+{
+  tracks.emplace_back();
+}
+
+void write_particle_track(const Particle& p)
+{
+  tracks.back().push_back({p.coord[0].xyz});
+}
+
+void finalize_particle_track(const Particle& p)
+{
+  std::stringstream filename;
+  filename << settings::path_output << "track_" << simulation::current_batch
+    << '_' << simulation::current_gen << '_' << p.id << ".h5";
+
+  // Determine number of coordinates for each particle
+  std::vector<int> n_coords;
+  for (auto& coords : tracks) {
+    n_coords.push_back(coords.size());
+  }
+
+  #pragma omp critical (FinalizeParticleTrack)
+  {
+    hid_t file_id = file_open(filename.str().c_str(), 'w');
+    write_attribute(file_id, "filetype", "track");
+    write_attribute(file_id, "version", VERSION_TRACK);
+    write_attribute(file_id, "n_particles", tracks.size());
+    write_attribute(file_id, "n_coords", n_coords);
+    for (int i = 1; i <= tracks.size(); ++i) {
+      const auto& t {tracks[i-1]};
+      size_t n = t.size();
+      xt::xtensor<double, 2> data({n,3});
+      for (int j = 0; j < n; ++j) {
+        data(j, 0) = t[j].x;
+        data(j, 1) = t[j].y;
+        data(j, 2) = t[j].z;
+      }
+      std::string name = "coordinates_" + std::to_string(i);
+      write_dataset(file_id, name.c_str(), data);
+    }
+    file_close(file_id);
+  }
+
+  // Clear particle tracks
+  tracks.clear();
+}
+
+} // namespace openmc