diff --git a/include/openmc/constants.h b/include/openmc/constants.h index df13da3707d..779243b2aba 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -324,8 +324,15 @@ enum TallyScore { }; // Global tally parameters -constexpr int N_GLOBAL_TALLIES {4}; -enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE }; +enum class GlobalTally { + K_COLLISION, + K_ABSORPTION, + K_TRACKLENGTH, + LEAKAGE, + K_TRACKLENGTH_SQ, // for calculation of stddev for generational k + SIZE +}; +constexpr int N_GLOBAL_TALLIES {static_cast(GlobalTally::SIZE)}; // Miscellaneous constexpr int C_NONE {-1}; diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b456fee21eb..d2bb462773e 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -21,7 +21,8 @@ namespace openmc { namespace simulation { -extern double keff_generation; //!< Single-generation k on each processor +extern array + keff_generation; //!< Single-generation k on each processor extern array k_sum; //!< Used to reduce sum and sum_sq extern vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 3e4e24e1d06..b81055799bf 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/array.h" #include "openmc/mesh.h" #include "openmc/particle.h" #include "openmc/vector.h" @@ -45,7 +46,7 @@ extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; -extern vector k_generation; +extern vector> k_generation; extern vector work_index; } // namespace simulation diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index fb1aaf7b985..00df7d4419e 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -49,6 +49,8 @@ void write_source_bank(hid_t group_id, span source_bank, void read_source_bank( hid_t group_id, vector& sites, bool distribute); +void write_global_tallies(hid_t file_id); +void read_global_tallies(hid_t file_id); void write_tally_results_nr(hid_t file_id); void restart_set_keff(); void write_unstructured_mesh_results(); diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 3beeb9d5ac5..fbd4f8ae1e6 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -225,6 +225,7 @@ extern "C" int32_t n_realizations; extern double global_tally_absorption; extern double global_tally_collision; extern double global_tally_tracklength; +extern double global_tally_tracklength_sq; extern double global_tally_leakage; //============================================================================== diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index c17b16597f9..fe8e83e3e2e 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -127,11 +127,11 @@ def global_tallies(): """ ptr = POINTER(c_double)() _dll.openmc_global_tallies(ptr) - array = as_array(ptr, (4, 3)) + array = as_array(ptr, (5, 3)) # Get sum, sum-of-squares, and number of realizations - sum_ = array[:, 1] - sum_sq = array[:, 2] + sum_ = array[:4, 1] + sum_sq = array[:4, 2] n = num_realizations() # Determine mean diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 715becf4882..e02d6ab60db 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -8,6 +8,7 @@ import numpy as np from pathlib import Path from uncertainties import ufloat +from uncertainties.unumpy import uarray import openmc import openmc.checkvalue as cv @@ -80,6 +81,8 @@ class StatePoint: .. versionadded:: 0.13.1 meshes : dict Dictionary whose keys are mesh IDs and whose values are MeshBase objects + multiplication : uncertainties.UFloat + Estimator for multiplication of initial source neutrons n_batches : int Number of batches n_inactive : int @@ -106,6 +109,8 @@ class StatePoint: 'E', 'wgt', 'delayed_group', 'surf_id', and 'particle', corresponding to the position, direction, energy, weight, delayed group, surface ID and particle type of the source site, respectively. + source_effectiveness : uncertainties.UFloat + Source effectiveness factor. source_present : bool Indicate whether source sites are present sparse : bool @@ -242,8 +247,8 @@ def global_tallies(self): ('mean', 'f8'), ('std_dev', 'f8')]) gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength', 'leakage'] - gt['sum'] = data[:,1] - gt['sum_sq'] = data[:,2] + gt['sum'] = data[:,0] + gt['sum_sq'] = data[:,1] # Calculate mean and sample standard deviation of mean n = self.n_realizations @@ -264,7 +269,13 @@ def k_cmfd(self): @property def k_generation(self): if self.run_mode == 'eigenvalue': - return self._f['k_generation'][()] + arr = self._f['k_generation'][()] + if arr.ndim==1: + return arr + elif arr.ndim==2: + return uarray(arr[:,0],arr[:,1]) + else: + raise ValueError(f'k_generation shape ({arr.shape}) must be either 1d or 2d') else: return None @@ -317,6 +328,17 @@ def meshes(self): self._meshes_read = True return self._meshes + + @property + def multiplication(self): + if self.run_mode == 'eigenvalue': + k_gen = self.k_generation + k_inactive = k_gen[:self.n_inactive] + cumprod = np.cumprod(k_inactive) + cumprod[-1] *= 1/(1-self.keff) + return 1.0+np.sum(cumprod) + else: + return None @property def n_batches(self): @@ -365,6 +387,15 @@ def stride(self): @property def source(self): return self._f['source_bank'][()] if self.source_present else None + + @property + def source_efficiency(self): + if self.run_mode == 'eigenvalue': + k_gen = self.k_generation + k_inactive = k_gen[:self.n_inactive] + return np.prod(k_inactive/self.keff) + else: + return None @property def source_present(self): diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index a2120a006d5..df2d1f896ab 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -36,7 +36,7 @@ namespace openmc { namespace simulation { -double keff_generation; +array keff_generation; array k_sum; vector entropy; xt::xtensor source_frac; @@ -52,16 +52,21 @@ void calculate_generation_keff() const auto& gt = simulation::global_tallies; // Get keff for this generation by subtracting off the starting value - simulation::keff_generation = + simulation::keff_generation[0] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - - simulation::keff_generation; + simulation::keff_generation[0]; + simulation::keff_generation[1] = + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) - + simulation::keff_generation[1]; - double keff_reduced; + array keff_reduced; #ifdef OPENMC_MPI if (settings::solver_type != SolverType::RANDOM_RAY) { // Combine values across all processors - MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE, - MPI_SUM, mpi::intracomm); + MPI_Allreduce(&simulation::keff_generation[0], &keff_reduced[0], 1, + MPI_DOUBLE, MPI_SUM, mpi::intracomm); + MPI_Allreduce(&simulation::keff_generation[1], &keff_reduced[1], 1, + MPI_DOUBLE, MPI_SUM, mpi::intracomm); } else { // If using random ray, MPI parallelism is provided by domain replication. // As such, all fluxes will be reduced at the end of each transport sweep, @@ -77,10 +82,13 @@ void calculate_generation_keff() // Normalize single batch estimate of k // TODO: This should be normalized by total_weight, not by n_particles if (settings::solver_type != SolverType::RANDOM_RAY) { - keff_reduced /= settings::n_particles; + keff_reduced[0] /= settings::n_particles; + keff_reduced[1] /= settings::n_particles; } - - simulation::k_generation.push_back(keff_reduced); + double k_mean = keff_reduced[0]; + double k_std = std::sqrt( + (keff_reduced[1] - std::pow(k_mean, 2)) / (settings::n_particles - 1)); + simulation::k_generation.push_back({k_mean, k_std}); } void synchronize_bank() @@ -379,11 +387,11 @@ void calculate_average_keff() if (n <= 0) { // For inactive generations, use current generation k as estimate for next // generation - simulation::keff = simulation::k_generation[i]; + simulation::keff = simulation::k_generation[i][0]; } else { // Sample mean of keff - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i][0]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i][0], 2); // Determine mean simulation::keff = simulation::k_sum[0] / n; @@ -660,7 +668,13 @@ void write_eigenvalue_hdf5(hid_t group) { write_dataset(group, "n_inactive", settings::n_inactive); write_dataset(group, "generations_per_batch", settings::gen_per_batch); - write_dataset(group, "k_generation", simulation::k_generation); + auto n = simulation::k_generation.size(); + xt::xtensor k_generation({n, 2}); + for (int i = 0; i < n; ++i) { + k_generation(i, 0) = simulation::k_generation[i][0]; + k_generation(i, 1) = simulation::k_generation[i][1]; + } + write_dataset(group, "k_generation", k_generation); if (settings::entropy_on) { write_dataset(group, "entropy", simulation::entropy); } @@ -675,9 +689,13 @@ void write_eigenvalue_hdf5(hid_t group) void read_eigenvalue_hdf5(hid_t group) { read_dataset(group, "generations_per_batch", settings::gen_per_batch); - int n = simulation::restart_batch * settings::gen_per_batch; + size_t n = simulation::restart_batch * settings::gen_per_batch; + xt::xtensor k_generation({n, 2}); + read_dataset(group, "k_generation", k_generation); simulation::k_generation.resize(n); - read_dataset(group, "k_generation", simulation::k_generation); + for (int i = 0; i < n; ++i) { + simulation::k_generation[i] = {k_generation(i, 0), k_generation(i, 1)}; + } if (settings::entropy_on) { read_dataset(group, "entropy", simulation::entropy); } diff --git a/src/output.cpp b/src/output.cpp index 0a14e8843d8..852509ec525 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -19,6 +19,7 @@ #endif #include "xtensor/xview.hpp" +#include "openmc/array.h" #include "openmc/capi.h" #include "openmc/cell.h" #include "openmc/constants.h" @@ -370,11 +371,12 @@ void print_build_info() void print_columns() { if (settings::entropy_on) { - fmt::print(" Bat./Gen. k Entropy Average k \n" - " ========= ======== ======== ====================\n"); + fmt::print( + " Bat./Gen. k Entropy Average k \n" + " ========= ==================== ======== ====================\n"); } else { - fmt::print(" Bat./Gen. k Average k\n" - " ========= ======== ====================\n"); + fmt::print(" Bat./Gen. k Average k\n" + " ========= ==================== ====================\n"); } } @@ -392,7 +394,10 @@ void print_generation() // write out batch/generation and generation k-effective auto batch_and_gen = std::to_string(simulation::current_batch) + "/" + std::to_string(simulation::current_gen); - fmt::print(" {:>9} {:8.5f}", batch_and_gen, simulation::k_generation[idx]); + array k_gen; + k_gen[0] = simulation::k_generation[idx][0]; + k_gen[1] = simulation::k_generation[idx][1]; + fmt::print(" {:>9} {:8.5f} +/-{:8.5f}", batch_and_gen, k_gen[0], k_gen[1]); // write out entropy info if (settings::entropy_on) { diff --git a/src/particle.cpp b/src/particle.cpp index f5ad45d80fe..9a724c62329 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -478,8 +478,11 @@ void Particle::event_death() global_tally_absorption += keff_tally_absorption(); #pragma omp atomic global_tally_collision += keff_tally_collision(); + double val = keff_tally_tracklength(); #pragma omp atomic - global_tally_tracklength += keff_tally_tracklength(); + global_tally_tracklength += val; +#pragma omp atomic + global_tally_tracklength_sq += val * val; #pragma omp atomic global_tally_leakage += keff_tally_leakage(); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b840..53c675c7113 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -498,6 +498,7 @@ void RandomRaySimulation::simulate() // Store random ray k-eff into OpenMC's native k-eff variable global_tally_tracklength = k_eff_; + global_tally_tracklength_sq = std::pow(k_eff_, 2); } // Execute all tallying tasks, if this is an active batch diff --git a/src/simulation.cpp b/src/simulation.cpp index b9986f47375..5ce33112945 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,5 +1,6 @@ #include "openmc/simulation.h" +#include "openmc/array.h" #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/container_util.h" @@ -318,7 +319,7 @@ int64_t work_per_rank; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; -vector k_generation; +vector> k_generation; vector work_index; } // namespace simulation @@ -506,8 +507,10 @@ void initialize_generation() ufs_count_sites(); // Store current value of tracklength k - simulation::keff_generation = simulation::global_tallies( - GlobalTally::K_TRACKLENGTH, TallyResult::VALUE); + auto& gt = simulation::global_tallies; + simulation::keff_generation = { + gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE), + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE)}; } } @@ -522,6 +525,8 @@ void finalize_generation() global_tally_absorption; gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength; + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) += + global_tally_tracklength_sq; } gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; @@ -530,6 +535,7 @@ void finalize_generation() global_tally_collision = 0.0; global_tally_absorption = 0.0; global_tally_tracklength = 0.0; + global_tally_tracklength_sq = 0.0; } global_tally_leakage = 0.0; diff --git a/src/state_point.cpp b/src/state_point.cpp index 8195c486507..da70d72604b 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -253,7 +253,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) if (settings::reduce_tallies) { // Write global tallies - write_dataset(file_id, "global_tallies", simulation::global_tallies); + write_global_tallies(file_id); // Write tallies if (model::active_tallies.size() > 0) { @@ -354,13 +354,13 @@ void restart_set_keff() { if (simulation::restart_batch > settings::n_inactive) { for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i][0]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i][0], 2); } int n = settings::gen_per_batch * simulation::n_realizations; simulation::keff = simulation::k_sum[0] / n; } else { - simulation::keff = simulation::k_generation.back(); + simulation::keff = simulation::k_generation.back()[0]; } } @@ -484,8 +484,7 @@ extern "C" int openmc_statepoint_load(const char* filename) if (mpi::master) { #endif // Read global tally data - read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE, H5S_ALL, - false, simulation::global_tallies.data()); + read_global_tallies(file_id); // Check if tally results are present bool present; @@ -917,6 +916,34 @@ void write_unstructured_mesh_results() } } } +void write_global_tallies(hid_t file_id) +{ + // Get global tallies + auto& gt = simulation::global_tallies; + auto write_view = xt::view(gt, + xt::range(static_cast(GlobalTally::K_COLLISION), + static_cast(GlobalTally::LEAKAGE) + 1), + xt::range(static_cast(TallyResult::SUM), + static_cast(TallyResult::SUM_SQ) + 1)); + auto gt_reduced = xt::empty_like(write_view); + gt_reduced = write_view; + write_dataset(file_id, "global_tallies", gt_reduced); +} + +void read_global_tallies(hid_t file_id) +{ + // Get global tallies + auto& gt = simulation::global_tallies; + auto gt_view = xt::view(gt, + xt::range(static_cast(GlobalTally::K_COLLISION), + static_cast(GlobalTally::LEAKAGE) + 1), + xt::range(static_cast(TallyResult::SUM), + static_cast(TallyResult::SUM_SQ) + 1)); + auto gt_reduced = xt::empty_like(gt_view); + read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE, H5S_ALL, + false, gt_reduced.data()); + gt_view = gt_reduced; +} void write_tally_results_nr(hid_t file_id) { @@ -951,7 +978,7 @@ void write_tally_results_nr(hid_t file_id) // Write out global tallies sum and sum_sq if (mpi::master) { - write_dataset(file_id, "global_tallies", gt); + write_global_tallies(file_id); } for (const auto& t : model::tallies) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ae0bffe6eb0..5da73e7565d 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -76,6 +76,7 @@ int32_t n_realizations {0}; double global_tally_absorption; double global_tally_collision; double global_tally_tracklength; +double global_tally_tracklength_sq; double global_tally_leakage; //============================================================================== diff --git a/tests/regression_tests/multiplication/__init__.py b/tests/regression_tests/multiplication/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/multiplication/inputs_true.dat b/tests/regression_tests/multiplication/inputs_true.dat new file mode 100644 index 00000000000..b98a95ebb32 --- /dev/null +++ b/tests/regression_tests/multiplication/inputs_true.dat @@ -0,0 +1,36 @@ + + + + mgxs.h5 + + + + + + + + + + + + eigenvalue + 100000 + 20 + 10 + + + 0 -1000 -1000 10 1000 1000 + + + 10000000.0 1.0 + + + + false + + multi-group + + false + + + diff --git a/tests/regression_tests/multiplication/results_true.dat b/tests/regression_tests/multiplication/results_true.dat new file mode 100644 index 00000000000..b8a2310dfd5 --- /dev/null +++ b/tests/regression_tests/multiplication/results_true.dat @@ -0,0 +1,10 @@ +k-combined: +8.271126E-01 4.868738E-04 +multiplication: +8.963514E+00 4.210376E-02 +k1: +1.378182E+00 4.087893E-03 +kf: +8.269380E-01 1.049151E-03 +g*: +1.663722E+00 1.856058E-02 diff --git a/tests/regression_tests/multiplication/test.py b/tests/regression_tests/multiplication/test.py new file mode 100644 index 00000000000..622783a8d52 --- /dev/null +++ b/tests/regression_tests/multiplication/test.py @@ -0,0 +1,104 @@ +"""This test is based on a simple 4-group slab model from +"MCNP Calculations of Subcritical Fixed and Fission Multiplication Factors", +LA-UR-10-00141 which can be found at https://mcnp.lanl.gov/pdf_files/TechReport_2010_LANL_LA-UR-10-00141_KiedrowskiBrown.pdf +""" +import openmc +from openmc.stats import delta_function +import numpy as np +import pytest +from openmc.examples import slab_mg +import os +import glob + +from tests.testing_harness import PyAPITestHarness + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + def _get_results(self, hash_output=False): + outstr = super()._get_results(hash_output=hash_output) + # Read the statepoint file. + statepoint = glob.glob(self._sp_name)[0] + with openmc.StatePoint(statepoint) as sp: + # Write out multiplication. + outstr += 'multiplication:\n' + form = '{0:12.6E} {1:12.6E}\n' + M = sp.multiplication + outstr += form.format(M.n, M.s) + + # Write out k1. + outstr += 'k1:\n' + form = '{0:12.6E} {1:12.6E}\n' + k1 = sp.k_generation[0] + outstr += form.format(k1.n, k1.s) + + # Write out kf. + outstr += 'kf:\n' + form = '{0:12.6E} {1:12.6E}\n' + kf = 1-k1/(M-1) + outstr += form.format(kf.n, kf.s) + + # Write out g*. + outstr += 'g*:\n' + form = '{0:12.6E} {1:12.6E}\n' + g_star = sp.source_efficiency + outstr += form.format(g_star.n, g_star.s) + return outstr + + +@pytest.fixture() +def slab_model(): + model = slab_mg(mgxslib_name='mgxs.h5') + right_boundary = model.geometry.get_all_surfaces()[2] + right_boundary.coefficients['x0'] = 10.0 + + cell = model.geometry.get_all_cells()[1] + + mat = model.geometry.get_all_materials()[1] + mat.set_density('macro', 0.01) + ########################################################################### + # Create multigroup data + + # Instantiate the energy group data + ebins = np.geomspace(1e-5, 20.0e6, 5) + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + nusigma_f = np.array([9.6,5.4,5.2,2.5]) + sigma_s = np.array([[0.5,0.5,0.5,0.5], + [0.0,1.0,0.5,0.5], + [0.0,0.0,1.5,0.5], + [0.0,0.0,0.0,2.0]]) + sigma_t = np.array([5.0,5.0,5.0,5.0]) + sigma_a = sigma_t - sigma_s.sum(axis=1) + chi = np.array([0.0,0.2,0.8,0.0]) + mat_data = openmc.XSdata('mat_1', groups) + mat_data.order = 0 + mat_data.set_total(sigma_t) + mat_data.set_nu_fission(nusigma_f) + mat_data.set_chi(chi) + mat_data.set_absorption(sigma_a) + mat_data.set_scatter_matrix(sigma_s[...,np.newaxis]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdata(mat_data) + mg_cross_sections_file.export_to_hdf5() + + # Settings + model.settings.particles = 100000 + model.settings.batches = 20 + model.settings.inactive = 10 + + space = openmc.stats.Box([0,-1000,-1000],[10,1000,1000]) + model.settings.source = openmc.IndependentSource( + space=space, energy = delta_function(10e6)) + + return model + + +def test_multiplication(slab_model): + harness = MGXSTestHarness("statepoint.20.h5", model=slab_model) + harness.main() diff --git a/tests/regression_tests/source_effectiveness/__init__.py b/tests/regression_tests/source_effectiveness/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/source_effectiveness/inputs_true.dat b/tests/regression_tests/source_effectiveness/inputs_true.dat new file mode 100644 index 00000000000..94e7fdcd963 --- /dev/null +++ b/tests/regression_tests/source_effectiveness/inputs_true.dat @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + eigenvalue + 100000 + 20 + 10 + + + 0.0 0.0 0.0 + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/source_effectiveness/results_true.dat b/tests/regression_tests/source_effectiveness/results_true.dat new file mode 100644 index 00000000000..b4980e9290d --- /dev/null +++ b/tests/regression_tests/source_effectiveness/results_true.dat @@ -0,0 +1,10 @@ +k-combined: +9.960995E-01 5.611810E-04 +multiplication: +3.646446E+02 5.066747E+01 +k1: +1.185417E+00 2.963267E-03 +kf: +9.967402E-01 4.542722E-04 +g*: +1.425744E+00 1.461587E-02 diff --git a/tests/regression_tests/source_effectiveness/test.py b/tests/regression_tests/source_effectiveness/test.py new file mode 100644 index 00000000000..28cd75d00c1 --- /dev/null +++ b/tests/regression_tests/source_effectiveness/test.py @@ -0,0 +1,95 @@ +"""This model is based on a numerical example at doi:10.1016/s0306-4549(98)00048-6 +""" +import openmc +from openmc.stats import spherical_uniform, Point, Mixture, Watt +import numpy as np +import pytest +import os +import glob + +from tests.testing_harness import PyAPITestHarness + + +class SourceTestHarness(PyAPITestHarness): + def _get_results(self, hash_output=False): + outstr = super()._get_results(hash_output=hash_output) + # Read the statepoint file. + statepoint = glob.glob(self._sp_name)[0] + with openmc.StatePoint(statepoint) as sp: + # Write out multiplication. + outstr += 'multiplication:\n' + form = '{0:12.6E} {1:12.6E}\n' + M = sp.multiplication + outstr += form.format(M.n, M.s) + + # Write out k1. + outstr += 'k1:\n' + form = '{0:12.6E} {1:12.6E}\n' + k1 = sp.k_generation[0] + outstr += form.format(k1.n, k1.s) + + # Write out kf. + outstr += 'kf:\n' + form = '{0:12.6E} {1:12.6E}\n' + kf = 1-k1/(M-1) + outstr += form.format(kf.n, kf.s) + + # Write out g*. + outstr += 'g*:\n' + form = '{0:12.6E} {1:12.6E}\n' + g_star = sp.source_efficiency + outstr += form.format(g_star.n, g_star.s) + return outstr + +@pytest.fixture() +def sources(): + point = Point((0.0,0.0,0.0)) + uniform = spherical_uniform(r_outer=8.85) + cf252_spec = Watt(1.0250e6,2.926e-6) + u235_spec = Watt(0.7747e6,4.852e-6) + u238_spec = Watt(0.6483e6,6.811e-6) + cf252 = openmc.IndependentSource(space=point, + energy=cf252_spec, + strength=100) + u235 = openmc.IndependentSource(space=uniform, + energy=u235_spec, + strength=0.5) + u238 = openmc.IndependentSource(space=uniform, + energy=u238_spec, + strength=57.3) + return [cf252,u235,u238] + + +@pytest.fixture() +def sphere_model(sources): + + mat = openmc.Material(name='heu') + mat.add_element('U',1,percent_type="ao",enrichment=92.2) + mat.set_density('g/cm3',18.6) + + materials = openmc.Materials([mat]) + + + + surface = openmc.Sphere(r=8.85,boundary_type='vacuum') + cell = openmc.Cell(fill=mat,region=-surface) + + root = openmc.Universe() + root.add_cell(cell) + + geometry = openmc.Geometry(root) + + model = openmc.Model(geometry=geometry,materials=materials) + + # Settings + model.settings.particles = 100000 + model.settings.batches = 20 + model.settings.inactive = 10 + model.settings.source = sources + + return model + + +def test_source_effectiveness(sphere_model): + harness = SourceTestHarness("statepoint.20.h5", model=sphere_model) + harness.main()