Merge pull request #1889 from nelsonag/vol_calc_mg

Fix bug associated with volume calc in MG mode

src/material.cpp

@@ -349,24 +349,26 @@ Material::~Material()
 void Material::finalize()
 {
   // Set fissionable if any nuclide is fissionable
-  for (const auto& i_nuc : nuclide_) {
-    if (data::nuclides[i_nuc]->fissionable_) {
-      fissionable_ = true;
-      break;
+  if (settings::run_CE) {
+    for (const auto& i_nuc : nuclide_) {
+      if (data::nuclides[i_nuc]->fissionable_) {
+        fissionable_ = true;
+        break;
+      }
     }
+
+    // Generate material bremsstrahlung data for electrons and positrons
+    if (settings::photon_transport &&
+        settings::electron_treatment == ElectronTreatment::TTB) {
+      this->init_bremsstrahlung();
+    }
+
+    // Assign thermal scattering tables
+    this->init_thermal();
   }
 
-  // Generate material bremsstrahlung data for electrons and positrons
-  if (settings::photon_transport &&
-      settings::electron_treatment == ElectronTreatment::TTB) {
-    this->init_bremsstrahlung();
-  }
-
-  // Assign thermal scattering tables
-  this->init_thermal();
-
-  // Normalize density
-  this->normalize_density();
+// Normalize density
+this->normalize_density();
 }
 
 void Material::normalize_density()

src/mgxs_interface.cpp

@@ -123,10 +123,14 @@ void MgxsInterface::create_macro_xs()
   // Therefore type(nuclides[mat->nuclide_[0]]) dictates type(macroxs).
   // At the same time, we will find the scattering type, as that will dictate
   // how we allocate the scatter object within macroxs.
+
   for (int i = 0; i < model::materials.size(); ++i) {
+    // First we have to normalize the densities as it has not been called yet
+    // for MG mode
+    auto& mat {model::materials[i]};
+    mat->finalize();
     if (kTs[i].size() > 0) {
       // Convert atom_densities to a vector
-      auto& mat {model::materials[i]};
       vector<double> atom_densities(
         mat->atom_density_.begin(), mat->atom_density_.end());
 

src/volume_calc.cpp

@@ -8,6 +8,7 @@
 #include "openmc/hdf5_interface.h"
 #include "openmc/material.h"
 #include "openmc/message_passing.h"
+#include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
 #include "openmc/random_lcg.h"
@@ -237,7 +238,8 @@ vector<VolumeCalculation::Result> VolumeCalculation::execute() const
       // Create 2D array to store atoms/uncertainty for each nuclide. Later this
       // is compressed into vectors storing only those nuclides that are
       // non-zero
-      auto n_nuc = data::nuclides.size();
+      auto n_nuc = settings::run_CE ? data::nuclides.size()
+                                    : data::mg.nuclides_.size();
       xt::xtensor<double, 2> atoms({n_nuc, 2}, 0.0);
 
 #ifdef OPENMC_MPI
@@ -442,7 +444,8 @@ void VolumeCalculation::to_hdf5(
 
     vector<std::string> nucnames;
     for (int i_nuc : result.nuclides) {
-      nucnames.push_back(data::nuclides[i_nuc]->name_);
+      nucnames.push_back(settings::run_CE ? data::nuclides[i_nuc]->name_
+                                          : data::mg.nuclides_[i_nuc].name);
     }
 
     // Create array of total # of atoms with uncertainty for each nuclide

tests/regression_tests/volume_calc/inputs_true_mg.dat

@@ -0,0 +1,76 @@
+<?xml version='1.0' encoding='utf-8'?>
+<geometry>
+  <cell id="1" material="2" region="-1 -3 5" universe="0" />
+  <cell id="2" material="1" region="-2 3" universe="0" />
+  <cell id="3" material="1" region="-4 -3" universe="0" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="1" type="z-cylinder" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 5.0 1.0" id="2" type="sphere" />
+  <surface coeffs="5.0" id="3" type="z-plane" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 -5.0 1.0" id="4" type="sphere" />
+  <surface coeffs="-5.0" id="5" type="z-plane" />
+</geometry>
+<?xml version='1.0' encoding='utf-8'?>
+<materials>
+  <cross_sections>mg_lib.h5</cross_sections>
+  <material id="1">
+    <density units="g/cc" value="1.0" />
+    <nuclide ao="2.0" name="H1" />
+    <nuclide ao="1.0" name="O16" />
+    <nuclide ao="0.0001" name="B10" />
+  </material>
+  <material depletable="true" id="2">
+    <density units="g/cc" value="4.5" />
+    <nuclide ao="1.0" name="U235" />
+    <nuclide ao="0.1" name="Mo99" />
+  </material>
+</materials>
+<?xml version='1.0' encoding='utf-8'?>
+<settings>
+  <run_mode>volume</run_mode>
+  <energy_mode>multi-group</energy_mode>
+  <volume_calc>
+    <domain_type>cell</domain_type>
+    <domain_ids>1 2 3</domain_ids>
+    <samples>100000</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+  </volume_calc>
+  <volume_calc>
+    <domain_type>material</domain_type>
+    <domain_ids>1 2</domain_ids>
+    <samples>100000</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+  </volume_calc>
+  <volume_calc>
+    <domain_type>universe</domain_type>
+    <domain_ids>0</domain_ids>
+    <samples>100000</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+  </volume_calc>
+  <volume_calc>
+    <domain_type>cell</domain_type>
+    <domain_ids>1 2 3</domain_ids>
+    <samples>100</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+    <threshold threshold="0.1" type="std_dev" />
+  </volume_calc>
+  <volume_calc>
+    <domain_type>material</domain_type>
+    <domain_ids>1 2</domain_ids>
+    <samples>100</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+    <threshold threshold="0.1" type="rel_err" />
+  </volume_calc>
+  <volume_calc>
+    <domain_type>cell</domain_type>
+    <domain_ids>1 2 3</domain_ids>
+    <samples>100</samples>
+    <lower_left>-1.0 -1.0 -6.0</lower_left>
+    <upper_right>1.0 1.0 6.0</upper_right>
+    <threshold threshold="0.05" type="variance" />
+  </volume_calc>
+</settings>

tests/regression_tests/volume_calc/test.py

@@ -1,6 +1,8 @@
 import os
 import glob
-import sys
+
+import numpy as np
+import pytest
 
 import openmc
 
@@ -9,12 +11,15 @@ from tests.testing_harness import PyAPITestHarness
 
 class VolumeTest(PyAPITestHarness):
 
-    def __init__(self, *args, **kwargs):
+    def __init__(self, is_ce, *args, **kwargs):
         super().__init__(*args, **kwargs)
 
         self.exp_std_dev = 1e-01
         self.exp_rel_err = 1e-01
         self.exp_variance = 5e-02
+        self.is_ce = is_ce
+        if not is_ce:
+            self.inputs_true = 'inputs_true_mg.dat'
 
     def _build_inputs(self):
         # Define materials
@@ -22,7 +27,8 @@ class VolumeTest(PyAPITestHarness):
         water.add_nuclide('H1', 2.0)
         water.add_nuclide('O16', 1.0)
         water.add_nuclide('B10', 0.0001)
-        water.add_s_alpha_beta('c_H_in_H2O')
+        if self.is_ce:
+            water.add_s_alpha_beta('c_H_in_H2O')
         water.set_density('g/cc', 1.0)
 
         fuel = openmc.Material(2)
@@ -31,6 +37,8 @@ class VolumeTest(PyAPITestHarness):
         fuel.set_density('g/cc', 4.5)
 
         materials = openmc.Materials((water, fuel))
+        if not self.is_ce:
+            materials.cross_sections = 'mg_lib.h5'
         materials.export_to_xml()
 
         cyl = openmc.ZCylinder(surface_id=1, r=1.0, boundary_type='vacuum')
@@ -68,9 +76,52 @@ class VolumeTest(PyAPITestHarness):
         # Define settings
         settings = openmc.Settings()
         settings.run_mode = 'volume'
+        if not self.is_ce:
+            settings.energy_mode = 'multi-group'
         settings.volume_calculations = vol_calcs
         settings.export_to_xml()
 
+        # Create the MGXS file if necessary
+        if not self.is_ce:
+            groups = openmc.mgxs.EnergyGroups(group_edges=[0., 20.e6])
+            mg_xs_file = openmc.MGXSLibrary(groups)
+
+            nu = [2.]
+            fiss = [1.]
+            capture = [1.]
+            absorption_fissile = np.add(fiss, capture)
+            absorption_other = capture
+            scatter = np.array([[[1.]]])
+            total_fissile = np.add(absorption_fissile,
+                                   np.sum(scatter[:, :, 0], axis=1))
+            total_other = np.add(absorption_other,
+                                 np.sum(scatter[:, :, 0], axis=1))
+            chi = [1.]
+
+            for iso in ['H1', 'O16', 'B10', 'Mo99', 'U235']:
+                mat = openmc.XSdata(iso, groups)
+                mat.order = 0
+                mat.atomic_weight_ratio = \
+                    openmc.data.atomic_mass(iso) / openmc.data.NEUTRON_MASS
+                mat.set_scatter_matrix(scatter)
+                if iso == 'U235':
+                    mat.set_nu_fission(np.multiply(nu, fiss))
+                    mat.set_absorption(absorption_fissile)
+                    mat.set_total(total_fissile)
+                    mat.set_chi(chi)
+                else:
+                    mat.set_absorption(absorption_other)
+                    mat.set_total(total_other)
+                mg_xs_file.add_xsdata(mat)
+            mg_xs_file.export_to_hdf5('mg_lib.h5')
+
+    def _cleanup(self):
+        super()._cleanup()
+        output = ['mg_lib.h5']
+        for f in output:
+            if os.path.exists(f):
+                os.remove(f)
+
     def _get_results(self):
         outstr = ''
         for i, filename in enumerate(sorted(glob.glob('volume_*.h5'))):
@@ -116,6 +167,8 @@ class VolumeTest(PyAPITestHarness):
     def _test_output_created(self):
         pass
 
-def test_volume_calc():
-    harness = VolumeTest('')
+
+@pytest.mark.parametrize('is_ce', [True, False])
+def test_volume_calc(is_ce):
+    harness = VolumeTest(is_ce, '')
     harness.main()