Added supportfor group-wise betas

include/openmc/xsdata.h

@@ -37,7 +37,7 @@ class XsData {
     //  the HDF5 file when beta is provided.
     void
     fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-         size_t energy_groups, size_t delayed_groups);
+         size_t energy_groups, size_t delayed_groups, bool is_isotropic);
 
     //! \brief Reads fission data formatted as chi and nu-fission vectors from
     //  the HDF5 file when beta is not provided.
@@ -55,7 +55,7 @@ class XsData {
     //  the HDF5 file when beta is provided.
     void
     fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-         size_t energy_groups, size_t delayed_groups);
+         size_t energy_groups, size_t delayed_groups, bool is_isotropic);
 
     //! \brief Reads fission data formatted as a nu-fission matrix from
     //  the HDF5 file when beta is not provided.

src/xsdata.cpp

@@ -4,7 +4,6 @@
 #include <cstdlib>
 #include <algorithm>
 #include <numeric>
-#include <iostream>
 
 #include "xtensor/xview.hpp"
 #include "xtensor/xindex_view.hpp"
@@ -124,7 +123,7 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
 
 void
 XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-     size_t energy_groups, size_t delayed_groups)
+     size_t energy_groups, size_t delayed_groups, bool is_isotropic)
 {
   // Data is provided as nu-fission and chi with a beta for delayed info
 
@@ -145,17 +144,35 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   xt::xtensor<double, 2> temp_nufiss({n_ang, energy_groups}, 0.);
   read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true);
 
-  // Get beta
-  xt::xtensor<double, 2> temp_beta({n_ang, delayed_groups}, 0.);
-  read_nd_vector(xsdata_grp, "beta", temp_beta, 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, delayed_groups}, 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 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());
+    // 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, delayed_groups, energy_groups},
+                                     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
@@ -219,7 +236,7 @@ XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang,
 
 void
 XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-     size_t energy_groups, size_t delayed_groups)
+     size_t energy_groups, size_t delayed_groups, bool is_isotropic)
 {
   // Data is provided as nu-fission and chi with a beta for delayed info
 
@@ -227,32 +244,65 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   xt::xtensor<double, 3> temp_matrix({n_ang, energy_groups, energy_groups}, 0.);
   read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);
 
-  // Get beta
-  xt::xtensor<double, 2> temp_beta({n_ang, delayed_groups}, 0.);
-  read_nd_vector(xsdata_grp, "beta", temp_beta, 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, delayed_groups}, 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});
+    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);
+    // 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);
 
-  // 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-prompt
+    chi_prompt = xt::view(1.0 - temp_beta_sum, xt::all(), xt::newaxis(),
+                          xt::newaxis()) * temp_matrix;
 
-  // 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());
+    // 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());
 
-  //Normalize both
+  } else if (beta_ndims == ndim_target + 1) {
+    xt::xtensor<double, 3> temp_beta({n_ang, delayed_groups, energy_groups}, 0.);
+    read_nd_vector(xsdata_grp, "beta", temp_beta, true);
+
+    xt::xtensor<double, 2> temp_beta_sum({n_ang, energy_groups}, 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 = chi_prompt / xt::view(xt::sum(chi_prompt, {2}), xt::all(),
                                      xt::all(), xt::newaxis());
 
@@ -335,7 +385,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups,
     } else {
       if (object_exists(xsdata_grp, "beta")) {
         fission_vector_beta_from_hdf5(xsdata_grp, n_ang, energy_groups,
-             delayed_groups);
+             delayed_groups, is_isotropic);
       } else {
         fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups,
              delayed_groups);
@@ -347,7 +397,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups,
     } else {
       if (object_exists(xsdata_grp, "beta")) {
         fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, energy_groups,
-             delayed_groups);
+             delayed_groups, is_isotropic);
       } else {
         fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups,
              delayed_groups);

tests/regression_tests/mg_basic_delayed/inputs_true.dat

@@ -4,11 +4,15 @@
   <cell id="2" material="2" region="2 -3" universe="0" />
   <cell id="3" material="3" region="3 -4" universe="0" />
   <cell id="4" material="4" region="4 -5" universe="0" />
+  <cell id="5" material="5" region="5 -6" universe="0" />
+  <cell id="6" material="6" region="6 -7" universe="0" />
   <surface boundary="reflective" coeffs="0.0" id="1" type="x-plane" />
-  <surface coeffs="232.3625" id="2" type="x-plane" />
-  <surface coeffs="464.725" id="3" type="x-plane" />
-  <surface coeffs="697.0875000000001" id="4" type="x-plane" />
-  <surface boundary="vacuum" coeffs="929.45" id="5" type="x-plane" />
+  <surface coeffs="154.90833333333333" id="2" type="x-plane" />
+  <surface coeffs="309.81666666666666" id="3" type="x-plane" />
+  <surface coeffs="464.725" id="4" type="x-plane" />
+  <surface coeffs="619.6333333333333" id="5" type="x-plane" />
+  <surface coeffs="774.5416666666666" id="6" type="x-plane" />
+  <surface boundary="vacuum" coeffs="929.45" id="7" type="x-plane" />
 </geometry>
 <?xml version='1.0' encoding='utf-8'?>
 <materials>
@@ -29,6 +33,14 @@
     <density units="macro" value="1.0" />
     <macroscopic name="mat_4" />
   </material>
+  <material id="5" name="vec group beta">
+    <density units="macro" value="1.0" />
+    <macroscopic name="mat_5" />
+  </material>
+  <material id="6" name="matrix group beta">
+    <density units="macro" value="1.0" />
+    <macroscopic name="mat_6" />
+  </material>
 </materials>
 <?xml version='1.0' encoding='utf-8'?>
 <settings>
@@ -38,7 +50,7 @@
   <inactive>5</inactive>
   <source strength="1.0">
     <space type="box">
-      <parameters>0.0 -1000.0 -1000.0 232.3625 1000.0 1000.0</parameters>
+      <parameters>0.0 -1000.0 -1000.0 154.90833333333333 1000.0 1000.0</parameters>
     </space>
   </source>
   <output>

tests/regression_tests/mg_basic_delayed/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-9.928704E-01 2.679667E-02
+1.003463E+00 2.173155E-02

tests/regression_tests/mg_basic_delayed/test.py

@@ -85,6 +85,29 @@ def create_library():
     mat_4.set_total(total)
     mg_cross_sections_file.add_xsdata(mat_4)
 
+    # Make the base data that uses chi & nu-fiss vectors with a group-wise beta
+    mat_5 = openmc.XSdata('mat_5', groups)
+    mat_5.order = 1
+    mat_5.num_delayed_groups = 2
+    mat_5.set_beta(np.stack([beta] * groups.num_groups))
+    mat_5.set_nu_fission(np.multiply(nu, fiss))
+    mat_5.set_absorption(absorption)
+    mat_5.set_scatter_matrix(scatter)
+    mat_5.set_total(total)
+    mat_5.set_chi(chi)
+    mg_cross_sections_file.add_xsdata(mat_5)
+
+    # Make a version that uses a nu-fission matrix with a group-wise beta
+    mat_6 = openmc.XSdata('mat_6', groups)
+    mat_6.order = 1
+    mat_6.num_delayed_groups = 2
+    mat_6.set_beta(np.stack([beta] * groups.num_groups))
+    mat_6.set_nu_fission(np.outer(np.multiply(nu, fiss), chi))
+    mat_6.set_absorption(absorption)
+    mat_6.set_scatter_matrix(scatter)
+    mat_6.set_total(total)
+    mg_cross_sections_file.add_xsdata(mat_6)
+
     # Write the file
     mg_cross_sections_file.export_to_hdf5('2g.h5')
 
@@ -99,8 +122,10 @@ class MGXSTestHarness(PyAPITestHarness):
 
 def test_mg_basic_delayed():
     create_library()
-    model = slab_mg(num_regions=4, mat_names=['vec beta', 'vec no beta',
-                                              'matrix beta', 'matrix no beta'])
+    model = slab_mg(num_regions=6, mat_names=['vec beta', 'vec no beta',
+                                              'matrix beta', 'matrix no beta',
+                                              'vec group beta',
+                                              'matrix group beta'])
 
     harness = PyAPITestHarness('statepoint.10.h5', model)
     harness.main()