From ed56a6821563fdea90d61a37e4679d5abfff4320 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 19 Jan 2016 20:39:23 -0500 Subject: [PATCH] Improve diff tally docs and comments --- docs/source/usersguide/input.rst | 21 ++++++++++++ src/input_xml.F90 | 4 +-- src/output.F90 | 2 +- src/relaxng/tallies.rnc | 8 ++++- src/relaxng/tallies.rng | 34 ++++++++++++++++++++ src/state_point.F90 | 2 +- src/tally.F90 | 55 ++++++++++++++++++++++---------- src/tally_header.F90 | 5 +-- src/tracking.F90 | 4 +-- 9 files changed, 109 insertions(+), 26 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 30e9ed07b9..2553e856bf 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1618,6 +1618,27 @@ The ```` element accepts the following sub-elements: *Default*: "all" + :derivative: + Tally the first-order derivative of the quantity of interest with respect to + some material perturbation. Differential tallies are currently only + implemented for collision and analog estimators. + + The ``derivative`` element has the following attributes/sub-elements: + + :variable: + The independent variable of the derivative. Accepted options are + "density" and "nuclide_density". A "density" derivative will give the + derivative with respect to the density of the material in [g / cm^3]. A + "nuclide_density" derivative will give the derivative with respect to + the density of a particular nuclide in units of [atom / b / cm]. + + :material: + The perturbed material. (Necessary for both "density" and + "nuclide_density") + + :nuclide: + The perturbed nuclide. (Necessary only for "nuclide_density") + ```` Element ------------------ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4cd803cb14..6f77fd96fe 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2996,11 +2996,11 @@ contains select case(temp_str) case("density") - t % deriv % dep_var = DIFF_DENSITY + t % deriv % variable = DIFF_DENSITY call get_node_value(node_deriv, "material", t % deriv % diff_material) case("nuclide_density") - t % deriv % dep_var = DIFF_NUCLIDE_DENSITY + t % deriv % variable = DIFF_NUCLIDE_DENSITY call get_node_value(node_deriv, "material", t % deriv % diff_material) call get_node_value(node_deriv, "nuclide", word) diff --git a/src/output.F90 b/src/output.F90 index 5f1b86ef6f..fde80f94f7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1019,7 +1019,7 @@ contains ! Write derivative information. if (allocated(t % deriv)) then - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_DENSITY) write(unit=unit_tally, fmt="(' Density derivative Material ',A)") & to_str(t % deriv % diff_material) diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 0f3672c6f2..11b5feec6e 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -41,7 +41,13 @@ element tallies { (element type { xsd:string } | attribute type { xsd:string }) & (element threshold { xsd:double} | attribute threshold { xsd:double }) & (element scores { list { xsd:string { maxLength = "20" }+ } } | attribute scores { list { xsd:string { maxLength = "20"}+ } } )? - }* + }* & + element derivative { + (element variable { ( "density" | "nuclide_density") } | + attribute variable { ( "density" | "nuclide_density") } ) & + (element material { xsd:int } | attribute material { xsd:int }) & + (element nuclide { xsd:string } | attribute nuclide { xsd:string })? + }? }* & element assume_separate { xsd:boolean }? diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index 36bd5cb85f..83c532b79a 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -260,5 +260,39 @@ + + + + + + density + nuclide_density + + + + + density + nuclide_density + + + + + + + + + + + + + + + + + + + + + diff --git a/src/state_point.F90 b/src/state_point.F90 index 5f33eaeb53..93eaaa69fe 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -311,7 +311,7 @@ contains if (allocated(tally % deriv)) then call write_dataset(tally_group, "derivative present", 1) deriv_group = create_group(tally_group, "derivative") - select case (tally % deriv % dep_var) + select case (tally % deriv % variable) case (DIFF_DENSITY) call write_dataset(deriv_group, "dependent variable", "density") call write_dataset(deriv_group, "material", & diff --git a/src/tally.F90 b/src/tally.F90 index bd6716f254..ba7581a609 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -772,7 +772,7 @@ contains end if case (ESTIMATOR_COLLISION) - if (t % deriv % dep_var == DIFF_NUCLIDE_DENSITY) then + if (t % deriv % variable == DIFF_NUCLIDE_DENSITY) then scoring_diff_nuclide = & (materials(p % material) % id == t % deriv % diff_material) & .and. (i_nuclide == t % deriv % diff_nuclide) @@ -801,7 +801,7 @@ contains score = score * t % deriv % flux_deriv case (SCORE_TOTAL) - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_NUCLIDE_DENSITY) if (i_nuclide == -1 .and. & @@ -819,7 +819,7 @@ contains end select case (SCORE_ABSORPTION) - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_NUCLIDE_DENSITY) if (i_nuclide == -1 .and. & @@ -837,7 +837,7 @@ contains end select case (SCORE_FISSION) - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_NUCLIDE_DENSITY) if (i_nuclide == -1 .and. & @@ -855,7 +855,7 @@ contains end select case (SCORE_NU_FISSION) - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_NUCLIDE_DENSITY) if (i_nuclide == -1 .and. & @@ -873,7 +873,7 @@ contains end select case (SCORE_KEFF) - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_DENSITY) score = score * t % deriv % flux_deriv @@ -2383,12 +2383,13 @@ contains end subroutine score_surface_current !=============================================================================== -! SCORE_TRACK_DERIVATIVE +! SCORE_TRACK_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron travelling through a perturbed material. !=============================================================================== subroutine score_track_derivative(p, distance) type(particle), intent(in) :: p - real(8), intent(in) :: distance + real(8), intent(in) :: distance ! Neutron flight distance integer :: i @@ -2396,13 +2397,16 @@ contains do i = 1, active_tallies % size() associate (t => tallies(active_tallies % get_item(i))) - if (.not. allocated(t % deriv)) cycle + if (.not. allocated(t % deriv)) cycle ! Ignore non-differential tallies - select case (t % deriv % dep_var) + select case (t % deriv % variable) case (DIFF_DENSITY) associate (mat => materials(p % material)) if (mat % id == t % deriv % diff_material) then + ! phi = e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist + ! (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist t % deriv % flux_deriv = t % deriv % flux_deriv & - distance * material_xs % total / mat % density_gpcc end if @@ -2411,6 +2415,9 @@ contains case (DIFF_NUCLIDE_DENSITY) associate (mat => materials(p % material)) if (mat % id == t % deriv % diff_material) then + ! phi = e^(-Sigma_tot * dist) + ! (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist + ! (1 / phi) * (d_phi / d_N) = - sigma_tot * dist t % deriv % flux_deriv = t % deriv % flux_deriv & - distance * micro_xs(t % deriv % diff_nuclide) % total end if @@ -2418,10 +2425,11 @@ contains end select end associate end do - end subroutine + end subroutine score_track_derivative !=============================================================================== -! SCORE_COLLISION_DERIVATIVE +! SCORE_COLLISION_DERIVATIVE Adjust flux derivatives on differential tallies to +! account for a neutron colliding in the perturbed material. !=============================================================================== subroutine score_collision_derivative(p) @@ -2433,12 +2441,15 @@ contains do i = 1, active_tallies % size() associate (t => tallies(active_tallies % get_item(i))) - if (.not. allocated(t % deriv)) cycle - select case (t % deriv % dep_var) + if (.not. allocated(t % deriv)) cycle ! Ignore non-differential tallies + select case (t % deriv % variable) case (DIFF_DENSITY) associate (mat => materials(p % material)) if (mat % id == t % deriv % diff_material) then + ! phi = Sigma_MT + ! (1 / phi) * (d_phi / d_rho) = (d_Sigma_MT / d_rho) / Sigma_MT + ! (1 / phi) * (d_phi / d_rho) = 1 / rho t % deriv % flux_deriv = t % deriv % flux_deriv & + ONE / mat % density_gpcc end if @@ -2448,12 +2459,18 @@ contains associate (mat => materials(p % material)) if (mat % id == t % deriv % diff_material & .and. p % event_nuclide == t % deriv % diff_nuclide) then + ! Find the index in this material for the diff_nuclide. do j = 1, mat % n_nuclides if (mat % nuclide(j) == t % deriv % diff_nuclide) exit end do + ! Make sure we found the nuclide. if (mat % nuclide(j) /= t % deriv % diff_nuclide) then call fatal_error("Couldn't find the right nuclide.") end if + ! phi = Sigma_MT + ! (1 / phi) * (d_phi / d_N) = (d_Sigma_MT / d_N) / Sigma_MT + ! (1 / phi) * (d_phi / d_N) = sigma_MT / Sigma_MT + ! (1 / phi) * (d_phi / d_N) = 1 / N t % deriv % flux_deriv = t % deriv % flux_deriv & + ONE / mat % atom_density(j) end if @@ -2461,16 +2478,20 @@ contains end select end associate end do - end subroutine + end subroutine score_collision_derivative - subroutine clear_diff_flux_derivs() +!=============================================================================== +! ZERO_FLUX_DERIVS Set the flux derivatives on differential tallies to zero. +!=============================================================================== + + subroutine zero_flux_derivs() integer :: i do i = 1, n_tallies if (allocated(tallies(i) % deriv)) then tallies(i) % deriv % flux_deriv = ZERO end if end do - end subroutine + end subroutine zero_flux_derivs !=============================================================================== ! GET_NEXT_BIN determines the next scoring bin for a particular filter variable diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 9acd05a5ec..a7b455d745 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -61,12 +61,13 @@ module tally_header !=============================================================================== -! TALLYDERIVATIVE +! TALLYDERIVATIVE describes a first-order derivative that can be applied to +! tallies. !=============================================================================== type TallyDerivative real(8) :: flux_deriv - integer :: dep_var + integer :: variable integer :: diff_material integer :: diff_nuclide end type TallyDerivative diff --git a/src/tracking.F90 b/src/tracking.F90 index 1ad379490c..587b678095 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -15,7 +15,7 @@ module tracking use tally, only: score_analog_tally, score_tracklength_tally, & score_collision_tally, score_surface_current, & score_track_derivative, & - score_collision_derivative, clear_diff_flux_derivs + score_collision_derivative, zero_flux_derivs use track_output, only: initialize_particle_track, write_particle_track, & add_particle_track, finalize_particle_track @@ -63,7 +63,7 @@ contains endif ! Every particle starts with no accumulated flux derivative. - call clear_diff_flux_derivs() + call zero_flux_derivs() EVENT_LOOP: do ! If the cell hasn't been determined based on the particle's location,