Fix random ray solver to correctly simulate fixed source problems with fissionable materials (#3106)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
John Tramm 2024-08-15 20:47:05 -05:00 committed by GitHub
parent e7bc9ba23c
commit b22656e57f
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
8 changed files with 772 additions and 21 deletions

View file

@ -804,7 +804,7 @@ def random_ray_lattice() -> openmc.Model:
azimuthal_cells.append(azimuthal_cell)
# Create a geometry with the azimuthal universes
pincell = openmc.Universe(cells=azimuthal_cells)
pincell = openmc.Universe(cells=azimuthal_cells, name='pincell')
########################################
# Define a moderator lattice universe

View file

@ -158,31 +158,31 @@ void FlatSourceDomain::update_neutron_source(double k_eff)
}
}
if (settings::run_mode == RunMode::EIGENVALUE) {
// Add fission source if in eigenvalue mode
// Add fission source
#pragma omp parallel for
for (int sr = 0; sr < n_source_regions_; sr++) {
int material = material_[sr];
for (int sr = 0; sr < n_source_regions_; sr++) {
int material = material_[sr];
for (int e_out = 0; e_out < negroups_; e_out++) {
float sigma_t = data::mg.macro_xs_[material].get_xs(
MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
float fission_source = 0.0f;
for (int e_out = 0; e_out < negroups_; e_out++) {
float sigma_t = data::mg.macro_xs_[material].get_xs(
MgxsType::TOTAL, e_out, nullptr, nullptr, nullptr, t, a);
float fission_source = 0.0f;
for (int e_in = 0; e_in < negroups_; e_in++) {
float scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
float nu_sigma_f = data::mg.macro_xs_[material].get_xs(
MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
float chi = data::mg.macro_xs_[material].get_xs(
MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
fission_source += nu_sigma_f * scalar_flux * chi;
}
source_[sr * negroups_ + e_out] +=
fission_source * inverse_k_eff / sigma_t;
for (int e_in = 0; e_in < negroups_; e_in++) {
float scalar_flux = scalar_flux_old_[sr * negroups_ + e_in];
float nu_sigma_f = data::mg.macro_xs_[material].get_xs(
MgxsType::NU_FISSION, e_in, nullptr, nullptr, nullptr, t, a);
float chi = data::mg.macro_xs_[material].get_xs(
MgxsType::CHI_PROMPT, e_in, &e_out, nullptr, nullptr, t, a);
fission_source += nu_sigma_f * scalar_flux * chi;
}
source_[sr * negroups_ + e_out] +=
fission_source * inverse_k_eff / sigma_t;
}
} else {
// Add external source if in fixed source mode
}
// Add external source if in fixed source mode
if (settings::run_mode == RunMode::FIXED_SOURCE) {
#pragma omp parallel for
for (int se = 0; se < n_source_elements_; se++) {
source_[se] += external_source_[se];

View file

@ -0,0 +1,140 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0"/>
<macroscopic name="UO2"/>
</material>
<material id="2" name="Water">
<density units="macro" value="1.0"/>
<macroscopic name="LWTR"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
<cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
<cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
<cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
<cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
<cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
<cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
<cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
<cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
<cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
<cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
<cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
<cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
<cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
<cell id="18" material="2" name="moderator infinite" universe="8"/>
<cell id="19" material="2" name="moderator infinite 2" universe="9"/>
<cell fill="10" id="20" universe="12"/>
<cell fill="11" id="21" universe="13"/>
<cell fill="14" id="22" name="assembly" region="18 -19 20 -21" universe="15"/>
<lattice id="10">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</lower_left>
<universes>
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8 </universes>
</lattice>
<lattice id="11">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</lower_left>
<universes>
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 9 </universes>
</lattice>
<lattice id="14">
<pitch>1.26 1.26</pitch>
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<universes>
2 12
12 13 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
<surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
<surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
<surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
<surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
<surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
<surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
<surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
<surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
<surface boundary="reflective" coeffs="-1.26" id="18" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="1.26" id="19" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.26" id="20" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="1.26" id="21" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>30</particles>
<batches>125</batches>
<inactive>100</inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="box">
<parameters>1.134 -1.26 -1.0 1.26 -1.134 1.0</parameters>
</space>
<energy type="discrete">
<parameters>2e-05 0.0735 20.0 200.0 2000.0 750000.0 2000000.0 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285</parameters>
</energy>
<constraints>
<domain_type>universe</domain_type>
<domain_ids>9</domain_ids>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>40.0</distance_active>
<distance_inactive>40.0</distance_inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
<source_shape>flat</source_shape>
</random_ray>
</settings>
<tallies>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<upper_right>1.26 1.26</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<tally id="1" name="Mesh tally">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

View file

@ -0,0 +1,169 @@
tally 1:
1.582116E+02
1.005480E+03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.897563E+01
1.396228E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.038901E+01
1.667203E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.465551E+01
2.438101E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.260230E+01
1.110141E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.438906E+01
3.579643E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.698287E+01
1.314208E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.002152E+02
1.608182E+03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.374897E+01
2.181176E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.636086E+01
2.787512E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
3.000395E+01
3.611990E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.824715E+01
1.360847E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.621963E+01
3.714397E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.954480E+01
1.436984E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.067903E+02
4.585529E+02
1.257560E+01
6.364831E+00
3.060650E+01
3.770132E+01
4.983480E+01
9.972910E+01
2.349541E+00
2.217300E-01
5.718312E+00
1.313392E+00
1.965462E+01
1.548888E+01
2.011494E-01
1.622334E-03
4.895573E-01
9.609705E-03
2.318671E+01
2.155841E+01
2.445032E-01
2.397352E-03
5.950720E-01
1.420043E-02
5.322867E+01
1.136138E+02
1.968847E-01
1.554358E-03
4.791838E-01
9.207284E-03
1.191175E+02
5.692803E+02
5.806982E-02
1.354810E-04
1.436897E-01
8.295235E-04
8.022963E+01
2.589171E+02
3.465139E-01
4.882327E-03
9.638109E-01
3.777193E-02
1.555305E+02
9.711601E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.836421E+01
1.367538E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.042089E+01
1.673401E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.477708E+01
2.463124E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.317851E+01
1.134205E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.634473E+01
3.724422E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.879813E+01
1.396454E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

View file

@ -0,0 +1,140 @@
<?xml version='1.0' encoding='utf-8'?>
<model>
<materials>
<cross_sections>mgxs.h5</cross_sections>
<material id="1" name="UO2 fuel">
<density units="macro" value="1.0"/>
<macroscopic name="UO2"/>
</material>
<material id="2" name="Water">
<density units="macro" value="1.0"/>
<macroscopic name="LWTR"/>
</material>
</materials>
<geometry>
<cell id="1" material="1" name="fuel inner a" region="-2" universe="1"/>
<cell id="2" material="1" name="fuel inner b" region="2 -3" universe="1"/>
<cell id="3" material="1" name="fuel inner c" region="3 -1" universe="1"/>
<cell id="4" material="2" name="moderator inner a" region="1 -4" universe="1"/>
<cell id="5" material="2" name="moderator outer b" region="4 -5" universe="1"/>
<cell id="6" material="2" name="moderator outer c" region="5" universe="1"/>
<cell fill="1" id="7" name="azimuthal_cell_0" region="6 -7" universe="2"/>
<cell fill="1" id="8" name="azimuthal_cell_1" region="7 -8" universe="2"/>
<cell fill="1" id="9" name="azimuthal_cell_2" region="8 -9" universe="2"/>
<cell fill="1" id="10" name="azimuthal_cell_3" region="9 -10" universe="2"/>
<cell fill="1" id="11" name="azimuthal_cell_4" region="10 -11" universe="2"/>
<cell fill="1" id="12" name="azimuthal_cell_5" region="11 -12" universe="2"/>
<cell fill="1" id="13" name="azimuthal_cell_6" region="12 -13" universe="2"/>
<cell fill="1" id="14" name="azimuthal_cell_7" region="13 -6" universe="2"/>
<cell id="18" material="2" name="moderator infinite" universe="8"/>
<cell id="19" material="2" name="moderator infinite 2" universe="9"/>
<cell fill="10" id="20" universe="12"/>
<cell fill="11" id="21" universe="13"/>
<cell fill="14" id="22" name="assembly" region="18 -19 20 -21" universe="15"/>
<lattice id="10">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</lower_left>
<universes>
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8 </universes>
</lattice>
<lattice id="11">
<pitch>0.126 0.126</pitch>
<dimension>10 10</dimension>
<lower_left>-0.63 -0.63</lower_left>
<universes>
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 8
8 8 8 8 8 8 8 8 8 9 </universes>
</lattice>
<lattice id="14">
<pitch>1.26 1.26</pitch>
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<universes>
2 12
12 13 </universes>
</lattice>
<surface coeffs="0.0 0.0 0.54" id="1" name="Fuel OR" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.33" id="2" name="inner ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.45" id="3" name="inner ring b" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.6" id="4" name="outer ring a" type="z-cylinder"/>
<surface coeffs="0.0 0.0 0.69" id="5" name="outer ring b" type="z-cylinder"/>
<surface coeffs="-0.0 1.0 0 0" id="6" type="plane"/>
<surface coeffs="-0.7071067811865475 0.7071067811865476 0 0" id="7" type="plane"/>
<surface coeffs="-1.0 6.123233995736766e-17 0 0" id="8" type="plane"/>
<surface coeffs="-0.7071067811865476 -0.7071067811865475 0 0" id="9" type="plane"/>
<surface coeffs="-1.2246467991473532e-16 -1.0 0 0" id="10" type="plane"/>
<surface coeffs="0.7071067811865475 -0.7071067811865477 0 0" id="11" type="plane"/>
<surface coeffs="1.0 -1.8369701987210297e-16 0 0" id="12" type="plane"/>
<surface coeffs="0.7071067811865477 0.7071067811865474 0 0" id="13" type="plane"/>
<surface boundary="reflective" coeffs="-1.26" id="18" name="minimum x" type="x-plane"/>
<surface boundary="reflective" coeffs="1.26" id="19" name="maximum x" type="x-plane"/>
<surface boundary="reflective" coeffs="-1.26" id="20" name="minimum y" type="y-plane"/>
<surface boundary="reflective" coeffs="1.26" id="21" name="maximum y" type="y-plane"/>
</geometry>
<settings>
<run_mode>fixed source</run_mode>
<particles>30</particles>
<batches>125</batches>
<inactive>100</inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="box">
<parameters>1.134 -1.26 -1.0 1.26 -1.134 1.0</parameters>
</space>
<energy type="discrete">
<parameters>2e-05 0.0735 20.0 200.0 2000.0 750000.0 2000000.0 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285 0.14285714285714285</parameters>
</energy>
<constraints>
<domain_type>universe</domain_type>
<domain_ids>9</domain_ids>
</constraints>
</source>
<energy_mode>multi-group</energy_mode>
<random_ray>
<distance_active>40.0</distance_active>
<distance_inactive>40.0</distance_inactive>
<source particle="neutron" strength="1.0" type="independent">
<space type="box">
<parameters>-1.26 -1.26 -1 1.26 1.26 1</parameters>
</space>
</source>
<volume_normalized_flux_tallies>False</volume_normalized_flux_tallies>
<source_shape>linear_xy</source_shape>
</random_ray>
</settings>
<tallies>
<mesh id="1">
<dimension>2 2</dimension>
<lower_left>-1.26 -1.26</lower_left>
<upper_right>1.26 1.26</upper_right>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="2" type="energy">
<bins>1e-05 0.0635 10.0 100.0 1000.0 500000.0 1000000.0 20000000.0</bins>
</filter>
<tally id="1" name="Mesh tally">
<filters>1 2</filters>
<scores>flux fission nu-fission</scores>
<estimator>analog</estimator>
</tally>
</tallies>
</model>

View file

@ -0,0 +1,169 @@
tally 1:
1.573422E+02
9.945961E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.868035E+01
1.382578E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.028997E+01
1.651549E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.453137E+01
2.414255E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.232506E+01
1.098683E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.385720E+01
3.539613E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.664229E+01
1.298333E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.001324E+02
1.606858E+03
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
7.358736E+01
2.171786E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.625416E+01
2.765491E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.987675E+01
3.582038E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.797311E+01
1.348247E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.573509E+01
3.677339E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.926615E+01
1.423707E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
1.060513E+02
4.521166E+02
1.247077E+01
6.257701E+00
3.035137E+01
3.706675E+01
4.962504E+01
9.887387E+01
2.338861E+00
2.196709E-01
5.692320E+00
1.301195E+00
1.957564E+01
1.536441E+01
2.003053E-01
1.608650E-03
4.875030E-01
9.528647E-03
2.309275E+01
2.138375E+01
2.434816E-01
2.377213E-03
5.925856E-01
1.408114E-02
5.300741E+01
1.126617E+02
1.960322E-01
1.540737E-03
4.771089E-01
9.126597E-03
1.184809E+02
5.630856E+02
5.774249E-02
1.339082E-04
1.428798E-01
8.198938E-04
7.971031E+01
2.554943E+02
3.441954E-01
4.814135E-03
9.573621E-01
3.724437E-02
1.549617E+02
9.639834E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.808547E+01
1.354629E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.029953E+01
1.654123E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
2.462331E+01
2.433425E+01
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.283851E+01
1.120036E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
9.572799E+01
3.677288E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00
5.842566E+01
1.378719E+02
0.000000E+00
0.000000E+00
0.000000E+00
0.000000E+00

View file

@ -0,0 +1,133 @@
import os
import openmc
from openmc.examples import random_ray_lattice
from openmc.utility_funcs import change_directory
import pytest
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)
@pytest.mark.parametrize("shape", ["flat", "linear_xy"])
def test_random_ray_source(shape):
with change_directory(shape):
openmc.reset_auto_ids()
# The general strategy is to reuse the random_ray_lattice model,
# but redfine some of the geometry to make it a good
# subcritical multiplication problem. We then also add in
# a fixed source term.
model = random_ray_lattice()
# Begin by updating the random ray settings for fixed source
settings = model.settings
settings.random_ray['source_shape'] = shape
settings.run_mode = 'fixed source'
settings.particles = 30
settings.random_ray['distance_active'] = 40.0
settings.random_ray['distance_inactive'] = 40.0
settings.random_ray['volume_normalized_flux_tallies'] = False
# This problem needs about 2k iterations to converge,
# but for regression testing we only need a few hundred
# to ensure things are working as expected. With
# only 100 inactive batches, tallies will still be off
# by 3x or more. For validation against MGMC, be sure
# to increase the batch counts.
settings.batches = 125
settings.inactive = 100
########################################
# Define the alternative geometry
pitch = 1.26
for material in model.materials:
if material.name == 'Water':
water = material
# The new geometry replaces two of the fuel pins with
# moderator, reducing k-eff to around 0.84. We also
# add a special universe in the corner of one of the moderator
# regions to use as a domain constraint for the source
moderator_infinite = openmc.Cell(fill=water, name='moderator infinite')
mu = openmc.Universe(cells=[moderator_infinite])
moderator_infinite2 = openmc.Cell(fill=water, name='moderator infinite 2')
mu2 = openmc.Universe(cells=[moderator_infinite2])
n_sub = 10
lattice = openmc.RectLattice()
lattice.lower_left = [-pitch/2.0, -pitch/2.0]
lattice.pitch = [pitch/n_sub, pitch/n_sub]
lattice.universes = [[mu] * n_sub for _ in range(n_sub)]
lattice2 = openmc.RectLattice()
lattice2.lower_left = [-pitch/2.0, -pitch/2.0]
lattice2.pitch = [pitch/n_sub, pitch/n_sub]
lattice2.universes = [[mu] * n_sub for _ in range(n_sub)]
lattice2.universes[n_sub-1][n_sub-1] = mu2
mod_lattice_cell = openmc.Cell(fill=lattice)
mod_lattice_uni = openmc.Universe(cells=[mod_lattice_cell])
mod_lattice_cell2 = openmc.Cell(fill=lattice2)
mod_lattice_uni2 = openmc.Universe(cells=[mod_lattice_cell2])
lattice2x2 = openmc.RectLattice()
lattice2x2.lower_left = [-pitch, -pitch]
lattice2x2.pitch = [pitch, pitch]
universes = model.geometry.get_all_universes()
for universe in universes.values():
if universe.name == 'pincell':
pincell = universe
lattice2x2.universes = [
[pincell, mod_lattice_uni],
[mod_lattice_uni, mod_lattice_uni2]
]
box = openmc.model.RectangularPrism(
pitch*2, pitch*2, boundary_type='reflective')
assembly = openmc.Cell(fill=lattice2x2, region=-box, name='assembly')
root = openmc.Universe(name='root universe', cells=[assembly])
model.geometry = openmc.Geometry(root)
########################################
# Define the fixed source term
s = 1.0 / 7.0
strengths = [s, s, s, s, s, s, s]
midpoints = [2.0e-5, 0.0735, 20.0, 2.0e2, 2.0e3, 0.75e6, 2.0e6]
energy_distribution = openmc.stats.Discrete(x=midpoints, p=strengths)
lower_left_src = [pitch - pitch/10.0, -pitch, -1.0]
upper_right_src = [pitch, -pitch + pitch/10.0, 1.0]
spatial_distribution = openmc.stats.Box(
lower_left_src, upper_right_src, only_fissionable=False)
settings.source = openmc.IndependentSource(
space=spatial_distribution,
energy=energy_distribution,
constraints={'domains': [mu2]},
strength=1.0
)
########################################
# Run test
harness = MGXSTestHarness('statepoint.125.h5', model)
harness.main()