xsdata.cpp

#include "openmc/xsdata.h"

#include <algorithm>
#include <cmath>
#include <cstdlib>
#include <numeric>

#include "xtensor/xbuilder.hpp"
#include "xtensor/xindex_view.hpp"
#include "xtensor/xmath.hpp"
#include "xtensor/xview.hpp"

#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/mgxs_interface.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"

namespace openmc {

//==============================================================================
// XsData class methods
//==============================================================================

XsData::XsData(bool fissionable, AngleDistributionType scatter_format,
  int n_pol, int n_azi, size_t n_groups, size_t n_d_groups)
  : n_g_(n_groups), n_dg_(n_d_groups)
{
  size_t n_ang = n_pol * n_azi;

  // check to make sure scatter format is OK before we allocate
  if (scatter_format != AngleDistributionType::HISTOGRAM &&
      scatter_format != AngleDistributionType::TABULAR &&
      scatter_format != AngleDistributionType::LEGENDRE) {
    fatal_error("Invalid scatter_format!");
  }
  // allocate all [temperature][angle][in group] quantities
  vector<size_t> shape {n_ang, n_g_};
  total = xt::zeros<double>(shape);
  absorption = xt::zeros<double>(shape);
  inverse_velocity = xt::zeros<double>(shape);
  if (fissionable) {
    fission = xt::zeros<double>(shape);
    nu_fission = xt::zeros<double>(shape);
    prompt_nu_fission = xt::zeros<double>(shape);
    kappa_fission = xt::zeros<double>(shape);
  }

  // allocate decay_rate; [temperature][angle][delayed group]
  shape[1] = n_dg_;
  decay_rate = xt::zeros<double>(shape);

  if (fissionable) {
    shape = {n_ang, n_dg_, n_g_};
    // allocate delayed_nu_fission; [temperature][angle][delay group][in group]
    delayed_nu_fission = xt::zeros<double>(shape);

    // chi_prompt; [temperature][angle][in group][out group]
    shape = {n_ang, n_g_, n_g_};
    chi_prompt = xt::zeros<double>(shape);

    // chi_delayed; [temperature][angle][delay group][in group][out group]
    shape = {n_ang, n_dg_, n_g_, n_g_};
    chi_delayed = xt::zeros<double>(shape);
  }

  for (int a = 0; a < n_ang; a++) {
    if (scatter_format == AngleDistributionType::HISTOGRAM) {
      scatter.emplace_back(new ScattDataHistogram);
    } else if (scatter_format == AngleDistributionType::TABULAR) {
      scatter.emplace_back(new ScattDataTabular);
    } else if (scatter_format == AngleDistributionType::LEGENDRE) {
      scatter.emplace_back(new ScattDataLegendre);
    }
  }
}

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

void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable,
  AngleDistributionType scatter_format,
  AngleDistributionType final_scatter_format, int order_data, bool is_isotropic,
  int n_pol, int n_azi)
{
  // Reconstruct the dimension information so it doesn't need to be passed
  size_t n_ang = n_pol * n_azi;
  size_t energy_groups = total.shape()[1];

  // Set the fissionable-specific data
  if (fissionable) {
    fission_from_hdf5(xsdata_grp, n_ang, is_isotropic);
  }
  // Get the non-fission-specific data
  read_nd_vector(xsdata_grp, "decay-rate", decay_rate);
  read_nd_vector(xsdata_grp, "absorption", absorption, true);
  read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);

  // Get scattering data
  scatter_from_hdf5(
    xsdata_grp, n_ang, scatter_format, final_scatter_format, order_data);

  // Check absorption to ensure it is not 0 since it is often the
  // denominator in tally methods
  xt::filtration(absorption, xt::equal(absorption, 0.)) = 1.e-10;

  // Get or calculate the total x/s
  if (object_exists(xsdata_grp, "total")) {
    read_nd_vector(xsdata_grp, "total", total);
  } else {
    for (size_t a = 0; a < n_ang; a++) {
      for (size_t gin = 0; gin < energy_groups; gin++) {
        total(a, gin) = absorption(a, gin) + scatter[a]->scattxs[gin];
      }
    }
  }

  // Fix if total is 0, since it is in the denominator when tallying
  xt::filtration(total, xt::equal(total, 0.)) = 1.e-10;
}

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

void XsData::fission_vector_beta_from_hdf5(
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
{
  // Data is provided as nu-fission and chi with a beta for delayed info

  // Get chi
  xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "chi", temp_chi, true);

  // Normalize chi by summing over the outgoing groups for each incoming angle
  temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());

  // Now every incoming group in prompt_chi and delayed_chi is the normalized
  // chi we just made
  chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
  chi_delayed =
    xt::view(temp_chi, xt::all(), xt::newaxis(), xt::newaxis(), xt::all());

  // Get nu-fission
  xt::xtensor<double, 2> temp_nufiss({n_ang, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true);

  // Get beta (strategy will depend upon the number of dimensions in beta)
  hid_t beta_dset = open_dataset(xsdata_grp, "beta");
  int beta_ndims = dataset_ndims(beta_dset);
  close_dataset(beta_dset);
  int ndim_target = 1;
  if (!is_isotropic)
    ndim_target += 2;
  if (beta_ndims == ndim_target) {
    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);

    // Set prompt_nu_fission = (1. - beta_total)*nu_fission
    prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1}));

    // Set delayed_nu_fission as beta * nu_fission
    delayed_nu_fission =
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) *
      xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all());
  } else if (beta_ndims == ndim_target + 1) {
    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);

    // Set prompt_nu_fission = (1. - beta_total)*nu_fission
    prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1}));

    // Set delayed_nu_fission as beta * nu_fission
    delayed_nu_fission =
      temp_beta * xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all());
  }
}

void XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
{
  // Data is provided separately as prompt + delayed nu-fission and chi

  // Get chi-prompt
  xt::xtensor<double, 2> temp_chi_p({n_ang, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true);

  // Normalize chi by summing over the outgoing groups for each incoming angle
  temp_chi_p /= xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis());

  // Get chi-delayed
  xt::xtensor<double, 3> temp_chi_d({n_ang, n_dg_, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true);

  // Normalize chi by summing over the outgoing groups for each incoming angle
  temp_chi_d /=
    xt::view(xt::sum(temp_chi_d, {2}), xt::all(), xt::all(), xt::newaxis());

  // Now assign the prompt and delayed chis by replicating for each incoming
  // group
  chi_prompt = xt::view(temp_chi_p, xt::all(), xt::newaxis(), xt::all());
  chi_delayed =
    xt::view(temp_chi_d, xt::all(), xt::all(), xt::newaxis(), xt::all());

  // Get prompt and delayed nu-fission directly
  read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, true);
  read_nd_vector(xsdata_grp, "delayed-nu-fission", delayed_nu_fission, true);
}

void XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
{
  // No beta is provided and there is no prompt/delay distinction.
  // Therefore, the code only considers the data as prompt.

  // Get chi
  xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "chi", temp_chi, true);

  // Normalize chi by summing over the outgoing groups for each incoming angle
  temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());

  // Now every incoming group in self.chi is the normalized chi we just made
  chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());

  // Get nu-fission directly
  read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission, true);
}

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

void XsData::fission_matrix_beta_from_hdf5(
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
{
  // Data is provided as nu-fission and chi with a beta for delayed info

  // Get nu-fission matrix
  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);

  // Get beta (strategy will depend upon the number of dimensions in beta)
  hid_t beta_dset = open_dataset(xsdata_grp, "beta");
  int beta_ndims = dataset_ndims(beta_dset);
  close_dataset(beta_dset);
  int ndim_target = 1;
  if (!is_isotropic)
    ndim_target += 2;
  if (beta_ndims == ndim_target) {
    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);

    xt::xtensor<double, 1> temp_beta_sum({n_ang}, 0.);
    temp_beta_sum = xt::sum(temp_beta, {1});

    // prompt_nu_fission is the sum of this matrix over outgoing groups and
    // multiplied by (1 - beta_sum)
    prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum);

    // Store chi-prompt
    chi_prompt =
      xt::view(1.0 - temp_beta_sum, xt::all(), xt::newaxis(), xt::newaxis()) *
      temp_matrix;

    // delayed_nu_fission is the sum of this matrix over outgoing groups and
    // multiplied by beta
    delayed_nu_fission =
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) *
      xt::view(xt::sum(temp_matrix, {2}), xt::all(), xt::newaxis(), xt::all());

    // Store chi-delayed
    chi_delayed =
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis(), xt::newaxis()) *
      xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());

  } else if (beta_ndims == ndim_target + 1) {
    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);

    xt::xtensor<double, 2> temp_beta_sum({n_ang, n_g_}, 0.);
    temp_beta_sum = xt::sum(temp_beta, {1});

    // prompt_nu_fission is the sum of this matrix over outgoing groups and
    // multiplied by (1 - beta_sum)
    prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum);

    // Store chi-prompt
    chi_prompt =
      xt::view(1.0 - temp_beta_sum, xt::all(), xt::all(), xt::newaxis()) *
      temp_matrix;

    // delayed_nu_fission is the sum of this matrix over outgoing groups and
    // multiplied by beta
    delayed_nu_fission = temp_beta * xt::view(xt::sum(temp_matrix, {2}),
                                       xt::all(), xt::newaxis(), xt::all());

    // Store chi-delayed
    chi_delayed =
      xt::view(temp_beta, xt::all(), xt::all(), xt::all(), xt::newaxis()) *
      xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
  }

  // Normalize both chis
  chi_prompt /=
    xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());

  chi_delayed /= xt::view(
    xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis());
}

void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
{
  // Data is provided separately as prompt + delayed nu-fission and chi

  // Get the prompt nu-fission matrix
  xt::xtensor<double, 3> temp_matrix_p({n_ang, n_g_, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true);

  // prompt_nu_fission is the sum over outgoing groups
  prompt_nu_fission = xt::sum(temp_matrix_p, {2});

  // chi_prompt is this matrix but normalized over outgoing groups, which we
  // have already stored in prompt_nu_fission
  chi_prompt = temp_matrix_p /
               xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());

  // Get the delayed nu-fission matrix
  xt::xtensor<double, 4> temp_matrix_d({n_ang, n_dg_, n_g_, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true);

  // delayed_nu_fission is the sum over outgoing groups
  delayed_nu_fission = xt::sum(temp_matrix_d, {3});

  // chi_prompt is this matrix but normalized over outgoing groups, which we
  // have already stored in prompt_nu_fission
  chi_delayed = temp_matrix_d / xt::view(delayed_nu_fission, xt::all(),
                                  xt::all(), xt::all(), xt::newaxis());
}

void XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
{
  // No beta is provided and there is no prompt/delay distinction.
  // Therefore, the code only considers the data as prompt.

  // Get nu-fission matrix
  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
  read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);

  // prompt_nu_fission is the sum over outgoing groups
  prompt_nu_fission = xt::sum(temp_matrix, {2});

  // chi_prompt is this matrix but normalized over outgoing groups, which we
  // have already stored in prompt_nu_fission
  chi_prompt = temp_matrix /
               xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
}

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

void XsData::fission_from_hdf5(
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
{
  // Get the fission and kappa_fission data xs; these are optional
  read_nd_vector(xsdata_grp, "fission", fission);
  read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);

  // Get the data; the strategy for doing so depends on if the data is provided
  // as a nu-fission matrix or a set of chi and nu-fission vectors
  if (object_exists(xsdata_grp, "chi") ||
      object_exists(xsdata_grp, "chi-prompt")) {
    if (n_dg_ == 0) {
      fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang);
    } else {
      if (object_exists(xsdata_grp, "beta")) {
        fission_vector_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic);
      } else {
        fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang);
      }
    }
  } else {
    if (n_dg_ == 0) {
      fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang);
    } else {
      if (object_exists(xsdata_grp, "beta")) {
        fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic);
      } else {
        fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang);
      }
    }
  }

  // Combine prompt_nu_fission and delayed_nu_fission into nu_fission
  if (n_dg_ == 0) {
    nu_fission = prompt_nu_fission;
  } else {
    nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1});
  }
}

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

void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
  AngleDistributionType scatter_format,
  AngleDistributionType final_scatter_format, int order_data)
{
  if (!object_exists(xsdata_grp, "scatter_data")) {
    fatal_error("Must provide scatter_data group!");
  }
  hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");

  // Get the outgoing group boundary indices
  xt::xtensor<int, 2> gmin({n_ang, n_g_}, 0.);
  read_nd_vector(scatt_grp, "g_min", gmin, true);
  xt::xtensor<int, 2> gmax({n_ang, n_g_}, 0.);
  read_nd_vector(scatt_grp, "g_max", gmax, true);

  // Make gmin and gmax start from 0 vice 1 as they do in the library
  gmin -= 1;
  gmax -= 1;

  // Now use this info to find the length of a vector to hold the flattened
  // data.
  size_t length = order_data * xt::sum(gmax - gmin + 1)();

  double_4dvec input_scatt(n_ang, double_3dvec(n_g_));
  xt::xtensor<double, 1> temp_arr({length}, 0.);
  read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);

  // Compare the number of orders given with the max order of the problem;
  // strip off the superfluous orders if needed
  int order_dim;
  if (scatter_format == AngleDistributionType::LEGENDRE) {
    order_dim = std::min(order_data - 1, settings::max_order) + 1;
  } else {
    order_dim = order_data;
  }

  // convert the flattened temp_arr to a jagged array for passing to
  // scatt data
  size_t temp_idx = 0;
  for (size_t a = 0; a < n_ang; a++) {
    for (size_t gin = 0; gin < n_g_; gin++) {
      input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
      for (size_t i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
        input_scatt[a][gin][i_gout].resize(order_dim);
        for (size_t l = 0; l < order_dim; l++) {
          input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++];
        }
        // Adjust index for the orders we didnt take
        temp_idx += (order_data - order_dim);
      }
    }
  }

  // Get multiplication matrix
  double_3dvec temp_mult(n_ang, double_2dvec(n_g_));
  if (object_exists(scatt_grp, "multiplicity_matrix")) {
    temp_arr.resize({length / order_data});
    read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);

    // convert the flat temp_arr to a jagged array for passing to scatt data
    size_t temp_idx = 0;
    for (size_t a = 0; a < n_ang; a++) {
      for (size_t gin = 0; gin < n_g_; gin++) {
        temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
        for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
          temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
        }
      }
    }
  } else {
    // Use a default: multiplicities are 1.0.
    for (size_t a = 0; a < n_ang; a++) {
      for (size_t gin = 0; gin < n_g_; gin++) {
        temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
        for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
          temp_mult[a][gin][i_gout] = 1.;
        }
      }
    }
  }
  close_group(scatt_grp);

  // Finally, convert the Legendre data to tabular, if needed
  if (scatter_format == AngleDistributionType::LEGENDRE &&
      final_scatter_format == AngleDistributionType::TABULAR) {
    for (size_t a = 0; a < n_ang; a++) {
      ScattDataLegendre legendre_scatt;
      xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
      xt::xtensor<int, 1> in_gmax = xt::view(gmax, a, xt::all());

      legendre_scatt.init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]);

      // Now create a tabular version of legendre_scatt
      convert_legendre_to_tabular(
        legendre_scatt, *static_cast<ScattDataTabular*>(scatter[a].get()));

      scatter_format = final_scatter_format;
    }
  } else {
    // We are sticking with the current representation
    // Initialize the ScattData object with this data
    for (size_t a = 0; a < n_ang; a++) {
      xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
      xt::xtensor<int, 1> in_gmax = xt::view(gmax, a, xt::all());
      scatter[a]->init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]);
    }
  }
}

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

void XsData::combine(
  const vector<XsData*>& those_xs, const vector<double>& scalars)
{
  // Combine the non-scattering data
  for (size_t i = 0; i < those_xs.size(); i++) {
    XsData* that = those_xs[i];
    if (!equiv(*that))
      fatal_error("Cannot combine the XsData objects!");
    double scalar = scalars[i];
    total += scalar * that->total;
    absorption += scalar * that->absorption;
    if (i == 0) {
      inverse_velocity = that->inverse_velocity;
    }
    if (that->prompt_nu_fission.shape()[0] > 0) {
      nu_fission += scalar * that->nu_fission;
      prompt_nu_fission += scalar * that->prompt_nu_fission;
      kappa_fission += scalar * that->kappa_fission;
      fission += scalar * that->fission;
      delayed_nu_fission += scalar * that->delayed_nu_fission;
      chi_prompt += scalar *
                    xt::view(xt::sum(that->prompt_nu_fission, {1}), xt::all(),
                      xt::newaxis(), xt::newaxis()) *
                    that->chi_prompt;
      chi_delayed += scalar *
                     xt::view(xt::sum(that->delayed_nu_fission, {2}), xt::all(),
                       xt::all(), xt::newaxis(), xt::newaxis()) *
                     that->chi_delayed;
    }
    decay_rate += scalar * that->decay_rate;
  }

  // Ensure the chi_prompt and chi_delayed are normalized to 1 for each
  // azimuthal angle and delayed group (for chi_delayed)
  chi_prompt /=
    xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
  chi_delayed /= xt::view(
    xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis());

  // Allow the ScattData object to combine itself
  for (size_t a = 0; a < total.shape()[0]; a++) {
    // Build vector of the scattering objects to incorporate
    vector<ScattData*> those_scatts(those_xs.size());
    for (size_t i = 0; i < those_xs.size(); i++) {
      those_scatts[i] = those_xs[i]->scatter[a].get();
    }

    // Now combine these guys
    scatter[a]->combine(those_scatts, scalars);
  }
}

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

bool XsData::equiv(const XsData& that)
{
  return (absorption.shape() == that.absorption.shape());
}

} // namespace openmc