diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 38238598e..e62b2941c 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -366,7 +366,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY }; enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID }; enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY }; -enum class RandomRaySampleMethod { PRNG, HALTON }; +enum class RandomRaySampleMethod { PRNG, HALTON, S2 }; //============================================================================== // Geometry Constants diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index 7fd68b83e..b61d2d67a 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -41,6 +41,7 @@ public: uint64_t transport_history_based_single_ray(); SourceSite sample_prng(); SourceSite sample_halton(); + SourceSite sample_s2(); //---------------------------------------------------------------------------- // Static data members diff --git a/openmc/settings.py b/openmc/settings.py index 3cf662e31..2ad1d06e7 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -207,7 +207,11 @@ class Settings: default is 'False'. :sample_method: Sampling method for the ray starting location and direction of - travel. Options are `prng` (default) or 'halton`. + travel. Options are `prng` (default), `halton`, or `s2`. `s2` + modifies the `prng` sampling method such that rays are sampled + with directions (-1, 0, 0) or (1, 0, 0). This is used for verification + against analytic transport benchmarks which are often derivied with + a reduced angular domain. :source_region_meshes: List of tuples where each tuple contains a mesh and a list of domains. Each domain is an instance of openmc.Material, openmc.Cell, @@ -1351,7 +1355,7 @@ class Settings: 'openmc.Material, openmc.Cell, or openmc.Universe.') elif key == 'sample_method': cv.check_value('sample method', value, - ('prng', 'halton')) + ('prng', 'halton', 's2')) elif key == 'diagonal_stabilization_rho': cv.check_type('diagonal stabilization rho', value, Real) cv.check_greater_than('diagonal stabilization rho', diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 44e161a8e..821c5d1ba 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -792,6 +792,9 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) case RandomRaySampleMethod::HALTON: site = sample_halton(); break; + case RandomRaySampleMethod::S2: + site = sample_s2(); + break; default: fatal_error("Unknown sample method for random ray transport."); } @@ -874,4 +877,27 @@ SourceSite RandomRay::sample_halton() return site; } +SourceSite RandomRay::sample_s2() +{ + // set random number seed + int64_t particle_seed = + (simulation::current_batch - 1) * settings::n_particles + id(); + init_particle_seeds(particle_seed, seeds()); + stream() = STREAM_TRACKING; + + // Get spatial component of the ray_source_ + SpatialDistribution* space = + dynamic_cast(RandomRay::ray_source_.get())->space(); + + SourceSite site; + + // Sample spatial distribution + site.r = space->sample(current_seed()).first; + + // Sample either left or right for S2 (flashlight) transport. + site.u = {prn(current_seed()) < 0.5 ? -1 : 1, 0.0, 0.0}; + + return site; +} + } // namespace openmc diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 5970dc65a..80cbfc3fe 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -576,9 +576,18 @@ void RandomRaySimulation::print_results_random_ray( fatal_error("Invalid random ray source shape"); } fmt::print(" Source Shape = {}\n", shape); - std::string sample_method = - (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "PRNG" - : "Halton"; + std::string sample_method; + switch (RandomRay::sample_method_) { + case RandomRaySampleMethod::PRNG: + sample_method = "PRNG"; + break; + case RandomRaySampleMethod::HALTON: + sample_method = "Halton"; + break; + case RandomRaySampleMethod::S2: + sample_method = "PRNG S2"; + break; + } fmt::print(" Sample Method = {}\n", sample_method); if (domain_->is_transport_stabilization_needed_) { diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468f..4b356d028 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -326,6 +326,8 @@ void get_run_parameters(pugi::xml_node node_base) RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } else if (temp_str == "halton") { RandomRay::sample_method_ = RandomRaySampleMethod::HALTON; + } else if (temp_str == "s2") { + RandomRay::sample_method_ = RandomRaySampleMethod::S2; } else { fatal_error("Unrecognized sample method: " + temp_str); } diff --git a/tests/regression_tests/random_ray_s2/__init__.py b/tests/regression_tests/random_ray_s2/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/random_ray_s2/inputs_true.dat b/tests/regression_tests/random_ray_s2/inputs_true.dat new file mode 100644 index 000000000..aad8ea28b --- /dev/null +++ b/tests/regression_tests/random_ray_s2/inputs_true.dat @@ -0,0 +1,69 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + fixed source + 100 + 30 + 10 + + + 0.0 -5.0 -5.0 40.0 5.0 5.0 + + + 1000.0 1.0 + + + material + 1 + + + multi-group + + 100.0 + 400.0 + + + 0.0 -5.0 -5.0 40.0 5.0 5.0 + + + flat + s2 + + + + + + + + 10 1 1 + 0.0 -5.0 -5.0 + 40.0 5.0 5.0 + + + + + 1 + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_s2/results_true.dat b/tests/regression_tests/random_ray_s2/results_true.dat new file mode 100644 index 000000000..5fb6f7904 --- /dev/null +++ b/tests/regression_tests/random_ray_s2/results_true.dat @@ -0,0 +1,21 @@ +tally 1: +1.153280E+01 +6.650271E+00 +1.413926E+01 +9.995935E+00 +1.579543E+01 +1.247479E+01 +1.676986E+01 +1.406141E+01 +1.722054E+01 +1.482734E+01 +1.722054E+01 +1.482734E+01 +1.676986E+01 +1.406141E+01 +1.579543E+01 +1.247479E+01 +1.413926E+01 +9.995935E+00 +1.153280E+01 +6.650271E+00 diff --git a/tests/regression_tests/random_ray_s2/test.py b/tests/regression_tests/random_ray_s2/test.py new file mode 100644 index 000000000..712d9c124 --- /dev/null +++ b/tests/regression_tests/random_ray_s2/test.py @@ -0,0 +1,82 @@ +import os + +import openmc +import openmc.model +import numpy as np + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_s2(): + NUM_SOURCE_REGIONS = 10 + L = 40.0 + + # Make the simple MGXS library and material. + groups = openmc.mgxs.EnergyGroups(group_edges=[1e-5, 20.0e6]) + + domain_mat_data = openmc.XSdata('domain', groups) + domain_mat_data.order = 0 + domain_mat_data.set_total([0.1]) + domain_mat_data.set_absorption([0.1]) + domain_mat_data.set_scatter_matrix(np.rollaxis(np.array([[[0.0]]]), 0, 3)) + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas([domain_mat_data]) + mg_cross_sections_file.export_to_hdf5() + + domain_data = openmc.Macroscopic('domain') + domain_mat = openmc.Material(name='domain') + domain_mat.set_density('macro', 1.0) + domain_mat.add_macroscopic(domain_data) + + container = openmc.model.RectangularPrism(width=10.0, height=10.0, axis='x', + origin=(0.0, 0.0), boundary_type='reflective') + left = openmc.XPlane(x0 = 0.0, boundary_type='vacuum') + right = openmc.XPlane(x0 = L, boundary_type='vacuum') + cell = [openmc.Cell(region = +left & -right & -container, fill = domain_mat)] + + model = openmc.model.Model() + model.geometry = openmc.Geometry(root=openmc.Universe(cells=cell)) + model.materials = openmc.Materials([domain_mat]) + model.materials.cross_sections = './mgxs.h5' + + mesh = openmc.RegularMesh() + mesh.dimension = (NUM_SOURCE_REGIONS, 1, 1) + mesh.lower_left = (0.0, -5.0, -5.0) + mesh.upper_right = (L, 5.0, 5.0) + + tally = openmc.Tally(name="LR") + tally.filters = [openmc.MeshFilter(mesh)] + tally.scores = ['flux'] + tally.estimator = 'tracklength' + model.tallies.append(tally) + + uniform_dist = openmc.stats.Box((0.0, -5.0, -5.0), (L, 5.0, 5.0)) + model.settings.source = [ + openmc.IndependentSource(space=uniform_dist, + energy=openmc.stats.Discrete(x = 1e3, p = 1.0), + constraints={'domains' : [domain_mat]}) + ] + model.settings.energy_mode = "multi-group" + model.settings.batches = 30 + model.settings.inactive = 10 + model.settings.particles = 100 + model.settings.run_mode = 'fixed source' + model.settings.random_ray['distance_inactive'] = 100.0 + model.settings.random_ray['distance_active'] = 400.0 + model.settings.random_ray['ray_source'] = openmc.IndependentSource(space=uniform_dist) + model.settings.random_ray['source_shape'] = 'flat' + model.settings.random_ray['sample_method'] = 's2' + model.settings.random_ray['source_region_meshes'] = [(mesh, [model.geometry.root_universe])] + + model.export_to_model_xml() + + harness = MGXSTestHarness('statepoint.30.h5', model) + harness.main()