From 7bb1df8c7402f7c637530d4ba4695625e898e80a Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 21 Nov 2013 16:55:54 -0500 Subject: [PATCH 01/45] added res scat variables to template files --- src/templates/settings.rnc | 13 +++++++++++++ src/templates/settings_t.xml | 10 ++++++++++ 2 files changed, 23 insertions(+) diff --git a/src/templates/settings.rnc b/src/templates/settings.rnc index 1c7cabafa..e2f5757ad 100644 --- a/src/templates/settings.rnc +++ b/src/templates/settings.rnc @@ -112,5 +112,18 @@ element settings { attribute lower_left { list { xsd:double+ } }) & (element upper_right { list { xsd:double+ } } | attribute upper_right { list { xsd:double+ } }) + }? & + + element resonance_scattering { + (element method { xsd:positiveInteger } | + attribute method { xsd:positiveInteger }) & + (element scatterer { xsd:positiveInteger } | + attribute scatterer { xsd:positiveInteger }) & + (element xs_label { xsd:positiveInteger } | + attribute xs_label { xsd:positiveInteger }) & + (element xs_label_0K { xsd:positiveInteger } | + attribute xs_label_0K { xsd:positiveInteger }) & + (element E_min { xsd:double } | attribute E_min { xsd:double }) & + (element E_max { xsd:double } | attribute E_max { xsd:double }) }? } diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index 02ee8a0f1..c00826f81 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -42,6 +42,15 @@ + + + + + + + + + @@ -63,6 +72,7 @@ + From 882d7c3cb78e3ce2e79d5649de9856b3c14c3eb8 Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 22 Nov 2013 09:23:51 -0500 Subject: [PATCH 02/45] added type for 0K nuclides --- src/ace_header.F90 | 57 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 5e2f39a00..0048d455b 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -143,6 +143,40 @@ module ace_header procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide +!=============================================================================== +! NUCLIDE_0K contains all the data for an ACE-format continuous-energy cross +! section that is required for treatment of resonance scattering. The ACE format +! (A Compact ENDF format) is used in MCNP and several other Monte Carlo codes. +!=============================================================================== + + type Nuclide_0K + character(10) :: name ! name of nuclide, e.g. 92235.03c + integer :: zaid ! Z and A identifier, e.g. 92235 + integer :: listing ! index in xs_listings + real(8) :: awr ! weight of nucleus in neutron masses + real(8) :: kT ! temperature in MeV (k*T) + + ! Energy grid information + integer :: n_grid ! # of nuclide grid points + integer, allocatable :: grid_index(:) ! pointers to union grid + real(8), allocatable :: energy(:) ! energy values corresponding to xs + + ! CDF of neutron velocity x cross section + real(8), allocatable :: xs_cdf(:) + + ! Microscopic elastic cross section + real(8), allocatable :: elastic(:) ! elastic scattering + + ! Unresolved resonance data + logical :: urr_present + integer :: urr_inelastic + type(UrrData), pointer :: urr_data => null() + + ! Type-Bound procedures + contains + procedure :: clear => nuclide_clear ! Deallocates Nuclide + end type Nuclide_0K + !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -361,4 +395,27 @@ module ace_header end subroutine nuclide_clear +!=============================================================================== +! NUCLIDE_0K_CLEAR resets and deallocates data in Nuclide_0K. +!=============================================================================== + + subroutine nuclide_0K_clear(this) + + class(Nuclide_0K), intent(inout) :: this ! The Nuclide object to clear + + integer :: i ! Loop counter + + if (allocated(this % grid_index)) & + deallocate(this % grid_index) + + if (allocated(this % energy)) & + deallocate(this % elastic, this % xs_cdf) + + if (associated(this % urr_data)) then + call this % urr_data % clear() + deallocate(this % urr_data) + end if + + end subroutine nuclide_0K_clear + end module ace_header From 46d23f60a28f60e445af5caed6fb34e4eb6ec574 Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 22 Nov 2013 09:59:33 -0500 Subject: [PATCH 03/45] added 0K xs global variables --- src/global.F90 | 18 ++++++++++++++---- 1 file changed, 14 insertions(+), 4 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 637ea315e..4b0a4b1d0 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS + MaterialMacroXS, Nuclide_0K use bank_header, only: Bank use cmfd_header use constants @@ -61,6 +61,7 @@ module global ! Cross section arrays type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections + type(Nuclide_0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings @@ -68,9 +69,10 @@ module global type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide type(MaterialMacroXS) :: material_xs ! Cache for current material - integer :: n_nuclides_total ! Number of nuclide cross section tables - integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables - integer :: n_listings ! Number of listings in cross_sections.xml + integer :: n_nuclides_total ! Number of nuclide cross section tables + integer :: n_nuclides_0K_total ! Number of 0K nuclide cross section tables + integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables + integer :: n_listings ! Number of listings in cross_sections.xml ! Dictionaries to look up cross sections and listings type(DictCharInt) :: nuclide_dict @@ -401,6 +403,14 @@ contains end do deallocate(nuclides) end if + ! Deallocate 0K cross section data, listings + if (allocated(nuclides_0K)) then + ! First call the clear routines + do i = 1, size(nuclides_0K) + call nuclides_0K(i) % clear() + end do + deallocate(nuclides_0K) + end if if (allocated(sab_tables)) deallocate(sab_tables) if (allocated(xs_listings)) deallocate(xs_listings) if (allocated(micro_xs)) deallocate(micro_xs) From f9311b5780ffb48a31cbe106a251c31e0f5df23b Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 22 Nov 2013 10:00:32 -0500 Subject: [PATCH 04/45] modified template files to allow for multiple scatterers --- src/templates/settings.rnc | 16 ++++++++-------- src/templates/settings_t.xml | 8 ++++---- 2 files changed, 12 insertions(+), 12 deletions(-) diff --git a/src/templates/settings.rnc b/src/templates/settings.rnc index e2f5757ad..f846e950d 100644 --- a/src/templates/settings.rnc +++ b/src/templates/settings.rnc @@ -115,14 +115,14 @@ element settings { }? & element resonance_scattering { - (element method { xsd:positiveInteger } | - attribute method { xsd:positiveInteger }) & - (element scatterer { xsd:positiveInteger } | - attribute scatterer { xsd:positiveInteger }) & - (element xs_label { xsd:positiveInteger } | - attribute xs_label { xsd:positiveInteger }) & - (element xs_label_0K { xsd:positiveInteger } | - attribute xs_label_0K { xsd:positiveInteger }) & + (element method { xsd:string { maxLength = "16" } } | + attribute method { xsd:string { maxLength = "16" } }) & + (element scatterer { list { xsd:string { maxLength = "12" }+ } } | + attribute scatterer { list { xsd:string { maxLength = "12" }+ } }) & + (element xs_label { list { xsd:string { maxLength = "12" }+ } } | + attribute xs_label { list { xsd:string { maxLength = "12" }+ } }) & + (element xs_label_0K { list { xsd:string { maxLength = "12" }+ } } | + attribute xs_label_0K { list { xsd:string { maxLength = "12" }+ } }) & (element E_min { xsd:double } | attribute E_min { xsd:double }) & (element E_max { xsd:double } | attribute E_max { xsd:double }) }? diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index c00826f81..7fbf624d2 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -43,10 +43,10 @@ - - - - + + + + From 66aff51ccbddc44a664cbcd3ee45c475ba58e3f2 Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 22 Nov 2013 10:20:23 -0500 Subject: [PATCH 05/45] added logical indicating resonant nuclides --- src/ace_header.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 0048d455b..1be53134a 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -108,6 +108,9 @@ module ace_header real(8), allocatable :: absorption(:) ! absorption (MT > 100) real(8), allocatable :: heating(:) ! heating + ! Resonance scattering info + logical :: resonant ! is nuclide to be treated as a resonant scatterer? + ! Fission information logical :: fissionable ! nuclide is fissionable? logical :: has_partial_fission ! nuclide has partial fission reactions? From e1cedd23b953290e4cf72bd8c259902f3ea3fd5f Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 22 Nov 2013 11:55:41 -0500 Subject: [PATCH 06/45] first crack at weight correction method (WCM) target sampling --- src/physics.F90 | 170 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 130 insertions(+), 40 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index fb1f89194..411722efe 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -627,14 +627,19 @@ contains !=============================================================================== ! SAMPLE_TARGET_VELOCITY samples the target velocity based on the free gas ! scattering formulation used by most Monte Carlo codes. Excellent documentation -! for this method can be found in FRA-TM-123. +! for this method can be found in FRA-TM-123. Methods for correctly accounting +! for the energy dependence of cross sections in treating resonance elastic +! scattering such as the DBRC, WCM, and a new, accelerated scheme are also +! implemented here. !=============================================================================== - subroutine sample_target_velocity(nuc, v_target, E, uvw) + subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt) type(Nuclide), pointer :: nuc real(8), intent(out) :: v_target(3) + real(8), intent(in) :: v_neut(3) real(8), intent(in) :: E + real(8), intent(inout) :: wgt real(8), intent(in) :: uvw(3) real(8) :: kT ! equilibrium temperature of target in MeV @@ -647,62 +652,147 @@ contains real(8) :: beta_vt ! beta * speed of target real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target + real(8) :: E_old ! tmp storage of current energy + logical :: reject ! resample if true + real(8) :: E_rel ! trial relative energy + real(8) :: xs_0K ! 0K xs at E_rel + real(8) :: wcf ! weight correction factor ! Determine equilibrium temperature in MeV kT = nuc % kT - ! Check if energy is above threshold - if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then - v_target = ZERO - return + ! Check if nuclide is a resonant scatterer and which sampling scheme + ! to use based on neutron energy + if (nuc % resonant) then + if (E > E_max_res_scat) then + v_target = ZERO + return + else if (E < E_min_res_scat) then + target_sampling = 'cxs' + end if + else + if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then + v_target = ZERO + return + else + target_sampling = 'cxs' end if - ! calculate beta - beta_vn = sqrt(nuc%awr * E / kT) + E_old = E - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) + ! reject unless criteria are satisfied + reject = .true. - do - ! Sample two random numbers - r1 = prn() - r2 = prn() + select case (target_sampling) + case ('svt') + ! calculate beta, alpha + beta_vn = sqrt(nuc%awr * E / kT) + alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler + do + ! Sample two random numbers + r1 = prn() + r2 = prn() - beta_vt_sq = -log(r1*r2) + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler + beta_vt_sq = -log(r1*r2) - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if + else + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) + c = cos(PI/TWO * prn()) + beta_vt_sq = -log(r1) - log(r2)*c*c + end if - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) + ! Sample cosine of angle between neutron and target velocity + mu = TWO*prn() - ONE - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do + ! Determine rejection probability + accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & + /(beta_vn + beta_vt) - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) exit + end do - ! determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) + ! determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/nuc % awr) + + ! determine velocity vector of target nucleus based on neutron's velocity + ! and the sampled angle between them + v_target = vt * rotate_angle(uvw, mu) + + case ('wcm') + ! calculate beta, alpha + beta_vn = sqrt(nuc%awr * E / kT) + alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) + + do + ! Sample two random numbers + r1 = prn() + r2 = prn() + + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler + + beta_vt_sq = -log(r1*r2) + + else + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler + + c = cos(PI/TWO * prn()) + beta_vt_sq = -log(r1) - log(r2)*c*c + end if + + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) + + ! Sample cosine of angle between neutron and target velocity + mu = TWO*prn() - ONE + + ! Determine rejection probability + accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & + /(beta_vn + beta_vt) + + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) exit + end do + + ! determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/nuc % awr) + + ! determine velocity vector of target nucleus based on neutron's velocity + ! and the sampled angle between them + v_target = vt * rotate_angle(uvw, mu) + + ! adjust particle weight + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + E = E_rel + if (grid_method == GRID_UNION) call find_energy_index(E) + call calculate_elastic_xs(nuc % i_nuclide_0K, NONE) + xs_0K = micro_xs(i_nuclide_0K) % elastic + wcf = xs_0K / micro_xs(i_nuclide) % elastic + wgt = wcf * wgt + + E = E_old + + case default + message = "Not a recognized resonance scattering treatment!" + call fatal_error() + end select end subroutine sample_target_velocity From 0a106a6eaee718b6fa52ce90cf9ad6887bb3ca67 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 09:49:22 -0800 Subject: [PATCH 07/45] allow for diff E ranges and methods for each scatterer --- src/templates/settings.rnc | 10 ++++++---- src/templates/settings_t.xml | 6 +++--- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/templates/settings.rnc b/src/templates/settings.rnc index f846e950d..735fc6c4f 100644 --- a/src/templates/settings.rnc +++ b/src/templates/settings.rnc @@ -115,15 +115,17 @@ element settings { }? & element resonance_scattering { - (element method { xsd:string { maxLength = "16" } } | - attribute method { xsd:string { maxLength = "16" } }) & + (element method { list { xsd:string { maxLength = "16" }+ } | + attribute method { list { xsd:string { maxLength = "16" }+ }) & (element scatterer { list { xsd:string { maxLength = "12" }+ } } | attribute scatterer { list { xsd:string { maxLength = "12" }+ } }) & (element xs_label { list { xsd:string { maxLength = "12" }+ } } | attribute xs_label { list { xsd:string { maxLength = "12" }+ } }) & (element xs_label_0K { list { xsd:string { maxLength = "12" }+ } } | attribute xs_label_0K { list { xsd:string { maxLength = "12" }+ } }) & - (element E_min { xsd:double } | attribute E_min { xsd:double }) & - (element E_max { xsd:double } | attribute E_max { xsd:double }) + (element E_min { list { xsd:double+ } } | + attribute E_min { list { xsd:double+ } }) & + (element E_max { list { xsd:double+ } } | + attribute E_max { list { xsd:double+ } }) }? } diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index 7fbf624d2..dbceb7550 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -43,12 +43,12 @@ - + - - + + From 9d6d66e33c4cdb4f4136801f72861fed76948946 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 09:51:09 -0800 Subject: [PATCH 08/45] added resonance scattering header --- src/res_scat_header.F90 | 60 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 src/res_scat_header.F90 diff --git a/src/res_scat_header.F90 b/src/res_scat_header.F90 new file mode 100644 index 000000000..3007acc80 --- /dev/null +++ b/src/res_scat_header.F90 @@ -0,0 +1,60 @@ +module res_scat_header + + implicit none + private + + type, public :: ResScatterer + + character(10) :: name ! name of nuclide, e.g. 92235.03c + integer :: zaid ! Z and A identifier, e.g. 92235 + integer :: listing ! index in xs_listings + real(8) :: awr ! weight of nucleus in neutron masses + real(8) :: kT ! temperature in MeV (k*T) + + ! Energy grid information + integer :: n_grid + integer, allocatable :: grid_index(:) ! pointers to union grid + real(8), allocatable :: energy(:) ! energy values corresponding to xs + + ! CDF of neutron velocity x cross section + real(8), allocatable :: xs_cdf(:) + + ! Microscopic elastic cross section + real(8), allocatable :: elastic(:) + + ! lower cutoff energy for resonance scattering + real(8) :: E_min + + ! upper cutoff energy for resonance scattering + real(8) :: E_max + + ! target velocity sampling scheme + character(16) :: scheme + + ! Type-Bound procedures + contains + procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer + + end type ResScatterer + +contains + +!=============================================================================== +! RES_SCATTERER_CLEAR resets and deallocates data in ResScatterer. +!=============================================================================== + + subroutine res_scatterer_clear(this) + + class(ResScatterer), intent(inout) :: this ! the ResScatterer object to clear + + integer :: i ! Loop counter + + if (allocated(this % grid_index)) & + deallocate(this % grid_index) + + if (allocated(this % energy)) & + deallocate(this % elastic, this % xs_cdf, this % energy) + + end subroutine res_scatterer_clear + +end module res_scat_header From 03edddff8a579051fee61dc60e9fae792b6040c8 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 09:52:52 -0800 Subject: [PATCH 09/45] cleaned up variables to switch to the ones defined in res_scat_header --- src/ace_header.F90 | 58 +--------- src/global.F90 | 35 ++++-- src/physics.F90 | 271 +++++++++++++++++++++++---------------------- 3 files changed, 162 insertions(+), 202 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 1be53134a..4f0bb6e5d 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -110,6 +110,7 @@ module ace_header ! Resonance scattering info logical :: resonant ! is nuclide to be treated as a resonant scatterer? + integer :: i_0K ! index into 0K nuclides (resonant scatterers) ! Fission information logical :: fissionable ! nuclide is fissionable? @@ -146,40 +147,6 @@ module ace_header procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide -!=============================================================================== -! NUCLIDE_0K contains all the data for an ACE-format continuous-energy cross -! section that is required for treatment of resonance scattering. The ACE format -! (A Compact ENDF format) is used in MCNP and several other Monte Carlo codes. -!=============================================================================== - - type Nuclide_0K - character(10) :: name ! name of nuclide, e.g. 92235.03c - integer :: zaid ! Z and A identifier, e.g. 92235 - integer :: listing ! index in xs_listings - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) - - ! Energy grid information - integer :: n_grid ! # of nuclide grid points - integer, allocatable :: grid_index(:) ! pointers to union grid - real(8), allocatable :: energy(:) ! energy values corresponding to xs - - ! CDF of neutron velocity x cross section - real(8), allocatable :: xs_cdf(:) - - ! Microscopic elastic cross section - real(8), allocatable :: elastic(:) ! elastic scattering - - ! Unresolved resonance data - logical :: urr_present - integer :: urr_inelastic - type(UrrData), pointer :: urr_data => null() - - ! Type-Bound procedures - contains - procedure :: clear => nuclide_clear ! Deallocates Nuclide - end type Nuclide_0K - !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -398,27 +365,4 @@ module ace_header end subroutine nuclide_clear -!=============================================================================== -! NUCLIDE_0K_CLEAR resets and deallocates data in Nuclide_0K. -!=============================================================================== - - subroutine nuclide_0K_clear(this) - - class(Nuclide_0K), intent(inout) :: this ! The Nuclide object to clear - - integer :: i ! Loop counter - - if (allocated(this % grid_index)) & - deallocate(this % grid_index) - - if (allocated(this % energy)) & - deallocate(this % elastic, this % xs_cdf) - - if (associated(this % urr_data)) then - call this % urr_data % clear() - deallocate(this % urr_data) - end if - - end subroutine nuclide_0K_clear - end module ace_header diff --git a/src/global.F90 b/src/global.F90 index 4b0a4b1d0..02f618fbb 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS, Nuclide_0K + MaterialMacroXS use bank_header, only: Bank use cmfd_header use constants @@ -11,6 +11,7 @@ module global use mesh_header, only: StructuredMesh use plot_header, only: ObjectPlot use set_header, only: SetInt + use res_scat_header, only: ResScatterer use source_header, only: ExtSource use tally_header, only: TallyObject, TallyMap, TallyResult use timer_header, only: Timer @@ -61,7 +62,6 @@ module global ! Cross section arrays type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections - type(Nuclide_0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides type(SAlphaBeta), allocatable, target :: sab_tables(:) ! S(a,b) tables type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings @@ -70,7 +70,6 @@ module global type(MaterialMacroXS) :: material_xs ! Cache for current material integer :: n_nuclides_total ! Number of nuclide cross section tables - integer :: n_nuclides_0K_total ! Number of 0K nuclide cross section tables integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables integer :: n_listings ! Number of listings in cross_sections.xml @@ -369,7 +368,19 @@ module global logical :: output_summary = .false. logical :: output_xs = .false. logical :: output_tallies = .true. + + ! ============================================================================ + ! RESONANCE SCATTERING VARIABLES + ! Is resonance scattering treated? + logical :: treat_res_scat = .false. + + ! Number of resonant scatterers + integer :: n_res_scatterers_total + + ! Main object + type(ResScatterer), allocatable, target :: res_scatterers(:) + !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & !$omp& trace, thread_id, current_work, matching_bins) @@ -403,18 +414,20 @@ contains end do deallocate(nuclides) end if - ! Deallocate 0K cross section data, listings - if (allocated(nuclides_0K)) then - ! First call the clear routines - do i = 1, size(nuclides_0K) - call nuclides_0K(i) % clear() - end do - deallocate(nuclides_0K) - end if + if (allocated(sab_tables)) deallocate(sab_tables) if (allocated(xs_listings)) deallocate(xs_listings) if (allocated(micro_xs)) deallocate(micro_xs) + ! Deallocate resonance scattering data + if (allocated(res_scatterers)) then + ! First call the clear routines + do i = 1, size(res_scatterers) + call res_scatterers(i) % clear() + end do + deallocate(res_scatterers) + end if + ! Deallocate external source if (allocated(external_source % params_space)) & deallocate(external_source % params_space) diff --git a/src/physics.F90 b/src/physics.F90 index 411722efe..abe1c3303 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,6 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants + use cross_section, only: elastic_0K_xs use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed @@ -14,6 +15,7 @@ module physics use particle_header, only: Particle use particle_restart_write, only: write_particle_restart use random_lcg, only: prn + use res_scat_header, only: ResScatterer use search, only: binary_search use string, only: to_str @@ -344,7 +346,7 @@ contains ! Perform collision physics for elastic scattering call elastic_scatter(i_nuclide, rxn, & - p % E, p % coord0 % uvw, p % mu) + p % E, p % coord0 % uvw, p % mu, p % wgt) end if p % event_MT = ELASTIC @@ -403,13 +405,14 @@ contains ! target. !=============================================================================== - subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab) + subroutine elastic_scatter(i_nuclide, rxn, E, uvw, mu_lab, wgt) integer, intent(in) :: i_nuclide type(Reaction), pointer :: rxn real(8), intent(inout) :: E real(8), intent(inout) :: uvw(3) real(8), intent(out) :: mu_lab + real(8), intent(inout) :: wgt real(8) :: awr ! atomic weight ratio of target real(8) :: mu_cm ! cosine of polar angle in center-of-mass @@ -432,7 +435,7 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then - call sample_target_velocity(nuc, v_t, E, uvw) + call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, micro_xs(i_nuclide) % elastic) else v_t = ZERO end if @@ -633,14 +636,16 @@ contains ! implemented here. !=============================================================================== - subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt) + subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - type(Nuclide), pointer :: nuc - real(8), intent(out) :: v_target(3) - real(8), intent(in) :: v_neut(3) - real(8), intent(in) :: E - real(8), intent(inout) :: wgt - real(8), intent(in) :: uvw(3) + type(Nuclide), pointer :: nuc ! target nuclide at temperature + type(ResScatterer), pointer :: nuc_0K => null() ! (resonant) target nuclide at 0K + + real(8), intent(out) :: v_target(3) + real(8), intent(in) :: v_neut(3) + real(8), intent(in) :: E + real(8), intent(in) :: uvw(3) + real(8), intent(inout) :: wgt real(8) :: kT ! equilibrium temperature of target in MeV real(8) :: alpha ! probability of sampling f2 over f1 @@ -652,148 +657,146 @@ contains real(8) :: beta_vt ! beta * speed of target real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target - real(8) :: E_old ! tmp storage of current energy - logical :: reject ! resample if true +! real(8) :: E_old ! tmp storage of current energy real(8) :: E_rel ! trial relative energy real(8) :: xs_0K ! 0K xs at E_rel + real(8) :: xs_eff ! effective elastic xs at temperature T real(8) :: wcf ! weight correction factor + logical :: reject ! resample if true + + nuc_0K => res_scatterers(nuc % i_0K) + ! Determine equilibrium temperature in MeV kT = nuc % kT ! Check if nuclide is a resonant scatterer and which sampling scheme ! to use based on neutron energy if (nuc % resonant) then - if (E > E_max_res_scat) then - v_target = ZERO - return - else if (E < E_min_res_scat) then - target_sampling = 'cxs' - end if + if (E > nuc_0K % E_max) then + v_target = ZERO + return + else if (E < nuc_0K % E_min) then + nuc_0K % scheme = 'cxs' + end if else - if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then - v_target = ZERO - return - else - target_sampling = 'cxs' + if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then + v_target = ZERO + return + else + nuc_0K % scheme = 'cxs' + end if end if - E_old = E - ! reject unless criteria are satisfied reject = .true. - - select case (target_sampling) - case ('svt') - ! calculate beta, alpha - beta_vn = sqrt(nuc%awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) - - ! determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - + + select case (nuc_0K % scheme) + case ('cxs') + ! calculate beta, alpha + beta_vn = sqrt(nuc%awr * E / kT) + alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) + + do + ! Sample two random numbers + r1 = prn() + r2 = prn() + + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler + + beta_vt_sq = -log(r1*r2) + + else + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler + + c = cos(PI/TWO * prn()) + beta_vt_sq = -log(r1) - log(r2)*c*c + end if + + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) + + ! Sample cosine of angle between neutron and target velocity + mu = TWO*prn() - ONE + + ! Determine rejection probability + accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & + /(beta_vn + beta_vt) + + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) exit + end do + + ! determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/nuc % awr) + + ! determine velocity vector of target nucleus based on neutron's velocity + ! and the sampled angle between them + v_target = vt * rotate_angle(uvw, mu) + case ('wcm') - ! calculate beta, alpha - beta_vn = sqrt(nuc%awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) - - ! determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - - ! adjust particle weight - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - E = E_rel - if (grid_method == GRID_UNION) call find_energy_index(E) - call calculate_elastic_xs(nuc % i_nuclide_0K, NONE) - xs_0K = micro_xs(i_nuclide_0K) % elastic - wcf = xs_0K / micro_xs(i_nuclide) % elastic - wgt = wcf * wgt - - E = E_old - + ! calculate beta, alpha + beta_vn = sqrt(nuc%awr * E / kT) + alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) + + do + ! Sample two random numbers + r1 = prn() + r2 = prn() + + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler + + beta_vt_sq = -log(r1*r2) + + else + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler + + c = cos(PI/TWO * prn()) + beta_vt_sq = -log(r1) - log(r2)*c*c + end if + + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) + + ! Sample cosine of angle between neutron and target velocity + mu = TWO*prn() - ONE + + ! Determine rejection probability + accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & + /(beta_vn + beta_vt) + + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) exit + end do + + ! determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/nuc % awr) + + ! determine velocity vector of target nucleus based on neutron's velocity + ! and the sampled angle between them + v_target = vt * rotate_angle(uvw, mu) + + ! adjust particle weight + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + xs_0K = elastic_0K_xs(E_rel, nuc_0K, nuc % i_0K) + wcf = xs_0K / xs_eff + wgt = wcf * wgt + case default - message = "Not a recognized resonance scattering treatment!" - call fatal_error() + message = "Not a recognized resonance scattering treatment!" + call fatal_error() end select - + end subroutine sample_target_velocity !=============================================================================== From 4c9dd3f094b9c238e2d09ccb7ff3c321479d1aa5 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 09:53:45 -0800 Subject: [PATCH 10/45] added required dependencies --- src/DEPENDENCIES | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index e404cd991..afb953272 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -158,6 +158,7 @@ global.o: hdf5_interface.o global.o: material_header.o global.o: mesh_header.o global.o: plot_header.o +global.o: res_scat_header.o global.o: set_header.o global.o: source_header.o global.o: tally_header.o @@ -289,6 +290,7 @@ particle_restart_write.o: string.o physics.o: ace_header.o physics.o: constants.o +physics.o: cross_section.o physics.o: endf.o physics.o: error.o physics.o: fission.o From bead663f806dd48a5b3e94a936e95223e9a74740 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 09:54:34 -0800 Subject: [PATCH 11/45] function for calculating 0K xs --- src/cross_section.F90 | 51 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 299033590..e05b18949 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -478,4 +478,55 @@ contains end subroutine find_energy_index +!=============================================================================== +! ELASTIC_0K_XS determines the microscopic elastic cross section for a +! nuclide of a given index in the nuclides array at the energy of the particle +!=============================================================================== + + function elastic_0K_xs(E, nuc_0K, i_0K) result(xs_out) + + type(ResScatterer), pointer :: nuc_0K + + integer, intent(in) :: i_0K ! index into 0K nuclides array + integer :: i_grid ! index on nuclide energy grid + + real(8) :: f ! interp factor on nuclide energy grid + real(8), intent(in) :: E ! trial energy + real(8) :: xs_out + + ! Determine index on nuclide energy grid + select case (grid_method) + case (GRID_UNION) + ! If we're using the unionized grid with pointers, finding the index on + ! the nuclide energy grid is as simple as looking up the pointer + + call find_energy_index(E) + i_grid = nuc_0K % grid_index(union_grid_index) + + case (GRID_NUCLIDE) + ! If we're not using the unionized grid, we have to do a binary search on + ! the nuclide energy grid in order to determine which points to + ! interpolate between + + if (E < nuc_0K % energy(1)) then + i_grid = 1 + elseif (E > nuc_0K % energy(nuc_0K % n_grid)) then + i_grid = nuc_0K % n_grid - 1 + else + i_grid = binary_search(nuc_0K % energy, nuc_0K % n_grid, E) + end if + + end select + + ! check for rare case where two energy points are the same + if (nuc_0K % energy(i_grid) == nuc_0K % energy(i_grid+1)) i_grid = i_grid + 1 + + ! calculate interpolation factor + f = (E - nuc_0K%energy(i_grid))/(nuc_0K%energy(i_grid+1) - nuc_0K%energy(i_grid)) + + ! Calculate microscopic nuclide elastic cross section + xs_out = (ONE - f) * nuc_0K % elastic(i_grid) + f * nuc_0K % elastic(i_grid+1) + + end function elastic_0K_xs + end module cross_section From 99249940e5ed150b63e9ec49e8f0d0aa296e0154 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 23 Nov 2013 15:12:22 -0800 Subject: [PATCH 12/45] moved contents of res_scat_header to ace_header --- src/DEPENDENCIES | 1 - src/Makefile | 4 +-- src/ace_header.F90 | 56 ++++++++++++++++++++++++++++++++++++++ src/global.F90 | 3 +-- src/input_xml.F90 | 10 +++++++ src/physics.F90 | 3 +-- src/res_scat_header.F90 | 60 ----------------------------------------- 7 files changed, 70 insertions(+), 67 deletions(-) delete mode 100644 src/res_scat_header.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index afb953272..016163ce8 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -158,7 +158,6 @@ global.o: hdf5_interface.o global.o: material_header.o global.o: mesh_header.o global.o: plot_header.o -global.o: res_scat_header.o global.o: set_header.o global.o: source_header.o global.o: tally_header.o diff --git a/src/Makefile b/src/Makefile index 4300d42e6..1a75b941c 100644 --- a/src/Makefile +++ b/src/Makefile @@ -16,7 +16,7 @@ COMPILER = gnu DEBUG = no PROFILE = no OPTIMIZE = no -MPI = no +MPI = yes OPENMP = no HDF5 = no PETSC = no @@ -25,7 +25,7 @@ PETSC = no # External Library Paths #=============================================================================== -MPI_DIR = /opt/mpich/3.0.4-$(COMPILER) +MPI_DIR = /opt/mpich/3.0.2-$(COMPILER) HDF5_DIR = /opt/hdf5/1.8.11-$(COMPILER) PHDF5_DIR = /opt/phdf5/1.8.11-$(COMPILER) PETSC_DIR = /opt/petsc/3.4.2-$(COMPILER) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 4f0bb6e5d..d8b3178e7 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -147,6 +147,44 @@ module ace_header procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide +!=============================================================================== +! RESSCATTERER contains all the data for a resonance scattering nuclide at 0K +!=============================================================================== + + type :: ResScatterer + + character(10) :: name ! name of nuclide, e.g. 92235.03c + integer :: zaid ! Z and A identifier, e.g. 92235 + integer :: listing ! index in xs_listings + real(8) :: awr ! weight of nucleus in neutron masses + real(8) :: kT ! temperature in MeV (k*T) + + ! Energy grid information + integer :: n_grid + integer, allocatable :: grid_index(:) ! pointers to union grid + real(8), allocatable :: energy(:) ! energy values corresponding to xs + + ! CDF of neutron velocity x cross section + real(8), allocatable :: xs_cdf(:) + + ! Microscopic elastic cross section + real(8), allocatable :: elastic(:) + + ! lower cutoff energy for resonance scattering + real(8) :: E_min + + ! upper cutoff energy for resonance scattering + real(8) :: E_max + + ! target velocity sampling scheme + character(16) :: scheme + + ! Type-Bound procedures + contains + procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer + + end type ResScatterer + !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -365,4 +403,22 @@ module ace_header end subroutine nuclide_clear +!=============================================================================== +! RES_SCATTERER_CLEAR resets and deallocates data in ResScatterer. +!=============================================================================== + + subroutine res_scatterer_clear(this) + + class(ResScatterer), intent(inout) :: this ! the ResScatterer object to clear + + integer :: i ! Loop counter + + if (allocated(this % grid_index)) & + deallocate(this % grid_index) + + if (allocated(this % energy)) & + deallocate(this % elastic, this % xs_cdf, this % energy) + + end subroutine res_scatterer_clear + end module ace_header diff --git a/src/global.F90 b/src/global.F90 index 02f618fbb..2b634db40 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS + MaterialMacroXS, ResScatterer use bank_header, only: Bank use cmfd_header use constants @@ -11,7 +11,6 @@ module global use mesh_header, only: StructuredMesh use plot_header, only: ObjectPlot use set_header, only: SetInt - use res_scat_header, only: ResScatterer use source_header, only: ExtSource use tally_header, only: TallyObject, TallyMap, TallyResult use timer_header, only: Timer diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9bedbc8b3..536628ba3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -602,6 +602,16 @@ contains #endif end if + ! set up resonance scattering treatment + call lower_case(resonance_scattering_ % method(1)) + if (resonance_scattering_ % method(1) /= '') then + treat_res_scat = .true. + n_res_scatterers_total = size(resonance_scattering_ % scatterer) + do i = 1, n_res_scatterers_total + + end do + end if + end subroutine read_settings_xml !=============================================================================== diff --git a/src/physics.F90 b/src/physics.F90 index abe1c3303..fdd4f50f8 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,6 +1,6 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy + use ace_header, only: Nuclide, Reaction, DistEnergy, ResScatterer use constants use cross_section, only: elastic_0K_xs use endf, only: reaction_name @@ -15,7 +15,6 @@ module physics use particle_header, only: Particle use particle_restart_write, only: write_particle_restart use random_lcg, only: prn - use res_scat_header, only: ResScatterer use search, only: binary_search use string, only: to_str diff --git a/src/res_scat_header.F90 b/src/res_scat_header.F90 deleted file mode 100644 index 3007acc80..000000000 --- a/src/res_scat_header.F90 +++ /dev/null @@ -1,60 +0,0 @@ -module res_scat_header - - implicit none - private - - type, public :: ResScatterer - - character(10) :: name ! name of nuclide, e.g. 92235.03c - integer :: zaid ! Z and A identifier, e.g. 92235 - integer :: listing ! index in xs_listings - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) - - ! Energy grid information - integer :: n_grid - integer, allocatable :: grid_index(:) ! pointers to union grid - real(8), allocatable :: energy(:) ! energy values corresponding to xs - - ! CDF of neutron velocity x cross section - real(8), allocatable :: xs_cdf(:) - - ! Microscopic elastic cross section - real(8), allocatable :: elastic(:) - - ! lower cutoff energy for resonance scattering - real(8) :: E_min - - ! upper cutoff energy for resonance scattering - real(8) :: E_max - - ! target velocity sampling scheme - character(16) :: scheme - - ! Type-Bound procedures - contains - procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer - - end type ResScatterer - -contains - -!=============================================================================== -! RES_SCATTERER_CLEAR resets and deallocates data in ResScatterer. -!=============================================================================== - - subroutine res_scatterer_clear(this) - - class(ResScatterer), intent(inout) :: this ! the ResScatterer object to clear - - integer :: i ! Loop counter - - if (allocated(this % grid_index)) & - deallocate(this % grid_index) - - if (allocated(this % energy)) & - deallocate(this % elastic, this % xs_cdf, this % energy) - - end subroutine res_scatterer_clear - -end module res_scat_header From 78746e1a928a07100e051626e123704065f11b5b Mon Sep 17 00:00:00 2001 From: walshjon Date: Sun, 24 Nov 2013 16:40:41 -0800 Subject: [PATCH 13/45] started setting up reading in 0K xs --- src/ace.F90 | 204 +++++++++++++++++++++++++++++++++++++++++++++ src/ace_header.F90 | 43 +++++----- src/input_xml.F90 | 8 +- 3 files changed, 233 insertions(+), 22 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index d54aa83c2..8ae447d4c 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -35,6 +35,8 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: k ! index over S(a,b) tables in material + integer :: l ! index over 0K nuclides + integer :: index_res_scatterers ! index into res_scatterers integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables @@ -58,6 +60,8 @@ contains ! ========================================================================== ! READ ALL ACE CROSS SECTION TABLES + index_res_scatterers = 0 + ! Loop over all files MATERIAL_LOOP: do i = 1, n_materials mat => materials(i) @@ -71,6 +75,15 @@ contains name = xs_listings(i_listing) % name alias = xs_listings(i_listing) % alias + ! check if nuclide's a resonant scatterer + do l = 1, n_res_scatterers_total + if (name == res_scatterers(l) % name) then + nuclides(i_nuclide) % resonant = .true. + index_res_scatterers = index_res_scatterers + 1 + nuclides(i_nuclide) % i_0k = index_res_scatterers + end if + end do + ! Keep track of what listing is associated with this nuclide nuc => nuclides(i_nuclide) nuc % listing = i_listing @@ -170,6 +183,26 @@ contains ! Avoid some valgrind leak errors call already_read % clear() + + ! check if nuclide's a resonant scatterer + do l = 1, n_res_scatterers_total + name = res_scatterers(l) % scatterer + i_listing = xs_listing_dict % get_key(name) + name = xs_listings(i_listing) % name + alias = xs_listings(i_listing) % alias + + nuc => nuclides(i_nuclide) + nuc % listing = i_listing + + ! Read the ACE table into the appropriate entry on the nuclides + ! array + call read_ace_table(i_nuclide, i_listing) + + ! Add name and alias to dictionary + call already_read % add(name) + call already_read % add(alias) + end if + end do end subroutine read_xs @@ -344,6 +377,177 @@ contains end subroutine read_ace_table +!=============================================================================== +! READ_0K_ACE_TABLE reads a single cross section table in either ASCII or binary +! format. This routine reads the header data for each table and then calls +! appropriate subroutines to parse the actual data. +!=============================================================================== + + subroutine read_0K_ace_table(i_table, i_listing) + + integer, intent(in) :: i_table ! index in nuclides/sab_tables + integer, intent(in) :: i_listing ! index in xs_listings + + integer :: i ! loop index for XSS records + integer :: j, j1, j2 ! indices in XSS + integer :: record_length ! Fortran record length + integer :: location ! location of ACE table + integer :: entries ! number of entries on each record + integer :: length ! length of ACE table + integer :: in = 7 ! file unit + integer :: zaids(16) ! list of ZAIDs (only used for S(a,b)) + integer :: filetype ! filetype (ASCII or BINARY) + real(8) :: kT ! temperature of table + real(8) :: awrs(16) ! list of atomic weight ratios (not used) + real(8) :: awr ! atomic weight ratio for table + logical :: file_exists ! does ACE library exist? + character(7) :: readable ! is ACE library readable? + character(10) :: name ! name of ACE table + character(10) :: date_ ! date ACE library was processed + character(10) :: mat ! material identifier + character(70) :: comment ! comment for ACE table + character(MAX_FILE_LEN) :: filename ! path to ACE cross section library + type(Nuclide), pointer :: nuc => null() + type(SAlphaBeta), pointer :: sab => null() + type(XsListing), pointer :: listing => null() + + ! determine path, record length, and location of table + listing => xs_listings(i_listing) + filename = listing % path + record_length = listing % recl + location = listing % location + entries = listing % entries + filetype = listing % filetype + + ! Check if ACE library exists and is readable + inquire(FILE=filename, EXIST=file_exists, READ=readable) + if (.not. file_exists) then + message = "ACE library '" // trim(filename) // "' does not exist!" + call fatal_error() + elseif (readable(1:3) == 'NO') then + message = "ACE library '" // trim(filename) // "' is not readable! & + &Change file permissions with chmod command." + call fatal_error() + end if + + ! display message + message = "Loading ACE cross section table: " // listing % name + call write_message(6) + + if (filetype == ASCII) then + ! ======================================================================= + ! READ ACE TABLE IN ASCII FORMAT + + ! Find location of table + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read') + rewind(UNIT=in) + do i = 1, location - 1 + read(UNIT=in, FMT=*) + end do + + ! Read first line of header + read(UNIT=in, FMT='(A10,2G12.0,1X,A10)') name, awr, kT, date_ + + ! Check that correct xs was found -- if cross_sections.xml is broken, the + ! location of the table may be wrong + if(adjustl(name) /= adjustl(listing % name)) then + message = "XS listing entry " // trim(listing % name) // " did not & + &match ACE data, " // trim(name) // " found instead." + call fatal_error() + end if + + ! Read more header and NXS and JXS + read(UNIT=in, FMT=100) comment, mat, & + (zaids(i), awrs(i), i=1,16), NXS, JXS +100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& + ,8I9/8I9/8I9/8I9/8I9/8I9) + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read XSS array + read(UNIT=in, FMT='(4G20.0)') XSS + + ! Close ACE file + close(UNIT=in) + + elseif (filetype == BINARY) then + ! ======================================================================= + ! READ ACE TABLE IN BINARY FORMAT + + ! Open ACE file + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read', & + ACCESS='direct', RECL=record_length) + + ! Read all header information + read(UNIT=in, REC=location) name, awr, kT, date_, & + comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read remaining records with XSS + do i = 1, (length + entries - 1)/entries + j1 = 1 + (i-1)*entries + j2 = min(length, j1 + entries - 1) + read(UNIT=IN, REC=location + i) (XSS(j), j=j1,j2) + end do + + ! Close ACE file + close(UNIT=in) + end if + + ! ========================================================================== + ! PARSE DATA BASED ON NXS, JXS, AND XSS ARRAYS + + select case(listing % type) + case (ACE_NEUTRON) + nuc => nuclides(i_table) + nuc % name = name + nuc % awr = awr + nuc % kT = kT + nuc % zaid = NXS(2) + + ! read all blocks + call read_esz(nuc) + call read_nu_data(nuc) + call read_reactions(nuc) + call read_angular_dist(nuc) + call read_energy_dist(nuc) + call read_unr_res(nuc) + + ! Currently subcritical fixed source calculations are not allowed. Thus, + ! if any fissionable material is found in a fixed source calculation, + ! abort the run. + if (run_mode == MODE_FIXEDSOURCE .and. nuc % fissionable) then + message = "Cannot have fissionable material in a fixed source run." + call fatal_error() + end if + + ! for fissionable nuclides, precalculate microscopic nu-fission cross + ! sections so that we don't need to call the nu_total function during + ! cross section lookups + + if (nuc % fissionable) call generate_nu_fission(nuc) + + case (ACE_THERMAL) + sab => sab_tables(i_table) + sab % name = name + sab % awr = awr + sab % kT = kT + sab % zaid = zaids(1) + + call read_thermal_data(sab) + end select + + deallocate(XSS) + if(associated(nuc)) nullify(nuc) + if(associated(sab)) nullify(sab) + + end subroutine read_0K_ace_table + !=============================================================================== ! READ_ESZ - reads through the ESZ block. This block contains the energy grid, ! total xs, absorption xs, elastic scattering xs, and heating numbers. diff --git a/src/ace_header.F90 b/src/ace_header.F90 index d8b3178e7..425d2a36c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -153,36 +153,37 @@ module ace_header type :: ResScatterer - character(10) :: name ! name of nuclide, e.g. 92235.03c - integer :: zaid ! Z and A identifier, e.g. 92235 - integer :: listing ! index in xs_listings - real(8) :: awr ! weight of nucleus in neutron masses - real(8) :: kT ! temperature in MeV (k*T) + character(10) :: name ! name of nuclide, e.g. 92235.70c + character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c +! integer :: zaid ! Z and A identifier, e.g. 92235 +! integer :: listing ! index in xs_listings +! real(8) :: awr ! weight of nucleus in neutron masses +! real(8) :: kT ! temperature in MeV (k*T) - ! Energy grid information - integer :: n_grid - integer, allocatable :: grid_index(:) ! pointers to union grid - real(8), allocatable :: energy(:) ! energy values corresponding to xs + ! Energy grid information +! integer :: n_grid +! integer, allocatable :: grid_index(:) ! pointers to union grid +! real(8), allocatable :: energy(:) ! energy values corresponding to xs - ! CDF of neutron velocity x cross section - real(8), allocatable :: xs_cdf(:) + ! CDF of neutron velocity x cross section +! real(8), allocatable :: xs_cdf(:) - ! Microscopic elastic cross section - real(8), allocatable :: elastic(:) + ! Microscopic elastic cross section +! real(8), allocatable :: elastic(:) - ! lower cutoff energy for resonance scattering + ! lower cutoff energy for resonance scattering real(8) :: E_min - ! upper cutoff energy for resonance scattering + ! upper cutoff energy for resonance scattering real(8) :: E_max - ! target velocity sampling scheme + ! target velocity sampling scheme character(16) :: scheme - - ! Type-Bound procedures - contains - procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer - + + ! Type-Bound procedures + contains + procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer + end type ResScatterer !=============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 536628ba3..32ea3e2f6 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -607,8 +607,13 @@ contains if (resonance_scattering_ % method(1) /= '') then treat_res_scat = .true. n_res_scatterers_total = size(resonance_scattering_ % scatterer) + allocate(res_scatterers(n_res_scatterers_total)) do i = 1, n_res_scatterers_total - + res_scatterers(i) % name = trim(resonance_scattering_ % xs_label(i)) + res_scatterers(i) % name_0K = trim(resonance_scattering_ % xs_0K_label(i)) + res_scatterers(i) % scheme = trim(resonance_scattering_ % method(i)) + res_scatterers(i) % E_min = resonance_scattering_ % E_min(i) + res_scatterers(i) % E_max = resonance_scattering_ % E_max(i) end do end if @@ -1052,6 +1057,7 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides + integer :: k ! loop index for resonant scatterers integer :: n ! number of nuclides integer :: n_sab ! number of sab tables for a material integer :: index_list ! index in xs_listings array From e77d2d51ff692b054d77e37e1c920c23ec5e5762 Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 26 Nov 2013 10:09:35 -0800 Subject: [PATCH 14/45] incorporated ResScatterer derived type into Nuclide --- src/ace_header.F90 | 71 +++++++--------------------------------------- 1 file changed, 10 insertions(+), 61 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 425d2a36c..b88003207 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -109,8 +109,13 @@ module ace_header real(8), allocatable :: heating(:) ! heating ! Resonance scattering info - logical :: resonant ! is nuclide to be treated as a resonant scatterer? - integer :: i_0K ! index into 0K nuclides (resonant scatterers) + logical :: resonant = .false. ! resonant scatterer? + character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c + character(16) :: scheme ! target velocity sampling scheme + real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section + real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section + real(8) :: E_min ! lower cutoff energy for res scattering + real(8) :: E_max ! upper cutoff energy for res scattering ! Fission information logical :: fissionable ! nuclide is fissionable? @@ -146,46 +151,7 @@ module ace_header contains procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide - -!=============================================================================== -! RESSCATTERER contains all the data for a resonance scattering nuclide at 0K -!=============================================================================== - - type :: ResScatterer - - character(10) :: name ! name of nuclide, e.g. 92235.70c - character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c -! integer :: zaid ! Z and A identifier, e.g. 92235 -! integer :: listing ! index in xs_listings -! real(8) :: awr ! weight of nucleus in neutron masses -! real(8) :: kT ! temperature in MeV (k*T) - - ! Energy grid information -! integer :: n_grid -! integer, allocatable :: grid_index(:) ! pointers to union grid -! real(8), allocatable :: energy(:) ! energy values corresponding to xs - - ! CDF of neutron velocity x cross section -! real(8), allocatable :: xs_cdf(:) - - ! Microscopic elastic cross section -! real(8), allocatable :: elastic(:) - - ! lower cutoff energy for resonance scattering - real(8) :: E_min - - ! upper cutoff energy for resonance scattering - real(8) :: E_max - - ! target velocity sampling scheme - character(16) :: scheme - ! Type-Bound procedures - contains - procedure :: clear => res_scatterer_clear ! Deallocates resonant scatterer - - end type ResScatterer - !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -363,8 +329,9 @@ module ace_header deallocate(this % grid_index) if (allocated(this % energy)) & - deallocate(this % total, this % elastic, this % fission, & - this % nu_fission, this % absorption) + deallocate(this % total, this % elastic, this % fission, & + this % nu_fission, this % absorption, & + this % elastic_0K, this % xs_cdf) if (allocated(this % heating)) & deallocate(this % heating) @@ -404,22 +371,4 @@ module ace_header end subroutine nuclide_clear -!=============================================================================== -! RES_SCATTERER_CLEAR resets and deallocates data in ResScatterer. -!=============================================================================== - - subroutine res_scatterer_clear(this) - - class(ResScatterer), intent(inout) :: this ! the ResScatterer object to clear - - integer :: i ! Loop counter - - if (allocated(this % grid_index)) & - deallocate(this % grid_index) - - if (allocated(this % energy)) & - deallocate(this % elastic, this % xs_cdf, this % energy) - - end subroutine res_scatterer_clear - end module ace_header From e69706b92188e8b15d6f7faf68228258bcac2342 Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 26 Nov 2013 10:12:51 -0800 Subject: [PATCH 15/45] removed separate 0K subroutines --- src/ace.F90 | 215 +++------------------------------------------------- 1 file changed, 11 insertions(+), 204 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 8ae447d4c..0cb97c647 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -35,11 +35,11 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: k ! index over S(a,b) tables in material - integer :: l ! index over 0K nuclides - integer :: index_res_scatterers ! index into res_scatterers + integer :: l ! index over res_scatterers integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables + integer :: i_res_scat = 0 ! index in res_scatterers integer :: m ! position for sorting integer :: temp_nuclide ! temporary value for sorting integer :: temp_table ! temporary value for sorting @@ -60,8 +60,6 @@ contains ! ========================================================================== ! READ ALL ACE CROSS SECTION TABLES - index_res_scatterers = 0 - ! Loop over all files MATERIAL_LOOP: do i = 1, n_materials mat => materials(i) @@ -75,19 +73,19 @@ contains name = xs_listings(i_listing) % name alias = xs_listings(i_listing) % alias - ! check if nuclide's a resonant scatterer - do l = 1, n_res_scatterers_total - if (name == res_scatterers(l) % name) then - nuclides(i_nuclide) % resonant = .true. - index_res_scatterers = index_res_scatterers + 1 - nuclides(i_nuclide) % i_0k = index_res_scatterers - end if - end do - ! Keep track of what listing is associated with this nuclide nuc => nuclides(i_nuclide) nuc % listing = i_listing + do l = 1, n_res_scatterers_total + if (name == res_scatterers(l) % name) then + i_res_scat = i_res_scat + 1 + nuclides(i_nuclide) % resonant = .true. + nuclides(i_nuclide) % i_0K = i_res_scat + exit + end if + end do + ! Read the ACE table into the appropriate entry on the nuclides ! array call read_ace_table(i_nuclide, i_listing) @@ -183,26 +181,6 @@ contains ! Avoid some valgrind leak errors call already_read % clear() - - ! check if nuclide's a resonant scatterer - do l = 1, n_res_scatterers_total - name = res_scatterers(l) % scatterer - i_listing = xs_listing_dict % get_key(name) - name = xs_listings(i_listing) % name - alias = xs_listings(i_listing) % alias - - nuc => nuclides(i_nuclide) - nuc % listing = i_listing - - ! Read the ACE table into the appropriate entry on the nuclides - ! array - call read_ace_table(i_nuclide, i_listing) - - ! Add name and alias to dictionary - call already_read % add(name) - call already_read % add(alias) - end if - end do end subroutine read_xs @@ -377,177 +355,6 @@ contains end subroutine read_ace_table -!=============================================================================== -! READ_0K_ACE_TABLE reads a single cross section table in either ASCII or binary -! format. This routine reads the header data for each table and then calls -! appropriate subroutines to parse the actual data. -!=============================================================================== - - subroutine read_0K_ace_table(i_table, i_listing) - - integer, intent(in) :: i_table ! index in nuclides/sab_tables - integer, intent(in) :: i_listing ! index in xs_listings - - integer :: i ! loop index for XSS records - integer :: j, j1, j2 ! indices in XSS - integer :: record_length ! Fortran record length - integer :: location ! location of ACE table - integer :: entries ! number of entries on each record - integer :: length ! length of ACE table - integer :: in = 7 ! file unit - integer :: zaids(16) ! list of ZAIDs (only used for S(a,b)) - integer :: filetype ! filetype (ASCII or BINARY) - real(8) :: kT ! temperature of table - real(8) :: awrs(16) ! list of atomic weight ratios (not used) - real(8) :: awr ! atomic weight ratio for table - logical :: file_exists ! does ACE library exist? - character(7) :: readable ! is ACE library readable? - character(10) :: name ! name of ACE table - character(10) :: date_ ! date ACE library was processed - character(10) :: mat ! material identifier - character(70) :: comment ! comment for ACE table - character(MAX_FILE_LEN) :: filename ! path to ACE cross section library - type(Nuclide), pointer :: nuc => null() - type(SAlphaBeta), pointer :: sab => null() - type(XsListing), pointer :: listing => null() - - ! determine path, record length, and location of table - listing => xs_listings(i_listing) - filename = listing % path - record_length = listing % recl - location = listing % location - entries = listing % entries - filetype = listing % filetype - - ! Check if ACE library exists and is readable - inquire(FILE=filename, EXIST=file_exists, READ=readable) - if (.not. file_exists) then - message = "ACE library '" // trim(filename) // "' does not exist!" - call fatal_error() - elseif (readable(1:3) == 'NO') then - message = "ACE library '" // trim(filename) // "' is not readable! & - &Change file permissions with chmod command." - call fatal_error() - end if - - ! display message - message = "Loading ACE cross section table: " // listing % name - call write_message(6) - - if (filetype == ASCII) then - ! ======================================================================= - ! READ ACE TABLE IN ASCII FORMAT - - ! Find location of table - open(UNIT=in, FILE=filename, STATUS='old', ACTION='read') - rewind(UNIT=in) - do i = 1, location - 1 - read(UNIT=in, FMT=*) - end do - - ! Read first line of header - read(UNIT=in, FMT='(A10,2G12.0,1X,A10)') name, awr, kT, date_ - - ! Check that correct xs was found -- if cross_sections.xml is broken, the - ! location of the table may be wrong - if(adjustl(name) /= adjustl(listing % name)) then - message = "XS listing entry " // trim(listing % name) // " did not & - &match ACE data, " // trim(name) // " found instead." - call fatal_error() - end if - - ! Read more header and NXS and JXS - read(UNIT=in, FMT=100) comment, mat, & - (zaids(i), awrs(i), i=1,16), NXS, JXS -100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& - ,8I9/8I9/8I9/8I9/8I9/8I9) - - ! determine table length - length = NXS(1) - allocate(XSS(length)) - - ! Read XSS array - read(UNIT=in, FMT='(4G20.0)') XSS - - ! Close ACE file - close(UNIT=in) - - elseif (filetype == BINARY) then - ! ======================================================================= - ! READ ACE TABLE IN BINARY FORMAT - - ! Open ACE file - open(UNIT=in, FILE=filename, STATUS='old', ACTION='read', & - ACCESS='direct', RECL=record_length) - - ! Read all header information - read(UNIT=in, REC=location) name, awr, kT, date_, & - comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS - - ! determine table length - length = NXS(1) - allocate(XSS(length)) - - ! Read remaining records with XSS - do i = 1, (length + entries - 1)/entries - j1 = 1 + (i-1)*entries - j2 = min(length, j1 + entries - 1) - read(UNIT=IN, REC=location + i) (XSS(j), j=j1,j2) - end do - - ! Close ACE file - close(UNIT=in) - end if - - ! ========================================================================== - ! PARSE DATA BASED ON NXS, JXS, AND XSS ARRAYS - - select case(listing % type) - case (ACE_NEUTRON) - nuc => nuclides(i_table) - nuc % name = name - nuc % awr = awr - nuc % kT = kT - nuc % zaid = NXS(2) - - ! read all blocks - call read_esz(nuc) - call read_nu_data(nuc) - call read_reactions(nuc) - call read_angular_dist(nuc) - call read_energy_dist(nuc) - call read_unr_res(nuc) - - ! Currently subcritical fixed source calculations are not allowed. Thus, - ! if any fissionable material is found in a fixed source calculation, - ! abort the run. - if (run_mode == MODE_FIXEDSOURCE .and. nuc % fissionable) then - message = "Cannot have fissionable material in a fixed source run." - call fatal_error() - end if - - ! for fissionable nuclides, precalculate microscopic nu-fission cross - ! sections so that we don't need to call the nu_total function during - ! cross section lookups - - if (nuc % fissionable) call generate_nu_fission(nuc) - - case (ACE_THERMAL) - sab => sab_tables(i_table) - sab % name = name - sab % awr = awr - sab % kT = kT - sab % zaid = zaids(1) - - call read_thermal_data(sab) - end select - - deallocate(XSS) - if(associated(nuc)) nullify(nuc) - if(associated(sab)) nullify(sab) - - end subroutine read_0K_ace_table - !=============================================================================== ! READ_ESZ - reads through the ESZ block. This block contains the energy grid, ! total xs, absorption xs, elastic scattering xs, and heating numbers. From 4c4e875262963f9f36efbc0b8af59390b6560b0d Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 26 Nov 2013 10:18:06 -0800 Subject: [PATCH 16/45] removed extraneous res scatter variables --- src/global.F90 | 14 +------------- 1 file changed, 1 insertion(+), 13 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 2b634db40..5f32374c5 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS, ResScatterer + MaterialMacroXS use bank_header, only: Bank use cmfd_header use constants @@ -368,18 +368,6 @@ module global logical :: output_xs = .false. logical :: output_tallies = .true. - ! ============================================================================ - ! RESONANCE SCATTERING VARIABLES - - ! Is resonance scattering treated? - logical :: treat_res_scat = .false. - - ! Number of resonant scatterers - integer :: n_res_scatterers_total - - ! Main object - type(ResScatterer), allocatable, target :: res_scatterers(:) - !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & !$omp& trace, thread_id, current_work, matching_bins) From 0a01092cc855861cac5e6f28f0139d9efa6c3ff7 Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 26 Nov 2013 10:42:32 -0800 Subject: [PATCH 17/45] removed dependence on ResScatterer type --- src/ace.F90 | 13 +------------ src/cross_section.F90 | 33 +++++++++++++++------------------ src/global.F90 | 9 --------- src/input_xml.F90 | 16 +--------------- src/physics.F90 | 26 +++++++++++++++----------- 5 files changed, 32 insertions(+), 65 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 0cb97c647..005a17d17 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -35,11 +35,9 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: k ! index over S(a,b) tables in material - integer :: l ! index over res_scatterers integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables - integer :: i_res_scat = 0 ! index in res_scatterers integer :: m ! position for sorting integer :: temp_nuclide ! temporary value for sorting integer :: temp_table ! temporary value for sorting @@ -66,7 +64,7 @@ contains NUCLIDE_LOOP: do j = 1, mat % n_nuclides name = mat % names(j) - + print*, name if (.not. already_read % contains(name)) then i_listing = xs_listing_dict % get_key(name) i_nuclide = nuclide_dict % get_key(name) @@ -77,15 +75,6 @@ contains nuc => nuclides(i_nuclide) nuc % listing = i_listing - do l = 1, n_res_scatterers_total - if (name == res_scatterers(l) % name) then - i_res_scat = i_res_scat + 1 - nuclides(i_nuclide) % resonant = .true. - nuclides(i_nuclide) % i_0K = i_res_scat - exit - end if - end do - ! Read the ACE table into the appropriate entry on the nuclides ! array call read_ace_table(i_nuclide, i_listing) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index e05b18949..086d758af 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -479,20 +479,17 @@ contains end subroutine find_energy_index !=============================================================================== -! ELASTIC_0K_XS determines the microscopic elastic cross section for a +! ELASTIC_0K_XS determines the microscopic 0K elastic cross section for a ! nuclide of a given index in the nuclides array at the energy of the particle !=============================================================================== - function elastic_0K_xs(E, nuc_0K, i_0K) result(xs_out) + function elastic_0K_xs(E, nuc) result(xs_out) - type(ResScatterer), pointer :: nuc_0K - - integer, intent(in) :: i_0K ! index into 0K nuclides array - integer :: i_grid ! index on nuclide energy grid - - real(8) :: f ! interp factor on nuclide energy grid - real(8), intent(in) :: E ! trial energy - real(8) :: xs_out + type(Nuclide), pointer :: nuc ! target nuclide at temperature + integer :: i_grid ! index on nuclide energy grid + real(8) :: f ! interp factor on nuclide energy grid + real(8), intent(in) :: E ! trial energy + real(8) :: xs_out ! Determine index on nuclide energy grid select case (grid_method) @@ -501,31 +498,31 @@ contains ! the nuclide energy grid is as simple as looking up the pointer call find_energy_index(E) - i_grid = nuc_0K % grid_index(union_grid_index) + i_grid = nuc % grid_index(union_grid_index) case (GRID_NUCLIDE) ! If we're not using the unionized grid, we have to do a binary search on ! the nuclide energy grid in order to determine which points to ! interpolate between - if (E < nuc_0K % energy(1)) then + if (E < nuc % energy(1)) then i_grid = 1 - elseif (E > nuc_0K % energy(nuc_0K % n_grid)) then - i_grid = nuc_0K % n_grid - 1 + elseif (E > nuc % energy(nuc % n_grid)) then + i_grid = nuc % n_grid - 1 else - i_grid = binary_search(nuc_0K % energy, nuc_0K % n_grid, E) + i_grid = binary_search(nuc % energy, nuc % n_grid, E) end if end select ! check for rare case where two energy points are the same - if (nuc_0K % energy(i_grid) == nuc_0K % energy(i_grid+1)) i_grid = i_grid + 1 + if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1 ! calculate interpolation factor - f = (E - nuc_0K%energy(i_grid))/(nuc_0K%energy(i_grid+1) - nuc_0K%energy(i_grid)) + f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid)) ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc_0K % elastic(i_grid) + f * nuc_0K % elastic(i_grid+1) + xs_out = (ONE - f) * nuc % elastic(i_grid) + f * nuc % elastic(i_grid+1) end function elastic_0K_xs diff --git a/src/global.F90 b/src/global.F90 index 5f32374c5..8b8e929a6 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -406,15 +406,6 @@ contains if (allocated(xs_listings)) deallocate(xs_listings) if (allocated(micro_xs)) deallocate(micro_xs) - ! Deallocate resonance scattering data - if (allocated(res_scatterers)) then - ! First call the clear routines - do i = 1, size(res_scatterers) - call res_scatterers(i) % clear() - end do - deallocate(res_scatterers) - end if - ! Deallocate external source if (allocated(external_source % params_space)) & deallocate(external_source % params_space) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 32ea3e2f6..5cfaa4bb2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -602,21 +602,6 @@ contains #endif end if - ! set up resonance scattering treatment - call lower_case(resonance_scattering_ % method(1)) - if (resonance_scattering_ % method(1) /= '') then - treat_res_scat = .true. - n_res_scatterers_total = size(resonance_scattering_ % scatterer) - allocate(res_scatterers(n_res_scatterers_total)) - do i = 1, n_res_scatterers_total - res_scatterers(i) % name = trim(resonance_scattering_ % xs_label(i)) - res_scatterers(i) % name_0K = trim(resonance_scattering_ % xs_0K_label(i)) - res_scatterers(i) % scheme = trim(resonance_scattering_ % method(i)) - res_scatterers(i) % E_min = resonance_scattering_ % E_min(i) - res_scatterers(i) % E_max = resonance_scattering_ % E_max(i) - end do - end if - end subroutine read_settings_xml !=============================================================================== @@ -1223,6 +1208,7 @@ contains else call list_density % append(-nuc % wo) end if + end do INDIVIDUAL_NUCLIDES ! ======================================================================= diff --git a/src/physics.F90 b/src/physics.F90 index fdd4f50f8..9acdc1393 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1,8 +1,8 @@ module physics - use ace_header, only: Nuclide, Reaction, DistEnergy, ResScatterer + use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: elastic_0K_xs + use cross_section, only: elastic_0K_xs use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed @@ -638,7 +638,6 @@ contains subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) type(Nuclide), pointer :: nuc ! target nuclide at temperature - type(ResScatterer), pointer :: nuc_0K => null() ! (resonant) target nuclide at 0K real(8), intent(out) :: v_target(3) real(8), intent(in) :: v_neut(3) @@ -664,7 +663,7 @@ contains logical :: reject ! resample if true - nuc_0K => res_scatterers(nuc % i_0K) + character(80) :: sampling_scheme ! Determine equilibrium temperature in MeV kT = nuc % kT @@ -672,25 +671,30 @@ contains ! Check if nuclide is a resonant scatterer and which sampling scheme ! to use based on neutron energy if (nuc % resonant) then - if (E > nuc_0K % E_max) then + print*, 'aaaaaaaaaa' + sampling_scheme = nuc % scheme + sampling_scheme = trim(sampling_scheme) + if (E > nuc % E_max) then v_target = ZERO return - else if (E < nuc_0K % E_min) then - nuc_0K % scheme = 'cxs' + else if (E < nuc % E_min) then + sampling_scheme = 'cxs' + sampling_scheme = trim(sampling_scheme) end if else if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then v_target = ZERO return else - nuc_0K % scheme = 'cxs' + sampling_scheme = 'cxs' + sampling_scheme = trim(sampling_scheme) end if end if ! reject unless criteria are satisfied reject = .true. - select case (nuc_0K % scheme) + select case (sampling_scheme) case ('cxs') ! calculate beta, alpha beta_vn = sqrt(nuc%awr * E / kT) @@ -784,10 +788,10 @@ contains ! determine velocity vector of target nucleus based on neutron's velocity ! and the sampled angle between them v_target = vt * rotate_angle(uvw, mu) - + ! adjust particle weight E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = elastic_0K_xs(E_rel, nuc_0K, nuc % i_0K) + xs_0K = elastic_0K_xs(E_rel, nuc) wcf = xs_0K / xs_eff wgt = wcf * wgt From 2c5d1ac4500fd901027bc9a4727d6f71cba32b46 Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 27 Nov 2013 00:18:39 -0800 Subject: [PATCH 18/45] wcm w/ arbitrary number of nuclides tentatively working --- src/ace.F90 | 116 +++++++++++++++++++++++++---------- src/ace_header.F90 | 33 ++++++++-- src/cross_section.F90 | 61 +++++++++--------- src/energy_grid.F90 | 22 ++++++- src/global.F90 | 15 ++++- src/input_xml.F90 | 17 ++++- src/physics.F90 | 5 +- src/templates/settings_t.xml | 12 ++-- 8 files changed, 201 insertions(+), 80 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 005a17d17..02eafa691 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -35,6 +35,7 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material integer :: k ! index over S(a,b) tables in material + integer :: n ! index over resonant scatterers integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables @@ -64,7 +65,7 @@ contains NUCLIDE_LOOP: do j = 1, mat % n_nuclides name = mat % names(j) - print*, name + if (.not. already_read % contains(name)) then i_listing = xs_listing_dict % get_key(name) i_nuclide = nuclide_dict % get_key(name) @@ -79,6 +80,25 @@ contains ! array call read_ace_table(i_nuclide, i_listing) + if (treat_res_scat) then + do n = 1, n_res_scatterers_total + if (name == nuclides_0K(n) % name) then + nuclides(i_nuclide) % resonant = .true. + nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K + nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % name_0K) + nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme + nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % scheme) + nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min + nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max + if (.not. already_read % contains(nuclides(i_nuclide) % name_0K)) then + i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % name_0K) + call read_ace_table(i_nuclide, i_listing) + end if + exit + end if + end do + end if + ! Add name and alias to dictionary call already_read % add(name) call already_read % add(alias) @@ -206,6 +226,7 @@ contains type(Nuclide), pointer :: nuc => null() type(SAlphaBeta), pointer :: sab => null() type(XsListing), pointer :: listing => null() + type(SetChar) :: already_read ! determine path, record length, and location of table listing => xs_listings(i_listing) @@ -301,18 +322,22 @@ contains select case(listing % type) case (ACE_NEUTRON) nuc => nuclides(i_table) - nuc % name = name - nuc % awr = awr - nuc % kT = kT - nuc % zaid = NXS(2) + if (trim(adjustl(name)) /= nuc % name_0K) then + nuc % name = name + nuc % awr = awr + nuc % kT = kT + nuc % zaid = NXS(2) + end if ! read all blocks call read_esz(nuc) - call read_nu_data(nuc) - call read_reactions(nuc) - call read_angular_dist(nuc) - call read_energy_dist(nuc) - call read_unr_res(nuc) + if (.not. allocated(nuc % energy_0K)) then + call read_nu_data(nuc) + call read_reactions(nuc) + call read_angular_dist(nuc) + call read_energy_dist(nuc) + call read_unr_res(nuc) + end if ! Currently subcritical fixed source calculations are not allowed. Thus, ! if any fissionable material is found in a fixed source calculation, @@ -326,7 +351,9 @@ contains ! sections so that we don't need to call the nu_total function during ! cross section lookups - if (nuc % fissionable) call generate_nu_fission(nuc) + if (nuc % fissionable .and. .not. allocated(nuc % energy_0K)) then + call generate_nu_fission(nuc) + end if case (ACE_THERMAL) sab => sab_tables(i_table) @@ -357,35 +384,58 @@ contains ! determine number of energy points NE = NXS(3) - nuc % n_grid = NE ! allocate storage for energy grid and cross section arrays - allocate(nuc % energy(NE)) - allocate(nuc % total(NE)) - allocate(nuc % elastic(NE)) - allocate(nuc % fission(NE)) - allocate(nuc % nu_fission(NE)) - allocate(nuc % absorption(NE)) + if (allocated(nuc % energy)) then + nuc % n_grid_0K = NE + allocate(nuc % energy_0K(NE)) + allocate(nuc % elastic_0K(NE)) - ! initialize cross sections - nuc % total = ZERO - nuc % elastic = ZERO - nuc % fission = ZERO - nuc % nu_fission = ZERO - nuc % absorption = ZERO + nuc % elastic_0K = ZERO - ! Read data from XSS -- only the energy grid, elastic scattering and heating - ! cross section values are actually read from here. The total and absorption - ! cross sections are reconstructed from the partial reaction data. + ! Read data from XSS -- only the energy grid, elastic scattering and heating + ! cross section values are actually read from here. The total and absorption + ! cross sections are reconstructed from the partial reaction data. - XSS_index = 1 - nuc % energy = get_real(NE) + XSS_index = 1 + nuc % energy_0K = get_real(NE) - ! Skip total and absorption - XSS_index = XSS_index + 2*NE + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic_0K = get_real(NE) - ! Continue reading elastic scattering and heating - nuc % elastic = get_real(NE) + else + nuc % n_grid = NE + allocate(nuc % energy(NE)) + allocate(nuc % total(NE)) + allocate(nuc % elastic(NE)) + allocate(nuc % fission(NE)) + allocate(nuc % nu_fission(NE)) + allocate(nuc % absorption(NE)) + + ! initialize cross sections + nuc % total = ZERO + nuc % elastic = ZERO + nuc % fission = ZERO + nuc % nu_fission = ZERO + nuc % absorption = ZERO + + ! Read data from XSS -- only the energy grid, elastic scattering and heating + ! cross section values are actually read from here. The total and absorption + ! cross sections are reconstructed from the partial reaction data. + + XSS_index = 1 + nuc % energy = get_real(NE) + + ! Skip total and absorption + XSS_index = XSS_index + 2*NE + + ! Continue reading elastic scattering and heating + nuc % elastic = get_real(NE) + + end if end subroutine read_esz diff --git a/src/ace_header.F90 b/src/ace_header.F90 index b88003207..8e41e77c4 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -112,8 +112,11 @@ module ace_header logical :: resonant = .false. ! resonant scatterer? character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme - real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section + integer :: n_grid_0K + real(8), allocatable :: grid_index_0K(:) ! pointers to union grid + real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section + real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section real(8) :: E_min ! lower cutoff energy for res scattering real(8) :: E_max ! upper cutoff energy for res scattering @@ -151,7 +154,23 @@ module ace_header contains procedure :: clear => nuclide_clear ! Deallocates Nuclide end type Nuclide - + +!=============================================================================== +! NUCLIDE_0K contains all 0K cross section data and other parameters needed to +! treat resonance scattering +!=============================================================================== + + type Nuclide0K + + character(10) :: nuclide ! name of nuclide, e.g. U-238 + character(16) :: scheme ! target velocity sampling scheme + character(10) :: name ! name of nuclide, e.g. 92235.03c + character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c + real(8) :: E_min ! lower cutoff energy for res scattering + real(8) :: E_max ! upper cutoff energy for res scattering + + end type Nuclide0K + !=============================================================================== ! SALPHABETA contains S(a,b) data for thermal neutron scattering, typically off ! of light isotopes such as water, graphite, Be, etc @@ -328,10 +347,16 @@ module ace_header if (allocated(this % grid_index)) & deallocate(this % grid_index) + if (allocated(this % grid_index_0K)) & + deallocate(this % grid_index_0K) + if (allocated(this % energy)) & deallocate(this % total, this % elastic, this % fission, & - this % nu_fission, this % absorption, & - this % elastic_0K, this % xs_cdf) + this % nu_fission, this % absorption) + + if (allocated(this % energy_0K)) & + deallocate(this % elastic_0K, this % xs_cdf) + if (allocated(this % heating)) & deallocate(this % heating) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 086d758af..7a30330a1 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -468,6 +468,7 @@ contains ! if particle's energy is outside of energy grid range, set to first or last ! index. Otherwise, do a binary search through the union energy grid. + if (E < e_grid(1)) then union_grid_index = 1 elseif (E > e_grid(n_grid)) then @@ -488,42 +489,42 @@ contains type(Nuclide), pointer :: nuc ! target nuclide at temperature integer :: i_grid ! index on nuclide energy grid real(8) :: f ! interp factor on nuclide energy grid - real(8), intent(in) :: E ! trial energy + real(8), intent(inout) :: E ! trial energy real(8) :: xs_out ! Determine index on nuclide energy grid select case (grid_method) - case (GRID_UNION) - ! If we're using the unionized grid with pointers, finding the index on - ! the nuclide energy grid is as simple as looking up the pointer - - call find_energy_index(E) - i_grid = nuc % grid_index(union_grid_index) - - case (GRID_NUCLIDE) - ! If we're not using the unionized grid, we have to do a binary search on - ! the nuclide energy grid in order to determine which points to - ! interpolate between - - if (E < nuc % energy(1)) then - i_grid = 1 - elseif (E > nuc % energy(nuc % n_grid)) then - i_grid = nuc % n_grid - 1 - else - i_grid = binary_search(nuc % energy, nuc % n_grid, E) - end if - - end select - - ! check for rare case where two energy points are the same - if (nuc % energy(i_grid) == nuc % energy(i_grid+1)) i_grid = i_grid + 1 + case (GRID_UNION) + ! If we're using the unionized grid with pointers, finding the index on + ! the nuclide energy grid is as simple as looking up the pointer - ! calculate interpolation factor - f = (E - nuc%energy(i_grid))/(nuc%energy(i_grid+1) - nuc%energy(i_grid)) + call find_energy_index(E) + i_grid = nuc % grid_index_0K(union_grid_index) - ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc % elastic(i_grid) + f * nuc % elastic(i_grid+1) + case (GRID_NUCLIDE) + ! If we're not using the unionized grid, we have to do a binary search on + ! the nuclide energy grid in order to determine which points to + ! interpolate between - end function elastic_0K_xs + if (E < nuc % energy(1)) then + i_grid = 1 + elseif (E > nuc % energy(nuc % n_grid)) then + i_grid = nuc % n_grid - 1 + else + i_grid = binary_search(nuc % energy, nuc % n_grid, E) + end if + + end select + + ! check for rare case where two energy points are the same + if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) i_grid = i_grid + 1 + + ! calculate interpolation factor + f = (E - nuc%energy_0K(i_grid))/(nuc%energy_0K(i_grid+1) - nuc%energy_0K(i_grid)) + + ! Calculate microscopic nuclide elastic cross section + xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) + + end function elastic_0K_xs end module cross_section diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 66b95be0d..41f79d1ad 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -28,6 +28,7 @@ contains do i = 1, n_nuclides_total nuc => nuclides(i) call add_grid_points(list, nuc % energy) + if (nuc % resonant) call add_grid_points(list, nuc % energy_0K) end do ! Set size of unionized energy grid @@ -132,10 +133,10 @@ contains do i = 1, n_nuclides_total nuc => nuclides(i) allocate(nuc % grid_index(n_grid)) - + index_e = 1 energy = nuc % energy(index_e) - + do j = 1, n_grid union_energy = e_grid(j) if (union_energy >= energy .and. index_e < nuc % n_grid) then @@ -144,6 +145,23 @@ contains end if nuc % grid_index(j) = index_e - 1 end do + + if (nuc % resonant) then + allocate(nuc % grid_index_0K(n_grid)) + + index_e = 1 + energy = nuc % energy_0K(index_e) + + do j = 1, n_grid + union_energy = e_grid(j) + if (union_energy >= energy .and. index_e < nuc % n_grid_0K) then + index_e = index_e + 1 + energy = nuc % energy_0K(index_e) + end if + nuc % grid_index_0K(j) = index_e - 1 + end do + end if + end do end subroutine grid_pointers diff --git a/src/global.F90 b/src/global.F90 index 8b8e929a6..dac10ee8a 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -1,7 +1,7 @@ module global use ace_header, only: Nuclide, SAlphaBeta, xsListing, NuclideMicroXS, & - MaterialMacroXS + MaterialMacroXS, Nuclide0K use bank_header, only: Bank use cmfd_header use constants @@ -367,7 +367,14 @@ module global logical :: output_summary = .false. logical :: output_xs = .false. logical :: output_tallies = .true. - + + ! ============================================================================ + ! RESONANCE SCATTERING VARIABLES + + logical :: treat_res_scat = .false. + integer :: n_res_scatterers_total = 0 + type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides + !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & !$omp& trace, thread_id, current_work, matching_bins) @@ -402,6 +409,10 @@ contains deallocate(nuclides) end if + if (allocated(nuclides_0K)) then + deallocate(nuclides_0K) + end if + if (allocated(sab_tables)) deallocate(sab_tables) if (allocated(xs_listings)) deallocate(xs_listings) if (allocated(micro_xs)) deallocate(micro_xs) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 5cfaa4bb2..af4d289fe 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -602,6 +602,21 @@ contains #endif end if + ! Resonance scattering parameters + if (lbound(resonance_scattering_ % scatterer, 1) > 0) then + treat_res_scat = .true. + n_res_scatterers_total = size(resonance_scattering_ % scatterer) + allocate(nuclides_0K(n_res_scatterers_total)) + do i = 1, n_res_scatterers_total + nuclides_0K(i) % nuclide = trim(resonance_scattering_ % scatterer(i)) + nuclides_0K(i) % scheme = trim(resonance_scattering_ % method(i)) + nuclides_0K(i) % name = trim(resonance_scattering_ % xs_label(i)) + nuclides_0K(i) % name_0K = trim(resonance_scattering_ % xs_label_0K(i)) + nuclides_0K(i) % E_min = resonance_scattering_ % E_min(i) + nuclides_0K(i) % E_max = resonance_scattering_ % E_max(i) + end do + end if + end subroutine read_settings_xml !=============================================================================== @@ -1301,7 +1316,7 @@ contains else mat % nuclide(j) = nuclide_dict % get_key(name) end if - + ! Copy name and atom/weight percent mat % names(j) = name mat % atom_density(j) = list_density % get_item(j) diff --git a/src/physics.F90 b/src/physics.F90 index 9acdc1393..45660f50c 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -655,7 +655,6 @@ contains real(8) :: beta_vt ! beta * speed of target real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target -! real(8) :: E_old ! tmp storage of current energy real(8) :: E_rel ! trial relative energy real(8) :: xs_0K ! 0K xs at E_rel real(8) :: xs_eff ! effective elastic xs at temperature T @@ -671,7 +670,6 @@ contains ! Check if nuclide is a resonant scatterer and which sampling scheme ! to use based on neutron energy if (nuc % resonant) then - print*, 'aaaaaaaaaa' sampling_scheme = nuc % scheme sampling_scheme = trim(sampling_scheme) if (E > nuc % E_max) then @@ -794,6 +792,9 @@ contains xs_0K = elastic_0K_xs(E_rel, nuc) wcf = xs_0K / xs_eff wgt = wcf * wgt + if (wgt < ZERO) then +! print*, wgt, wcf, xs_0K, xs_eff + end if case default message = "Not a recognized resonance scattering treatment!" diff --git a/src/templates/settings_t.xml b/src/templates/settings_t.xml index dbceb7550..57710206b 100644 --- a/src/templates/settings_t.xml +++ b/src/templates/settings_t.xml @@ -43,12 +43,12 @@ - - - - - - + + + + + + From b3b02b8370d4ed02e382f3f1fc79307c7c58fd2e Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 27 Nov 2013 07:36:22 -0800 Subject: [PATCH 19/45] DBRC tentatively working --- src/physics.F90 | 124 ++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 119 insertions(+), 5 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index 45660f50c..012fbd2bd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,7 +2,8 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: elastic_0K_xs + use cross_section, only: elastic_0K_xs, find_energy_index, & + union_grid_index use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed @@ -641,7 +642,7 @@ contains real(8), intent(out) :: v_target(3) real(8), intent(in) :: v_neut(3) - real(8), intent(in) :: E + real(8), intent(inout) :: E real(8), intent(in) :: uvw(3) real(8), intent(inout) :: wgt @@ -659,6 +660,18 @@ contains real(8) :: xs_0K ! 0K xs at E_rel real(8) :: xs_eff ! effective elastic xs at temperature T real(8) :: wcf ! weight correction factor + real(8) :: E_red ! reduced energy + real(8) :: E_low ! lowest practical relative energy + real(8) :: E_up ! highest practical relative energy + real(8) :: E_old ! tmp storage of current energy + real(8) :: awr ! target/neutron mass ratio + integer :: i_E_low ! 0K index to lowest practical relative energy + integer :: i_E_up ! 0K index to highest practical relative energy + real(8) :: xs_max ! max 0K xs over practical relative energies + real(8) :: xs_low ! 0K xs at lowest practical relative energy + real(8) :: xs_up ! 0K xs at highest practical relative energy + real(8) :: m ! slope for interpolation + real(8) :: R_dbrc ! DBRC rejection criterion logical :: reject ! resample if true @@ -666,6 +679,7 @@ contains ! Determine equilibrium temperature in MeV kT = nuc % kT + awr = nuc % awr ! Check if nuclide is a resonant scatterer and which sampling scheme ! to use based on neutron energy @@ -792,9 +806,109 @@ contains xs_0K = elastic_0K_xs(E_rel, nuc) wcf = xs_0K / xs_eff wgt = wcf * wgt - if (wgt < ZERO) then -! print*, wgt, wcf, xs_0K, xs_eff - end if + + case ('dbrc') + ! reduced neutron energy + E_red = sqrt((awr * E) / kT) + + ! lower limit for range that max xs is determined over + E_low = (((E_red - 4.0_8)**2) * kT) / awr + + ! upper limit for range that max xs is determined over + E_up = (((E_red + 4.0_8)**2) * kT) / awr + + ! incident neutron energy + E_old = E + + ! find index and calculate 0K xs at lower energy bound + E = E_low + call find_energy_index(E) + ! nuclide energy grid index + i_E_low = nuc % grid_index(union_grid_index) + + ! find index and calculate 0K xs at upper energy bound + E = E_up + call find_energy_index(E) + ! nuclide energy grid index + i_E_up = nuc % grid_index(union_grid_index) + + ! reset particle energy to incident value + ! (it was only changed to calculate 0K xs at different energies) + E = E_old + + ! xs at lower bounding index + xs_low = nuc % elastic_0K(i_E_low) + + ! slope + m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + + ! actual lower bound 0K xs + xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) + + ! xs at upper bounding index + xs_up = nuc % elastic_0K(i_E_up) + + ! slope + m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + + ! actual upper bound 0K xs + xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) + + ! get max 0K xs value + xs_max = max(xs_low, maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) + + ! calculate beta + beta_vn = sqrt(awr * E / kT) + alpha = ONE / (ONE + sqrt(pi) * beta_vn / TWO) + + do + ! Sample two random numbers + r1 = prn() + r2 = prn() + + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler + beta_vt_sq = -log(r1 * r2) + else + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler + c = cos(PI / TWO * prn()) + beta_vt_sq = -log(r1) - log(r2) * c**2 + end if + + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) + + ! Sample cosine of angle between neutron and target velocity + mu = TWO * prn() - ONE + + ! Determine rejection probability + accept_prob = sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) & + / (beta_vn + beta_vt) + + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) then + vt = sqrt(beta_vt_sq * kT / awr) + v_target = vt * rotate_angle(uvw, mu) + + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + + xs_0K = elastic_0K_xs(E_rel, nuc) + + R_dbrc = xs_0K / xs_max + + if (prn() < R_dbrc) then + reject = .false. + end if + end if + + if (.not. reject) exit + end do case default message = "Not a recognized resonance scattering treatment!" From 9b03eedc90d9d4b0147866b0fda0949f6626376b Mon Sep 17 00:00:00 2001 From: walshjon Date: Sun, 1 Dec 2013 18:03:11 -0800 Subject: [PATCH 20/45] fixed calculation of 0K xs --- src/ace_header.F90 | 5 ++++- src/physics.F90 | 9 ++++----- 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 8e41e77c4..a19795d8e 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -355,7 +355,10 @@ module ace_header this % nu_fission, this % absorption) if (allocated(this % energy_0K)) & - deallocate(this % elastic_0K, this % xs_cdf) + deallocate(this % elastic_0K) + + if (allocated(this % xs_cdf)) & + deallocate(this % xs_cdf) if (allocated(this % heating)) & deallocate(this % heating) diff --git a/src/physics.F90 b/src/physics.F90 index 012fbd2bd..2c1640202 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -824,13 +824,13 @@ contains E = E_low call find_energy_index(E) ! nuclide energy grid index - i_E_low = nuc % grid_index(union_grid_index) + i_E_low = nuc % grid_index_0K(union_grid_index) ! find index and calculate 0K xs at upper energy bound E = E_up call find_energy_index(E) ! nuclide energy grid index - i_E_up = nuc % grid_index(union_grid_index) + i_E_up = nuc % grid_index_0K(union_grid_index) ! reset particle energy to incident value ! (it was only changed to calculate 0K xs at different energies) @@ -894,6 +894,7 @@ contains ! Perform rejection sampling on vt and mu if (prn() < accept_prob) then vt = sqrt(beta_vt_sq * kT / awr) + v_target = vt * rotate_angle(uvw, mu) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) @@ -902,9 +903,7 @@ contains R_dbrc = xs_0K / xs_max - if (prn() < R_dbrc) then - reject = .false. - end if + if (prn() < R_dbrc) reject = .false. end if if (.not. reject) exit From 3c6781facf213f4d09d0a870e6941237d5932b20 Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 6 Dec 2013 19:09:53 -0800 Subject: [PATCH 21/45] new free gas subroutine to remove duplicate code in DBRC, CXS, WCM --- src/Makefile | 4 +- src/ace_header.F90 | 2 +- src/cross_section.F90 | 9 +- src/energy_grid.F90 | 4 +- src/global.F90 | 6 +- src/input_xml.F90 | 3 +- src/physics.F90 | 451 ++++++++++++++++++------------------------ 7 files changed, 202 insertions(+), 277 deletions(-) diff --git a/src/Makefile b/src/Makefile index 94da0f528..03c370e07 100644 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,7 @@ COMPILER = gnu DEBUG = no PROFILE = no OPTIMIZE = no -MPI = yes +MPI = no OPENMP = no HDF5 = no PETSC = no @@ -22,7 +22,7 @@ PETSC = no # External Library Paths #=============================================================================== -MPI_DIR = /opt/mpich/3.0.2-$(COMPILER) +MPI_DIR = /opt/mpich/3.0.4-$(COMPILER) HDF5_DIR = /opt/hdf5/1.8.11-$(COMPILER) PHDF5_DIR = /opt/phdf5/1.8.11-$(COMPILER) PETSC_DIR = /opt/petsc/3.4.2-$(COMPILER) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index a19795d8e..de9d115b3 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -112,7 +112,7 @@ module ace_header logical :: resonant = .false. ! resonant scatterer? character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme - integer :: n_grid_0K + integer :: n_grid_0K ! number of 0K energy grid points real(8), allocatable :: grid_index_0K(:) ! pointers to union grid real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 7a30330a1..a8be4f5c8 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -468,7 +468,6 @@ contains ! if particle's energy is outside of energy grid range, set to first or last ! index. Otherwise, do a binary search through the union energy grid. - if (E < e_grid(1)) then union_grid_index = 1 elseif (E > e_grid(n_grid)) then @@ -480,11 +479,11 @@ contains end subroutine find_energy_index !=============================================================================== -! ELASTIC_0K_XS determines the microscopic 0K elastic cross section for a -! nuclide of a given index in the nuclides array at the energy of the particle +! CALCULATE_0K_ELASTIC_XS determines the microscopic 0K elastic cross section +! for a given nuclide at the trial relative energy used in resonance scattering !=============================================================================== - function elastic_0K_xs(E, nuc) result(xs_out) + function calculate_0K_elastic_xs(E, nuc) result(xs_out) type(Nuclide), pointer :: nuc ! target nuclide at temperature integer :: i_grid ! index on nuclide energy grid @@ -525,6 +524,6 @@ contains ! Calculate microscopic nuclide elastic cross section xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) - end function elastic_0K_xs + end function calculate_0K_elastic_xs end module cross_section diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 41f79d1ad..097b74ff5 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -133,10 +133,10 @@ contains do i = 1, n_nuclides_total nuc => nuclides(i) allocate(nuc % grid_index(n_grid)) - + index_e = 1 energy = nuc % energy(index_e) - + do j = 1, n_grid union_energy = e_grid(j) if (union_energy >= energy .and. index_e < nuc % n_grid) then diff --git a/src/global.F90 b/src/global.F90 index dac10ee8a..9f646ba04 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -68,9 +68,9 @@ module global type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide type(MaterialMacroXS) :: material_xs ! Cache for current material - integer :: n_nuclides_total ! Number of nuclide cross section tables - integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables - integer :: n_listings ! Number of listings in cross_sections.xml + integer :: n_nuclides_total ! Number of nuclide cross section tables + integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables + integer :: n_listings ! Number of listings in cross_sections.xml ! Dictionaries to look up cross sections and listings type(DictCharInt) :: nuclide_dict diff --git a/src/input_xml.F90 b/src/input_xml.F90 index cc55a242b..95c6ab5f8 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1457,7 +1457,6 @@ contains call get_node_value(node_nuc, "wo", temp_dble) call list_density % append(-temp_dble) end if - end do INDIVIDUAL_NUCLIDES ! ======================================================================= @@ -1557,7 +1556,7 @@ contains else mat % nuclide(j) = nuclide_dict % get_key(name) end if - + ! Copy name and atom/weight percent mat % names(j) = name mat % atom_density(j) = list_density % get_item(j) diff --git a/src/physics.F90 b/src/physics.F90 index 2c1640202..722678cd9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,7 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: elastic_0K_xs, find_energy_index, & + use cross_section, only: calculate_0K_elastic_xs, find_energy_index, & union_grid_index use endf, only: reaction_name use error, only: fatal_error, warning @@ -435,7 +435,8 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then - call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, micro_xs(i_nuclide) % elastic) + call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & + & micro_xs(i_nuclide) % elastic) else v_t = ZERO end if @@ -628,9 +629,8 @@ contains end subroutine sab_scatter !=============================================================================== -! SAMPLE_TARGET_VELOCITY samples the target velocity based on the free gas -! scattering formulation used by most Monte Carlo codes. Excellent documentation -! for this method can be found in FRA-TM-123. Methods for correctly accounting +! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section +! free gas model is the default method. Methods for correctly accounting ! for the energy dependence of cross sections in treating resonance elastic ! scattering such as the DBRC, WCM, and a new, accelerated scheme are also ! implemented here. @@ -638,15 +638,151 @@ contains subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - type(Nuclide), pointer :: nuc ! target nuclide at temperature + type(Nuclide), pointer :: nuc ! target nuclide at temperature - real(8), intent(out) :: v_target(3) - real(8), intent(in) :: v_neut(3) - real(8), intent(inout) :: E - real(8), intent(in) :: uvw(3) - real(8), intent(inout) :: wgt + real(8), intent(out) :: v_target(3) ! target velocity + real(8), intent(in) :: v_neut(3) ! neutron velocity + real(8), intent(in) :: E ! particle energy + real(8), intent(in) :: uvw(3) ! direction cosines + real(8), intent(inout) :: wgt ! particle weight + + real(8) :: awr ! target/neutron mass ratio + real(8) :: kT ! equilibrium temperature of target in MeV + real(8) :: E_rel ! trial relative energy + real(8) :: xs_0K ! 0K xs at E_rel + real(8) :: xs_eff ! effective elastic xs at temperature T + real(8) :: wcf ! weight correction factor + real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) + real(8) :: E_low ! lowest practical relative energy + real(8) :: E_up ! highest practical relative energy + real(8) :: xs_max ! max 0K xs over practical relative energies + real(8) :: xs_low ! 0K xs at lowest practical relative energy + real(8) :: xs_up ! 0K xs at highest practical relative energy + real(8) :: m ! slope for interpolation + real(8) :: R_dbrc ! DBRC rejection criterion + + integer :: i_E_low ! 0K index to lowest practical relative energy + integer :: i_E_up ! 0K index to highest practical relative energy + + logical :: reject ! resample if true + + character(80) :: sampling_scheme ! method of target velocity sampling + + kT = nuc % kT + awr = nuc % awr + + ! Check if nuclide is a resonant scatterer + if (nuc % resonant) then + + ! sampling scheme to use + sampling_scheme = nuc % scheme + sampling_scheme = trim(sampling_scheme) + + ! upper resonance scattering energy bound (target is at rest above this E) + if (E > nuc % E_max) then + v_target = ZERO + return + + ! lower resonance scattering energy bound (should be no resonances below) + else if (E < nuc % E_min) then + sampling_scheme = 'cxs' + sampling_scheme = trim(sampling_scheme) + end if + + ! Otherwise, use free gas model + else + if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then + v_target = ZERO + return + else + sampling_scheme = 'cxs' + sampling_scheme = trim(sampling_scheme) + end if + end if + + ! Use appropriate target velocity sampling method + select case (sampling_scheme) + + case ('cxs') + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + case ('wcm') + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + ! Adjust weight as prescribed by the weight correction method (wcm) + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + xs_0K = calculate_0K_elastic_xs(E_rel, nuc) + wcf = xs_0K / xs_eff + wgt = wcf * wgt + + case ('dbrc') + E_red = sqrt((awr * E) / kT) + E_low = (((E_red - 4.0_8)**2) * kT) / awr + E_up = (((E_red + 4.0_8)**2) * kT) / awr + + ! find lower and upper energy bound indices + call find_energy_index(E_low) + i_E_low = nuc % grid_index_0K(union_grid_index) + call find_energy_index(E_up) + i_E_up = nuc % grid_index_0K(union_grid_index) + + ! interpolate xs since we're not exactly at the energy indices + xs_low = nuc % elastic_0K(i_E_low) + m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) + xs_up = nuc % elastic_0K(i_E_up) + m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) + + ! get max 0K xs value over range of practical relative energies + xs_max = max(xs_low, & + & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) + + reject = .true. + + ! sample target velocities until one is accepted by the DBRC + do + + ! sample target velocity with the constant cross section (cxs) approx. + call sample_cxs_target_velocity(nuc, v_target, E, uvw) + + ! perform Doppler broadening rejection correction (dbrc) + E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) + xs_0K = calculate_0K_elastic_xs(E_rel, nuc) + R_dbrc = xs_0K / xs_max + if (prn() < R_dbrc) reject = .false. + if (.not. reject) exit + end do + + case default + message = "Not a recognized resonance scattering treatment!" + call fatal_error() + end select + + end subroutine sample_target_velocity + +!=============================================================================== +! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas +! scattering formulation, used by most Monte Carlo codes, in which cross section +! is assumed to be constant in energy. Excellent documentation for this method +! can be found in FRA-TM-123. +!=============================================================================== + + subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw) + + type(Nuclide), pointer :: nuc ! target nuclide at temperature + real(8), intent(out) :: v_target(3) + real(8), intent(in) :: E + real(8), intent(in) :: uvw(3) real(8) :: kT ! equilibrium temperature of target in MeV + real(8) :: awr ! target/neutron mass ratio real(8) :: alpha ! probability of sampling f2 over f1 real(8) :: mu ! cosine of angle between neutron and target vel real(8) :: r1, r2 ! pseudo-random numbers @@ -656,265 +792,56 @@ contains real(8) :: beta_vt ! beta * speed of target real(8) :: beta_vt_sq ! (beta * speed of target)^2 real(8) :: vt ! speed of target - real(8) :: E_rel ! trial relative energy - real(8) :: xs_0K ! 0K xs at E_rel - real(8) :: xs_eff ! effective elastic xs at temperature T - real(8) :: wcf ! weight correction factor - real(8) :: E_red ! reduced energy - real(8) :: E_low ! lowest practical relative energy - real(8) :: E_up ! highest practical relative energy - real(8) :: E_old ! tmp storage of current energy - real(8) :: awr ! target/neutron mass ratio - integer :: i_E_low ! 0K index to lowest practical relative energy - integer :: i_E_up ! 0K index to highest practical relative energy - real(8) :: xs_max ! max 0K xs over practical relative energies - real(8) :: xs_low ! 0K xs at lowest practical relative energy - real(8) :: xs_up ! 0K xs at highest practical relative energy - real(8) :: m ! slope for interpolation - real(8) :: R_dbrc ! DBRC rejection criterion - logical :: reject ! resample if true - - character(80) :: sampling_scheme - - ! Determine equilibrium temperature in MeV kT = nuc % kT awr = nuc % awr - ! Check if nuclide is a resonant scatterer and which sampling scheme - ! to use based on neutron energy - if (nuc % resonant) then - sampling_scheme = nuc % scheme - sampling_scheme = trim(sampling_scheme) - if (E > nuc % E_max) then - v_target = ZERO - return - else if (E < nuc % E_min) then - sampling_scheme = 'cxs' - sampling_scheme = trim(sampling_scheme) - end if - else - if (E >= FREE_GAS_THRESHOLD * kT .and. nuc % awr > ONE) then - v_target = ZERO - return + beta_vn = sqrt(awr * E / kT) + alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) + + do + ! Sample two random numbers + r1 = prn() + r2 = prn() + + if (prn() < alpha) then + ! With probability alpha, we sample the distribution p(y) = + ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte + ! Carlo sampler + + beta_vt_sq = -log(r1*r2) + else - sampling_scheme = 'cxs' - sampling_scheme = trim(sampling_scheme) + ! With probability 1-alpha, we sample the distribution p(y) = y^2 * + ! e^(-y^2). This can be done with sampling scheme C61 from the Monte + ! Carlo sampler + + c = cos(PI/TWO * prn()) + beta_vt_sq = -log(r1) - log(r2)*c*c end if - end if - - ! reject unless criteria are satisfied - reject = .true. + + ! Determine beta * vt + beta_vt = sqrt(beta_vt_sq) + + ! Sample cosine of angle between neutron and target velocity + mu = TWO*prn() - ONE + + ! Determine rejection probability + accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & + /(beta_vn + beta_vt) + + ! Perform rejection sampling on vt and mu + if (prn() < accept_prob) exit + end do - select case (sampling_scheme) - case ('cxs') - ! calculate beta, alpha - beta_vn = sqrt(nuc%awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) - - ! determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - - case ('wcm') - ! calculate beta, alpha - beta_vn = sqrt(nuc%awr * E / kT) - alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - - beta_vt_sq = -log(r1*r2) - - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - - c = cos(PI/TWO * prn()) - beta_vt_sq = -log(r1) - log(r2)*c*c - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO*prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) exit - end do - - ! determine speed of target nucleus - vt = sqrt(beta_vt_sq*kT/nuc % awr) - - ! determine velocity vector of target nucleus based on neutron's velocity - ! and the sampled angle between them - v_target = vt * rotate_angle(uvw, mu) - - ! adjust particle weight - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = elastic_0K_xs(E_rel, nuc) - wcf = xs_0K / xs_eff - wgt = wcf * wgt - - case ('dbrc') - ! reduced neutron energy - E_red = sqrt((awr * E) / kT) - - ! lower limit for range that max xs is determined over - E_low = (((E_red - 4.0_8)**2) * kT) / awr - - ! upper limit for range that max xs is determined over - E_up = (((E_red + 4.0_8)**2) * kT) / awr - - ! incident neutron energy - E_old = E - - ! find index and calculate 0K xs at lower energy bound - E = E_low - call find_energy_index(E) - ! nuclide energy grid index - i_E_low = nuc % grid_index_0K(union_grid_index) - - ! find index and calculate 0K xs at upper energy bound - E = E_up - call find_energy_index(E) - ! nuclide energy grid index - i_E_up = nuc % grid_index_0K(union_grid_index) - - ! reset particle energy to incident value - ! (it was only changed to calculate 0K xs at different energies) - E = E_old - - ! xs at lower bounding index - xs_low = nuc % elastic_0K(i_E_low) - - ! slope - m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & - & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - - ! actual lower bound 0K xs - xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) - - ! xs at upper bounding index - xs_up = nuc % elastic_0K(i_E_up) - - ! slope - m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & - & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - - ! actual upper bound 0K xs - xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) - - ! get max 0K xs value - xs_max = max(xs_low, maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) - - ! calculate beta - beta_vn = sqrt(awr * E / kT) - alpha = ONE / (ONE + sqrt(pi) * beta_vn / TWO) - - do - ! Sample two random numbers - r1 = prn() - r2 = prn() - - if (prn() < alpha) then - ! With probability alpha, we sample the distribution p(y) = - ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte - ! Carlo sampler - beta_vt_sq = -log(r1 * r2) - else - ! With probability 1-alpha, we sample the distribution p(y) = y^2 * - ! e^(-y^2). This can be done with sampling scheme C61 from the Monte - ! Carlo sampler - c = cos(PI / TWO * prn()) - beta_vt_sq = -log(r1) - log(r2) * c**2 - end if - - ! Determine beta * vt - beta_vt = sqrt(beta_vt_sq) - - ! Sample cosine of angle between neutron and target velocity - mu = TWO * prn() - ONE - - ! Determine rejection probability - accept_prob = sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) & - / (beta_vn + beta_vt) - - ! Perform rejection sampling on vt and mu - if (prn() < accept_prob) then - vt = sqrt(beta_vt_sq * kT / awr) - - v_target = vt * rotate_angle(uvw, mu) - - E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - - xs_0K = elastic_0K_xs(E_rel, nuc) - - R_dbrc = xs_0K / xs_max - - if (prn() < R_dbrc) reject = .false. - end if - - if (.not. reject) exit - end do - - case default - message = "Not a recognized resonance scattering treatment!" - call fatal_error() - end select + ! determine speed of target nucleus + vt = sqrt(beta_vt_sq*kT/awr) - end subroutine sample_target_velocity + ! determine velocity vector of target nucleus based on neutron's velocity + ! and the sampled angle between them + v_target = vt * rotate_angle(uvw, mu) + + end subroutine sample_cxs_target_velocity !=============================================================================== ! CREATE_FISSION_SITES determines the average total, prompt, and delayed From e59f05e3c1cd826af2cd9b3d6ede52f3aa41a85e Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 7 Dec 2013 17:31:32 -0800 Subject: [PATCH 22/45] cleaned up error handling, comments --- src/Makefile | 4 +- src/ace.F90 | 31 ++++++----- src/ace_header.F90 | 16 +++--- src/cross_section.F90 | 6 +-- src/energy_grid.F90 | 3 ++ src/global.F90 | 6 +-- src/input_xml.F90 | 112 +++++++++++++++++++++++++++++---------- src/physics.F90 | 38 ++++++------- src/relaxng/settings.rnc | 26 ++++----- 9 files changed, 154 insertions(+), 88 deletions(-) diff --git a/src/Makefile b/src/Makefile index 03c370e07..94da0f528 100644 --- a/src/Makefile +++ b/src/Makefile @@ -13,7 +13,7 @@ COMPILER = gnu DEBUG = no PROFILE = no OPTIMIZE = no -MPI = no +MPI = yes OPENMP = no HDF5 = no PETSC = no @@ -22,7 +22,7 @@ PETSC = no # External Library Paths #=============================================================================== -MPI_DIR = /opt/mpich/3.0.4-$(COMPILER) +MPI_DIR = /opt/mpich/3.0.2-$(COMPILER) HDF5_DIR = /opt/hdf5/1.8.11-$(COMPILER) PHDF5_DIR = /opt/phdf5/1.8.11-$(COMPILER) PETSC_DIR = /opt/petsc/3.4.2-$(COMPILER) diff --git a/src/ace.F90 b/src/ace.F90 index 02eafa691..894d7fcdc 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -80,18 +80,23 @@ contains ! array call read_ace_table(i_nuclide, i_listing) + ! 0K resonant scatterer information, if treating resonance scattering if (treat_res_scat) then do n = 1, n_res_scatterers_total if (name == nuclides_0K(n) % name) then nuclides(i_nuclide) % resonant = .true. nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K - nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % name_0K) + nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % & + & name_0K) nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme - nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % scheme) + nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % & + & scheme) nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max - if (.not. already_read % contains(nuclides(i_nuclide) % name_0K)) then - i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % name_0K) + if (.not. already_read % contains(nuclides(i_nuclide) % & + & name_0K)) then + i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % & + & name_0K) call read_ace_table(i_nuclide, i_listing) end if exit @@ -226,7 +231,6 @@ contains type(Nuclide), pointer :: nuc => null() type(SAlphaBeta), pointer :: sab => null() type(XsListing), pointer :: listing => null() - type(SetChar) :: already_read ! determine path, record length, and location of table listing => xs_listings(i_listing) @@ -321,6 +325,8 @@ contains select case(listing % type) case (ACE_NEUTRON) + + ! only read in a resonant scatterers info once nuc => nuclides(i_table) if (trim(adjustl(name)) /= nuc % name_0K) then nuc % name = name @@ -331,6 +337,8 @@ contains ! read all blocks call read_esz(nuc) + + ! don't read unnecessary 0K data for resonant scatterers if (.not. allocated(nuc % energy_0K)) then call read_nu_data(nuc) call read_reactions(nuc) @@ -349,7 +357,8 @@ contains ! for fissionable nuclides, precalculate microscopic nu-fission cross ! sections so that we don't need to call the nu_total function during - ! cross section lookups + ! cross section lookups (except if we're dealing w/ 0K data for resonant + ! scatterers) if (nuc % fissionable .and. .not. allocated(nuc % energy_0K)) then call generate_nu_fission(nuc) @@ -386,17 +395,13 @@ contains NE = NXS(3) ! allocate storage for energy grid and cross section arrays + + ! read in 0K data if we've already read in non-0K data if (allocated(nuc % energy)) then nuc % n_grid_0K = NE allocate(nuc % energy_0K(NE)) allocate(nuc % elastic_0K(NE)) - nuc % elastic_0K = ZERO - - ! Read data from XSS -- only the energy grid, elastic scattering and heating - ! cross section values are actually read from here. The total and absorption - ! cross sections are reconstructed from the partial reaction data. - XSS_index = 1 nuc % energy_0K = get_real(NE) @@ -406,7 +411,7 @@ contains ! Continue reading elastic scattering and heating nuc % elastic_0K = get_real(NE) - else + else ! read in non-0K data nuc % n_grid = NE allocate(nuc % energy(NE)) allocate(nuc % total(NE)) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index de9d115b3..14cc98492 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -156,18 +156,18 @@ module ace_header end type Nuclide !=============================================================================== -! NUCLIDE_0K contains all 0K cross section data and other parameters needed to -! treat resonance scattering +! NUCLIDE_0K temporarily contains all 0K cross section data and other parameters +! needed to treat resonance scattering before transferring them to NUCLIDE !=============================================================================== type Nuclide0K - character(10) :: nuclide ! name of nuclide, e.g. U-238 - character(16) :: scheme ! target velocity sampling scheme - character(10) :: name ! name of nuclide, e.g. 92235.03c - character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c - real(8) :: E_min ! lower cutoff energy for res scattering - real(8) :: E_max ! upper cutoff energy for res scattering + character(10) :: nuclide ! name of nuclide, e.g. U-238 + character(16) :: scheme = 'dbrc' ! target velocity sampling scheme + character(10) :: name ! name of nuclide, e.g. 92235.03c + character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c + real(8) :: E_min = 0.1e-6 ! lower cutoff energy for res scattering + real(8) :: E_max = 1000.0e-6 ! upper cutoff energy for res scattering end type Nuclide0K diff --git a/src/cross_section.F90 b/src/cross_section.F90 index a8be4f5c8..9b5073235 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -487,9 +487,9 @@ contains type(Nuclide), pointer :: nuc ! target nuclide at temperature integer :: i_grid ! index on nuclide energy grid - real(8) :: f ! interp factor on nuclide energy grid - real(8), intent(inout) :: E ! trial energy - real(8) :: xs_out + real(8) :: f ! interp factor on nuclide energy grid + real(8), intent(inout) :: E ! trial energy + real(8) :: xs_out ! 0K xs at trial energy ! Determine index on nuclide energy grid select case (grid_method) diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 097b74ff5..642b86b10 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -28,6 +28,8 @@ contains do i = 1, n_nuclides_total nuc => nuclides(i) call add_grid_points(list, nuc % energy) + + ! add 0K points to the grid if (nuc % resonant) call add_grid_points(list, nuc % energy_0K) end do @@ -146,6 +148,7 @@ contains nuc % grid_index(j) = index_e - 1 end do + ! set pointers for 0K energy grid to the unionized grid if (nuc % resonant) then allocate(nuc % grid_index_0K(n_grid)) diff --git a/src/global.F90 b/src/global.F90 index 9f646ba04..4da06301b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -371,9 +371,9 @@ module global ! ============================================================================ ! RESONANCE SCATTERING VARIABLES - logical :: treat_res_scat = .false. - integer :: n_res_scatterers_total = 0 - type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides + logical :: treat_res_scat = .false. ! is resonance scattering treated? + integer :: n_res_scatterers_total = 0 ! total number of resonant scatterers + type(Nuclide0K), allocatable, target :: nuclides_0K(:) ! 0K nuclides info !$omp threadprivate(micro_xs, material_xs, fission_bank, n_bank, message, & !$omp& trace, thread_id, current_work, matching_bins) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 95c6ab5f8..6fcbb1c7b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -61,17 +61,19 @@ contains character(MAX_FILE_LEN) :: env_variable character(MAX_WORD_LEN) :: type character(MAX_LINE_LEN) :: filename - type(Node), pointer :: doc => null() - type(Node), pointer :: node_mode => null() - type(Node), pointer :: node_source => null() - type(Node), pointer :: node_dist => null() - type(Node), pointer :: node_cutoff => null() - type(Node), pointer :: node_entropy => null() - type(Node), pointer :: node_ufs => null() - type(Node), pointer :: node_sp => null() - type(Node), pointer :: node_output => null() - type(Node), pointer :: node_verb => null() - type(Node), pointer :: node_res_scat => null() + type(Node), pointer :: doc => null() + type(Node), pointer :: node_mode => null() + type(Node), pointer :: node_source => null() + type(Node), pointer :: node_dist => null() + type(Node), pointer :: node_cutoff => null() + type(Node), pointer :: node_entropy => null() + type(Node), pointer :: node_ufs => null() + type(Node), pointer :: node_sp => null() + type(Node), pointer :: node_output => null() + type(Node), pointer :: node_verb => null() + type(Node), pointer :: node_res_scat => null() + type(Node), pointer :: node_scatterer => null() + type(NodeList), pointer :: node_scat_list => null() ! Display output message message = "Reading settings XML file..." @@ -712,26 +714,72 @@ contains ! Resonance scattering parameters if (check_for_node(doc, "resonance_scattering")) then - ! Get pointer to output node call get_node_ptr(doc, "resonance_scattering", node_res_scat) - ! Check if any resonant scatterers specified - if (get_arraysize_string(node_res_scat, "scatterer") >= 1) then + call get_node_list(node_res_scat, "scatterer", node_scat_list) + + ! check that a nuclide is specified + if (get_list_size(node_scat_list) >= 1) then treat_res_scat = .true. - n_res_scatterers_total = get_arraysize_string(node_res_scat, "scatterer") + n_res_scatterers_total = get_list_size(node_scat_list) + + ! store 0K info for resonant scatterers allocate(nuclides_0K(n_res_scatterers_total)) - call get_node_array(node_res_scat, "scatterer", nuclides_0K % nuclide) - call get_node_array(node_res_scat, "method", nuclides_0K % scheme) - call get_node_array(node_res_scat, "xs_label", nuclides_0K % name) - call get_node_array(node_res_scat, "xs_label_0K", nuclides_0K % name_0K) - call get_node_array(node_res_scat, "E_min", nuclides_0K % E_min) - call get_node_array(node_res_scat, "E_max", nuclides_0K % E_max) + do i = 1, n_res_scatterers_total + call get_list_item(node_scat_list, i, node_scatterer) + + ! check to make sure a nuclide is specified + if (.not. check_for_node(node_scatterer, "nuclide")) then + message = "No nuclide specified for scatterer " // trim(to_str(i)) & + // " in settings.xml file!" + call fatal_error() + end if + call get_node_value(node_scatterer, "nuclide", & + nuclides_0K(i) % nuclide) + + if (check_for_node(node_scatterer, "method")) then + call get_node_value(node_scatterer, "method", & + nuclides_0K(i) % scheme) + end if + + ! check to make sure xs name for which method is applied is given + if (.not. check_for_node(node_scatterer, "xs_label")) then + message = "Must specify the temperature dependent name of " // '' & + //"scatterer " // trim(to_str(i)) // " given in cross_sections.xml" + call fatal_error() + end if + call get_node_value(node_scatterer, "xs_label", & + nuclides_0K(i) % name) + + ! check to make sure 0K xs name for which method is applied is given + if (.not. check_for_node(node_scatterer, "xs_label_0K")) then + message = "Must specify the 0K name of " // '' & + //"scatterer "// trim(to_str(i)) // " given in cross_sections.xml" + call fatal_error() + end if + call get_node_value(node_scatterer, "xs_label_0K", & + nuclides_0K(i) % name_0K) + + if (check_for_node(node_scatterer, "E_min")) then + call get_node_value(node_scatterer, "E_min", & + nuclides_0K(i) % E_min) + end if + + if (check_for_node(node_scatterer, "E_max")) then + call get_node_value(node_scatterer, "E_max", & + nuclides_0K(i) % E_max) + end if + + nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide) + nuclides_0K(i) % scheme = trim(nuclides_0K(i) % scheme) + call lower_case(nuclides_0K(i) % scheme) + nuclides_0K(i) % name = trim(nuclides_0K(i) % name) + nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K) + end do + else + message = "No resonant scatterers are specified within the " // "" & + // "resonance_scattering element in settings.xml" + call fatal_error() end if - do i = 1, n_res_scatterers_total - nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide) - nuclides_0K(i) % scheme = trim(nuclides_0K(i) % scheme) - nuclides_0K(i) % name = trim(nuclides_0K(i) % name) - nuclides_0K(i) % name_0K = trim(nuclides_0K(i) % name_0K) - end do end if ! Close settings XML file @@ -1245,7 +1293,6 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides - integer :: k ! loop index for resonant scatterers integer :: n ! number of nuclides integer :: n_sab ! number of sab tables for a material integer :: index_list ! index in xs_listings array @@ -3002,6 +3049,15 @@ contains end if end do + ! Check that 0K nuclides are listed in the cross_sections.xml file + do i = 1, n_res_scatterers_total + if (.not. xs_listing_dict % has_key(trim(nuclides_0K(i) % name_0K))) then + message = "Could not find nuclide " // trim(nuclides_0K(i) % name_0K) // & + " in cross_sections.xml file!" + call fatal_error() + end if + end do + ! Close cross sections XML file call close_xmldoc(doc) diff --git a/src/physics.F90 b/src/physics.F90 index 722678cd9..dec8db1ad 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -638,7 +638,7 @@ contains subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff) - type(Nuclide), pointer :: nuc ! target nuclide at temperature + type(Nuclide), pointer :: nuc ! target nuclide at temperature T real(8), intent(out) :: v_target(3) ! target velocity real(8), intent(in) :: v_neut(3) ! neutron velocity @@ -671,7 +671,7 @@ contains kT = nuc % kT awr = nuc % awr - ! Check if nuclide is a resonant scatterer + ! check if nuclide is a resonant scatterer if (nuc % resonant) then ! sampling scheme to use @@ -689,7 +689,7 @@ contains sampling_scheme = trim(sampling_scheme) end if - ! Otherwise, use free gas model + ! otherwise, use free gas model else if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then v_target = ZERO @@ -700,20 +700,20 @@ contains end if end if - ! Use appropriate target velocity sampling method + ! use appropriate target velocity sampling method select case (sampling_scheme) case ('cxs') ! sample target velocity with the constant cross section (cxs) approx. call sample_cxs_target_velocity(nuc, v_target, E, uvw) - + case ('wcm') ! sample target velocity with the constant cross section (cxs) approx. call sample_cxs_target_velocity(nuc, v_target, E, uvw) - ! Adjust weight as prescribed by the weight correction method (wcm) + ! adjust weight as prescribed by the weight correction method (wcm) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) xs_0K = calculate_0K_elastic_xs(E_rel, nuc) wcf = xs_0K / xs_eff @@ -798,46 +798,46 @@ contains beta_vn = sqrt(awr * E / kT) alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO) - + do ! Sample two random numbers r1 = prn() r2 = prn() - + if (prn() < alpha) then ! With probability alpha, we sample the distribution p(y) = ! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte ! Carlo sampler - + beta_vt_sq = -log(r1*r2) - + else ! With probability 1-alpha, we sample the distribution p(y) = y^2 * ! e^(-y^2). This can be done with sampling scheme C61 from the Monte ! Carlo sampler - + c = cos(PI/TWO * prn()) beta_vt_sq = -log(r1) - log(r2)*c*c end if - + ! Determine beta * vt beta_vt = sqrt(beta_vt_sq) - + ! Sample cosine of angle between neutron and target velocity mu = TWO*prn() - ONE - + ! Determine rejection probability accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & /(beta_vn + beta_vt) - + ! Perform rejection sampling on vt and mu if (prn() < accept_prob) exit end do - - ! determine speed of target nucleus + + ! Determine speed of target nucleus vt = sqrt(beta_vt_sq*kT/awr) - - ! determine velocity vector of target nucleus based on neutron's velocity + + ! Determine velocity vector of target nucleus based on neutron's velocity ! and the sampled angle between them v_target = vt * rotate_angle(uvw, mu) diff --git a/src/relaxng/settings.rnc b/src/relaxng/settings.rnc index 609815f41..eb78483a3 100644 --- a/src/relaxng/settings.rnc +++ b/src/relaxng/settings.rnc @@ -119,17 +119,19 @@ element settings { }? & element resonance_scattering { - (element method { list { xsd:string { maxLength = "16" }+ } | - attribute method { list { xsd:string { maxLength = "16" }+ }) & - (element scatterer { list { xsd:string { maxLength = "12" }+ } } | - attribute scatterer { list { xsd:string { maxLength = "12" }+ } }) & - (element xs_label { list { xsd:string { maxLength = "12" }+ } } | - attribute xs_label { list { xsd:string { maxLength = "12" }+ } }) & - (element xs_label_0K { list { xsd:string { maxLength = "12" }+ } } | - attribute xs_label_0K { list { xsd:string { maxLength = "12" }+ } }) & - (element E_min { list { xsd:double+ } } | - attribute E_min { list { xsd:double+ } }) & - (element E_max { list { xsd:double+ } } | - attribute E_max { list { xsd:double+ } }) + element scatterer { + (element nuclide { xsd:string { maxLength = "12" } } | + attribute nuclide { xsd:string { maxLength = "12" } }) & + (element method { xsd:string { maxLength = "16" } } | + attribute method { xsd:string { maxLength = "16" } }) & + (element xs_label { xsd:string { maxLength = "12" } } | + attribute xs_label { xsd:string { maxLength = "12" } }) & + (element xs_label_0K { xsd:string { maxLength = "12" } } | + attribute xs_label_0K { xsd:string { maxLength = "12" } }) & + (element E_min { xsd:double } | + attribute E_min { xsd:double }) & + (element E_max { xsd:double } | + attribute E_max { xsd:double })? + }* }? } From 649f3b109a1bb0c0fb906ef987c1a99e7c0258ae Mon Sep 17 00:00:00 2001 From: walshjon Date: Sun, 8 Dec 2013 20:51:30 -0800 Subject: [PATCH 23/45] changed 0K grid pointers to integers --- src/ace_header.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 14cc98492..12808dd9f 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -113,7 +113,7 @@ module ace_header character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme integer :: n_grid_0K ! number of 0K energy grid points - real(8), allocatable :: grid_index_0K(:) ! pointers to union grid + integer, allocatable :: grid_index_0K(:) ! pointers to union grid real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section From 771670f87303cb392b6f40fdced2e866005ede09 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 23 Dec 2013 17:01:24 -0800 Subject: [PATCH 24/45] handle negative 0K elastic xs values by setting to ZERO --- src/ace.F90 | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/src/ace.F90 b/src/ace.F90 index e902b516d..6aec1702d 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -398,6 +398,7 @@ contains type(Nuclide), pointer :: nuc integer :: NE ! number of energy points for total and elastic cross sections + integer :: i ! index in 0K elastic xs array for this nuclide ! determine number of energy points NE = NXS(3) @@ -419,6 +420,12 @@ contains ! Continue reading elastic scattering and heating nuc % elastic_0K = get_real(NE) + ! Negative cross sections result in a CDF that is not monotonically + ! increasing. Set all negative xs values to ZERO. + do i = 1, nuc % n_grid_0K + if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO + end do + else ! read in non-0K data nuc % n_grid = NE allocate(nuc % energy(NE)) From cf739db807adc1b57e05b5e0cd649f730ead817a Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 10 Mar 2014 16:08:53 -0700 Subject: [PATCH 25/45] added accelerated sampling scheme --- src/ace.F90 | 17 ++++++-- src/physics.F90 | 105 +++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 104 insertions(+), 18 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index 6aec1702d..41596e2fe 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -400,6 +400,8 @@ contains integer :: NE ! number of energy points for total and elastic cross sections integer :: i ! index in 0K elastic xs array for this nuclide + real(8) :: xs_cdf_sum = ZERO ! xs cdf value + ! determine number of energy points NE = NXS(3) @@ -410,7 +412,9 @@ contains nuc % n_grid_0K = NE allocate(nuc % energy_0K(NE)) allocate(nuc % elastic_0K(NE)) + allocate(nuc % xs_cdf(NE)) nuc % elastic_0K = ZERO + nuc % xs_cdf = ZERO XSS_index = 1 nuc % energy_0K = get_real(NE) @@ -420,10 +424,17 @@ contains ! Continue reading elastic scattering and heating nuc % elastic_0K = get_real(NE) - ! Negative cross sections result in a CDF that is not monotonically - ! increasing. Set all negative xs values to ZERO. - do i = 1, nuc % n_grid_0K + do i = 1, nuc % n_grid_0K - 1 + + ! Negative cross sections result in a CDF that is not monotonically + ! increasing. Set all negative xs values to ZERO. if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO + + ! build xs cdf + xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) & + & + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO & + & * (nuc % energy_0K(i+1) - nuc % energy_0K(i)) + nuc % xs_cdf(i) = xs_cdf_sum end do else ! read in non-0K data diff --git a/src/physics.F90 b/src/physics.F90 index adbfa2871..e16c9d306 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -745,25 +745,36 @@ contains real(8), intent(in) :: uvw(3) ! direction cosines real(8), intent(inout) :: wgt ! particle weight - real(8) :: awr ! target/neutron mass ratio - real(8) :: kT ! equilibrium temperature of target in MeV - real(8) :: E_rel ! trial relative energy - real(8) :: xs_0K ! 0K xs at E_rel - real(8) :: xs_eff ! effective elastic xs at temperature T - real(8) :: wcf ! weight correction factor - real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) - real(8) :: E_low ! lowest practical relative energy - real(8) :: E_up ! highest practical relative energy - real(8) :: xs_max ! max 0K xs over practical relative energies - real(8) :: xs_low ! 0K xs at lowest practical relative energy - real(8) :: xs_up ! 0K xs at highest practical relative energy - real(8) :: m ! slope for interpolation - real(8) :: R_dbrc ! DBRC rejection criterion + real(8) :: awr ! target/neutron mass ratio + real(8) :: kT ! equilibrium temperature of target in MeV + real(8) :: E_rel ! trial relative energy + real(8) :: xs_0K ! 0K xs at E_rel + real(8) :: xs_eff ! effective elastic xs at temperature T + real(8) :: wcf ! weight correction factor + real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1) + real(8) :: E_low ! lowest practical relative energy + real(8) :: E_up ! highest practical relative energy + real(8) :: E_mode ! most probable Maxwellian energy + real(8) :: E_t_max ! highest practical target energy + real(8) :: E_t ! trial target energy + real(8) :: xs_max ! max 0K xs over practical relative energies + real(8) :: xs_low ! 0K xs at lowest practical relative energy + real(8) :: xs_up ! 0K xs at highest practical relative energy + real(8) :: m ! slope for interpolation + real(8) :: R_dbrc ! DBRC rejection criterion + real(8) :: R_speed ! target speed rejection criterion + real(8) :: cdf_low ! xs cdf at lowest practical relative energy + real(8) :: cdf_up ! xs cdf at highest practical relative energy + real(8) :: cdf_rel ! trial xs cdf value + real(8) :: p_mode ! probability at most probable energy + real(8) :: p_t ! probability at trial target energy + real(8) :: mu ! cosine between neutron and target velocities integer :: i_E_low ! 0K index to lowest practical relative energy integer :: i_E_up ! 0K index to highest practical relative energy + integer :: i_E_rel ! index to trial relative energy - logical :: reject ! resample if true + logical :: reject ! resample if true character(80) :: sampling_scheme ! method of target velocity sampling @@ -859,6 +870,70 @@ contains if (.not. reject) exit end do + case ('arts') + E_red = sqrt((awr * E) / kT) + E_low = (((E_red - 4.0_8)**2) * kT) / awr + E_up = (((E_red + 4.0_8)**2) * kT) / awr + + ! find lower and upper energy bound indices + call find_energy_index(E_low) + i_E_low = nuc % grid_index_0K(union_grid_index) + call find_energy_index(E_up) + i_E_up = nuc % grid_index_0K(union_grid_index) + + ! interpolate xs CDF since we're not exactly at the energy indices + m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + cdf_low = nuc % xs_cdf(i_E_low - 1) + m * (E_low - nuc % energy_0K(i_E_low)) + m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + cdf_up = nuc % xs_cdf(i_E_up - 1) + m * (E_up - nuc % energy_0K(i_E_up)) + + ! values used to sample the Maxwellian + E_mode = kT + p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) & + & * exp(-E_mode / kT) + E_t_max = 16.0_8 * E_mode + + reject = .true. + + do + + ! perform Maxwellian rejection sampling + E_t = E_t_max * prn()**2 + p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) & + & * exp(-E_t / kT) + R_speed = p_t / p_mode + + if (prn() < R_speed) then + + ! sample a relative energy using the xs cdf + cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) + i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & + & i_E_up - i_E_low + 2, cdf_rel) + E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) + m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & + & - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & + & / (nuc % energy_0K(i_E_low + i_E_rel) & + & - nuc % energy_0K(i_E_low + i_E_rel - 1)) + E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m + + ! perform rejection sampling on cosine between + ! neutron and target velocities + mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t)) + + if (abs(mu) < ONE) then + + ! set and accept target velocity + E_t = E_t / awr + v_target = sqrt(E_t) * rotate_angle(uvw, mu) + reject = .false. + end if + end if + + if (.not. reject) exit + end do + case default message = "Not a recognized resonance scattering treatment!" call fatal_error() From 9d28e6817019964279ffb81aacd40c5bb32df053 Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 19 Mar 2014 17:23:51 -0400 Subject: [PATCH 26/45] Update publications.rst --- docs/source/publications.rst | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/publications.rst b/docs/source/publications.rst index df7ca639b..70ad1a6b6 100644 --- a/docs/source/publications.rst +++ b/docs/source/publications.rst @@ -4,6 +4,10 @@ Publications ============ +- Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated sampling + of the free gas resonance elastic scattering kernel," *Ann. Nucl. Energy*, + **69**, 116--124 (2014). ``_ + - Benoit Forget, Sheng Xu, and Kord Smith, "Direct Doppler broadening in Monte Carlo simulations using the multipole representation," *Ann. Nucl. Energy*, **64**, 78--85 (2014). ``_ From 67a1e50b964039e09ef9f918e6f969a8fe9117d9 Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 19 Mar 2014 18:30:49 -0400 Subject: [PATCH 27/45] updated input.rst w/ resonance scattering info --- docs/source/usersguide/input.rst | 56 ++++++++++++++++++++++++++++++++ 1 file changed, 56 insertions(+) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 83eb2a3d8..5cc8f46b5 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -237,6 +237,62 @@ or sub-elements and can be set to either "false" or "true". *Default*: true +```` Element +---------------------- + +The ``resonance_scattering`` element can contain one or more of the following +attributes or sub-elements: + + :scatterer: + An element with attributes/sub-elements called ``nuclide``, ``method``, + ``xs_label``, ``xs_label_0K``, ``E_min``, and ``E_max``. The ``nuclide`` + attribute is the name, as given by the ``name`` attribute within the + ``nuclide`` sub-element of the ``material`` element in ``materials.xml``, + of the nuclide to which a resonance scattering treatment is to be applied. + The ``method`` attribute gives the type of resonance scattering treatment + that is to be applied to the ``nuclide``. Acceptable inputs for the + ``method`` attribute are ``ARTS``, ``CXS``, ``WCM``, and ``DBRC``. + Descriptions of each of these methods are documented here_. The + ``xs_label`` attribute gives the label for the cross section data of the + ``nuclide`` at a given temperature. The ``xs_label_0K`` gives the label + for the 0 K cross section data for the ``nuclide``. The ``E_min`` attribute + gives the minimum energy above which the ``method`` is applied. The + ``E_max`` attribute gives the maximum energy below which the ``method`` is + applied. One example would be as follows: + + .. _here: http://dx.doi.org/10.1016/j.anucene.2014.01.017 + + .. code-block:: xml + + + + U-238 + ARTS + 92238.72c + 92238.00c + 5.0e-6 + 40.0e-6 + + + Pu-239 + dbrc + 94239.72c + 94239.00c + 0.01e-6 + 210.0e-6 + + + + .. note:: The free gas, constant cross section (``cxs``) scattering model, + which has historically been used by Monte Carlo codes to sample + target velocities, is the default ``method``. Below ``E_min``, + the ``cxs`` default method is applied and above ``E_max``, the + target-at-rest (asymptotic) kernel is used. An arbitrary number of + ``scatterer`` elements may be specified, each corresponding to a + single nuclide at a single material temperature. + + *Default*: None + ```` Element ---------------------- From 436b4fd95afb5be4deeb443bba9b98a38f342c79 Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 19 Mar 2014 19:04:21 -0400 Subject: [PATCH 28/45] added reference to resonance scattering paper --- docs/source/methods/physics.rst | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index 54dc91345..55051d597 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -997,9 +997,14 @@ Carlo code, it must be done in such a way that preserves the thermally-averaged reaction rate as per equation :eq:`doppler-broaden`. The method by which most Monte Carlo codes sample the target velocity for use in -elastic scattering kinematics is outlined in detail by [Gelbard]_. The -derivation here largely follows that of Gelbard. Let us first write the reaction -rate as a function of the velocity of the target nucleus: +elastic scattering kinematics is outlined in detail by [Gelbard]_. This is the +default method for sampling target velocities in OpenMC. Alternate methods for +treating thermal target motion which reproduce the exact elastic scattering kernel +have been implemented in OpenMC and are described by [Walsh]_. + +Here, the derivation of the default method largely follows that of Gelbard. Let +us first write the reaction rate as a function of the velocity of the target +nucleus: .. math:: :label: reaction-rate @@ -1559,6 +1564,10 @@ References .. [Gelbard] Ely M. Gelbard, "Epithermal Scattering in VIM," FRA-TM-123, Argonne National Laboratory (1979). + +.. [Walsh] Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated + sampling of the free gas resonance elastic scattering kernel," *Ann. Nucl. + Energy*, **69**, 116-124 (2014). http://dx.doi.org/10.1016/j.anucene.2014.01.017 .. [Levitt] Leo B. Levitt, "The Probability Table Method for Treating Unresolved Neutron Resonances in Monte Carlo Calculations," *Nucl. Sci. Eng.*, **49**, From 10aa60dc7843d16d06d2354b59dc83dc622fe428 Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 21 Mar 2014 12:09:30 -0700 Subject: [PATCH 29/45] Revert "added reference to resonance scattering paper" This reverts commit 436b4fd95afb5be4deeb443bba9b98a38f342c79. --- docs/source/methods/physics.rst | 15 +++------------ 1 file changed, 3 insertions(+), 12 deletions(-) diff --git a/docs/source/methods/physics.rst b/docs/source/methods/physics.rst index 55051d597..54dc91345 100644 --- a/docs/source/methods/physics.rst +++ b/docs/source/methods/physics.rst @@ -997,14 +997,9 @@ Carlo code, it must be done in such a way that preserves the thermally-averaged reaction rate as per equation :eq:`doppler-broaden`. The method by which most Monte Carlo codes sample the target velocity for use in -elastic scattering kinematics is outlined in detail by [Gelbard]_. This is the -default method for sampling target velocities in OpenMC. Alternate methods for -treating thermal target motion which reproduce the exact elastic scattering kernel -have been implemented in OpenMC and are described by [Walsh]_. - -Here, the derivation of the default method largely follows that of Gelbard. Let -us first write the reaction rate as a function of the velocity of the target -nucleus: +elastic scattering kinematics is outlined in detail by [Gelbard]_. The +derivation here largely follows that of Gelbard. Let us first write the reaction +rate as a function of the velocity of the target nucleus: .. math:: :label: reaction-rate @@ -1564,10 +1559,6 @@ References .. [Gelbard] Ely M. Gelbard, "Epithermal Scattering in VIM," FRA-TM-123, Argonne National Laboratory (1979). - -.. [Walsh] Jonathan A. Walsh, Benoit Forget, and Kord S. Smith, "Accelerated - sampling of the free gas resonance elastic scattering kernel," *Ann. Nucl. - Energy*, **69**, 116-124 (2014). http://dx.doi.org/10.1016/j.anucene.2014.01.017 .. [Levitt] Leo B. Levitt, "The Probability Table Method for Treating Unresolved Neutron Resonances in Monte Carlo Calculations," *Nucl. Sci. Eng.*, **49**, From 23929d49b9586dad42ae3acc3de8e406a8e9d0f3 Mon Sep 17 00:00:00 2001 From: walshjon Date: Sat, 22 Mar 2014 15:56:45 -0400 Subject: [PATCH 30/45] add resonance scattering test case --- tests/test_resonance_scattering/geometry.xml | 8 ++++ tests/test_resonance_scattering/materials.xml | 9 ++++ tests/test_resonance_scattering/results.py | 25 +++++++++++ .../results_true.dat | 2 + tests/test_resonance_scattering/settings.xml | 27 ++++++++++++ .../test_resonance_scattering.py | 44 +++++++++++++++++++ 6 files changed, 115 insertions(+) create mode 100644 tests/test_resonance_scattering/geometry.xml create mode 100644 tests/test_resonance_scattering/materials.xml create mode 100644 tests/test_resonance_scattering/results.py create mode 100644 tests/test_resonance_scattering/results_true.dat create mode 100644 tests/test_resonance_scattering/settings.xml create mode 100644 tests/test_resonance_scattering/test_resonance_scattering.py diff --git a/tests/test_resonance_scattering/geometry.xml b/tests/test_resonance_scattering/geometry.xml new file mode 100644 index 000000000..612e46132 --- /dev/null +++ b/tests/test_resonance_scattering/geometry.xml @@ -0,0 +1,8 @@ + + + + + + + + diff --git a/tests/test_resonance_scattering/materials.xml b/tests/test_resonance_scattering/materials.xml new file mode 100644 index 000000000..f1f941f6b --- /dev/null +++ b/tests/test_resonance_scattering/materials.xml @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/test_resonance_scattering/results.py b/tests/test_resonance_scattering/results.py new file mode 100644 index 000000000..8ff10971c --- /dev/null +++ b/tests/test_resonance_scattering/results.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python + +import sys + +# import statepoint +sys.path.append('../../src/utils') +import statepoint + +# read in statepoint file +if len(sys.argv) > 1: + sp = statepoint.StatePoint(sys.argv[1]) +else: + sp = statepoint.StatePoint('statepoint.10.binary') +sp.read_results() + +# set up output string +outstr = '' + +# write out k-combined +outstr += 'k-combined:\n' +outstr += "{0:12.6E} {1:12.6E}\n".format(sp.k_combined[0], sp.k_combined[1]) + +# write results to file +with open('results_test.dat','w') as fh: + fh.write(outstr) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat new file mode 100644 index 000000000..caf9941fe --- /dev/null +++ b/tests/test_resonance_scattering/results_true.dat @@ -0,0 +1,2 @@ +k-combined: +6.712684E-02 1.123945E-03 diff --git a/tests/test_resonance_scattering/settings.xml b/tests/test_resonance_scattering/settings.xml new file mode 100644 index 000000000..2ac47de97 --- /dev/null +++ b/tests/test_resonance_scattering/settings.xml @@ -0,0 +1,27 @@ + + + + + + U-238 + cxs + 92238.72c + 92238.72c + 5.0e-6 + 40.0e-6 + + + + + 10 + 5 + 1000 + + + + + -4 -4 -4 4 4 4 + + + + diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py new file mode 100644 index 000000000..64d20b611 --- /dev/null +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import os +from subprocess import Popen, STDOUT, PIPE, call +import filecmp +from nose_mpi import NoseMPI +import glob + +pwd = os.path.dirname(__file__) + +def setup(): + os.putenv('PWD', pwd) + os.chdir(pwd) + +def test_run(): + openmc_path = pwd + '/../../src/openmc' + if int(NoseMPI.mpi_np) > 0: + proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path], + stderr=STDOUT, stdout=PIPE) + else: + proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE) + print(proc.communicate()[0]) + returncode = proc.returncode + assert returncode == 0 + +def test_created_statepoint(): + statepoint = glob.glob(pwd + '/statepoint.10.*') + assert len(statepoint) == 1 + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5') + +def test_results(): + statepoint = glob.glob(pwd + '/statepoint.10.*') + call(['python', 'results.py', statepoint[0]]) + compare = filecmp.cmp('results_test.dat', 'results_true.dat') + if not compare: + os.rename('results_test.dat', 'results_error.dat') + assert compare + +def teardown(): + output = glob.glob(pwd + '/statepoint.10.*') + output.append(pwd + '/results_test.dat') + for f in output: + if os.path.exists(f): + os.remove(f) From b5ddfa0355c81730c2be5ba1fe334fdbe6cfb41d Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 16:06:15 -0400 Subject: [PATCH 31/45] removed extraneous character trimming --- src/physics.F90 | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index bc5cf9bc3..118e1c281 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -786,7 +786,6 @@ contains ! sampling scheme to use sampling_scheme = nuc % scheme - sampling_scheme = trim(sampling_scheme) ! upper resonance scattering energy bound (target is at rest above this E) if (E > nuc % E_max) then @@ -796,7 +795,6 @@ contains ! lower resonance scattering energy bound (should be no resonances below) else if (E < nuc % E_min) then sampling_scheme = 'cxs' - sampling_scheme = trim(sampling_scheme) end if ! otherwise, use free gas model @@ -806,7 +804,6 @@ contains return else sampling_scheme = 'cxs' - sampling_scheme = trim(sampling_scheme) end if end if From ca561f4d6c67695a419fd22069e42cb07b604d4a Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 16:14:48 -0400 Subject: [PATCH 32/45] changed calculate_0K_elastic_xs to 0K_elastic_xs --- src/cross_section.F90 | 6 +++--- src/physics.F90 | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 9b5073235..82e01c879 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -479,11 +479,11 @@ contains end subroutine find_energy_index !=============================================================================== -! CALCULATE_0K_ELASTIC_XS determines the microscopic 0K elastic cross section +! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section ! for a given nuclide at the trial relative energy used in resonance scattering !=============================================================================== - function calculate_0K_elastic_xs(E, nuc) result(xs_out) + function 0K_elastic_xs(E, nuc) result(xs_out) type(Nuclide), pointer :: nuc ! target nuclide at temperature integer :: i_grid ! index on nuclide energy grid @@ -524,6 +524,6 @@ contains ! Calculate microscopic nuclide elastic cross section xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) - end function calculate_0K_elastic_xs + end function 0K_elastic_xs end module cross_section diff --git a/src/physics.F90 b/src/physics.F90 index 118e1c281..0d27f1fc7 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,7 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: calculate_0K_elastic_xs, find_energy_index, & + use cross_section, only: 0K_elastic_xs, find_energy_index, & union_grid_index use endf, only: reaction_name use error, only: fatal_error, warning @@ -822,7 +822,7 @@ contains ! adjust weight as prescribed by the weight correction method (wcm) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = calculate_0K_elastic_xs(E_rel, nuc) + xs_0K = 0K_elastic_xs(E_rel, nuc) wcf = xs_0K / xs_eff wgt = wcf * wgt @@ -861,7 +861,7 @@ contains ! perform Doppler broadening rejection correction (dbrc) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = calculate_0K_elastic_xs(E_rel, nuc) + xs_0K = 0K_elastic_xs(E_rel, nuc) R_dbrc = xs_0K / xs_max if (prn() < R_dbrc) reject = .false. if (.not. reject) exit From daaf718870d22fc7349a05c60536d6977a8c6f5d Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 16:15:41 -0400 Subject: [PATCH 33/45] deallocate 0 K energy grid --- src/ace_header.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 2f2c4c3f3..70958259f 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -376,6 +376,9 @@ module ace_header this % nu_fission, this % absorption) if (allocated(this % energy_0K)) & + deallocate(this % energy_0K) + + if (allocated(this % elastic_0K)) & deallocate(this % elastic_0K) if (allocated(this % xs_cdf)) & From 423136518f7c79c4af065e0e00848f19f493d215 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 16:23:01 -0400 Subject: [PATCH 34/45] changed default resonance scattering treatment to arts --- src/ace_header.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 70958259f..bbf1c6fb9 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -163,10 +163,10 @@ module ace_header type Nuclide0K character(10) :: nuclide ! name of nuclide, e.g. U-238 - character(16) :: scheme = 'dbrc' ! target velocity sampling scheme + character(16) :: scheme = 'arts' ! target velocity sampling scheme character(10) :: name ! name of nuclide, e.g. 92235.03c character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c - real(8) :: E_min = 0.1e-6 ! lower cutoff energy for res scattering + real(8) :: E_min = 0.01e-6 ! lower cutoff energy for res scattering real(8) :: E_max = 1000.0e-6 ! upper cutoff energy for res scattering end type Nuclide0K From a38c3e393750f4e393e3fa13f5f010cb5caa7fad Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 17:37:22 -0400 Subject: [PATCH 35/45] changed 0K_elastic_xs to elastic_xs_0K --- src/cross_section.F90 | 4 ++-- src/physics.F90 | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 82e01c879..efb38fe0e 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -483,7 +483,7 @@ contains ! for a given nuclide at the trial relative energy used in resonance scattering !=============================================================================== - function 0K_elastic_xs(E, nuc) result(xs_out) + function elastic_xs_0K(E, nuc) result(xs_out) type(Nuclide), pointer :: nuc ! target nuclide at temperature integer :: i_grid ! index on nuclide energy grid @@ -524,6 +524,6 @@ contains ! Calculate microscopic nuclide elastic cross section xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) - end function 0K_elastic_xs + end function elastic_xs_0K end module cross_section diff --git a/src/physics.F90 b/src/physics.F90 index 0d27f1fc7..4584416bd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,7 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: 0K_elastic_xs, find_energy_index, & + use cross_section, only: elastic_xs_0K, find_energy_index, & union_grid_index use endf, only: reaction_name use error, only: fatal_error, warning @@ -822,7 +822,7 @@ contains ! adjust weight as prescribed by the weight correction method (wcm) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = 0K_elastic_xs(E_rel, nuc) + xs_0K = elastic_xs_0K(E_rel, nuc) wcf = xs_0K / xs_eff wgt = wcf * wgt @@ -861,7 +861,7 @@ contains ! perform Doppler broadening rejection correction (dbrc) E_rel = dot_product((v_neut - v_target), (v_neut - v_target)) - xs_0K = 0K_elastic_xs(E_rel, nuc) + xs_0K = elastic_xs_0K(E_rel, nuc) R_dbrc = xs_0K / xs_max if (prn() < R_dbrc) reject = .false. if (.not. reject) exit From 979b0b13bc70c9c562589e5492e5506847820911 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 18:02:39 -0400 Subject: [PATCH 36/45] cleaned up reading in of 0 K data --- src/ace.F90 | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index c5afb2a21..da3631c70 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -222,6 +222,7 @@ contains real(8) :: awrs(16) ! list of atomic weight ratios (not used) real(8) :: awr ! atomic weight ratio for table logical :: file_exists ! does ACE library exist? + logical :: data_0K ! are we reading 0K data? character(7) :: readable ! is ACE library readable? character(10) :: name ! name of ACE table character(10) :: date_ ! date ACE library was processed @@ -328,7 +329,10 @@ contains ! only read in a resonant scatterers info once nuc => nuclides(i_table) - if (trim(adjustl(name)) /= nuc % name_0K) then + data_0K = .false. + if (trim(adjustl(name)) == nuc % name_0K) then + data_0K = .true. + else nuc % name = name nuc % awr = awr nuc % kT = kT @@ -336,10 +340,12 @@ contains end if ! read all blocks - call read_esz(nuc) + call read_esz(nuc, data_0K) ! don't read unnecessary 0K data for resonant scatterers - if (.not. allocated(nuc % energy_0K)) then + if (data_0K) then + continue + else call read_nu_data(nuc) call read_reactions(nuc) call read_angular_dist(nuc) @@ -360,7 +366,7 @@ contains ! cross section lookups (except if we're dealing w/ 0K data for resonant ! scatterers) - if (nuc % fissionable .and. .not. allocated(nuc % energy_0K)) then + if (nuc % fissionable .and. .not. data_0K) then call generate_nu_fission(nuc) end if @@ -393,10 +399,12 @@ contains ! total xs, absorption xs, elastic scattering xs, and heating numbers. !=============================================================================== - subroutine read_esz(nuc) + subroutine read_esz(nuc, data_0K) type(Nuclide), pointer :: nuc + logical :: data_0K ! are we reading 0K data? + integer :: NE ! number of energy points for total and elastic cross sections integer :: i ! index in 0K elastic xs array for this nuclide @@ -408,7 +416,7 @@ contains ! allocate storage for energy grid and cross section arrays ! read in 0K data if we've already read in non-0K data - if (allocated(nuc % energy)) then + if (data_0K) then nuc % n_grid_0K = NE allocate(nuc % energy_0K(NE)) allocate(nuc % elastic_0K(NE)) From 044cad1f073fd51737c3ac18aa78e7c349bc5df3 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 7 Apr 2014 18:28:28 -0400 Subject: [PATCH 37/45] reflect changes to resonance_scattering defaults --- docs/source/usersguide/input.rst | 37 +++++++++++++++++--------------- 1 file changed, 20 insertions(+), 17 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 2ba8df75b..59209b52e 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -250,15 +250,16 @@ attributes or sub-elements: ``nuclide`` sub-element of the ``material`` element in ``materials.xml``, of the nuclide to which a resonance scattering treatment is to be applied. The ``method`` attribute gives the type of resonance scattering treatment - that is to be applied to the ``nuclide``. Acceptable inputs for the - ``method`` attribute are ``ARTS``, ``CXS``, ``WCM``, and ``DBRC``. - Descriptions of each of these methods are documented here_. The - ``xs_label`` attribute gives the label for the cross section data of the - ``nuclide`` at a given temperature. The ``xs_label_0K`` gives the label - for the 0 K cross section data for the ``nuclide``. The ``E_min`` attribute - gives the minimum energy above which the ``method`` is applied. The - ``E_max`` attribute gives the maximum energy below which the ``method`` is - applied. One example would be as follows: + that is to be applied to the ``nuclide``. Acceptable inputs - none of + which are case-sensitive - for the ``method`` attribute are ``ARTS``, + ``CXS``, ``WCM``, and ``DBRC``. Descriptions of each of these methods + are documented here_. The ``xs_label`` attribute gives the label for the + cross section data of the ``nuclide`` at a given temperature. The + ``xs_label_0K`` gives the label for the 0 K cross section data for the + ``nuclide``. The ``E_min`` attribute gives the minimum energy above + which the ``method`` is applied. The ``E_max`` attribute gives the + maximum energy below which the ``method`` is applied. One example would + be as follows: .. _here: http://dx.doi.org/10.1016/j.anucene.2014.01.017 @@ -283,15 +284,17 @@ attributes or sub-elements: - .. note:: The free gas, constant cross section (``cxs``) scattering model, - which has historically been used by Monte Carlo codes to sample - target velocities, is the default ``method``. Below ``E_min``, - the ``cxs`` default method is applied and above ``E_max``, the - target-at-rest (asymptotic) kernel is used. An arbitrary number of - ``scatterer`` elements may be specified, each corresponding to a - single nuclide at a single material temperature. + .. note:: If the ``resonance_scattering`` element is not given, the free gas, + constant cross section (``cxs``) scattering model, which has + historically been used by Monte Carlo codes to sample target + velocities, is used to treat the target motion of all nuclides. If + ``resonance_scattering`` is present, the ``cxs`` method is applied + below ``E_min`` and the target-at-rest (asymptotic) kernel is used + above ``E_max``. An arbitrary number of ``scatterer`` elements may + be specified, each corresponding to a single nuclide at a single + temperature. - *Default*: None + *Defaults*: None (scatterer), ARTS (method), 0.01 eV (E_min), 1.0 keV (E_max) ```` Element ---------------------- From f8c0d77c40868ebc988f01cee209c9ea9823691b Mon Sep 17 00:00:00 2001 From: walshjon Date: Fri, 25 Apr 2014 07:16:51 -0400 Subject: [PATCH 38/45] modified resonance scattering test --- .../test_resonance_scattering.py | 49 ++++++++++++------- 1 file changed, 31 insertions(+), 18 deletions(-) diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index 64d20b611..bcc3f8a81 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -3,42 +3,55 @@ import os from subprocess import Popen, STDOUT, PIPE, call import filecmp -from nose_mpi import NoseMPI import glob +from optparse import OptionParser -pwd = os.path.dirname(__file__) - -def setup(): - os.putenv('PWD', pwd) - os.chdir(pwd) +parser = OptionParser() +parser.add_option('--mpi_exec', dest='mpi_exec', default='') +parser.add_option('--mpi_np', dest='mpi_np', default='3') +parser.add_option('--exe', dest='exe') +(opts, args) = parser.parse_args() +cwd = os.getcwd() def test_run(): - openmc_path = pwd + '/../../src/openmc' - if int(NoseMPI.mpi_np) > 0: - proc = Popen([NoseMPI.mpi_exec, '-np', NoseMPI.mpi_np, openmc_path], + if opts.mpi_exec != '': + proc = Popen([opts.mpi_exec, '-np', opts.mpi_np, opts.exe, cwd], stderr=STDOUT, stdout=PIPE) else: - proc = Popen([openmc_path], stderr=STDOUT, stdout=PIPE) + proc = Popen([opts.exe, cwd], stderr=STDOUT, stdout=PIPE) print(proc.communicate()[0]) returncode = proc.returncode - assert returncode == 0 + assert returncode == 0, 'OpenMC did not exit successfully.' def test_created_statepoint(): - statepoint = glob.glob(pwd + '/statepoint.10.*') - assert len(statepoint) == 1 - assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5') + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + assert len(statepoint) == 1, 'Either multiple or no statepoint files exist.' + assert statepoint[0].endswith('binary') or statepoint[0].endswith('h5'),\ + 'Statepoint file is not a binary or hdf5 file.' def test_results(): - statepoint = glob.glob(pwd + '/statepoint.10.*') + statepoint = glob.glob(os.path.join(cwd, 'statepoint.10.*')) call(['python', 'results.py', statepoint[0]]) compare = filecmp.cmp('results_test.dat', 'results_true.dat') if not compare: os.rename('results_test.dat', 'results_error.dat') - assert compare + assert compare, 'Results do not agree.' def teardown(): - output = glob.glob(pwd + '/statepoint.10.*') - output.append(pwd + '/results_test.dat') + output = glob.glob(os.path.join(cwd, 'statepoint.10.*')) + output.append(os.path.join(cwd, 'results_test.dat')) for f in output: if os.path.exists(f): os.remove(f) + +if __name__ == '__main__': + + # test for openmc executable + if opts.exe is None: + raise Exception('Must specify OpenMC executable from command line with --exe.') + + # run tests + test_run() + test_created_statepoint() + test_results() + teardown() From a35c705248bb7c759e638effdbf56b67ee67af12 Mon Sep 17 00:00:00 2001 From: walshjon Date: Mon, 12 May 2014 19:36:18 -0700 Subject: [PATCH 39/45] set lower cdf bound to zero if lower energy bound is below available data --- src/physics.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/physics.F90 b/src/physics.F90 index 4584416bd..baf649b66 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -882,6 +882,7 @@ contains m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) cdf_low = nuc % xs_cdf(i_E_low - 1) + m * (E_low - nuc % energy_0K(i_E_low)) + if (i_E_low == 1) cdf_low = ZERO m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) cdf_up = nuc % xs_cdf(i_E_up - 1) + m * (E_up - nuc % energy_0K(i_E_up)) From 8367a9727cc42c122ed83bb270cc82327fb7d3f7 Mon Sep 17 00:00:00 2001 From: walshjon Date: Thu, 19 Jun 2014 10:09:43 -0400 Subject: [PATCH 40/45] modified res scat test to use nndc data and at a provided temperature --- tests/test_resonance_scattering/materials.xml | 2 +- tests/test_resonance_scattering/results_true.dat | 2 +- tests/test_resonance_scattering/settings.xml | 4 ++-- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/tests/test_resonance_scattering/materials.xml b/tests/test_resonance_scattering/materials.xml index f1f941f6b..52a8c04be 100644 --- a/tests/test_resonance_scattering/materials.xml +++ b/tests/test_resonance_scattering/materials.xml @@ -3,7 +3,7 @@ - + diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index caf9941fe..2649c4414 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -6.712684E-02 1.123945E-03 +6.822440E-02 6.306868E-04 diff --git a/tests/test_resonance_scattering/settings.xml b/tests/test_resonance_scattering/settings.xml index 2ac47de97..7ce4f23ac 100644 --- a/tests/test_resonance_scattering/settings.xml +++ b/tests/test_resonance_scattering/settings.xml @@ -5,8 +5,8 @@ U-238 cxs - 92238.72c - 92238.72c + 92238.71c + 92238.71c 5.0e-6 40.0e-6 From 39e272cd63b6c6785bb8c9e130f09186770830d4 Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 9 Sep 2014 13:18:43 -0400 Subject: [PATCH 41/45] added some checks on resonance scattering energy bounds and in test problem python script --- src/input_xml.F90 | 14 +++++++++++++- .../test_resonance_scattering.py | 10 ++++++---- 2 files changed, 19 insertions(+), 5 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7132d18e6..9b8b1a3c1 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -840,12 +840,24 @@ contains call get_node_value(node_scatterer, "E_min", & nuclides_0K(i) % E_min) end if - + + ! check that E_min is non-negative + if (E_min < ZERO) then + message = "Lower resonance scattering energy bound is negative" + call fatal_error() + end if + if (check_for_node(node_scatterer, "E_max")) then call get_node_value(node_scatterer, "E_max", & nuclides_0K(i) % E_max) end if + ! check that E_max is not less than E_min + if (E_max < E_min) then + message = "Lower resonance scattering energy bound exceeds upper" + call fatal_error() + end if + nuclides_0K(i) % nuclide = trim(nuclides_0K(i) % nuclide) nuclides_0K(i) % scheme = trim(nuclides_0K(i) % scheme) call lower_case(nuclides_0K(i) % scheme) diff --git a/tests/test_resonance_scattering/test_resonance_scattering.py b/tests/test_resonance_scattering/test_resonance_scattering.py index bcc3f8a81..6fdbf8745 100644 --- a/tests/test_resonance_scattering/test_resonance_scattering.py +++ b/tests/test_resonance_scattering/test_resonance_scattering.py @@ -51,7 +51,9 @@ if __name__ == '__main__': raise Exception('Must specify OpenMC executable from command line with --exe.') # run tests - test_run() - test_created_statepoint() - test_results() - teardown() + try: + test_run() + test_created_statepoint() + test_results() + finally: + teardown() From 20cffb0c4a850ab97c8262c86734d2f566549b03 Mon Sep 17 00:00:00 2001 From: walshjon Date: Tue, 9 Sep 2014 18:47:44 -0400 Subject: [PATCH 42/45] handle case with neutron being in the lowest energy grid bin or below --- src/physics.F90 | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/src/physics.F90 b/src/physics.F90 index ba317f5fd..32d317567 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -884,13 +884,24 @@ contains i_E_up = nuc % grid_index_0K(union_grid_index) ! interpolate xs CDF since we're not exactly at the energy indices - m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & - & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) - cdf_low = nuc % xs_cdf(i_E_low - 1) + m * (E_low - nuc % energy_0K(i_E_low)) - if (i_E_low == 1) cdf_low = ZERO + ! cdf value at lower bound attainable energy + if (i_E_low > 1) then + m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + cdf_low = nuc % xs_cdf(i_E_low - 1) & + & + m * (E_low - nuc % energy_0K(i_E_low)) + else + m = nuc % xs_cdf(i_E_low) & + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + cdf_low = m * (E_low - nuc % energy_0K(i_E_low)) + if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO + end if + + ! cdf value at upper bound attainable energy m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) - cdf_up = nuc % xs_cdf(i_E_up - 1) + m * (E_up - nuc % energy_0K(i_E_up)) + cdf_up = nuc % xs_cdf(i_E_up - 1) & + & + m * (E_up - nuc % energy_0K(i_E_up)) ! values used to sample the Maxwellian E_mode = kT From 4f8634680e04ab0144679c0a89d40d4ea8afb3f8 Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 10 Sep 2014 12:03:38 -0400 Subject: [PATCH 43/45] only use nuclide grids for resonance scattering --- src/ace_header.F90 | 4 ---- src/cross_section.F90 | 39 +++++++++++++----------------------- src/energy_grid.F90 | 18 ----------------- src/input_xml.F90 | 4 ++-- src/physics.F90 | 46 +++++++++++++++++++++++++++++++++---------- 5 files changed, 52 insertions(+), 59 deletions(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index c881a559d..b0fec99b6 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -117,7 +117,6 @@ module ace_header character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c character(16) :: scheme ! target velocity sampling scheme integer :: n_grid_0K ! number of 0K energy grid points - integer, allocatable :: grid_index_0K(:) ! pointers to union grid real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section @@ -372,9 +371,6 @@ module ace_header if (allocated(this % grid_index)) & deallocate(this % grid_index) - - if (allocated(this % grid_index_0K)) & - deallocate(this % grid_index_0K) if (allocated(this % energy)) & deallocate(this % energy, this % total, this % elastic, & diff --git a/src/cross_section.F90 b/src/cross_section.F90 index 1c6ed8ae8..0620cf11f 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -527,37 +527,26 @@ contains real(8) :: xs_out ! 0K xs at trial energy ! Determine index on nuclide energy grid - select case (grid_method) - case (GRID_UNION) - ! If we're using the unionized grid with pointers, finding the index on - ! the nuclide energy grid is as simple as looking up the pointer - - call find_energy_index(E) - i_grid = nuc % grid_index_0K(union_grid_index) - - case (GRID_NUCLIDE) - ! If we're not using the unionized grid, we have to do a binary search on - ! the nuclide energy grid in order to determine which points to - ! interpolate between - - if (E < nuc % energy(1)) then - i_grid = 1 - elseif (E > nuc % energy(nuc % n_grid)) then - i_grid = nuc % n_grid - 1 - else - i_grid = binary_search(nuc % energy, nuc % n_grid, E) - end if - - end select + if (E < nuc % energy_0K(1)) then + i_grid = 1 + elseif (E > nuc % energy_0K(nuc % n_grid_0K)) then + i_grid = nuc % n_grid_0K - 1 + else + i_grid = binary_search(nuc % energy_0K, nuc % n_grid_0K, E) + end if ! check for rare case where two energy points are the same - if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) i_grid = i_grid + 1 + if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then + i_grid = i_grid + 1 + end if ! calculate interpolation factor - f = (E - nuc%energy_0K(i_grid))/(nuc%energy_0K(i_grid+1) - nuc%energy_0K(i_grid)) + f = (E - nuc % energy_0K(i_grid)) & + & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) ! Calculate microscopic nuclide elastic cross section - xs_out = (ONE - f) * nuc % elastic_0K(i_grid) + f * nuc % elastic_0K(i_grid+1) + xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & + & + f * nuc % elastic_0K(i_grid + 1) end function elastic_xs_0K diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index 642b86b10..a38a305d3 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -147,24 +147,6 @@ contains end if nuc % grid_index(j) = index_e - 1 end do - - ! set pointers for 0K energy grid to the unionized grid - if (nuc % resonant) then - allocate(nuc % grid_index_0K(n_grid)) - - index_e = 1 - energy = nuc % energy_0K(index_e) - - do j = 1, n_grid - union_energy = e_grid(j) - if (union_energy >= energy .and. index_e < nuc % n_grid_0K) then - index_e = index_e + 1 - energy = nuc % energy_0K(index_e) - end if - nuc % grid_index_0K(j) = index_e - 1 - end do - end if - end do end subroutine grid_pointers diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b8b1a3c1..27995597c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -842,7 +842,7 @@ contains end if ! check that E_min is non-negative - if (E_min < ZERO) then + if (nuclides_0K(i) % E_min < ZERO) then message = "Lower resonance scattering energy bound is negative" call fatal_error() end if @@ -853,7 +853,7 @@ contains end if ! check that E_max is not less than E_min - if (E_max < E_min) then + if (nuclides_0K(i) % E_max < nuclides_0K(i) % E_min) then message = "Lower resonance scattering energy bound exceeds upper" call fatal_error() end if diff --git a/src/physics.F90 b/src/physics.F90 index 32d317567..e0fbf80d2 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -835,13 +835,26 @@ contains E_red = sqrt((awr * E) / kT) E_low = (((E_red - 4.0_8)**2) * kT) / awr E_up = (((E_red + 4.0_8)**2) * kT) / awr - + ! find lower and upper energy bound indices - call find_energy_index(E_low) - i_E_low = nuc % grid_index_0K(union_grid_index) - call find_energy_index(E_up) - i_E_up = nuc % grid_index_0K(union_grid_index) - + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if + ! interpolate xs since we're not exactly at the energy indices xs_low = nuc % elastic_0K(i_E_low) m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & @@ -878,10 +891,23 @@ contains E_up = (((E_red + 4.0_8)**2) * kT) / awr ! find lower and upper energy bound indices - call find_energy_index(E_low) - i_E_low = nuc % grid_index_0K(union_grid_index) - call find_energy_index(E_up) - i_E_up = nuc % grid_index_0K(union_grid_index) + ! lower index + if (E_low < nuc % energy_0K(1)) then + i_E_low = 1 + elseif (E_low > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_low = nuc % n_grid_0K - 1 + else + i_E_low = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_low) + end if + + ! upper index + if (E_up < nuc % energy_0K(1)) then + i_E_up = 1 + elseif (E_up > nuc % energy_0K(nuc % n_grid_0K)) then + i_E_up = nuc % n_grid_0K - 1 + else + i_E_up = binary_search(nuc % energy_0K, nuc % n_grid_0K, E_up) + end if ! interpolate xs CDF since we're not exactly at the energy indices ! cdf value at lower bound attainable energy From 4d10704a8771d25537e000c1e6e7195a061699ca Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 10 Sep 2014 11:35:09 -0700 Subject: [PATCH 44/45] updated resonance scattering test --- tests/test_resonance_scattering/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_resonance_scattering/results_true.dat b/tests/test_resonance_scattering/results_true.dat index 2649c4414..62b0ac223 100644 --- a/tests/test_resonance_scattering/results_true.dat +++ b/tests/test_resonance_scattering/results_true.dat @@ -1,2 +1,2 @@ k-combined: -6.822440E-02 6.306868E-04 +6.452021E-02 1.738736E-03 From 3a9f2268def5f2354069f8d5ea9d78a4f76491ca Mon Sep 17 00:00:00 2001 From: walshjon Date: Wed, 10 Sep 2014 12:20:58 -0700 Subject: [PATCH 45/45] do not add 0 K energy points to union grid --- docs/source/usersguide/input.rst | 6 +++--- src/ace_header.F90 | 4 ++-- src/energy_grid.F90 | 3 --- src/physics.F90 | 5 ++--- 4 files changed, 7 insertions(+), 11 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index f50c7ff21..c22008b79 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -274,7 +274,7 @@ attributes or sub-elements: of the nuclide to which a resonance scattering treatment is to be applied. The ``method`` attribute gives the type of resonance scattering treatment that is to be applied to the ``nuclide``. Acceptable inputs - none of - which are case-sensitive - for the ``method`` attribute are ``ARTS``, + which are case-sensitive - for the ``method`` attribute are ``ARES``, ``CXS``, ``WCM``, and ``DBRC``. Descriptions of each of these methods are documented here_. The ``xs_label`` attribute gives the label for the cross section data of the ``nuclide`` at a given temperature. The @@ -291,7 +291,7 @@ attributes or sub-elements: U-238 - ARTS + ARES 92238.72c 92238.00c 5.0e-6 @@ -317,7 +317,7 @@ attributes or sub-elements: be specified, each corresponding to a single nuclide at a single temperature. - *Defaults*: None (scatterer), ARTS (method), 0.01 eV (E_min), 1.0 keV (E_max) + *Defaults*: None (scatterer), ARES (method), 0.01 eV (E_min), 1.0 keV (E_max) ```` Element ---------------------- diff --git a/src/ace_header.F90 b/src/ace_header.F90 index b0fec99b6..6d339891c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -159,14 +159,14 @@ module ace_header end type Nuclide !=============================================================================== -! NUCLIDE_0K temporarily contains all 0K cross section data and other parameters +! NUCLIDE0K temporarily contains all 0K cross section data and other parameters ! needed to treat resonance scattering before transferring them to NUCLIDE !=============================================================================== type Nuclide0K character(10) :: nuclide ! name of nuclide, e.g. U-238 - character(16) :: scheme = 'arts' ! target velocity sampling scheme + character(16) :: scheme = 'ares' ! target velocity sampling scheme character(10) :: name ! name of nuclide, e.g. 92235.03c character(10) :: name_0K ! name of 0K nuclide, e.g. 92235.00c real(8) :: E_min = 0.01e-6 ! lower cutoff energy for res scattering diff --git a/src/energy_grid.F90 b/src/energy_grid.F90 index a38a305d3..66b95be0d 100644 --- a/src/energy_grid.F90 +++ b/src/energy_grid.F90 @@ -28,9 +28,6 @@ contains do i = 1, n_nuclides_total nuc => nuclides(i) call add_grid_points(list, nuc % energy) - - ! add 0K points to the grid - if (nuc % resonant) call add_grid_points(list, nuc % energy_0K) end do ! Set size of unionized energy grid diff --git a/src/physics.F90 b/src/physics.F90 index e0fbf80d2..04cc66afa 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -2,8 +2,7 @@ module physics use ace_header, only: Nuclide, Reaction, DistEnergy use constants - use cross_section, only: elastic_xs_0K, find_energy_index, & - union_grid_index + use cross_section, only: elastic_xs_0K use endf, only: reaction_name use error, only: fatal_error, warning use fission, only: nu_total, nu_delayed @@ -885,7 +884,7 @@ contains if (.not. reject) exit end do - case ('arts') + case ('ares') E_red = sqrt((awr * E) / kT) E_low = (((E_red - 4.0_8)**2) * kT) / awr E_up = (((E_red + 4.0_8)**2) * kT) / awr