Fixing issue with INF error values. Adding test case.

This commit is contained in:
Patrick Shriwise 2019-10-07 00:08:47 -05:00
parent 1d27972d50
commit f3fa1f2578
5 changed files with 76 additions and 26 deletions

View file

@ -336,7 +336,7 @@ class VolumeCalculation(object):
ur_elem = ET.SubElement(element, "upper_right")
ur_elem.text = ' '.join(str(x) for x in self.upper_right)
if self.threshold:
threshold_elem = ET.SubElement(element, "threshold")
threshold_elem.text = str(self.threshold)
threshold_elem.set("type", self.trigger_type)
trigger_elem = ET.SubElement(element, "threshold")
trigger_elem.set("type", self.trigger_type)
trigger_elem.set("threshold", str(self.threshold))
return element

View file

@ -69,19 +69,20 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
}
if (check_for_node(node, "threshold")) {
threshold_ = std::stod(get_node_value(node, "threshold"));
pugi::xml_node threshold_node = node.child("threshold");
threshold_ = std::stod(get_node_value(threshold_node, "threshold"));
if (threshold_ <= 0.0) {
std::stringstream msg;
msg << "Invalid error threshold " << threshold_ << " provided for a volume calculation.";
fatal_error(msg);
}
pugi::xml_node threshold_node = node.child("threshold");
std::string tmp = get_node_value(threshold_node, "type");
if (tmp == "variance") {
trigger_type_ = ThresholdType::VARIANCE;
} else if (tmp == "std_dev") {
trigger_type_ = ThresholdType::STD_DEV;
trigger_type_ = ThresholdType::STD_DEV;
} else if ( tmp == "rel_err") {
trigger_type_ = ThresholdType::REL_ERR;
} else {
@ -116,11 +117,11 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const {
while (true) {
// check maximum error value for all domains
max_val = -INFTY;
for (const auto& result : results) {
for (const auto& result : results) {
double val;
switch (trigger_type_) {
case ThresholdType::STD_DEV:
val = result.volume[1];
val = result.volume[1];
break;
case ThresholdType::REL_ERR:
val = result.volume[1] / result.volume[0];
@ -129,8 +130,9 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const {
val = result.volume[1] * result.volume[1];
break;
}
// update max
max_val = std::max(max_val, val);
// update max if entry is valid
if (val > 0.0) { max_val = std::max(max_val, val); }
}
// exit once we're below our error limit
@ -306,7 +308,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::_execute(size_t seed_o
for (int j = 0; j < master_indices[i_domain].size(); ++j) {
total_hits += master_hits[i_domain][j];
double f = static_cast<double>(master_hits[i_domain][j]) / n_samples_;
double var_f = f*(1.0 - f)/n_samples_;
double var_f = f*(1.0 - f) / n_samples_;
int i_material = master_indices[i_domain][j];
if (i_material == MATERIAL_VOID) continue;
@ -340,7 +342,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::_execute(size_t seed_o
} else {
result.nuclides.push_back(j);
result.atoms.push_back(0.0);
result.uncertainty.push_back(INFTY);
result.uncertainty.push_back(0.0);
}
}
}
@ -384,7 +386,7 @@ void VolumeCalculation::to_hdf5(const std::string& filename,
write_attribute(file_id, "trigger_type", "rel_err");
break;
}
}
if (domain_type_ == FILTER_CELL) {

View file

@ -48,4 +48,12 @@
<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.0001" type="std_dev" />
</volume_calc>
</settings>

View file

@ -2,15 +2,22 @@ Volume calculation 0
Domain 1: 31.47+/-0.07 cm^3
Domain 2: 2.093+/-0.031 cm^3
Domain 3: 2.049+/-0.031 cm^3
Cell Nuclide Atoms
0 1 U235 (3.482+/-0.008)e+23
1 1 Mo99 (3.482+/-0.008)e+22
2 2 H1 (1.400+/-0.021)e+23
3 2 O16 (7.00+/-0.10)e+22
4 2 B10 (7.00+/-0.10)e+18
5 3 H1 (1.370+/-0.021)e+23
6 3 O16 (6.85+/-0.10)e+22
7 3 B10 (6.85+/-0.10)e+18
Cell Nuclide Atoms
0 1 H1 0.0+/-0
1 1 O16 0.0+/-0
2 1 B10 0.0+/-0
3 1 U235 (3.482+/-0.008)e+23
4 1 Mo99 (3.482+/-0.008)e+22
5 2 H1 (1.400+/-0.021)e+23
6 2 O16 (7.00+/-0.10)e+22
7 2 B10 (7.00+/-0.10)e+18
8 2 U235 0.0+/-0
9 2 Mo99 0.0+/-0
10 3 H1 (1.370+/-0.021)e+23
11 3 O16 (6.85+/-0.10)e+22
12 3 B10 (6.85+/-0.10)e+18
13 3 U235 0.0+/-0
14 3 Mo99 0.0+/-0
Volume calculation 1
Domain 1: 4.14+/-0.04 cm^3
Domain 2: 31.47+/-0.07 cm^3
@ -18,8 +25,13 @@ Domain 2: 31.47+/-0.07 cm^3
0 1 H1 (2.770+/-0.029)e+23
1 1 O16 (1.385+/-0.014)e+23
2 1 B10 (1.385+/-0.014)e+19
3 2 U235 (3.482+/-0.008)e+23
4 2 Mo99 (3.482+/-0.008)e+22
3 1 U235 0.0+/-0
4 1 Mo99 0.0+/-0
5 2 H1 0.0+/-0
6 2 O16 0.0+/-0
7 2 B10 0.0+/-0
8 2 U235 (3.482+/-0.008)e+23
9 2 Mo99 (3.482+/-0.008)e+22
Volume calculation 2
Domain 0: 35.61+/-0.07 cm^3
Universe Nuclide Atoms
@ -28,3 +40,23 @@ Domain 0: 35.61+/-0.07 cm^3
2 0 B10 (1.385+/-0.014)e+19
3 0 U235 (3.482+/-0.008)e+23
4 0 Mo99 (3.482+/-0.008)e+22
Volume calculation 3
Domain 1: 31.35476+/-0.00010 cm^3
Domain 2: 2.09147+/-0.00005 cm^3
Domain 3: 2.110402+/-0.000033 cm^3
Cell Nuclide Atoms
0 1 H1 0.0+/-0
1 1 O16 0.0+/-0
2 1 B10 0.0+/-0
3 1 U235 (3.474110+/-0.000011)e+23
4 1 Mo99 (3.474110+/-0.000011)e+22
5 2 H1 (1.398833+/-0.000031)e+23
6 2 O16 (6.99416+/-0.00015)e+22
7 2 B10 (6.99416+/-0.00015)e+18
8 2 U235 0.0+/-0
9 2 Mo99 0.0+/-0
10 3 H1 (1.401947+/-0.000022)e+23
11 3 O16 (7.00974+/-0.00011)e+22
12 3 B10 (7.00974+/-0.00011)e+18
13 3 U235 0.0+/-0
14 3 Mo99 0.0+/-0

View file

@ -45,9 +45,12 @@ class VolumeTest(PyAPITestHarness):
vol_calcs = [
openmc.VolumeCalculation(list(root.cells.values()), 100000),
openmc.VolumeCalculation([water, fuel], 100000, ll, ur),
openmc.VolumeCalculation([root], 100000, ll, ur)
openmc.VolumeCalculation([root], 100000, ll, ur),
openmc.VolumeCalculation(list(root.cells.values()), 100)
]
vol_calcs[-1].set_trigger(1e-04, 'std_dev')
# Define settings
settings = openmc.Settings()
settings.run_mode = 'volume'
@ -62,6 +65,11 @@ class VolumeTest(PyAPITestHarness):
# Read volume calculation results
volume_calc = openmc.VolumeCalculation.from_hdf5(filename)
if volume_calc.samples == 100:
assert(volume_calc.trigger_type == 'std_dev')
assert(volume_calc.threshold == 1e-04)
# Write cell volumes and total # of atoms for each nuclide
for uid, volume in sorted(volume_calc.volumes.items()):
outstr += 'Domain {}: {} cm^3\n'.format(uid, volume)
@ -74,4 +82,4 @@ class VolumeTest(PyAPITestHarness):
def test_volume_calc():
harness = VolumeTest('')
harness._build_inputs()
harness.main()