particle_restart.cpp

#include "openmc/particle_restart.h"

#include "openmc/array.h"
#include "openmc/bank.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/photon.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/tally.h"
#include "openmc/track_output.h"

#include <algorithm> // for copy
#include <stdexcept>
#include <string>

namespace openmc {

void read_particle_restart(Particle& p, RunMode& previous_run_mode)
{
  // Write meessage
  write_message(
    5, "Loading particle restart file {}", settings::path_particle_restart);

  // Open file
  hid_t file_id = file_open(settings::path_particle_restart, 'r');

  // Read data from file
  bool legacy_particle_codes = true;
  if (attribute_exists(file_id, "version")) {
    array<int, 2> version;
    read_attribute(file_id, "version", version);
    if (version[0] > VERSION_PARTICLE_RESTART[0] ||
        (version[0] == VERSION_PARTICLE_RESTART[0] && version[1] >= 1)) {
      legacy_particle_codes = false;
    }
  }

  read_dataset(file_id, "current_batch", simulation::current_batch);
  read_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
  read_dataset(file_id, "current_generation", simulation::current_gen);
  read_dataset(file_id, "n_particles", settings::n_particles);
  std::string mode;
  read_dataset(file_id, "run_mode", mode);
  if (mode == "eigenvalue") {
    previous_run_mode = RunMode::EIGENVALUE;
  } else if (mode == "fixed source") {
    previous_run_mode = RunMode::FIXED_SOURCE;
  }
  read_dataset(file_id, "id", p.id());
  int type;
  read_dataset(file_id, "type", type);
  p.type() = legacy_particle_codes ? legacy_particle_index_to_type(type)
                                   : ParticleType {type};
  read_dataset(file_id, "weight", p.wgt());
  read_dataset(file_id, "energy", p.E());
  read_dataset(file_id, "xyz", p.r());
  read_dataset(file_id, "uvw", p.u());
  read_dataset(file_id, "time", p.time());

  // Set energy group and average energy in multi-group mode
  if (!settings::run_CE) {
    p.g() = p.E();
    p.E() = data::mg.energy_bin_avg_[p.g()];
  }

  // Set particle last attributes
  p.wgt_last() = p.wgt();
  p.r_last_current() = p.r();
  p.r_last() = p.r();
  p.u_last() = p.u();
  p.E_last() = p.E();
  p.g_last() = p.g();
  p.time_last() = p.time();

  // Close hdf5 file
  file_close(file_id);
}

void run_particle_restart()
{
  // Set verbosity high
  settings::verbosity = 10;

  // Initialize nuclear data (energy limits, log grid, etc.)
  initialize_data();

  // Initialize the particle to be tracked
  Particle p;

  // Read in the restart information
  RunMode previous_run_mode;
  read_particle_restart(p, previous_run_mode);

  // write track if that was requested on command line
  if (settings::write_all_tracks) {
    open_track_file();
    p.write_track() = true;
  }

  // Set all tallies to 0 for now (just tracking errors)
  model::tallies.clear();

  // Allocate progeny_per_particle if needed for shared secondary mode
  // (event_death() writes to this array). Set current_work to 0 since we
  // only have one particle being restarted.
  if (settings::use_shared_secondary_bank) {
    p.current_work() = 0;
    simulation::progeny_per_particle.resize(1, 0);
  }

  // Compute random number seed
  int64_t particle_seed = compute_transport_seed(p.id());
  init_particle_seeds(particle_seed, p.seeds());

  // Force calculation of cross-sections by setting last energy to zero
  if (settings::run_CE) {
    p.invalidate_neutron_xs();
  }

  // Prepare to write out particle track.
  if (p.write_track())
    add_particle_track(p);

  // Transport neutron
  transport_history_based_single_particle(p);

  // Write output if particle made it
  print_particle(p);

  if (settings::write_all_tracks) {
    close_track_file();
  }
}

// Function for automatically creating tracks file for lost particles
void run_lost_particle_track(Particle& lost)
{
  if (in_lost_track)
    return;
  in_lost_track = true;

#pragma omp critical(TrackFile)
  {
    if (track_file < 0 && !lost_particle_track_file_open) {
      open_track_file();
      lost_particle_track_file_open = true;
    }
  }

  Particle p;
  p.id() = lost.id();

  int64_t i = lost.current_work();
  SourceSite site;
  if (settings::run_mode == RunMode::EIGENVALUE) {
    site = simulation::source_bank[i];
  } else if (settings::run_mode == RunMode::FIXED_SOURCE &&
             settings::use_shared_secondary_bank &&
             i < simulation::shared_secondary_bank_read.size()) {
    site = simulation::shared_secondary_bank_read[i];
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
    int64_t id = compute_transport_seed(compute_particle_id(i + 1));
    uint64_t seed = init_seed(id, STREAM_SOURCE);
    site = sample_external_source(&seed);
  }

  p.wgt() = site.wgt;
  p.E() = site.E;
  p.r() = site.r;
  p.u() = site.u;
  p.time() = site.time;
  p.type() = lost.type();

  int64_t particle_seed = compute_transport_seed(p.id());
  init_particle_seeds(particle_seed, p.seeds());

  if (settings::run_CE) {
    p.invalidate_neutron_xs();
  }

  p.write_track() = true;
  add_particle_track(p);
  transport_history_based_single_particle(p);
  // finalize_particle_track called internally by transport, don't call again

  in_lost_track = false;
}

} // namespace openmc