Merge pull request #1949 from ojschumann/glibcxx_debug

Fixes for various bugs

include/openmc/hdf5_interface.h

@@ -247,11 +247,11 @@ inline void read_dataset(
 {
   // Create buffer to read data into
   auto n = dataset_typesize(obj_id, name);
-  char* buffer = new char[n];
+  std::vector<std::string::value_type> buffer(n, '\0');
 
   // Read attribute and set string
-  read_string(obj_id, name, n, buffer, indep);
-  str = std::string {buffer, n};
+  read_string(obj_id, name, n, buffer.data(), indep);
+  str = std::string {buffer.begin(), buffer.end()};
 }
 
 // array version

include/openmc/search.h

@@ -16,16 +16,14 @@ typename std::iterator_traits<It>::difference_type lower_bound_index(
 {
   if (*first == value)
     return 0;
-  It index = std::lower_bound(first, last, value) - 1;
-  return (index == last) ? -1 : index - first;
+  return std::lower_bound(first, last, value) - first - 1;
 }
 
 template<class It, class T>
 typename std::iterator_traits<It>::difference_type upper_bound_index(
   It first, It last, const T& value)
 {
-  It index = std::upper_bound(first, last, value) - 1;
-  return (index == last) ? -1 : index - first;
+  return std::upper_bound(first, last, value) - first - 1;
 }
 
 } // namespace openmc

src/distribution_energy.cpp

@@ -181,23 +181,9 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
     l = r > prn(seed) ? i + 1 : i;
   }
 
-  // Interpolation for energy E1 and EK
-  int n_energy_out = distribution_[i].e_out.size();
-  int n_discrete = distribution_[i].n_discrete;
-  double E_i_1 = distribution_[i].e_out[n_discrete];
-  double E_i_K = distribution_[i].e_out[n_energy_out - 1];
-
-  n_energy_out = distribution_[i + 1].e_out.size();
-  n_discrete = distribution_[i + 1].n_discrete;
-  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
-  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
-
-  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
-  double E_K = E_i_K + r * (E_i1_K - E_i_K);
-
   // Determine outgoing energy bin
-  n_energy_out = distribution_[l].e_out.size();
-  n_discrete = distribution_[l].n_discrete;
+  int n_energy_out = distribution_[l].e_out.size();
+  int n_discrete = distribution_[l].n_discrete;
   double r1 = prn(seed);
   double c_k = distribution_[l].c[0];
   int k = 0;
@@ -257,6 +243,20 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const
 
   // Now interpolate between incident energy bins i and i + 1
   if (!histogram_interp && n_energy_out > 1 && k >= n_discrete) {
+    // Interpolation for energy E1 and EK
+    n_energy_out = distribution_[i].e_out.size();
+    n_discrete = distribution_[i].n_discrete;
+    const double E_i_1 = distribution_[i].e_out[n_discrete];
+    const double E_i_K = distribution_[i].e_out[n_energy_out - 1];
+
+    n_energy_out = distribution_[i + 1].e_out.size();
+    n_discrete = distribution_[i + 1].n_discrete;
+    const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
+    const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];
+
+    const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
+    const double E_K = E_i_K + r * (E_i1_K - E_i_K);
+
     if (l == i) {
       return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
     } else {

src/mesh.cpp

@@ -158,7 +158,9 @@ std::string StructuredMesh::bin_label(int bin) const
 
 xt::xtensor<int, 1> StructuredMesh::get_x_shape() const
 {
-  return xt::adapt(shape_, {n_dimension_});
+  // because method is const, shape_ is const as well and can't be adapted
+  auto tmp_shape = shape_;
+  return xt::adapt(tmp_shape, {n_dimension_});
 }
 
 //==============================================================================

src/nuclide.cpp

@@ -455,7 +455,7 @@ void Nuclide::create_derived(
         xs_cdf_sum +=
           (std::sqrt(E[i]) * xs[i] + std::sqrt(E[i + 1]) * xs[i + 1]) / 2.0 *
           (E[i + 1] - E[i]);
-        xs_cdf_[i] = xs_cdf_sum;
+        xs_cdf_[i+1] = xs_cdf_sum;
       }
     }
   }

src/photon.cpp

@@ -220,7 +220,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   // Create Compton profile CDF
   auto n_profile = data::compton_profile_pz.size();
   profile_cdf_ = xt::empty<double>({n_shell, n_profile});
-  for (int i = 0; i < n_shell; ++i) {
+  for (int i = 0; i < profile_pdf_.shape(0); ++i) {
     double c = 0.0;
     profile_cdf_(i, 0) = 0.0;
     for (int j = 0; j < n_profile - 1; ++j) {

src/physics.cpp

@@ -936,18 +936,18 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
     } else if (sampling_method == ResScatMethod::rvs) {
       // interpolate xs CDF since we're not exactly at the energy indices
       // cdf value at lower bound attainable energy
-      double m = (nuc.xs_cdf_[i_E_low] - nuc.xs_cdf_[i_E_low - 1]) /
-                 (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]);
-      double cdf_low =
-        nuc.xs_cdf_[i_E_low - 1] + m * (E_low - nuc.energy_0K_[i_E_low]);
-      if (E_low <= nuc.energy_0K_.front())
-        cdf_low = 0.0;
+      double cdf_low = 0.0;
+      if (E_low > nuc.energy_0K_.front()) {
+        double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) /
+                   (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]);
+        cdf_low = nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]);
+      }
 
       // cdf value at upper bound attainable energy
-      m = (nuc.xs_cdf_[i_E_up] - nuc.xs_cdf_[i_E_up - 1]) /
+      double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) /
           (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
       double cdf_up =
-        nuc.xs_cdf_[i_E_up - 1] + m * (E_up - nuc.energy_0K_[i_E_up]);
+        nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
 
       while (true) {
         // directly sample Maxwellian
@@ -955,14 +955,15 @@ Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
 
         // sample a relative energy using the xs cdf
         double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low);
-        int i_E_rel = lower_bound_index(
-          &nuc.xs_cdf_[i_E_low - 1], &nuc.xs_cdf_[i_E_up + 1], cdf_rel);
+        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low, 
+                                        nuc.xs_cdf_.begin() + i_E_up+2, 
+                                        cdf_rel);
         double E_rel = nuc.energy_0K_[i_E_low + i_E_rel];
-        double m = (nuc.xs_cdf_[i_E_low + i_E_rel] -
-                     nuc.xs_cdf_[i_E_low + i_E_rel - 1]) /
+        double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] -
+                     nuc.xs_cdf_[i_E_low + i_E_rel]) /
                    (nuc.energy_0K_[i_E_low + i_E_rel + 1] -
                      nuc.energy_0K_[i_E_low + i_E_rel]);
-        E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel - 1]) / m;
+        E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m;
 
         // perform rejection sampling on cosine between
         // neutron and target velocities

src/plot.cpp

@@ -745,7 +745,7 @@ void output_png(Plot const& pl, const ImageData& data)
   // Clean up data structures
   std::fclose(fp);
   png_free_data(png_ptr, info_ptr, PNG_FREE_ALL, -1);
-  png_destroy_write_struct(&png_ptr, nullptr);
+  png_destroy_write_struct(&png_ptr, &info_ptr);
 }
 #endif
 

src/secondary_thermal.cpp

@@ -4,6 +4,8 @@
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
+#include <gsl/gsl-lite.hpp>
+
 #include "xtensor/xview.hpp"
 
 #include <cmath> // for log, exp
@@ -19,7 +21,8 @@ void get_energy_index(
   f = 0.0;
   if (E >= energies.front()) {
     i = lower_bound_index(energies.begin(), energies.end(), E);
-    f = (E - energies[i]) / (energies[i + 1] - energies[i]);
+    if (i + 1 < energies.size())
+      f = (E - energies[i]) / (energies[i + 1] - energies[i]);
   }
 }
 
@@ -32,25 +35,25 @@ CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs} {}
 void CoherentElasticAE::sample(
   double E_in, double& E_out, double& mu, uint64_t* seed) const
 {
-  // Get index and interpolation factor for elastic grid
-  int i;
-  double f;
+  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
+  E_out = E_in;
+
   const auto& energies {xs_.bragg_edges()};
-  get_energy_index(energies, E_in, i, f);
+
+  Expects(E_in >= energies.front());
+
+  const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
 
   // Sample a Bragg edge between 1 and i
+  // E[0] < E_in < E[i+1] -> can scatter in bragg edges 0..i
   const auto& factors = xs_.factors();
-  double prob = prn(seed) * factors[i + 1];
-  int k = 0;
-  if (prob >= factors.front()) {
-    k = lower_bound_index(factors.begin(), factors.begin() + (i + 1), prob);
-  }
+  const double prob = prn(seed) * factors[i];
+
+  const int k = std::lower_bound(factors.begin(), factors.begin() + i, prob) -
+                factors.begin();
 
   // Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
   mu = 1.0 - 2.0 * energies[k] / E_in;
-
-  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
-  E_out = E_in;
 }
 
 //==============================================================================

src/state_point.cpp

@@ -720,7 +720,9 @@ std::string dtype_member_names(hid_t dtype_id)
   int nmembers = H5Tget_nmembers(dtype_id);
   std::string names;
   for (int i = 0; i < nmembers; i++) {
-    names = names.append(H5Tget_member_name(dtype_id, i));
+    char* name = H5Tget_member_name(dtype_id, i);
+    names = names.append(name);
+    H5free_memory(name);
     if (i < nmembers - 1)
       names += ", ";
   }

src/tallies/tally_scoring.cpp

@@ -1509,7 +1509,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index,
   }
 
   // For shorthand, assign pointers to the material and nuclide xs set
-  auto& nuc_xs = data::mg.nuclides_[i_nuclide];
+  auto& nuc_xs = (i_nuclide >= 0) ? data::mg.nuclides_[i_nuclide] : data::mg.macro_xs_[p.material()];
   auto& macro_xs = data::mg.macro_xs_[p.material()];
 
   // Find the temperature and angle indices of interest

tests/regression_tests/salphabeta/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-8.214162E-01 2.262777E-02
+8.214164E-01 2.262776E-02

tests/regression_tests/triso/results_true.dat

@@ -1,2 +1,2 @@
 k-combined:
-1.653488E+00 2.128233E-02
+1.406055E+00 9.581396E-02