From 480d2ca6e578616793ffe5d592f40cbd36c7a2fc Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Thu, 13 Oct 2016 16:54:47 -0400 Subject: [PATCH] fixed issue in mgxs_header.F90 where nu cannot be calculated if a fission xs has not been given --- .../python/pincell_multigroup/build-xml.py | 6 +- examples/xml/pincell_multigroup/geometry.xml | 22 ++- examples/xml/pincell_multigroup/materials.xml | 58 ++------ examples/xml/pincell_multigroup/mgxs.h5 | Bin 16496 -> 16744 bytes examples/xml/pincell_multigroup/settings.xml | 50 ++----- examples/xml/pincell_multigroup/tallies.xml | 2 +- src/mgxs_header.F90 | 135 +++++------------- .../results_true.dat | 2 +- .../test_mgxs_library_ce_to_mg.py | 3 +- 9 files changed, 74 insertions(+), 204 deletions(-) diff --git a/examples/python/pincell_multigroup/build-xml.py b/examples/python/pincell_multigroup/build-xml.py index 31a6c52f7..216c56539 100644 --- a/examples/python/pincell_multigroup/build-xml.py +++ b/examples/python/pincell_multigroup/build-xml.py @@ -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() diff --git a/examples/xml/pincell_multigroup/geometry.xml b/examples/xml/pincell_multigroup/geometry.xml index 94c70cb34..7f2fb36e8 100644 --- a/examples/xml/pincell_multigroup/geometry.xml +++ b/examples/xml/pincell_multigroup/geometry.xml @@ -1,16 +1,10 @@ + - - - - - - - - - - - - - - + + + + + + + diff --git a/examples/xml/pincell_multigroup/materials.xml b/examples/xml/pincell_multigroup/materials.xml index f75d7ef6f..9b9ce2180 100644 --- a/examples/xml/pincell_multigroup/materials.xml +++ b/examples/xml/pincell_multigroup/materials.xml @@ -1,51 +1,11 @@ - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + diff --git a/examples/xml/pincell_multigroup/mgxs.h5 b/examples/xml/pincell_multigroup/mgxs.h5 index cbe338c1c7de692a9bed69901f1da186bc1d3348..debdf0257b3a70162986d343d179c444d985fccd 100644 GIT binary patch delta 1317 zcmZuwO-vI(6yBXSWp`<5AqHwoV3AUgf)td3mAC~BG)4^^jD(m#kxDfRl^+9*A=GHV z14-*UNc5r{v=9$$HYR!_XtcotjVF%vLP*r4df}{{7DMIjlpaGgRYVynh6;`aCq`!nM<>Ro z!;=(T);WM{;Ed$2k>KsXY;a@{{5qay;HK^~9fD?@CnN>*Wr$9~tYN$>aa=AT*Co4@ z8fRPz&L}_B5UGk6;YZe@8Yc{KaUXac%V>bq%=FM~idD{KAj;O#8`3IUhLFVwn}uVY zjM|`<%dN(}Hsy|1Goj|t&Si5-ca42Cy#XIgRR!Ch4z6N4adAnSDdK4%=52g|WUCgQc!MoMC*ONc^gqlSY zD3_Hyn6$G5A|^kw7Q;uA#-ErtXT4_4dR(rAOlYS~8aRmwg*`CQnnZkYR@XsK9Py`R z!?ZE=Zf`!*el?!y$Y+}LkZBT&^ygBU?NMTTE_X?;O`XjuxZHsiOtP_JCEMO1_4<8( VtLGOq^)elrdNU>eQ*ZrW^bfgsQp*4U delta 1109 zcmZ8fJ#5oJ6ux^Y>Ddiw=#W}fQm#K~LR*ng)DTn(hJ}$bQ56Z91)&TjY$=n4glGjr z+_waYp=_uNwY*drXxIW%Vd|I(Rq8-kP&TlA=Q!b_=-qes{N4B7_fB8Zj!A1dsOZ*w ziJv2s5Tdz%ZzbHHw|K-4q)TvK%hIU5UxRpo47sP243ijP1ZxqFm57Rs$YbHuX$;KBI~TE&t0ROIigwv&(B)ZvM} zsMXZ5nM}~9&iCX!s)pS7O_`%L=e{-!V+TS2S(EW(1~-N#)w--$r&ykzR8JPgb@wth zid8V<_%_YG*YGIwIu?}%5OR{T>`5%7813TY)J2s`*3-0iZl)M8?(UX!d9F=~O;DFE zu`wgggZj-0%ictuF2iB?^r!S0bOs>J;hVJI1#~Z|(L=wc-l-wCS3lGzG#+b?8X>qd4Eh2R|nfvl6ZZ2 zuuW_K&2MO%922M|Fcudj<(EFdVZpYG6f4u38glyWX%{^WXcQS#m$6)IK|QSGD6@m@ zYYcKF&rB{vkka)*m4?LT%Am^QWznxP9EX&iL*WO$1hy(1+FQZh8!gK4;i~ZVC}V+f YD==g2X1^KTm3}ii7Y5A8jrGBuf6M0rmH+?% diff --git a/examples/xml/pincell_multigroup/settings.xml b/examples/xml/pincell_multigroup/settings.xml index ab716d9e4..6b9080b5a 100644 --- a/examples/xml/pincell_multigroup/settings.xml +++ b/examples/xml/pincell_multigroup/settings.xml @@ -1,39 +1,15 @@ - + - - multi-group - - - - 100 - 10 - 1000 - - - - - - - -0.63 -0.63 -1E50 - 0.63 0.63 1E50 - - - - - - true - true - - - false - - ./mgxs.h5 - + + 1000 + 100 + 10 + + + + -0.63 -0.63 -1 0.63 0.63 1 + + + ./mgxs.h5 + multi-group diff --git a/examples/xml/pincell_multigroup/tallies.xml b/examples/xml/pincell_multigroup/tallies.xml index df65b461d..0b8b8ef92 100644 --- a/examples/xml/pincell_multigroup/tallies.xml +++ b/examples/xml/pincell_multigroup/tallies.xml @@ -1,4 +1,4 @@ - + 100 100 1 diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index a26453a8f..7c5784649 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -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 diff --git a/tests/test_mgxs_library_ce_to_mg/results_true.dat b/tests/test_mgxs_library_ce_to_mg/results_true.dat index 0b73724c9..27db7f3f6 100644 --- a/tests/test_mgxs_library_ce_to_mg/results_true.dat +++ b/tests/test_mgxs_library_ce_to_mg/results_true.dat @@ -1,2 +1,2 @@ k-combined: -1.140724E+00 2.952379E-02 +1.133362E+00 1.746521E-02 diff --git a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py index eee528e88..74a424597 100644 --- a/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py +++ b/tests/test_mgxs_library_ce_to_mg/test_mgxs_library_ce_to_mg.py @@ -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)