crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 13:45:36 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

fixed issue in mgxs_header.F90 where nu cannot be calculated if a fission xs has not been given

Browse source
This commit is contained in:
Sam Shaner 2016-10-13 16:54:47 -04:00
parent 89149ea7eb
commit 480d2ca6e5
9 changed files with 74 additions and 204 deletions
Download patch file Download diff file Expand all files Collapse all files

6
examples/python/pincell_multigroup/build-xml.py
View file

@ -19,7 +19,7 @@ groups = openmc.mgxs.EnergyGroups(group_edges=[1E-11, 0.0635E-6, 10.0E-6,
1.0E-4, 1.0E-3, 0.5, 1.0, 20.0])
# Instantiate the 7-group (C5G7) cross section data
uo2_xsdata = openmc.XSdata('UO2', groups, delayed_groups=8)
uo2_xsdata = openmc.XSdata('UO2', groups)
uo2_xsdata.order = 0
uo2_xsdata.set_total(
[0.1779492, 0.3298048, 0.4803882, 0.5543674, 0.3118013, 0.3951678,
@ -42,8 +42,6 @@ uo2_xsdata.set_nu_fission([2.005998E-02, 2.027303E-03, 1.570599E-02,
5.257105E-01])
uo2_xsdata.set_chi([5.8791E-01, 4.1176E-01, 3.3906E-04, 1.1761E-07, 0.0000E+00,
0.0000E+00, 0.0000E+00])
uo2_xsdata.set_beta([2.13333e-4, 1.04514e-3, 6.03969e-4, 1.33963e-3,
2.29386e-3, 7.05174e-4, 6.00381e-4, 2.07736e-4])
h2o_xsdata = openmc.XSdata('LWTR', groups)
h2o_xsdata.order = 0
@ -61,7 +59,7 @@ h2o_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
mg_cross_sections_file = openmc.MGXSLibrary(groups, delayed_groups=8)
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])
mg_cross_sections_file.export_to_hdf5()

22
examples/xml/pincell_multigroup/geometry.xml
View file

@ -1,16 +1,10 @@
<?xml version='1.0' encoding='utf-8'?>
<geometry>
<surface coeffs="0. 0. 0.540" id="1" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="20" type="x-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="21" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="22" type="y-plane" />
<surface boundary="reflective" coeffs=" 0.63" id="23" type="y-plane" />
<cell id="1" material="1" region=" -1" />
<cell id="2" material="7" region="1 20 -21 22 -23" />
<cell id="1" material="1" name="cell 1" region="-1" universe="0" />
<cell id="2" material="2" name="cell 2" region="((((1 4) -5) 6) -7)" universe="0" />
<surface coeffs="0 0 0.54" id="1" name="Fuel OR" type="z-cylinder" />
<surface boundary="reflective" coeffs="-0.63" id="4" name="left" type="x-plane" />
<surface boundary="reflective" coeffs="0.63" id="5" name="right" type="x-plane" />
<surface boundary="reflective" coeffs="-0.63" id="6" name="bottom" type="y-plane" />
<surface boundary="reflective" coeffs="0.63" id="7" name="top" type="y-plane" />
</geometry>

58
examples/xml/pincell_multigroup/materials.xml
View file

@ -1,51 +1,11 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<materials>
<!-- UO2 -->
<material id="1">
<density units="macro" value="1.0" />
<macroscopic name="UO2"/>
</material>
<!-- 4.3% MOX -->
<material id="2">
<density units="macro" value="1.0" />
<macroscopic name="MOX1"/>
</material>
<!-- 7.0 MOX -->
<material id="3">
<density units="macro" value="1.0" />
<macroscopic name="MOX2"/>
</material>
<!-- 8.0% MOX -->
<material id="4">
<density units="macro" value="1.0" />
<macroscopic name="MOX3"/>
</material>
<!-- Fission Chamber -->
<material id="5">
<density units="macro" value="1.0" />
<macroscopic name="FC"/>
</material>
<!-- Guide Tube -->
<material id="6">
<density units="macro" value="1.0" />
<macroscopic name="GT"/>
</material>
<!-- Water -->
<material id="7">
<density units="macro" value="1.0" />
<macroscopic name="LWTR"/>
</material>
<!-- Control Rod -->
<material id="8">
<density units="macro" value="1.0" />
<macroscopic name="CR"/>
</material>
<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>

BIN
examples/xml/pincell_multigroup/mgxs.h5
View file

Binary file not shown.

50
examples/xml/pincell_multigroup/settings.xml
View file

@ -1,39 +1,15 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<settings>
<energy_mode>multi-group</energy_mode>
<!--
Define how many particles to run and for how many batches
in an eigenvalue calculation mode
-->
<eigenvalue>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
</eigenvalue>
<!--
Start with uniformally distributed neutron source
with the default energy spectrum of a Maxwellian
and isotropic distribution.
-->
<source>
<space type="box">
<parameters>
-0.63 -0.63 -1E50
0.63 0.63 1E50
</parameters>
</space>
</source>
<output>
<summary>true</summary>
<tallies>true</tallies>
</output>
<survival_biasing>false</survival_biasing>
<cross_sections>./mgxs.h5</cross_sections>
<eigenvalue>
<particles>1000</particles>
<batches>100</batches>
<inactive>10</inactive>
</eigenvalue>
<source strength="1.0">
<space type="box">
<parameters>-0.63 -0.63 -1 0.63 0.63 1</parameters>
</space>
</source>
<cross_sections>./mgxs.h5</cross_sections>
<energy_mode>multi-group</energy_mode>
</settings>

2
examples/xml/pincell_multigroup/tallies.xml
View file

@ -1,4 +1,4 @@
<?xml version="1.0"?>
<?xml version='1.0' encoding='utf-8'?>
<tallies>
<mesh id="1" type="regular">
<dimension>100 100 1</dimension>

135
src/mgxs_header.F90
View file

@ -426,7 +426,7 @@ module mgxs_header
integer(HSIZE_T) :: dims(2)
real(8), allocatable :: temp_arr(:), temp_2d(:, :)
real(8), allocatable :: temp_beta(:, :)
real(8) :: dmu, mu, norm
real(8) :: dmu, mu, norm, chi_sum
integer :: order, order_dim, gin, gout, l, imu, length
type(VectorInt) :: temps_to_read
integer :: t, dg
@ -535,12 +535,12 @@ module mgxs_header
end do
! Normalize chi_prompt so its CDF goes to 1
if (sum(xs % chi_prompt(:, gin)) == ZERO) then
chi_sum =sum(xs % chi_prompt(:, gin))
if (chi_sum == ZERO) then
call fatal_error("Encountered chi for a group that sums to &
&zero")
else
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / &
sum(xs % chi_prompt(:, gin))
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum
end if
end do
@ -625,12 +625,12 @@ module mgxs_header
! Normalize chi so its CDF goes to 1
do gin = 1, energy_groups
if (sum(xs % chi_prompt(:, gin)) == ZERO) then
chi_sum = sum(xs % chi_prompt(:, gin))
if (chi_sum == ZERO) then
call fatal_error("Encountered chi for a group that sums to &
&zero")
else
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / &
sum(xs % chi_prompt(:, gin))
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum
end if
end do
@ -642,6 +642,11 @@ module mgxs_header
call fatal_error("nu-fission must be provided as a 1D or 2D &
&array")
end if
print *, xs % prompt_nu_fission
print *, xs % delayed_nu_fission
print *, xs % chi_prompt
print *, xs % chi_delayed
end if
! If chi_prompt provided, set chi_prompt
@ -659,12 +664,12 @@ module mgxs_header
end do
! Normalize chi so its CDF goes to 1
if (sum(xs % chi_prompt(:, gin)) == ZERO) then
chi_sum = sum(xs % chi_prompt(:, gin))
if (chi_sum == ZERO) then
call fatal_error("Encountered chi prompt for a group that &
&sums to zero")
else
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / &
sum(xs % chi_prompt(:, gin))
xs % chi_prompt(:, gin) = xs % chi_prompt(:, gin) / chi_sum
end if
end do
@ -695,13 +700,13 @@ module mgxs_header
end do
! Normalize chi so its CDF goes to 1
if (sum(xs % chi_delayed(:, gin, dg)) == ZERO) then
chi_sum = sum(xs % chi_delayed(:, gin, dg))
if (chi_sum == ZERO) then
call fatal_error("Encountered chi delayed for a group &
&that sums to zero")
else
xs % chi_delayed(:, gin, dg) = &
xs % chi_delayed(:, gin, dg) / &
sum(xs % chi_delayed(:, gin, dg))
xs % chi_delayed(:, gin, dg) / chi_sum
end if
end do
end do
@ -726,13 +731,13 @@ module mgxs_header
end do
! Normalize chi so its CDF goes to 1
if (sum(xs % chi_delayed(:, gin, dg)) == ZERO) then
chi_sum = sum(xs % chi_delayed(:, gin, dg))
if (chi_sum == ZERO) then
call fatal_error("Encountered chi delayed for a group &
&that sums to zero")
else
xs % chi_delayed(:, gin, dg) = &
xs % chi_delayed(:, gin, dg) / &
sum(xs % chi_delayed(:, gin, dg))
xs % chi_delayed(:, gin, dg) / chi_sum
end if
end do
end do
@ -1101,7 +1106,7 @@ module mgxs_header
integer(HID_T) :: xsdata, xsdata_grp, scatt_grp, ndims
integer(HSIZE_T) :: dims(4)
integer, allocatable :: int_arr(:)
real(8), allocatable :: temp_1d(:), temp_2d(:, :), temp_3d(:, :, :)
real(8), allocatable :: temp_1d(:), temp_3d(:, :, :)
real(8), allocatable :: temp_4d(:, :, :, :), temp_beta(:, :, :, :)
real(8) :: dmu, mu, norm
integer :: order, order_dim, gin, gout, l, imu, dg
@ -2713,36 +2718,6 @@ module mgxs_header
xs = this % xs(t) % prompt_nu_fission(gin) + &
sum(this % xs(t) % delayed_nu_fission(gin, :))
case('nu')
if (this % xs(t) % fission(gin) > ZERO) then
xs = (sum(this % xs(t) % delayed_nu_fission(gin, :)) + &
this % xs(t) % prompt_nu_fission(gin)) / &
this % xs(t) % fission(gin)
else
xs = ZERO
end if
case('nu_prompt')
if (this % xs(t) % fission(gin) > ZERO) then
xs = this % xs(t) % prompt_nu_fission(gin) / &
this % xs(t) % fission(gin)
else
xs = ZERO
end if
case('nu_delayed')
if (this % xs(t) % fission(gin) > ZERO) then
if (present(dg)) then
xs = this % xs(t) % delayed_nu_fission(gin, dg) &
/ this % xs(t) % fission(gin)
else
xs = sum(this % xs(t) % delayed_nu_fission(gin, :)) &
/ this % xs(t) % fission(gin)
end if
else
xs = ZERO
end if
case('chi_prompt')
if (present(gout)) then
xs = this % xs(t) % chi_prompt(gout,gin)
@ -2868,36 +2843,6 @@ module mgxs_header
xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) + &
sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol))
case('nu')
if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then
xs = (sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) + &
this % xs(t) % prompt_nu_fission(gin, iazi, ipol)) / &
this % xs(t) % fission(gin, iazi, ipol)
else
xs = ZERO
end if
case('nu_prompt')
if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then
xs = this % xs(t) % prompt_nu_fission(gin, iazi, ipol) / &
this % xs(t) % fission(gin, iazi, ipol)
else
xs = ZERO
end if
case('nu_delayed')
if (this % xs(t) % fission(gin, iazi, ipol) > ZERO) then
if (present(dg)) then
xs = this % xs(t) % delayed_nu_fission(gin, dg, iazi, ipol) / &
this % xs(t) % fission(gin, iazi, ipol)
else
xs = sum(this % xs(t) % delayed_nu_fission(gin, :, iazi, ipol)) &
/ this % xs(t) % fission(gin, iazi, ipol)
end if
else
xs = ZERO
end if
case('chi_prompt')
if (present(gout)) then
xs = this % xs(t) % chi_prompt(gout, gin, iazi, ipol)
@ -3006,20 +2951,20 @@ module mgxs_header
integer, intent(out) :: gout ! Sampled outgoing group
real(8) :: xi_pd ! Our random number for prompt/delayed
real(8) :: xi_gout ! Our random number for gout
real(8) :: prob_pd ! Running probability for prompt/delayed
real(8) :: prob_gout ! Running probability for gout
! Get nu and nu_prompt
real(8) :: nu, nu_prompt
nu = this % get_xs('nu', gin)
nu_prompt = this % get_xs('nu_prompt', gin)
real(8) :: prob_prompt
prob_prompt = this % get_xs('prompt_nu_fission', gin) / &
this % get_xs('nu_fission', gin)
! Sample random numbers
xi_pd = prn()
xi_gout = prn()
! Neutron is born prompt
if (xi_pd < nu_prompt / nu) then
if (xi_pd <= prob_prompt) then
! set the delayed group for the particle born from fission to 0
dg = 0
@ -3038,13 +2983,12 @@ module mgxs_header
! Get the delayed group
dg = 0
prob_pd = nu_prompt / nu
do while (xi_pd >= prob_pd)
do while (xi_pd >= prob_prompt)
dg = dg + 1
prob_pd = prob_pd + &
this % get_xs('nu_delayed', gin, dg=dg) &
/ nu
prob_prompt = prob_prompt + &
this % get_xs('delayed_nu_fision', gin, dg=dg) &
/ this % get_xs('nu_fission', gin)
end do
! Adjust dg in case of round off error
@ -3071,20 +3015,19 @@ module mgxs_header
integer, intent(out) :: gout ! Sampled outgoing group
real(8) :: xi_pd ! Our random number for prompt/delayed
real(8) :: xi_gout ! Our random number for gout
real(8) :: prob_pd ! Running probability for prompt/delayed
real(8) :: prob_gout ! Running probability for gout
real(8) :: nu, nu_prompt
real(8) :: prob_prompt
! Get nu and nu_prompt
nu = this % get_xs('nu', gin, uvw=uvw)
nu_prompt = this % get_xs('nu_prompt', gin, uvw=uvw)
prob_prompt = this % get_xs('prompt_nu_fission', gin, uvw=uvw) / &
this % get_xs('nu_fission', gin, uvw=uvw)
! Sample random numbers
xi_pd = prn()
xi_gout = prn()
! Neutron is born prompt
if (xi_pd < nu_prompt / nu) then
if (xi_pd <= prob_prompt) then
! set the delayed group for the particle born from fission to 0
dg = 0
@ -3103,12 +3046,12 @@ module mgxs_header
! Get the delayed group
dg = 0
prob_pd = nu_prompt / nu
do while (xi_pd < prob_pd)
do while (xi_pd < prob_prompt)
dg = dg + 1
prob_pd = prob_pd + &
this % get_xs('nu_delayed', gin, uvw=uvw, dg=dg) / nu
prob_prompt = prob_prompt + &
this % get_xs('delayed_nu_fission', gin, uvw=uvw, dg=dg) / &
this % get_xs('nu_fission', gin, uvw=uvw)
end do
! Adjust dg in case of round off error

2
tests/test_mgxs_library_ce_to_mg/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
1.140724E+00 2.952379E-02
1.133362E+00 1.746521E-02

3
tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py
View file

@ -20,8 +20,7 @@ class MGXSTestHarness(PyAPITestHarness):
super(MGXSTestHarness, self)._build_inputs()
# Initialize a two-group structure
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625e-6,
20.])
energy_groups = openmc.mgxs.EnergyGroups(group_edges=[0, 0.625e-6, 20.])
# Initialize MGXS Library for a few cross section types
self.mgxs_lib = openmc.mgxs.Library(self._input_set.geometry)