From 349c2a804bb43239f9806653442fa58ed7267247 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 20 Sep 2018 14:00:32 -0400 Subject: [PATCH 01/38] Add CMFD helper functions to access entropy_on, is_master, current_batch, implement vector_write --- openmc/capi/core.py | 20 ++++++++++++++++++++ openmc/capi/settings.py | 4 +++- src/cmfd_data.F90 | 2 +- src/cmfd_solver.F90 | 11 ++++++++--- src/tallies/tally_header.F90 | 2 ++ src/vector_header.F90 | 17 ++++++++++++++++- 6 files changed, 50 insertions(+), 6 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 53c3c32bd6..7e4d29afc3 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -62,6 +62,16 @@ def calculate_volumes(): _dll.openmc_calculate_volumes() +def current_batch(): + """Return the current batch of the simulation. + Returns + ------- + int + Current batch of the simulation + """ + return c_int.in_dll(_dll, 'openmc_current_batch').value + + def finalize(): """Finalize simulation and free memory""" _dll.openmc_finalize() @@ -214,6 +224,16 @@ def keff(): return (mean, std_dev) +def master(): + """Return whether processor is master processor or not. + Returns + ------- + bool + Whether is master processor or not + """ + return c_bool.in_dll(_dll, 'openmc_master').value + + def next_batch(): """Run next batch. diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index 1063d6463e..af6dd7d7a2 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -1,4 +1,5 @@ -from ctypes import c_int, c_int32, c_int64, c_double, c_char_p, POINTER +from ctypes import (c_int, c_int32, c_int64, c_double, c_char_p, c_bool, + POINTER) from . import _dll from .core import _DLLGlobal @@ -17,6 +18,7 @@ _dll.openmc_get_seed.restype = c_int64 class _Settings(object): # Attributes that are accessed through a descriptor batches = _DLLGlobal(c_int32, 'n_batches') + entropy_on = _DLLGlobal(c_bool, 'entropy_on') generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') particles = _DLLGlobal(c_int64, 'n_particles') diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ec07a4d04c..32a85cabfd 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -254,7 +254,7 @@ contains ! Set the energy bin if needed if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 + filter_matches(i_filter_ein) % bins % data(1) = 12*(ng - h) + 1 end if score_index = 0 diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index 99d482480a..01404d9b0e 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -146,8 +146,13 @@ contains call loss % assemble() call prod % assemble() if (cmfd_write_matrices) then - call loss % write('loss.dat') - call prod % write('prod.dat') + if (adjoint) then + call loss % write('adj_loss.dat') + call prod % write('adj_prod.dat') + else + call loss % write('loss.dat') + call prod % write('prod.dat') + end if end if ! Set norms to 0 @@ -740,7 +745,7 @@ contains else filename = 'fluxvec.dat' end if - ! TODO: call phi_n % write(filename) + call phi_n % write(filename) end if end subroutine extract_results diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 367214a930..db84056812 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -280,6 +280,7 @@ contains end subroutine tally_allocate_results + function tally_set_filters(this, filter_indices) result(err) class(TallyObject), intent(inout) :: this integer(C_INT32_T), intent(in) :: filter_indices(:) @@ -1025,6 +1026,7 @@ contains end if end function openmc_tally_set_scores + function openmc_tally_set_type(index, type) result(err) bind(C) ! Update the type of a tally that is already allocated integer(C_INT32_T), value, intent(in) :: index diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 1ada0fba74..dadaaa48ca 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -14,7 +14,7 @@ module vector_header procedure :: destroy => vector_destroy procedure :: add_value => vector_add_value procedure :: copy => vector_copy - ! TODO: procedure :: write => vector_write + procedure :: write => vector_write end type Vector contains @@ -88,4 +88,19 @@ contains end subroutine vector_copy +!=============================================================================== +! VECTOR_WRITE writes a vector to file +!=============================================================================== + subroutine vector_write(self, filename) + class(Vector), target, intent(inout) :: self ! vector instance + character(*), intent(in) :: filename ! filename to output to + integer :: unit_ + integer :: i + open(newunit=unit_, file=filename) + do i = 1, self % n + write(unit_,*) i, self % data(i) + end do + close(unit_) + end subroutine vector_write + end module vector_header From c818b5a50b945d51425e9cab32161a45b97129f0 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 20 Sep 2018 14:06:13 -0400 Subject: [PATCH 02/38] Update spacing --- openmc/capi/core.py | 12 ++++++++---- src/vector_header.F90 | 21 ++++++++++++++------- 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/openmc/capi/core.py b/openmc/capi/core.py index 7e4d29afc3..fa6517488d 100644 --- a/openmc/capi/core.py +++ b/openmc/capi/core.py @@ -64,11 +64,13 @@ def calculate_volumes(): def current_batch(): """Return the current batch of the simulation. - Returns + + Returns ------- int Current batch of the simulation - """ + + """ return c_int.in_dll(_dll, 'openmc_current_batch').value @@ -226,11 +228,13 @@ def keff(): def master(): """Return whether processor is master processor or not. - Returns + + Returns ------- bool Whether is master processor or not - """ + + """ return c_bool.in_dll(_dll, 'openmc_master').value diff --git a/src/vector_header.F90 b/src/vector_header.F90 index dadaaa48ca..27d773f888 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -91,16 +91,23 @@ contains !=============================================================================== ! VECTOR_WRITE writes a vector to file !=============================================================================== - subroutine vector_write(self, filename) - class(Vector), target, intent(inout) :: self ! vector instance + + subroutine vector_write(self, filename) + + class(Vector), target, intent(inout) :: self ! vector instance character(*), intent(in) :: filename ! filename to output to - integer :: unit_ + + integer :: unit_ integer :: i - open(newunit=unit_, file=filename) - do i = 1, self % n + + open(newunit=unit_, file=filename) + + do i = 1, self % n write(unit_,*) i, self % data(i) end do - close(unit_) - end subroutine vector_write + + close(unit_) + + end subroutine vector_write end module vector_header From 8a551b807450b60aa8b48363028b3eede4ef3dda Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 20 Sep 2018 14:33:13 -0400 Subject: [PATCH 03/38] Expose run_CE througb C API --- openmc/capi/settings.py | 1 + 1 file changed, 1 insertion(+) diff --git a/openmc/capi/settings.py b/openmc/capi/settings.py index af6dd7d7a2..115962395d 100644 --- a/openmc/capi/settings.py +++ b/openmc/capi/settings.py @@ -22,6 +22,7 @@ class _Settings(object): generations_per_batch = _DLLGlobal(c_int32, 'gen_per_batch') inactive = _DLLGlobal(c_int32, 'n_inactive') particles = _DLLGlobal(c_int64, 'n_particles') + run_CE = _DLLGlobal(c_bool, 'run_CE') verbosity = _DLLGlobal(c_int, 'verbosity') @property From 852cbb0109f3db6cd5f79d41902a66a6a617788d Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 25 Sep 2018 19:53:28 -0400 Subject: [PATCH 04/38] Initial pass at fixing issue #1053 --- src/simulation.F90 | 21 +++++++++++++++------ 1 file changed, 15 insertions(+), 6 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index 176dfb921b..481e9203b7 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -204,7 +204,8 @@ contains ! Reset total starting particle weight used for normalizing tallies total_weight = ZERO - if (n_inactive > 0 .and. current_batch == 1) then + if ((n_inactive > 0 .and. current_batch == 1) .or. & + (restart_run .and. restart_batch <= n_inactive .and. current_batch == restart_batch)) then ! Turn on inactive timer call time_inactive % start() elseif (current_batch == n_inactive + 1) then @@ -394,7 +395,7 @@ contains subroutine replay_batch_history ! Write message at beginning - if (current_batch == 1) then + if (n_realizations == 0) then call write_message("Replaying history from state point...", 6) end if @@ -414,10 +415,12 @@ contains end if ! Increment n_realizations as would ordinarily be done in finalize_batch - if (reduce_tallies) then - n_realizations = n_realizations + 1 - else - n_realizations = n_realizations + n_procs + if (current_batch > n_inactive) then + if (reduce_tallies) then + n_realizations = n_realizations + 1 + else + n_realizations = n_realizations + n_procs + end if end if ! Write message at end @@ -503,6 +506,12 @@ contains else current_batch = restart_batch - n_realizations*n_procs end if + ! If simulation restarted from inactive batch, decrement current_batch + ! by one to replay an additional batch so that keff is set before resuming + ! simulation + if (restart_batch <= n_inactive) then + current_batch = current_batch - 1 + end if n_realizations = 0 end if From 47cd40849c76bb51eb75671a83f2231a69dccda4 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Wed, 26 Sep 2018 17:21:15 -0400 Subject: [PATCH 05/38] Remove whitespace, add more comments --- src/simulation.F90 | 4 ++-- src/vector_header.F90 | 14 +++++++------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index 481e9203b7..f9fcc62bd5 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -205,7 +205,7 @@ contains total_weight = ZERO if ((n_inactive > 0 .and. current_batch == 1) .or. & - (restart_run .and. restart_batch <= n_inactive .and. current_batch == restart_batch)) then + (restart_run .and. restart_batch <= n_inactive .and. current_batch == restart_batch)) then ! Turn on inactive timer call time_inactive % start() elseif (current_batch == n_inactive + 1) then @@ -414,7 +414,7 @@ contains end do end if - ! Increment n_realizations as would ordinarily be done in finalize_batch + ! Increment n_realizations if in active batches as would ordinarily be done in finalize_batch if (current_batch > n_inactive) then if (reduce_tallies) then n_realizations = n_realizations + 1 diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 27d773f888..fd7726c79e 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -91,23 +91,23 @@ contains !=============================================================================== ! VECTOR_WRITE writes a vector to file !=============================================================================== - + subroutine vector_write(self, filename) - + class(Vector), target, intent(inout) :: self ! vector instance character(*), intent(in) :: filename ! filename to output to - + integer :: unit_ integer :: i - + open(newunit=unit_, file=filename) - + do i = 1, self % n write(unit_,*) i, self % data(i) end do - + close(unit_) - + end subroutine vector_write end module vector_header From 3ec09e8124116d3ce6d1ca867a9809a875bbab0b Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 2 Oct 2018 18:15:55 -0400 Subject: [PATCH 06/38] Redefine score_index to reflect correct order of bins --- src/cmfd_data.F90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 32a85cabfd..bf33bbd8eb 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -254,7 +254,7 @@ contains ! Set the energy bin if needed if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = 12*(ng - h) + 1 + filter_matches(i_filter_ein) % bins % data(1) = ng - h end if score_index = 0 @@ -265,39 +265,39 @@ contains ! Left surface cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_LEFT) + score_index + ng*OUT_LEFT) cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_LEFT) + score_index + ng*IN_LEFT) ! Right surface cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_RIGHT) + score_index + ng*IN_RIGHT) cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_RIGHT) + score_index + ng*OUT_RIGHT) ! Back surface cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_BACK) + score_index + ng*OUT_BACK) cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_BACK) + score_index + ng*IN_BACK) ! Front surface cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_FRONT) + score_index + ng*IN_FRONT) cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_FRONT) + score_index + ng*OUT_FRONT) ! Left surface cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_BOTTOM) + score_index + ng*OUT_BOTTOM) cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_BOTTOM) + score_index + ng*IN_BOTTOM) ! Right surface cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + IN_TOP) + score_index + ng*IN_TOP) cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + OUT_TOP) + score_index + ng*OUT_TOP) else if (ital == 4) then From f216ee3fb43da3a18062d90e00eb315f5c135416 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 2 Oct 2018 22:02:54 -0400 Subject: [PATCH 07/38] First pass at removing replay_batch_history --- src/simulation.F90 | 71 ++------------------------------------------- src/state_point.F90 | 16 +++++++++- 2 files changed, 17 insertions(+), 70 deletions(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index f9fcc62bd5..b51962a17e 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -87,13 +87,6 @@ contains call initialize_batch() - ! Handle restart runs - if (restart_run .and. current_batch <= restart_batch) then - call replay_batch_history() - status = STATUS_EXIT_NORMAL - return - end if - ! ======================================================================= ! LOOP OVER GENERATIONS GENERATION_LOOP: do current_gen = 1, gen_per_batch @@ -205,7 +198,7 @@ contains total_weight = ZERO if ((n_inactive > 0 .and. current_batch == 1) .or. & - (restart_run .and. restart_batch <= n_inactive .and. current_batch == restart_batch)) then + (restart_run .and. restart_batch < n_inactive .and. current_batch == restart_batch + 1)) then ! Turn on inactive timer call time_inactive % start() elseif (current_batch == n_inactive + 1) then @@ -387,48 +380,6 @@ contains end subroutine finalize_batch -!=============================================================================== -! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a -! batch using data read from a state point file -!=============================================================================== - - subroutine replay_batch_history - - ! Write message at beginning - if (n_realizations == 0) then - call write_message("Replaying history from state point...", 6) - end if - - if (run_mode == MODE_EIGENVALUE) then - do current_gen = 1, gen_per_batch - call calculate_average_keff() - - ! print out batch keff - if (verbosity >= 7) then - if (current_gen < gen_per_batch) then - if (master) call print_generation() - else - if (master) call print_batch_keff() - end if - end if - end do - end if - - ! Increment n_realizations if in active batches as would ordinarily be done in finalize_batch - if (current_batch > n_inactive) then - if (reduce_tallies) then - n_realizations = n_realizations + 1 - else - n_realizations = n_realizations + n_procs - end if - end if - - ! Write message at end - if (current_batch == restart_batch) then - call write_message("Resuming simulation...", 6) - end if - - end subroutine replay_batch_history !=============================================================================== ! INITIALIZE_SIMULATION @@ -492,29 +443,11 @@ contains ! file if (restart_run) then call load_state_point() + call write_message("Resuming simulation...", 6) else call initialize_source() end if - ! In restart, set the batch to begin from in order to reproduce the correct - ! average keff (used in sampling the fission bank). Use n_realizations from - ! the statepoint rather than n_inactive in case openmc_reset was called in - ! the previous run. - if (restart_run) then - if (reduce_tallies) then - current_batch = restart_batch - n_realizations - else - current_batch = restart_batch - n_realizations*n_procs - end if - ! If simulation restarted from inactive batch, decrement current_batch - ! by one to replay an additional batch so that keff is set before resuming - ! simulation - if (restart_batch <= n_inactive) then - current_batch = current_batch - 1 - end if - n_realizations = 0 - end if - ! Display header if (master) then if (run_mode == MODE_FIXEDSOURCE) then diff --git a/src/state_point.F90 b/src/state_point.F90 index 526066c529..91d0f10908 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -16,7 +16,7 @@ module state_point use bank_header, only: Bank use cmfd_header use constants - use eigenvalue, only: openmc_get_keff + use eigenvalue, only: openmc_get_keff, k_sum use endf, only: reaction_name use error, only: fatal_error, warning, write_message use hdf5_interface @@ -749,6 +749,20 @@ contains ! Read number of realizations for global tallies call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) + ! Set k_sum, keff, and current_batch based on whether restart file is part + ! of active cycle or inactive cycle + if (restart_batch > n_inactive) then + do i = n_inactive + 1, restart_batch + k_sum(1) = k_sum(1) + k_generation % data(i) + k_sum(2) = k_sum(2) + k_generation % data(i)**2 + end do + n = gen_per_batch*n_realizations + keff = k_sum(1) / n + else + keff = k_generation % data(n) + end if + current_batch = restart_batch + ! Check to make sure source bank is present if (path_source_point == path_state_point .and. .not. source_present) then call fatal_error("Source bank must be contained in statepoint restart & From 6ad57ae11dbba47dd335e3ec1f9e08f19b598433 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Tue, 2 Oct 2018 22:13:52 -0400 Subject: [PATCH 08/38] Turn on active timer for restart run in actives --- src/simulation.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/simulation.F90 b/src/simulation.F90 index b51962a17e..b1087e21d2 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -201,7 +201,8 @@ contains (restart_run .and. restart_batch < n_inactive .and. current_batch == restart_batch + 1)) then ! Turn on inactive timer call time_inactive % start() - elseif (current_batch == n_inactive + 1) then + elseif ((current_batch == n_inactive + 1) .or. & + (restart_run .and. restart_batch > n_inactive .and. current_batch == restart_batch + 1)) then ! Switch from inactive batch timer to active batch timer call time_inactive % stop() call time_active % start() From 8feefcd13a5b14fa174da562b29f9d0653bbc9d4 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Wed, 3 Oct 2018 11:38:07 -0400 Subject: [PATCH 09/38] Fix travis error --- tests/unit_tests/test_capi.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 9d3bc106d2..ae2d927191 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -375,13 +375,13 @@ def test_restart(capi_init): openmc.capi.next_batch() keff0 = openmc.capi.keff() - # Restart the simulation from the statepoint and the 5 active batches. + # Restart the simulation from the statepoint and the 3 remaining active batches. openmc.capi.simulation_finalize() openmc.capi.hard_reset() openmc.capi.finalize() openmc.capi.init(args=('-r', 'restart_test.h5')) openmc.capi.simulation_init() - for i in range(5): + for i in range(3): openmc.capi.next_batch() keff1 = openmc.capi.keff() openmc.capi.simulation_finalize() From 794ec19c5a7c6b6eef2a1b5fbcc64a278df26012 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 4 Oct 2018 19:29:59 -0500 Subject: [PATCH 10/38] Photon physics documentation --- docs/source/methods/photon_physics.rst | 679 ++++++++++++++++++++++++- 1 file changed, 660 insertions(+), 19 deletions(-) diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 4d58ee0db2..e006c1dd60 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -85,26 +85,26 @@ momentum transfer is traditionally expressed as .. math:: :label: momentum-transfer - x = \kappa \alpha \sqrt{1 - \mu} + x = a \kappa \sqrt{1 - \mu} -where :math:`\alpha` is the ratio of the photon energy to the electron rest -mass, and the coefficient :math:`\kappa` can be shown to be +where :math:`\kappa` is the ratio of the photon energy to the electron rest +mass, and the coefficient :math:`a` can be shown to be .. math:: - :label: kappa + :label: omega - \kappa = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, + a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`, -we have :math:`\mu = 1 - [x/(\kappa\alpha)]^2` and :math:`d\mu/dx^2 = --1/(\kappa\alpha)^2`. The probability density in :math:`x^2` is +we have :math:`\mu = 1 - [x/(a\kappa)]^2` and :math:`d\mu/dx^2 = +-1/(a\kappa)^2`. The probability density in :math:`x^2` is .. math:: :label: coherent-pdf-x2 p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi - r_e^2 A(\bar{x}^2,Z)}{(\kappa\alpha)^2 \sigma(E)} \left ( + r_e^2 A(\bar{x}^2,Z)}{(a\kappa)^2 \sigma(E)} \left ( \frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2 where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for @@ -113,7 +113,7 @@ where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for .. math:: :label: xmax - \bar{x} = \kappa \alpha \sqrt{2} = \frac{m_e c^2}{hc} \alpha, + \bar{x} = a \kappa \sqrt{2} = \frac{m_e c^2}{hc} \kappa, and :math:`A(x^2, Z)` is the integral of the square of the form factor: @@ -154,6 +154,8 @@ section. The complete algorithm is as follows: 6. If :math:`\xi_2 < (1 + \mu^2)/2`, accept :math:`\mu`. Otherwise, repeat the sampling at step 3. +.. _incoherent-sampling: + Incoherent (Compton) Scattering ------------------------------- @@ -166,11 +168,11 @@ the two authors who discovered it: .. math:: :label: klein-nishina - \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\alpha'}{\alpha} \right - )^2 \left [ \frac{\alpha'}{\alpha} + \frac{\alpha}{\alpha'} + \mu^2 - 1 + \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\kappa'}{\kappa} \right + )^2 \left [ \frac{\kappa'}{\kappa} + \frac{\kappa}{\kappa'} + \mu^2 - 1 \right ] -where :math:`\alpha` and :math:`\alpha'` are the ratios of the incoming and +where :math:`\kappa` and :math:`\kappa'` are the ratios of the incoming and exiting photon energies to the electron rest mass energy equivalent (0.511 MeV), respectively. Although it appears that the outgoing energy and angle are separate, there is actually a one-to-one relationship between them such that @@ -179,9 +181,9 @@ only one needs to be sampled: .. math:: :label: compton-energy-angle - \alpha' = \frac{\alpha}{1 + \alpha(1 - \mu)}. + \kappa' = \frac{\kappa}{1 + \kappa(1 - \mu)}. -Note that when :math:`\alpha'/\alpha` goes to one, i.e., scattering is elastic, +Note that when :math:`\kappa'/\kappa` goes to one, i.e., scattering is elastic, the Klein-Nishina cross section becomes identical to the Thomson cross section. In general though, the scattering is inelastic and is known as Compton scattering. When a photon interacts with a bound electron in an atom, the @@ -193,16 +195,16 @@ differential cross section for incoherent scattering is given by :label: incoherent-xs \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left ( - \frac{\alpha'}{\alpha} \right )^2 \left [ \frac{\alpha'}{\alpha} + - \frac{\alpha}{\alpha'} + \mu^2 - 1 \right ] S(x,Z) + \frac{\kappa'}{\kappa} \right )^2 \left [ \frac{\kappa'}{\kappa} + + \frac{\kappa}{\kappa'} + \mu^2 - 1 \right ] S(x,Z) where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first sample the Klein-Nishina cross section and then perform rejection sampling on the form factor. As in other codes, `Kahn's rejection method`_ is used for -:math:`\alpha < 3` and a direct method by Koblinger_ is used for :math:`\alpha +:math:`\kappa < 3` and a direct method by Koblinger_ is used for :math:`\kappa \ge 3`. The complete algorithm is as follows: -1. If :math:`\alpha < 3`, sample :math:`\mu` from the Klein-Nishina cross +1. If :math:`\kappa < 3`, sample :math:`\mu` from the Klein-Nishina cross section using Kahn's rejection method. Otherwise, use Koblinger's direct method. @@ -215,19 +217,367 @@ the form factor. As in other codes, `Kahn's rejection method`_ is used for Doppler Energy Broadening +++++++++++++++++++++++++ -LA-UR-04-0487_ and LA-UR-04-0488_ +Bound electrons are not at rest but have a momentum distribution that will +cause the energy of the scattered photon to be Doppler broadened. More tightly +bound electrons have a wider momentum distribution, so the energy spectrum of +photons scattering off inner shell electrons will be broadened the most. +In addition, scattering from bound electrons places a limit on the maximum +scattered photon energy: + +.. math:: + :label: max-energy-out + + E'_{\text{max}} = E - E_{b,i}, + +where :math:`E_{b,i}` is the binding energy of the :math:`i`-th subshell. + +Compton profiles :math:`J_i(p_z)` are used to account for the binding effects. +The quantity :math:`p_z = {\bf p} \cdot {\bf q}/q` is the projection of the +initial electron momentum on :math:`{\bf q}`, where the scattering vector +:math:`{\bf q} = {\bf p} - {\bf p'}` is the momentum gained by the photon, +:math:`{\bf p}` is the initial momentum of the electron, and :math:`{\bf p'}` +is the momentum of the scattered electron. Applying the conservation of energy +and momentum, :math:`p_z` can be written in terms of the photon energy and +scattering angle: + +.. math:: + :label: pz + + p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 - 2EE'\mu}}, + +where :math:`\alpha` is the fine structure constant. The maximum momentum +transferred, :math:`p_{z,\text{max}}`, can be calculated from :eq:`pz` using +:math:`E' = E'_{\text{max}}`. The Compton profile of the :math:`i`-th electron +subshell is defined as + +.. math:: + :label: compton-profile + + J_i(p_z) = \int \int \rho_i({\bf p}) dp_x dp_y, + +where :math:`\rho_i({\bf p})` is the initial electron momentum distribution. +:math:`J_i(p_z)` can be interpreted as the probability density function of +:math:`p_z`. + +The Doppler broadened energy of the Compton-scattered can be sampled by +selecting an electron shell, sampling a value of :math:`p_z` using the Compton +profile, and calculating the scattered photon energy. The methods and algorithm +used to do this are described in LA-UR-04-0487_ and LA-UR-04-0488_. The +sampling algorithm is summarized below: + +1. Sample :math:`\mu` from :eq:`incoherent-xs` using the algorithm described in + :ref:`incoherent-sampling`. + +2. Sample the electron subshell :math:`i` using the number of electrons per + shell as the PDF. + +3. Sample :math:`p_z` using :math:`J_i(p_z)` as the PDF. + +4. Calculate :math:`E'` by solving :eq:`pz` for :math:`E'` using the sampled + value of :math:`p_z`. + +5. If :math:`p_z < p_{z,\text{max}}` for shell :math:`i`, accept :math:`E'`. + Otherwise repeat from step 2. Compton Electrons +++++++++++++++++ +Because the Compton-scattered photons can transfer a large fraction of their +energy to the kinetic energy of the recoil electron, which may in turn go on to +lose its energy as bremsstrahlung radiation, it is necessary to accurately +model the angular and energy distributions of Compton electrons. The energy of +the recoil electron ejected from the :math:`i`-th subshell is simply given by + +.. math:: + :label: compton-electron-energy + + E_{-} = E - E' - E_{b,i}. + +The direction of the electron is assumed to be in the direction of the momentum +transfer, with the cosine of the polar angle given by + +.. math:: + :label: compton-electron-mu + + \mu_{-} = \frac{E - E'\mu}{\sqrt{E^2 +E'^2 - 2EE'\mu}} + +and the azimuthal angle :math:`\phi_{-} = \phi + \pi`, where :math:`\phi` is +the azimuthal angle of the photon. The vacancy left by the ejected electron is +filled through atomic relaxation. Photoelectric Effect -------------------- +In the photoelectric effect, the incident photon is absorbed by an atomic +electron, which is then emitted from the :math:`i`-th shell with kinetic energy + +.. math:: + :label: photoelectron-kinetic-energy + + E_{-} = E - E_{b,i}. + +Photoelectric emission is only possible when the photon energy exceeds the +binding energy of the shell. These binding energies are often referred to as +edge energies because the otherwise continuously decreasing cross section has +discontinuities at these points, creating the characteristic sawtooth shape. +The photoelectric effect dominates at low energies and is more important for +heavier elements. + +When simulating the photoelectric effect, the first step is to sample the +electron shell. The shell :math:`i` where the ionization occurs can be +considered a discrete random variable with probability density function + +.. math:: + :label: photoelectron-shell-pdf + + p_i = \frac{\sigma_{\text{pe},i}}{\sigma_{\text{pe}}}, + +where :math:`\sigma_{\text{pe},i}` is the cross section of the :math:`i`-th +shell, and the total photoelectric cross section of the atom, +:math:`\sigma_{\text{pe}}`, is the sum over the shell cross sections. Once the +shell has been sampled, the energy of the photoelectron is calculated using +:eq:`photoelectron-kinetic-energy`. + +To determine the direction of the photoelectron, we implement the method +described in [Kaltiaisenaho]_, which models the angular distribution of the +photoelectrons using the K-shell cross section derived by Sauter (K-shell +electrons are the most tightly bound, and they contribute the most to +:math:`\sigma_{\text{pe}}`). The non-relativistic Sauter distribution for +unpolarized photons can be approximated as + +.. math:: + :label: sauter + + \frac{d\sigma_{\text{pe}}}{d\mu_{-}} \propto + \frac{1 - \mu_{-}^2}{(1 - \beta_{-} \mu_{-})^4}, + +where :math:`\beta_{-}` is the ratio of the velocity of the electron to the +speed of light, + +.. math:: + :label: beta-2 + + \beta_{-} = \frac{\sqrt{(E_{-}(E_{-} + 2m_e c^2)}}{E_{-} + m_e c^2}. + +To sample :math:`\mu_{-}` from the Sauter distribution, we first express +:eq:`sauter` in the form: + +.. math:: + :label: photoelectron-mu-pdf + + f(\mu_{-}) = \frac{3}{2} \psi(\mu_{-}) g(\mu_{-}), + +where + +.. math:: + :label: mu-pdf-factors + + \psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 - \beta_{-}\mu_{-})^2}, \\ + g(\mu_{-}) &= \frac{1 - \beta_{-}^2}{2 (1 - \beta_{-}\mu_{-})^2}. + +In the interval :math:`(-1, 1)`, :math:`g(\mu_{-})` is a normalized PDF and +:math:`\psi(\mu_{-})` satisfies the condition :math:`0 < \psi(\mu_{-}) < 1`. +The following algorithm can now be used to sample :math:`\mu_{-}`: + +1. Using the inverse transform method, sample :math:`\mu_{-}` from + :math:`g(\mu_{-})` using the sampling formula + + .. math:: + + \mu_{-} = \frac{2\xi_1 + \beta_{-} - 1}{2\beta_{-}\xi_1 - \beta_{-} + 1}. + +2. If :math:`\xi_2 \le \psi(\mu_{-})`, accept :math:`\mu_{-}`. Otherwise, + repeat the sampling from step 1. + +The azimuthal angle is sampled uniformly on :math:`[0, 2\pi)`. + +The atom is left in an excited state with a vacancy in the :math:`i`-th shell +and decays to its ground state through a cascade of transitions that produce +fluorescent photons and Auger electrons. Pair Production --------------- +In electron-positron pair production, a photon is absorbed in the vicinity of +an atomic nucleus or an electron and an electron and positron are created. Pair +production is the dominant interaction with matter at high photon energies and +is more important for high-Z elements. When it takes place in the field of a +nucleus, energy is essentially conserved among the incident photon and the +resulting charged particles. Therefore, in order for pair production to occur, +the photon energy must be greater than the sum of the rest mass energies of the +electron and positron, i.e., :math:`E_{\text{threshold,pp}} = 2 m_e c^2 = +1.022` MeV. + +The photon can also interact in the field of an atomic electron. This process +is referred to as "triplet production" because the target electron is ejected +from the atom and three charged particles emerge from the interaction. In this +case, the recoiling electron also absorbs some energy, so the energy threshold +for triplet production is greater than that of pair production from atomic +nuclei, with :math:`E_{\text{threshold,tp}} = 4 m_e c^2 = 2.044` MeV. The ratio +of the triplet production cross section to the pair production cross section is +approximately 1/Z, so triplet production becomes increasingly unimportant for +high-Z elements. Though it can be significant in lighter elements, the momentum +of the recoil electron becomes negligible in the energy regime where pair +production dominates. For our purposes, it is a good approximation to treat +triplet production as pair production and only simulate the electron-positron +pair. + +Accurately modeling the creation of electron-positron pair is important because +the charged particles can go on to lose much of their energy as bremsstrahlung +radiation, and the subsequent annihilation of the positron with an electron +produces two additional photons. We sample the energy and direction of the +charged particles using a semiempirical model described in [Salvat]_. The +Bethe-Heitler differential cross section, given by + +.. math:: + :label: bethe-heitler + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = \alpha r_e^2 Z^2 + \left[ (\epsilon^2 + (1-\epsilon)^2) (\Phi_1 - 4f_C) + + \frac{2}{3}\epsilon(1-\epsilon)(\Phi_2 - 4f_C) \right], + +is used as a starting point, where :math:`\alpha` is the fine structure +constant, :math:`f_C` is the Coulomb correction function, :math:`\Phi_1` and +:math:`\Phi_2` are screening functions, and :math:`\epsilon = (E_{-} + m_e +c^2)/E` is the electron reduced energy (i.e., the fraction of the photon energy +given to the electron). :math:`\epsilon` can take values between +:math:`\epsilon_{\text{min}} = \kappa^{-1}` (when the kinetic energy of the +electron is zero) and :math:`\epsilon_{\text{max}} = 1 - \kappa^{-1}` (when the +kinetic energy of the positron is zero). + +The Coulomb correction, given by + +.. math:: + :label: coulomb-correction + + f_C = \alpha^{2}Z^{2} \big[&(1 + \alpha^{2}Z^{2})^{-1} + 0.202059 + - 0.03693\alpha^{2}Z^{2} + 0.00835\alpha^{4}Z^{4} \\ + &- 0.00201\alpha^{6}Z^{6} + 0.00049\alpha^{8}Z^{8} + - 0.00012\alpha^{10}Z^{10} + 0.00003\alpha^{12}Z^{12}\big] + +is introduced to correct for the fact that the Bethe-Heitler differential cross +section was derived using the Born approximation, which treats the Coulomb +interaction as a small perturbation. + +The screening functions :math:`\Phi_1` and :math:`\Phi_2` account for the +screening of the Coulomb field of the atomic nucleus by outer electrons. Since +they are given by integrals which include the atomic form factor, they must be +computed numerically for a realistic form factor. However, by assuming +exponential screening and using a simplified form factor, analytical +approximations of the screening functions can be derived: + +.. math:: + :label: screening-functions + + \Phi_1 &= 2 - 2\ln(1 + b^2) - 4b\arctan(b^{-1}) + 4\ln(Rm_{e}c/\hbar) \\ + \Phi_2 &= \frac{4}{3} - 2\ln(1 + b^2) + 2b^2 \left[ 4 - 4b\arctan(b^{-1}) + - 3\ln(1 + b^{-2}) \right] + 4\ln(Rm_{e}c/\hbar) + +where + +.. math:: + :label: b + + b = \frac{Rm_{e}c}{2\kappa\epsilon(1 - \epsilon)\hbar}. + +and :math:`R` is the screening radius. + +The differential cross section in :eq:`bethe-heitler` with the approximations +described above will not be accurate at low energies: the lower boundary of +:math:`\epsilon` will be shifted above :math:`\epsilon_{\text{min}}` and the +upper boundary of :math:`\epsilon` will be shifted below +:math:`\epsilon_{\text{max}}`. To offset this behavior, a correcting factor +:math:`F_0(\kappa, Z)` is used: + +.. math:: + :label: correcting-factor + + F_0(\kappa, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/\kappa)^{1/2} \\ + &+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/\kappa) \\ + &- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/\kappa)^{3/2} \\ + &+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/\kappa)^{2}. + +To aid sampling, the differential cross section used to sample :math:`\epsilon` +(minus the normalization constant) can now be expressed in the form + +.. math:: + :label: pp-pdf + + \frac{d\sigma_{\text{pp}}}{d\epsilon} = + u_1 \frac{\phi_1(\epsilon)}{\phi_1(1/2)} \pi_1(\epsilon) + + u_2 \frac{\phi_2(\epsilon)}{\phi_2(1/2)} \pi_2(\epsilon) + +where + +.. math:: + :label: u + + u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^2 \phi_1(1/2), \\ + u_2 &= \phi_2(1/2), + +.. math:: + :label: phi + + \phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(\kappa, Z), \\ + \phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(\kappa, Z), + +and + +.. math:: + :label: pi + + \pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^{-3} + \left(\frac{1}{2} - \epsilon\right)^2, \\ + \pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^{-1}. + +The functions in :eq:`phi` are non-negative and maximum at :math:`\epsilon = +1/2`. In the interval :math:`(\epsilon_{\text{min}}, \epsilon_{\text{max}})`, +the functions in :eq:`pi` are normalized PDFs and +:math:`\phi_i(\epsilon)/\phi_i(1/2)` satisfies the condition :math:`0 < +\phi_i(\epsilon)/\phi_i(1/2) < 1`. The following algorithm can now be used to +sample the reduced electron energy :math:`\epsilon`: + +1. Sample :math:`i` according to the point probabilities + :math:`p(i=1) = u_1/(u_1 + u_2)` and :math:`p(i=2) = u_2/(u_1 + u_2)`. + +2. Using the inverse transform method, sample :math:`\epsilon` from + :math:`\pi_i(\epsilon)` using the sampling formula + + .. math:: + + \epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{\kappa}\right) + (2\xi_1 - 1)^{1/3} ~~~~&\text{if}~~ i = 1 \\ + \epsilon &= \frac{1}{\kappa} + \left(\frac{1}{2} - + \frac{1}{\kappa}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2. + +3. If :math:`\xi_2 \le \phi_i(\epsilon)/\phi_i(1/2)`, accept + :math:`\epsilon`. Otherwise, repeat the sampling from step 1. + +Because charged particles have a much smaller range than the mean free path of +photons and because they immediately undergo multiple scattering events which +randomize their direction, it is sufficient to use a simplified model to sample +the direction of the electron and positron. The cosines of the polar angles are +sampled using the leading order term of the Sauter–Gluckstern–Hull +distribution, + +.. math:: + :label: sauter–gluckstern–hull + + \frac{d\sigma_{pp}}{d\Omega_{\pm}} = C(1 - \beta_{\pm}\mu_{\pm})^{-2}, + +where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the ratio of +the velocity of the charged particle to the speed of light given in :eq:`beta-2`. + +The inverse transform method is used to sample :math:`\mu_{-}` and :math:`\mu_{+}` +from :eq:`sauter–gluckstern–hull`, using the sampling formula + +.. math:: + :label: sample-mu + + \mu_{\pm} = \frac{2\xi - 1 + \beta_{\pm}}{(2\xi - 1)\beta_{\pm} + 1}. + +The azimuthal angles for the electron and positron are sampled independently +and uniformly on :math:`[0, 2\pi)`. ------------------- Secondary Processes @@ -248,19 +598,301 @@ additional photons. Atomic Relaxation ----------------- +When an electron is ejected from an atom and a vacancy is left in an inner +shell, an electron from a higher energy level will fill the vacancy. This +results in either a radiative transition, in which a photon with a +characteristic energy (fluorescence photon) is emitted, or non-radiative +transition, in which an electron from a shell that is farther out (Auger +electron) is emitted. If a non-radiative transition occurs, the new vacancy is +filled in the same manner, and as the process repeats a shower of photons and +electrons can be produced. + +The energy of a fluorescence photon is the equal to the energy difference +between the transition states, i.e., + +.. math:: + :label: fluorescence-photon-energy + + E = E_{b,v} - E_{b,i}, + +where :math:`E_{b,v}` is the binding energy of the vacancy shell and +:math:`E_{b,i}` is the binding energy of the shell from which the electron +transitioned. The energy of an Auger electron is given by + +.. math:: + :label: auger-electron-energy + + E_{-} = E_{b,v} - E_{b,i} - E_{b,a}, + +where :math:`E_{b,a}` is the binding energy of the shell from which the Auger +electron is emitted. While Auger electrons are low-energy so their range and +bremsstrahlung yield is small, fluorescence photons can travel far before +depositing their energy, so the relaxation process should be modeled in detail. + +Transition energies and probabilities are needed for each subshell to simulate +atomic relaxation. Starting with the initial shell vacancy, the following +recursive algorithm is used to fill vacancies and create fluorescence photons +and Auger electrons: + +1. If there are no transitions for the vacancy shell, create a fluorescence + photon assuming it is from a captured free electron and terminate. + +2. Sample a transition using the transition probabilities for the vacancy + shell as the PDF. + +3. Create either a fluorescence photon or Auger electron, sampling the + direction of the particle isotropically. + +4. If a non-radiative transition occurred, repeat from step 1 for the vacancy + left by the emitted Auger electron. + +5. Repeat from step 1 for vacancy left by the transition electron. Electron-Positron Annihilation ------------------------------ +When a positron collides with an electron, both particles are annihilated and +generally two photons with equal energy are created. If the kinetic energy of +the positron is high enough, the two photons can have different energies, and +the higher-energy photon is emitted preferentially in the direction of flight +of the positron. It is also possible to produce a single photon if the +interaction occurs with a bound electron, and in some cases three (or, rarely, +even more) photons can be emitted. However, the annihilation cross section is +largest for low-energy positrons, and as the positron energy decreases, the +angular distribution of the emitted photons becomes isotropic. + +In OpenMC, we assume the most likely case in which a low-energy positron (which +has already lost most of its energy to bremsstrahlung radiation) interacts with +an electron which is free and at rest. Two photons with energy equal to the +electron rest mass energy :math:`m_e c^2 = 0.511` MeV are emitted isotropically +in opposite directions. Bremsstrahlung -------------- +When a charged particle is decelerated in the field of an atom, some of its +kinetic energy is converted into electromagnetic radiation known as +bremsstrahlung, or 'braking radiation'. In each event, an electron or positron +with kinetic energy :math:`T` generates a photon with an energy :math:`E` +between :math:`0` and :math:`T`. Bremsstrahlung is described by a cross section +that is differential in photon energy, in the direction of the emitted photon, +and in the final direction of the charged particle. However, in Monte Carlo +simulations it is typical to integrate over the angular variables to obtain a +single differential cross section with respect to photon energy, which is often +expressed in the form + +.. math:: + :label: bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{Z^2}{\beta^2} \frac{1}{E} + \chi(Z, T, \kappa), + +where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, +\kappa)` is the scaled bremsstrahlung cross section, + +.. math:: + :label: scaled-bremsstrahlung-dcs + + \chi(Z, T, \kappa) = \frac{\beta^2}{Z^2} E \frac{d\sigma_{\text{br}}}{dE}. + +Because electrons are attracted to atomic nuclei whereas positrons are +repulsed, the cross section for positrons is smaller, though it approaches that +of electrons in the high energy limit. To obtain the positron cross section, we +multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used +in [Salvat]_, + +.. math:: + :label: positron-factor + + F_{\text{p}}(Z,T) = + & 1 - \text{exp}(-1.2359\times 10^{-1}t + 6.1274\times 10^{-2}t^2 - 3.1516\times 10^{-2}t^3 \\ + & + 7.7446\times 10^{-3}t^4 - 1.0595\times 10^{-3}t^5 + 7.0568\times 10^{-5}t^6 \\ + & - 1.8080\times 10^{-6}t^7), + +where + +.. math:: + :label: positron-factor-t + + t = \ln\left(1 + \frac{10^6}{Z^2}\frac{T}{\text{m}_\text{e}c^2} \right). + +:math:`F_{\text{p}}(Z,T)` is the ratio of the radiative stopping powers for +positrons and electrons. Stopping power describes the average energy loss per +unit path length of a charged particle as it passes through matter: + +.. math:: + :label: stopping-power + + -\frac{dT}{ds} = n \int E \frac{d\sigma}{dE} dE \equiv S(T), + +where :math:`n` is the number density of the material and :math:`d\sigma/dE` is +the cross section differential in energy loss. The total stopping power +:math:`S(T)` can be separated into two components: the radiative stopping +power :math:`S_{\text{rad}}(T)`, which refers to energy loss due to +bremsstrahlung, and the collision stopping power :math:`S_{\text{col}}(T)`, +which refers to the energy loss due to inelastic collisions with bound +electrons in the material that result in ionization and excitation. To obtain +the radiative stopping power for positrons, the radiative stopping power for +electrons is multiplied by :eq:`positron-factor`. Currently, the collision +stopping power for electrons is also used for positrons. + +While the models for photon interactions with matter described above can safely +assume interactions occur with free atoms, sampling the target atom based on +the macroscopic cross sections, molecular effects cannot necessarily be +disregarded for charged particle treatment. For compounds and mixtures, the +bremsstrahlung cross section is calculated using Bragg's additivity rule as + +.. math:: + :label: material-bremsstrahlung-dcs + + \frac{d\sigma_{\text{br}}}{dE} = \frac{1}{\beta^2 E} \sum_i \gamma_i Z^2_i + \chi(Z_i, T, \kappa), + +where the sum is over the constituent elements and :math:`\gamma_i` is the +atomic fraction of the :math:`i`-th element. Similarly, the radiative stopping +power is calculated using Bragg's additivity rule as + +.. math:: + :label: material-radiative-stopping-power + + S_{\text{rad}}(T) = \sum_i w_i S_{\text{rad},i}(T), + +where :math:`w_i` is the mass fraction of the :math:`i`-th element. The +collision stopping power, however, is a function of certain quantities such as +the mean excitation energy :math:`I` and the density effect correction +:math:`\delta_F` that depend on molecular properties. These quantities cannot +simply be summed over constituent elements in a compound, but should instead be +calculated for the material. Currently, we use Bragg's additivity rule to +calculate the collision stopping power as well, but this is not a good +approximation and should be fixed in the future. + .. _ttb: Thick-Target Bremsstrahlung Approximation +++++++++++++++++++++++++++++++++++++++++ +Since charged particles lose their energy on a much shorter distance scale than +neutral particles, not much error should be introduced by neglecting to +transport electrons. However, the bremsstrahlung emitted from high energy +electrons and positrons can travel far from the interaction site. Thus, even +without a full electron transport mode it is necessary to model bremsstrahlung. +We use a thick-target bremsstrahlung (TTB) approximation based on the models in +[Salvat]_ and [Kaltiaisenaho]_ for generating bremsstrahlung photons, which +assumes the charged particle loses all its energy in a single homogeneous +material region. + +To model bremsstrahlung using the TTB approximation, we need to know the number +of photons emitted by the charged particle and the energy distribution of the +photons. These quantities can be calculated using the continuous slowing down +approximation (CSDA). The CSDA assumes charged particles lose energy +continuously along their trajectory with a rate of energy loss equal to the +total stopping power, ignoring fluctuations in the energy loss. The +approximation is useful for expressing average quantities that describe how +charged particles slow down in matter. For example, the CSDA range approximates +the average path length a charged particle travels as it slows to rest: + +.. math:: + :label: csda-range + + R(T) = \int^T_0 \frac{dT'}{S(T')}. + +Actual path lengths will fluctuate around :math:`R(T)`. The average number of +photons emitted per unit path length is given by the inverse bremsstrahlung +mean free path: + +.. math:: + :label: inverse-bremsstrahlung-mfp + + \lambda_{\text{br}}^{-1}(T,E_{\text{cut}}) + = n\int_{E_{\text{cut}}}^T\frac{d\sigma_{\text{br}}}{dE}dE + = n\frac{Z^2}{\beta^2}\int_{\kappa_{\text{cut}}}^1\frac{1}{\kappa} + \chi(Z,T,\kappa)d\kappa. + +The lower limit of the integral in :eq:`inverse-bremsstrahlung-mfp` is non-zero +because the bremsstrahlung differential cross section diverges for small photon +energies but is finite for photon energies above some cutoff energy +:math:`E_{\text{cut}}`. The mean free path +:math:`\lambda_{\text{br}}^{-1}(T,E_{\text{cut}})` is used to calculate the +photon number yield, defined as the average number of photons emitted with +energy greater than :math:`E_{\text{cut}}` as the charged particle slows down +from energy :math:`T` to :math:`E_{\text{cut}}`. The photon number yield is +given by + +.. math:: + :label: photon-number-yield + + Y(T,E_{\text{cut}}) = \int^{R(T)}_{R(E_{\text{cut}})} + \lambda_{\text{br}}^{-1}(T',E_{\text{cut}})ds = \int_{E_{\text{cut}}}^T + \frac{\lambda_{\text{br}}^{-1}(T',E_{\text{cut}})}{S(T')}dT'. + +:math:`Y(T,E_{\text{cut}})` can be used to construct the energy spectrum of +bremsstrahlung photons: the number of photons created with energy between +:math:`E_1` and :math:`E_2` by a charged particle with initial kinetic energy +:math:`T` as it comes to rest is given by :math:`Y(T,E_1) - Y(T,E_2)`. + +To simulate the emission of bremsstrahlung photons, the total stopping power +and bremsstrahlung differential cross section for positrons and electrons must +be calculated for a given material using :eq:`material-bremsstrahlung-dcs` and +:eq:`material-radiative-stopping-power`. These quantities are used to build the +tabulated bremsstrahlung energy PDF and CDF for that material for each incident +energy :math:`T_k` on the energy grid. The following algorithm is then applied +to sample the photon energies: + +1. For an incident charged particle with energy :math:`T`, sample the number of + emitted photons as + + .. math:: + + N = \lfloor Y(T,E_{\text{cut}}) + \xi_1 \rfloor. + +2. Rather than interpolate the PDF between indices :math:`k` and :math:`k+1` + for which :math:`T_k < T < T_{k+1}`, which is computationally expensive, use + the composition method and sample from the PDF at either :math:`k` or + :math:`k+1`. Using linear interpolation on a logarithmic scale, the PDF can + be expressed as + + .. math:: + + p_{\text{br}}(T,E) = \pi_k p_{\text{br}}(T_k,E) + \pi_{k+1} + p_{\text{br}}(T_{k+1},E), + + where the interpolation weights are + + .. math:: + + \pi_k = \frac{\ln T_{k+1} - \ln T}{\ln T_{k+1} - \ln T_k},~~~ + \pi_{k+1} = \frac{\ln T - \ln T_k}{\ln T_{k+1} - \ln T_k}. + + Sample either the index :math:`i = k` or :math:`i = k+1` according to the + point probabilities :math:`\pi_{k}` and :math:`\pi_{k+1}`. + +3. Determine the maximum value of the CDF :math:`P_{\text{br,max}}`. + +3. Sample the photon energies using the inverse transform method with the + tabulated CDF :math:`P_{\text{br}}(T_i, E)` i.e., + + .. math:: + + E = E_j \left[ (1 + a_j) \frac{\xi_2 P_{\text{br,max}} - + P_{\text{br}}(T_i, E_j)} {E_j p_{\text{br}}(T_i, E_j)} + 1 + \right]^{\frac{1}{1 + a_j}} + + where the interpolation factor :math:`a_j` is given by + + .. math:: + + a_j = \frac{\ln p_{\text{br}}(T_i,E_{j+1}) - \ln p_{\text{br}}(T_i,E_j)} + {\ln E_{j+1} - \ln E_j} + + and :math:`P_{\text{br}}(T_i, E_j) \le \xi_2 P_{\text{br,max}} \le + P_{\text{br}}(T_i, E_{j+1})`. + +We ignore the range of the electron or positron, i.e., the bremsstrahlung +photons are produced in the same location that the charged particle was +created. The direction of the photons is assumed to be the same as the +direction of the incident charged particle, which is a reasonable approximation +at higher energies when the bremsstrahlung radiation is emitted at small +angles. .. _Koblinger: https://doi.org/10.13182/NSE75-A26663 @@ -273,3 +905,12 @@ Thick-Target Bremsstrahlung Approximation .. _LA-UR-04-0487: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0487.pdf .. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf + + .. rubric:: References + +.. [Kaltiaisenaho] T. Kaltiaisenaho, "Implementing a photon physics model in + Serpent 2." M.Sc. Thesis, Aalto University, 2016. + +.. [Salvat] F. Salvat, J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: + A Code System for Monte Carlo Simulation of Electron and Photon Transport," + OECD-NEA, Issy-les-Moulineaux, France (2011). From 9e0703c86d67d21e71f415e9338d6adbfc8fc948 Mon Sep 17 00:00:00 2001 From: amandalund Date: Thu, 4 Oct 2018 19:56:29 -0500 Subject: [PATCH 11/38] Proofread photon physics doc --- docs/source/methods/photon_physics.rst | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index e006c1dd60..0508225961 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -259,11 +259,11 @@ where :math:`\rho_i({\bf p})` is the initial electron momentum distribution. :math:`J_i(p_z)` can be interpreted as the probability density function of :math:`p_z`. -The Doppler broadened energy of the Compton-scattered can be sampled by +The Doppler broadened energy of the Compton-scattered photon can be sampled by selecting an electron shell, sampling a value of :math:`p_z` using the Compton -profile, and calculating the scattered photon energy. The methods and algorithm -used to do this are described in LA-UR-04-0487_ and LA-UR-04-0488_. The -sampling algorithm is summarized below: +profile, and calculating the scattered photon energy. The theory and methods +used to do this are described in detail in LA-UR-04-0487_ and LA-UR-04-0488_. +The sampling algorithm is summarized below: 1. Sample :math:`\mu` from :eq:`incoherent-xs` using the algorithm described in :ref:`incoherent-sampling`. @@ -286,7 +286,7 @@ Because the Compton-scattered photons can transfer a large fraction of their energy to the kinetic energy of the recoil electron, which may in turn go on to lose its energy as bremsstrahlung radiation, it is necessary to accurately model the angular and energy distributions of Compton electrons. The energy of -the recoil electron ejected from the :math:`i`-th subshell is simply given by +the recoil electron ejected from the :math:`i`-th subshell is given by .. math:: :label: compton-electron-energy From 4db16feec7cd417adcad77d692613f000d8f1081 Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 4 Oct 2018 21:27:50 -0400 Subject: [PATCH 12/38] Fix cmfd tally indexing --- src/cmfd_data.F90 | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index bf33bbd8eb..ef85949b9c 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -254,7 +254,7 @@ contains ! Set the energy bin if needed if (energy_filters) then - filter_matches(i_filter_ein) % bins % data(1) = ng - h + filter_matches(i_filter_ein) % bins % data(1) = ng - h + 1 end if score_index = 0 @@ -265,39 +265,39 @@ contains ! Left surface cmfd % current(1,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_LEFT) + score_index + 1 + ng*(OUT_LEFT - 1)) cmfd % current(2,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_LEFT) + score_index + 1 + ng*(IN_LEFT - 1)) ! Right surface cmfd % current(3,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_RIGHT) + score_index + 1 + ng*(IN_RIGHT - 1)) cmfd % current(4,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_RIGHT) + score_index + 1 + ng*(OUT_RIGHT - 1)) ! Back surface cmfd % current(5,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_BACK) + score_index + 1 + ng*(OUT_BACK - 1)) cmfd % current(6,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_BACK) + score_index + 1 + ng*(IN_BACK - 1)) ! Front surface cmfd % current(7,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_FRONT) + score_index + 1 + ng*(IN_FRONT - 1)) cmfd % current(8,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_FRONT) + score_index + 1 + ng*(OUT_FRONT - 1)) ! Left surface cmfd % current(9,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_BOTTOM) + score_index + 1 + ng*(OUT_BOTTOM - 1)) cmfd % current(10,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_BOTTOM) + score_index + 1 + ng*(IN_BOTTOM - 1)) ! Right surface cmfd % current(11,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*IN_TOP) + score_index + 1 + ng*(IN_TOP - 1)) cmfd % current(12,h,i,j,k) = t % results(RESULT_SUM, 1, & - score_index + ng*OUT_TOP) + score_index + 1 + ng*(OUT_TOP - 1)) else if (ital == 4) then From c9a3900842d9fb1450c00c0ae6d9d02acea4c4db Mon Sep 17 00:00:00 2001 From: Shikhar Kumar Date: Thu, 4 Oct 2018 22:07:01 -0400 Subject: [PATCH 13/38] Remove print statement in cmfd_execute.cpp --- src/cmfd_execute.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp index 7774cfe1d2..740c1f51f0 100644 --- a/src/cmfd_execute.cpp +++ b/src/cmfd_execute.cpp @@ -25,8 +25,6 @@ cmfd_populate_sourcecounts(int n_energy, const double* energies, auto& m = meshes.at(index_cmfd_mesh); xt::xarray counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside); - std::cout << counts << "\n"; - // Copy data from the xarray into the source counts array std::copy(counts.begin(), counts.end(), source_counts); } From f756d14233b4252e00e9063c16f9decab12ecdd7 Mon Sep 17 00:00:00 2001 From: amandalund Date: Fri, 5 Oct 2018 17:56:02 -0500 Subject: [PATCH 14/38] Address @paulromano comments on #1085 --- docs/source/methods/photon_physics.rst | 140 ++++++++++++------------- 1 file changed, 66 insertions(+), 74 deletions(-) diff --git a/docs/source/methods/photon_physics.rst b/docs/source/methods/photon_physics.rst index 0508225961..028a71de72 100644 --- a/docs/source/methods/photon_physics.rst +++ b/docs/source/methods/photon_physics.rst @@ -85,26 +85,26 @@ momentum transfer is traditionally expressed as .. math:: :label: momentum-transfer - x = a \kappa \sqrt{1 - \mu} + x = a k \sqrt{1 - \mu} -where :math:`\kappa` is the ratio of the photon energy to the electron rest +where :math:`k` is the ratio of the photon energy to the electron rest mass, and the coefficient :math:`a` can be shown to be .. math:: :label: omega - a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329, + a = \frac{m_e c^2}{\sqrt{2}hc} \approx 29.14329~\unicode{x212B}, where :math:`m_e` is the mass of the electron, :math:`c` is the speed of light in a vacuum, and :math:`h` is Planck's constant. Using :eq:`momentum-transfer`, -we have :math:`\mu = 1 - [x/(a\kappa)]^2` and :math:`d\mu/dx^2 = --1/(a\kappa)^2`. The probability density in :math:`x^2` is +we have :math:`\mu = 1 - [x/(ak)]^2` and :math:`d\mu/dx^2 = +-1/(ak)^2`. The probability density in :math:`x^2` is .. math:: :label: coherent-pdf-x2 p(x^2) dx^2 = p(\mu) \left | \frac{d\mu}{dx^2} \right | dx^2 = \frac{2\pi - r_e^2 A(\bar{x}^2,Z)}{(a\kappa)^2 \sigma(E)} \left ( + r_e^2 A(\bar{x}^2,Z)}{(ak)^2 \sigma(E)} \left ( \frac{1 + \mu^2}{2} \right ) \left ( \frac{F(x, Z)^2}{A(\bar{x}^2, Z)} \right ) dx^2 where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for @@ -113,7 +113,7 @@ where :math:`\bar{x}` is the maximum value of :math:`x` that occurs for .. math:: :label: xmax - \bar{x} = a \kappa \sqrt{2} = \frac{m_e c^2}{hc} \kappa, + \bar{x} = a k \sqrt{2} = \frac{m_e c^2}{hc} k, and :math:`A(x^2, Z)` is the integral of the square of the form factor: @@ -168,12 +168,11 @@ the two authors who discovered it: .. math:: :label: klein-nishina - \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{\kappa'}{\kappa} \right - )^2 \left [ \frac{\kappa'}{\kappa} + \frac{\kappa}{\kappa'} + \mu^2 - 1 - \right ] + \frac{d\sigma_{KN}}{d\mu} = \pi r_e^2 \left ( \frac{k'}{k} \right)^2 \left + [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 \right ] -where :math:`\kappa` and :math:`\kappa'` are the ratios of the incoming and -exiting photon energies to the electron rest mass energy equivalent (0.511 MeV), +where :math:`k` and :math:`k'` are the ratios of the incoming and exiting +photon energies to the electron rest mass energy equivalent (0.511 MeV), respectively. Although it appears that the outgoing energy and angle are separate, there is actually a one-to-one relationship between them such that only one needs to be sampled: @@ -181,32 +180,31 @@ only one needs to be sampled: .. math:: :label: compton-energy-angle - \kappa' = \frac{\kappa}{1 + \kappa(1 - \mu)}. + k' = \frac{k}{1 + k(1 - \mu)}. -Note that when :math:`\kappa'/\kappa` goes to one, i.e., scattering is elastic, -the Klein-Nishina cross section becomes identical to the Thomson cross -section. In general though, the scattering is inelastic and is known as Compton -scattering. When a photon interacts with a bound electron in an atom, the -Klein-Nishina formula must be modified to account for the binding effects. As in -the case of coherent scattering, this is done by means of a form factor. The -differential cross section for incoherent scattering is given by +Note that when :math:`k'/k` goes to one, i.e., scattering is elastic, the +Klein-Nishina cross section becomes identical to the Thomson cross section. In +general though, the scattering is inelastic and is known as Compton scattering. +When a photon interacts with a bound electron in an atom, the Klein-Nishina +formula must be modified to account for the binding effects. As in the case of +coherent scattering, this is done by means of a form factor. The differential +cross section for incoherent scattering is given by .. math:: :label: incoherent-xs \frac{d\sigma}{d\mu} = \frac{d\sigma_{KN}}{d\mu} S(x,Z) = \pi r_e^2 \left ( - \frac{\kappa'}{\kappa} \right )^2 \left [ \frac{\kappa'}{\kappa} + - \frac{\kappa}{\kappa'} + \mu^2 - 1 \right ] S(x,Z) + \frac{k'}{k} \right )^2 \left [ \frac{k'}{k} + \frac{k}{k'} + \mu^2 - 1 + \right ] S(x,Z) where :math:`S(x,Z)` is the form factor. The approach in OpenMC is to first sample the Klein-Nishina cross section and then perform rejection sampling on the form factor. As in other codes, `Kahn's rejection method`_ is used for -:math:`\kappa < 3` and a direct method by Koblinger_ is used for :math:`\kappa -\ge 3`. The complete algorithm is as follows: +:math:`k < 3` and a direct method by Koblinger_ is used for :math:`k \ge 3`. +The complete algorithm is as follows: -1. If :math:`\kappa < 3`, sample :math:`\mu` from the Klein-Nishina cross - section using Kahn's rejection method. Otherwise, use Koblinger's direct - method. +1. If :math:`k < 3`, sample :math:`\mu` from the Klein-Nishina cross section + using Kahn's rejection method. Otherwise, use Koblinger's direct method. 2. Calculate :math:`x` and :math:`\bar{x}` using :eq:`momentum-transfer` and :eq:`xmax`, respectively. @@ -243,7 +241,8 @@ scattering angle: .. math:: :label: pz - p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 - 2EE'\mu}}, + p_z = \frac{E - E' - EE'(1 - \mu)/(m_e c^2)}{-\alpha \sqrt{E^2 + E'^2 - + 2EE'\mu}}, where :math:`\alpha` is the fine structure constant. The maximum momentum transferred, :math:`p_{z,\text{max}}`, can be calculated from :eq:`pz` using @@ -269,7 +268,7 @@ The sampling algorithm is summarized below: :ref:`incoherent-sampling`. 2. Sample the electron subshell :math:`i` using the number of electrons per - shell as the PDF. + shell as the probability mass function. 3. Sample :math:`p_z` using :math:`J_i(p_z)` as the PDF. @@ -325,7 +324,7 @@ heavier elements. When simulating the photoelectric effect, the first step is to sample the electron shell. The shell :math:`i` where the ionization occurs can be -considered a discrete random variable with probability density function +considered a discrete random variable with probability mass function .. math:: :label: photoelectron-shell-pdf @@ -339,7 +338,7 @@ shell has been sampled, the energy of the photoelectron is calculated using :eq:`photoelectron-kinetic-energy`. To determine the direction of the photoelectron, we implement the method -described in [Kaltiaisenaho]_, which models the angular distribution of the +described in Kaltiaisenaho_, which models the angular distribution of the photoelectrons using the K-shell cross section derived by Sauter (K-shell electrons are the most tightly bound, and they contribute the most to :math:`\sigma_{\text{pe}}`). The non-relativistic Sauter distribution for @@ -372,10 +371,11 @@ where .. math:: :label: mu-pdf-factors - \psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 - \beta_{-}\mu_{-})^2}, \\ + \psi(\mu_{-}) &= \frac{(1 - \beta_{-}^2)(1 - \mu_{-}^2)}{(1 - + \beta_{-}\mu_{-})^2}, \\ g(\mu_{-}) &= \frac{1 - \beta_{-}^2}{2 (1 - \beta_{-}\mu_{-})^2}. -In the interval :math:`(-1, 1)`, :math:`g(\mu_{-})` is a normalized PDF and +In the interval :math:`[-1, 1]`, :math:`g(\mu_{-})` is a normalized PDF and :math:`\psi(\mu_{-})` satisfies the condition :math:`0 < \psi(\mu_{-}) < 1`. The following algorithm can now be used to sample :math:`\mu_{-}`: @@ -426,7 +426,7 @@ Accurately modeling the creation of electron-positron pair is important because the charged particles can go on to lose much of their energy as bremsstrahlung radiation, and the subsequent annihilation of the positron with an electron produces two additional photons. We sample the energy and direction of the -charged particles using a semiempirical model described in [Salvat]_. The +charged particles using a semiempirical model described in Salvat_. The Bethe-Heitler differential cross section, given by .. math:: @@ -441,9 +441,9 @@ constant, :math:`f_C` is the Coulomb correction function, :math:`\Phi_1` and :math:`\Phi_2` are screening functions, and :math:`\epsilon = (E_{-} + m_e c^2)/E` is the electron reduced energy (i.e., the fraction of the photon energy given to the electron). :math:`\epsilon` can take values between -:math:`\epsilon_{\text{min}} = \kappa^{-1}` (when the kinetic energy of the -electron is zero) and :math:`\epsilon_{\text{max}} = 1 - \kappa^{-1}` (when the -kinetic energy of the positron is zero). +:math:`\epsilon_{\text{min}} = k^{-1}` (when the kinetic energy of the electron +is zero) and :math:`\epsilon_{\text{max}} = 1 - k^{-1}` (when the kinetic +energy of the positron is zero). The Coulomb correction, given by @@ -478,7 +478,7 @@ where .. math:: :label: b - b = \frac{Rm_{e}c}{2\kappa\epsilon(1 - \epsilon)\hbar}. + b = \frac{Rm_{e}c}{2k\epsilon(1 - \epsilon)\hbar}. and :math:`R` is the screening radius. @@ -487,15 +487,15 @@ described above will not be accurate at low energies: the lower boundary of :math:`\epsilon` will be shifted above :math:`\epsilon_{\text{min}}` and the upper boundary of :math:`\epsilon` will be shifted below :math:`\epsilon_{\text{max}}`. To offset this behavior, a correcting factor -:math:`F_0(\kappa, Z)` is used: +:math:`F_0(k, Z)` is used: .. math:: :label: correcting-factor - F_0(\kappa, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/\kappa)^{1/2} \\ - &+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/\kappa) \\ - &- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/\kappa)^{3/2} \\ - &+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/\kappa)^{2}. + F_0(k, Z) =~& (0.1774 + 12.10\alpha Z - 11.18\alpha^{2}Z^{2})(2/k)^{1/2} \\ + &+ (8.523 + 73.26\alpha Z - 44.41\alpha^{2}Z^{2})(2/k) \\ + &- (13.52 + 121.1\alpha Z - 96.41\alpha^{2}Z^{2})(2/k)^{3/2} \\ + &+ (8.946 + 62.05\alpha Z - 63.41\alpha^{2}Z^{2})(2/k)^{2}. To aid sampling, the differential cross section used to sample :math:`\epsilon` (minus the normalization constant) can now be expressed in the form @@ -512,23 +512,23 @@ where .. math:: :label: u - u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^2 \phi_1(1/2), \\ + u_1 &= \frac{2}{3} \left(\frac{1}{2} - \frac{1}{k}\right)^2 \phi_1(1/2), \\ u_2 &= \phi_2(1/2), .. math:: :label: phi - \phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(\kappa, Z), \\ - \phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(\kappa, Z), + \phi_1(\epsilon) &= \frac{1}{2}(3\Phi_1 - \Phi_2) - 4f_{C}(Z) + F_0(k, Z), \\ + \phi_2(\epsilon) &= \frac{1}{4}(3\Phi_1 + \Phi_2) - 4f_{C}(Z) + F_0(k, Z), and .. math:: :label: pi - \pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^{-3} + \pi_1(\epsilon) &= \frac{3}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-3} \left(\frac{1}{2} - \epsilon\right)^2, \\ - \pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{\kappa}\right)^{-1}. + \pi_2(\epsilon) &= \frac{1}{2} \left(\frac{1}{2} - \frac{1}{k}\right)^{-1}. The functions in :eq:`phi` are non-negative and maximum at :math:`\epsilon = 1/2`. In the interval :math:`(\epsilon_{\text{min}}, \epsilon_{\text{max}})`, @@ -545,10 +545,10 @@ sample the reduced electron energy :math:`\epsilon`: .. math:: - \epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{\kappa}\right) + \epsilon &= \frac{1}{2} + \left(\frac{1}{2} - \frac{1}{k}\right) (2\xi_1 - 1)^{1/3} ~~~~&\text{if}~~ i = 1 \\ - \epsilon &= \frac{1}{\kappa} + \left(\frac{1}{2} - - \frac{1}{\kappa}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2. + \epsilon &= \frac{1}{k} + \left(\frac{1}{2} - + \frac{1}{k}\right) 2\xi_1 ~~~~&\text{if}~~ i = 2. 3. If :math:`\xi_2 \le \phi_i(\epsilon)/\phi_i(1/2)`, accept :math:`\epsilon`. Otherwise, repeat the sampling from step 1. @@ -563,13 +563,14 @@ distribution, .. math:: :label: sauter–gluckstern–hull - \frac{d\sigma_{pp}}{d\Omega_{\pm}} = C(1 - \beta_{\pm}\mu_{\pm})^{-2}, + p(\mu_{\pm}) = C(1 - \beta_{\pm}\mu_{\pm})^{-2}, -where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the ratio of -the velocity of the charged particle to the speed of light given in :eq:`beta-2`. +where :math:`C` is a normalization constant and :math:`\beta_{\pm}` is the +ratio of the velocity of the charged particle to the speed of light given in +:eq:`beta-2`. -The inverse transform method is used to sample :math:`\mu_{-}` and :math:`\mu_{+}` -from :eq:`sauter–gluckstern–hull`, using the sampling formula +The inverse transform method is used to sample :math:`\mu_{-}` and +:math:`\mu_{+}` from :eq:`sauter–gluckstern–hull`, using the sampling formula .. math:: :label: sample-mu @@ -638,7 +639,7 @@ and Auger electrons: photon assuming it is from a captured free electron and terminate. 2. Sample a transition using the transition probabilities for the vacancy - shell as the PDF. + shell as the probability mass function. 3. Create either a fluorescence photon or Auger electron, sampling the direction of the particle isotropically. @@ -688,18 +689,14 @@ expressed in the form \chi(Z, T, \kappa), where :math:`\kappa = E/T` is the reduced photon energy and :math:`\chi(Z, T, -\kappa)` is the scaled bremsstrahlung cross section, - -.. math:: - :label: scaled-bremsstrahlung-dcs - - \chi(Z, T, \kappa) = \frac{\beta^2}{Z^2} E \frac{d\sigma_{\text{br}}}{dE}. +\kappa)` is the scaled bremsstrahlung cross section, which is experimentally +measured. Because electrons are attracted to atomic nuclei whereas positrons are repulsed, the cross section for positrons is smaller, though it approaches that of electrons in the high energy limit. To obtain the positron cross section, we multiply :eq:`bremsstrahlung-dcs` by the :math:`\kappa`-independent factor used -in [Salvat]_, +in Salvat_, .. math:: :label: positron-factor @@ -777,9 +774,9 @@ transport electrons. However, the bremsstrahlung emitted from high energy electrons and positrons can travel far from the interaction site. Thus, even without a full electron transport mode it is necessary to model bremsstrahlung. We use a thick-target bremsstrahlung (TTB) approximation based on the models in -[Salvat]_ and [Kaltiaisenaho]_ for generating bremsstrahlung photons, which -assumes the charged particle loses all its energy in a single homogeneous -material region. +Salvat_ and Kaltiaisenaho_ for generating bremsstrahlung photons, which assumes +the charged particle loses all its energy in a single homogeneous material +region. To model bremsstrahlung using the TTB approximation, we need to know the number of photons emitted by the charged particle and the energy distribution of the @@ -906,11 +903,6 @@ angles. .. _LA-UR-04-0488: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-ur-04-0488.pdf - .. rubric:: References +.. _Kaltiaisenaho: https://aaltodoc.aalto.fi/bitstream/handle/123456789/21004/master_Kaltiaisenaho_Toni_2016.pdf -.. [Kaltiaisenaho] T. Kaltiaisenaho, "Implementing a photon physics model in - Serpent 2." M.Sc. Thesis, Aalto University, 2016. - -.. [Salvat] F. Salvat, J. M. Fernández-Varea, and J. Sempau, "PENELOPE-2011: - A Code System for Monte Carlo Simulation of Electron and Photon Transport," - OECD-NEA, Issy-les-Moulineaux, France (2011). +.. _Salvat: http://www.oecd-nea.org/globalsearch/download.php?doc=77434 From 08a8d3288f949bfbadc12549296fe745eda33daa Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 5 Oct 2018 21:07:44 -0400 Subject: [PATCH 15/38] Fix HDF5/xtensor mismatched free/delete error --- include/openmc/hdf5_interface.h | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 468c25980c..ea4284cb82 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -162,13 +162,14 @@ void read_attribute(hid_t obj_id, const char* name, xt::xarray& arr) std::size_t size = 1; for (const auto x : shape) size *= x; - T* buffer = new T[size]; + std::vector buffer; + buffer.resize(size); // Read data from attribute - read_attr(obj_id, name, H5TypeMap::type_id, buffer); + read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); // Adapt array into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape); + arr = xt::adapt(buffer, shape); } // overload for std::string @@ -244,13 +245,14 @@ void read_dataset(hid_t dset, xt::xarray& arr, bool indep=false) std::size_t size = 1; for (const auto x : shape) size *= x; - T* buffer = new T[size]; + std::vector buffer; + buffer.resize(size); // Read data from attribute - read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); // Adapt into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape); + arr = xt::adapt(buffer, shape); } template @@ -273,13 +275,14 @@ void read_dataset_as_shape(hid_t obj_id, const char* name, std::size_t size = 1; for (const auto x : arr.shape()) size *= x; - T* buffer = new T[size]; + std::vector buffer; + buffer.resize(size); // Read data from attribute - read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); // Adapt into xarray - arr = xt::adapt(buffer, size, xt::acquire_ownership(), arr.shape()); + arr = xt::adapt(buffer, arr.shape()); close_dataset(dset); } From 0792781dc577afe8c3ca41d3bdfa21e9a3bba2da Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 5 Oct 2018 22:49:29 -0400 Subject: [PATCH 16/38] Fix indexing below lower bound due to bin. search --- include/openmc/search.h | 1 + 1 file changed, 1 insertion(+) diff --git a/include/openmc/search.h b/include/openmc/search.h index 81370a3de0..c91ad18fe1 100644 --- a/include/openmc/search.h +++ b/include/openmc/search.h @@ -14,6 +14,7 @@ template typename std::iterator_traits::difference_type lower_bound_index(It first, It last, const T& value) { + if (*first == value) return 0; It index = std::lower_bound(first, last, value) - 1; return (index == last) ? -1 : index - first; } From a4b627a9e9fa1733a847971b37c2065e3b26697d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 10:37:53 -0500 Subject: [PATCH 17/38] Removing silo option from the voxel to vtk script. --- scripts/openmc-voxel-to-silovtk | 79 +++++++++++---------------------- 1 file changed, 26 insertions(+), 53 deletions(-) diff --git a/scripts/openmc-voxel-to-silovtk b/scripts/openmc-voxel-to-silovtk index e0cfc0a5ba..6e3e525e82 100755 --- a/scripts/openmc-voxel-to-silovtk +++ b/scripts/openmc-voxel-to-silovtk @@ -6,6 +6,7 @@ from argparse import ArgumentParser import numpy as np import h5py +import vtk def main(): @@ -14,8 +15,6 @@ def main(): parser.add_argument('voxel_file', help='Path to voxel file') parser.add_argument('-o', '--output', action='store', default='plot', help='Path to output SILO or VTK file.') - parser.add_argument('-s', '--silo', action='store_true', - default=False, help='Flag to convert to SILO instead of VTK.') args = parser.parse_args() # Read data from voxel file @@ -27,59 +26,33 @@ def main(): nx, ny, nz = dimension upper_right = lower_left + width*dimension + + grid = vtk.vtkImageData() + grid.SetDimensions(nx+1, ny+1, nz+1) + grid.SetOrigin(*lower_left) + grid.SetSpacing(*width) + + data = vtk.vtkDoubleArray() + data.SetName("id") + data.SetNumberOfTuples(nx*ny*nz) + for x in range(nx): + sys.stdout.write(" {}%\r".format(int(x/nx*100))) + sys.stdout.flush() + for y in range(ny): + for z in range(nz): + i = z*nx*ny + y*nx + x + data.SetValue(i, voxel_data[x, y, z]) + grid.GetCellData().AddArray(data) - if not args.silo: - import vtk - - grid = vtk.vtkImageData() - grid.SetDimensions(nx+1, ny+1, nz+1) - grid.SetOrigin(*lower_left) - grid.SetSpacing(*width) - - data = vtk.vtkDoubleArray() - data.SetName("id") - data.SetNumberOfTuples(nx*ny*nz) - for x in range(nx): - sys.stdout.write(" {}%\r".format(int(x/nx*100))) - sys.stdout.flush() - for y in range(ny): - for z in range(nz): - i = z*nx*ny + y*nx + x - data.SetValue(i, voxel_data[x, y, z]) - grid.GetCellData().AddArray(data) - - writer = vtk.vtkXMLImageDataWriter() - if vtk.vtkVersion.GetVTKMajorVersion() > 5: - writer.SetInputData(grid) - else: - writer.SetInput(grid) - if not args.output.endswith(".vti"): - args.output += ".vti" - writer.SetFileName(args.output) - writer.Write() - + writer = vtk.vtkXMLImageDataWriter() + if vtk.vtkVersion.GetVTKMajorVersion() > 5: + writer.SetInputData(grid) else: - import silomesh - - if not args.output.endswith(".silo"): - args.output += ".silo" - silomesh.init_silo(args.output) - meshparams = list(map(int, dimension)) + list(map(float, lower_left)) + \ - list(map(float, upper_right)) - silomesh.init_mesh('plot', *meshparams) - silomesh.init_var("id") - for x in range(nx): - sys.stdout.write(" {}%\r".format(int(x/nx*100))) - sys.stdout.flush() - for y in range(ny): - for z in range(nz): - silomesh.set_value(float(voxel_data[x, y, z]), - x + 1, y + 1, z + 1) - print() - silomesh.finalize_var() - silomesh.finalize_mesh() - silomesh.finalize_silo() - + writer.SetInput(grid) + if not args.output.endswith(".vti"): + args.output += ".vti" + writer.SetFileName(args.output) + writer.Write() if __name__ == '__main__': main() From a6d70cf88b9466db89169c9b41e878f449053848 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 10:38:19 -0500 Subject: [PATCH 18/38] Renaming voxel conversion script. --- scripts/{openmc-voxel-to-silovtk => openmc-voxel-to-vtk} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename scripts/{openmc-voxel-to-silovtk => openmc-voxel-to-vtk} (100%) diff --git a/scripts/openmc-voxel-to-silovtk b/scripts/openmc-voxel-to-vtk similarity index 100% rename from scripts/openmc-voxel-to-silovtk rename to scripts/openmc-voxel-to-vtk From c67d1e3d5950b3c3ca8639c037e4f0d3444e9637 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 12:25:06 -0500 Subject: [PATCH 19/38] Adding VTK image format to ignore. --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 7eb8f286ac..3832322966 100644 --- a/.gitignore +++ b/.gitignore @@ -69,6 +69,7 @@ scripts/G4EMLOW*/ # Images *.ppm *.voxel +*.vti # PyCharm project configuration files .idea From e9672ca1dd7dd6694849817afcd4bec93534684e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 12:25:23 -0500 Subject: [PATCH 20/38] Adding call to to the voxel conversion script. --- tests/regression_tests/plot/test.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index bfaef019c1..2fec279d22 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -1,7 +1,7 @@ import glob import hashlib import os - +from subprocess import call import h5py import openmc @@ -52,6 +52,10 @@ class PlotTestHarness(TestHarness): sha512.update(outstr) outstr = sha512.hexdigest() + # test the voxel to vtk conversion script + call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5')) + return outstr From 271bd799b2c5ef9dd72ac77cf0dcd948d82b6bf1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 12:25:44 -0500 Subject: [PATCH 21/38] Decreasing pixel size of the voxel image plot. --- tests/regression_tests/plot/plots.xml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/plot/plots.xml b/tests/regression_tests/plot/plots.xml index ecfe69125b..2cec871ea3 100644 --- a/tests/regression_tests/plot/plots.xml +++ b/tests/regression_tests/plot/plots.xml @@ -24,7 +24,7 @@ - 100 100 10 + 50 50 10 0. 0. 0. 20 20 10 From 4c6798b3e06ab61a7a45cd9f98a487edc2ef48ef Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 12:25:59 -0500 Subject: [PATCH 22/38] Updating results. --- tests/regression_tests/plot/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index ba8a49226a..f72a4ecd76 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -01ecda0f3820a49c8a41d8dc47d1e5c58767a04301621c2437231fcc04401ddea47b67d0529ca56a32d4d97b4f1416a2e0b6120d3bdc87d74a7e9889758a8808 \ No newline at end of file +3337721f8f6d3777dd36c404ed117fe5bc07b1d9691a0c91131d48e7e73e21191e72401cf9bb621e3250ca51345514f90ea612e0fb6236bec33ce78655988dec \ No newline at end of file From 755d0a5a5fa1a4e11564b495e159acf7b6dd474b Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 4 Oct 2018 12:43:32 -0500 Subject: [PATCH 23/38] Adding vtk install to tests. --- tests/regression_tests/plot/test.py | 8 ++++++-- tools/ci/travis-install.sh | 5 +++++ 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index 2fec279d22..e568fe86a9 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -53,8 +53,12 @@ class PlotTestHarness(TestHarness): outstr = sha512.hexdigest() # test the voxel to vtk conversion script - call(['../../../scripts/openmc-voxel-to-vtk'] + - glob.glob('plot_4.h5')) + try: + import vtk + call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5')) + except: + pass return outstr diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index c10c222793..db8d1ca5a2 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -26,5 +26,10 @@ python tools/ci/travis-install.py # Install Python API in editable mode pip install -e .[test] +a=$(dpkg --compare-versions $(python --version | cut -d" " -f 2) lt 3.7.0) +if [ $a -eq 0 ]; then + pip install -e .[vtk] +fi + # For uploading to coveralls pip install coveralls From b8a845073bd164ea245e632643431ccef90a5ed6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 7 Oct 2018 17:21:48 -0500 Subject: [PATCH 24/38] Using travis environment variable to check Python version. --- tools/ci/travis-install.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index db8d1ca5a2..cc707d861b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -26,8 +26,8 @@ python tools/ci/travis-install.py # Install Python API in editable mode pip install -e .[test] -a=$(dpkg --compare-versions $(python --version | cut -d" " -f 2) lt 3.7.0) -if [ $a -eq 0 ]; then +# conditionally install vtk +if [ $TRAVIS_PYTHON_VERSION -ne 3.7 ]; then pip install -e .[vtk] fi From a5b179525ccbdece30f80910bf15c8d521daa037 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 7 Oct 2018 18:35:40 -0500 Subject: [PATCH 25/38] Removing optional silo dependency. --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index cbf564480f..c0168d8805 100755 --- a/setup.py +++ b/setup.py @@ -64,7 +64,7 @@ kwargs = { # Optional dependencies 'extras_require': { 'test': ['pytest', 'pytest-cov'], - 'vtk': ['vtk', 'silomesh'], + 'vtk': ['vtk'], }, } From ad3f4ccf0c49f670f3256f8316c4648a3c1e3e94 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 7 Oct 2018 19:04:07 -0500 Subject: [PATCH 26/38] Adding plot.vti to list of expected files if vtk is present. --- tests/regression_tests/plot/test.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index e568fe86a9..9a6fb3baf2 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -18,6 +18,13 @@ class PlotTestHarness(TestHarness): def _run_openmc(self): openmc.plot_geometry(openmc_exec=config['exe']) + try: + import vtk + call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5')) + except: + pass + def _test_output_created(self): """Make sure *.ppm has been created.""" for fname in self._plot_names: @@ -52,18 +59,16 @@ class PlotTestHarness(TestHarness): sha512.update(outstr) outstr = sha512.hexdigest() - # test the voxel to vtk conversion script - try: - import vtk - call(['../../../scripts/openmc-voxel-to-vtk'] + - glob.glob('plot_4.h5')) - except: - pass - return outstr - def test_plot(): - harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', - 'plot_4.h5')) + expected_plots = ['plot_1.ppm', 'plot_2.ppm', + 'plot_3.ppm', 'plot_4.h5'] + try: + import vtk + expected_plots.append('plot.vti') + except: + pass + + harness = PlotTestHarness(expected_plots) harness.main() From 76aa15c6e14246a3fa30eef012e0110a9c128bf6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Oct 2018 10:00:53 -0500 Subject: [PATCH 27/38] Adding comments to plot test for clarity. --- tests/regression_tests/plot/test.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index 9a6fb3baf2..a16b3c43cd 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -18,6 +18,8 @@ class PlotTestHarness(TestHarness): def _run_openmc(self): openmc.plot_geometry(openmc_exec=config['exe']) + # TEMP: this should always be checked once vtk is added + # to Python3.7 try: import vtk call(['../../../scripts/openmc-voxel-to-vtk'] + @@ -64,6 +66,8 @@ class PlotTestHarness(TestHarness): def test_plot(): expected_plots = ['plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', 'plot_4.h5'] + # TEMP: this should always be checked once vtk is added + # to Python3.7 try: import vtk expected_plots.append('plot.vti') From eb14dbd9cdce11dd9341616349066ef5fae6826b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Oct 2018 13:19:47 -0400 Subject: [PATCH 28/38] Set vector size with constructor --- include/openmc/hdf5_interface.h | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index ea4284cb82..f0b3a245c7 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -162,8 +162,7 @@ void read_attribute(hid_t obj_id, const char* name, xt::xarray& arr) std::size_t size = 1; for (const auto x : shape) size *= x; - std::vector buffer; - buffer.resize(size); + std::vector buffer(size); // Read data from attribute read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); @@ -245,8 +244,7 @@ void read_dataset(hid_t dset, xt::xarray& arr, bool indep=false) std::size_t size = 1; for (const auto x : shape) size *= x; - std::vector buffer; - buffer.resize(size); + std::vector buffer(size); // Read data from attribute read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); @@ -275,8 +273,7 @@ void read_dataset_as_shape(hid_t obj_id, const char* name, std::size_t size = 1; for (const auto x : arr.shape()) size *= x; - std::vector buffer; - buffer.resize(size); + std::vector buffer(size); // Read data from attribute read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); From 2ef7ea8fe06816e92ee10eb893d1e5c93379f06a Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Oct 2018 14:09:14 -0500 Subject: [PATCH 29/38] PEP8 updates. --- tests/regression_tests/plot/test.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index a16b3c43cd..90e04f913f 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -2,6 +2,7 @@ import glob import hashlib import os from subprocess import call + import h5py import openmc @@ -63,9 +64,9 @@ class PlotTestHarness(TestHarness): return outstr + def test_plot(): - expected_plots = ['plot_1.ppm', 'plot_2.ppm', - 'plot_3.ppm', 'plot_4.h5'] + expected_plots = ['plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', 'plot_4.h5'] # TEMP: this should always be checked once vtk is added # to Python3.7 try: From cba6c5ad63b16303915cf4acb17de14ebfb589ed Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Oct 2018 14:45:15 -0500 Subject: [PATCH 30/38] Creating a separate test for the voxel script. --- tests/regression_tests/plot/plots.xml | 2 +- tests/regression_tests/plot/results_true.dat | 2 +- tests/regression_tests/plot/test.py | 22 +------ tests/regression_tests/plot_voxel/__init__.py | 0 .../regression_tests/plot_voxel/geometry.xml | 13 ++++ .../regression_tests/plot_voxel/materials.xml | 19 ++++++ tests/regression_tests/plot_voxel/plots.xml | 10 +++ .../plot_voxel/results_true.dat | 1 + .../regression_tests/plot_voxel/settings.xml | 13 ++++ tests/regression_tests/plot_voxel/test.py | 64 +++++++++++++++++++ 10 files changed, 124 insertions(+), 22 deletions(-) create mode 100644 tests/regression_tests/plot_voxel/__init__.py create mode 100644 tests/regression_tests/plot_voxel/geometry.xml create mode 100644 tests/regression_tests/plot_voxel/materials.xml create mode 100644 tests/regression_tests/plot_voxel/plots.xml create mode 100644 tests/regression_tests/plot_voxel/results_true.dat create mode 100644 tests/regression_tests/plot_voxel/settings.xml create mode 100644 tests/regression_tests/plot_voxel/test.py diff --git a/tests/regression_tests/plot/plots.xml b/tests/regression_tests/plot/plots.xml index 2cec871ea3..ecfe69125b 100644 --- a/tests/regression_tests/plot/plots.xml +++ b/tests/regression_tests/plot/plots.xml @@ -24,7 +24,7 @@ - 50 50 10 + 100 100 10 0. 0. 0. 20 20 10 diff --git a/tests/regression_tests/plot/results_true.dat b/tests/regression_tests/plot/results_true.dat index f72a4ecd76..ba8a49226a 100644 --- a/tests/regression_tests/plot/results_true.dat +++ b/tests/regression_tests/plot/results_true.dat @@ -1 +1 @@ -3337721f8f6d3777dd36c404ed117fe5bc07b1d9691a0c91131d48e7e73e21191e72401cf9bb621e3250ca51345514f90ea612e0fb6236bec33ce78655988dec \ No newline at end of file +01ecda0f3820a49c8a41d8dc47d1e5c58767a04301621c2437231fcc04401ddea47b67d0529ca56a32d4d97b4f1416a2e0b6120d3bdc87d74a7e9889758a8808 \ No newline at end of file diff --git a/tests/regression_tests/plot/test.py b/tests/regression_tests/plot/test.py index 90e04f913f..bfaef019c1 100644 --- a/tests/regression_tests/plot/test.py +++ b/tests/regression_tests/plot/test.py @@ -1,7 +1,6 @@ import glob import hashlib import os -from subprocess import call import h5py import openmc @@ -19,15 +18,6 @@ class PlotTestHarness(TestHarness): def _run_openmc(self): openmc.plot_geometry(openmc_exec=config['exe']) - # TEMP: this should always be checked once vtk is added - # to Python3.7 - try: - import vtk - call(['../../../scripts/openmc-voxel-to-vtk'] + - glob.glob('plot_4.h5')) - except: - pass - def _test_output_created(self): """Make sure *.ppm has been created.""" for fname in self._plot_names: @@ -66,14 +56,6 @@ class PlotTestHarness(TestHarness): def test_plot(): - expected_plots = ['plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', 'plot_4.h5'] - # TEMP: this should always be checked once vtk is added - # to Python3.7 - try: - import vtk - expected_plots.append('plot.vti') - except: - pass - - harness = PlotTestHarness(expected_plots) + harness = PlotTestHarness(('plot_1.ppm', 'plot_2.ppm', 'plot_3.ppm', + 'plot_4.h5')) harness.main() diff --git a/tests/regression_tests/plot_voxel/__init__.py b/tests/regression_tests/plot_voxel/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/tests/regression_tests/plot_voxel/geometry.xml b/tests/regression_tests/plot_voxel/geometry.xml new file mode 100644 index 0000000000..83619d9f78 --- /dev/null +++ b/tests/regression_tests/plot_voxel/geometry.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_voxel/materials.xml b/tests/regression_tests/plot_voxel/materials.xml new file mode 100644 index 0000000000..90b3542675 --- /dev/null +++ b/tests/regression_tests/plot_voxel/materials.xml @@ -0,0 +1,19 @@ + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/plot_voxel/plots.xml b/tests/regression_tests/plot_voxel/plots.xml new file mode 100644 index 0000000000..833329b427 --- /dev/null +++ b/tests/regression_tests/plot_voxel/plots.xml @@ -0,0 +1,10 @@ + + + + + 50 50 10 + 0. 0. 0. + 20 20 10 + + + diff --git a/tests/regression_tests/plot_voxel/results_true.dat b/tests/regression_tests/plot_voxel/results_true.dat new file mode 100644 index 0000000000..41193c9bc3 --- /dev/null +++ b/tests/regression_tests/plot_voxel/results_true.dat @@ -0,0 +1 @@ +20a0c3598bb698efa4bba65c5e55d71f4385e5b1bcf0067964e45001c12f5c0a729935a96409c760dbd3ec439ae6c676fce77d6e578b41f8f3e0ac68c9277acb \ No newline at end of file diff --git a/tests/regression_tests/plot_voxel/settings.xml b/tests/regression_tests/plot_voxel/settings.xml new file mode 100644 index 0000000000..adf256d2d4 --- /dev/null +++ b/tests/regression_tests/plot_voxel/settings.xml @@ -0,0 +1,13 @@ + + + + plot + + + 5 4 3 + -10 -10 -10 + 10 10 10 + + 1 + + diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py new file mode 100644 index 0000000000..1e57e557d6 --- /dev/null +++ b/tests/regression_tests/plot_voxel/test.py @@ -0,0 +1,64 @@ +import glob +import hashlib +import os +from subprocess import run +import importlib +import h5py +import openmc +import pytest + +from tests.testing_harness import TestHarness +from tests.regression_tests import config + + +pytestmark = pytest.mark.skipif( + importlib.util.find_spec('vtk') is None, + reason="vtk module is not installed.") +class PlotVoxelTestHarness(TestHarness): + """Specialized TestHarness for running OpenMC voxel plot tests.""" + def __init__(self, plot_names): + super().__init__(None) + self._plot_names = plot_names + + def _run_openmc(self): + openmc.plot_geometry(openmc_exec=config['exe']) + + run(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5'), check=True) + + def _test_output_created(self): + """Make sure *.ppm has been created.""" + for fname in self._plot_names: + assert os.path.exists(fname), 'Plot output file does not exist.' + + def _cleanup(self): + super()._cleanup() + for fname in self._plot_names: + if os.path.exists(fname): + os.remove(fname) + + def _get_results(self): + """Return a string hash of the plot files.""" + outstr = bytes() + + for fname in self._plot_names: + if fname.endswith('.h5'): + # Add voxel data to results + with h5py.File(fname, 'r') as fh: + outstr += fh.attrs['filetype'] + outstr += fh.attrs['num_voxels'].tostring() + outstr += fh.attrs['lower_left'].tostring() + outstr += fh.attrs['voxel_width'].tostring() + outstr += fh['data'].value.tostring() + + # Hash the information and return. + sha512 = hashlib.sha512() + sha512.update(outstr) + outstr = sha512.hexdigest() + + return outstr + + +def test_plot(): + harness = PlotVoxelTestHarness(('plot_4.h5', 'plot.vti')) + harness.main() From fd92aac49c54b1014426c15883234fbc27321e97 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Oct 2018 15:23:14 -0500 Subject: [PATCH 31/38] Removing mention of SILO in help string. --- scripts/openmc-voxel-to-vtk | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/openmc-voxel-to-vtk b/scripts/openmc-voxel-to-vtk index 6e3e525e82..c0d8e9a147 100755 --- a/scripts/openmc-voxel-to-vtk +++ b/scripts/openmc-voxel-to-vtk @@ -14,7 +14,7 @@ def main(): parser = ArgumentParser() parser.add_argument('voxel_file', help='Path to voxel file') parser.add_argument('-o', '--output', action='store', - default='plot', help='Path to output SILO or VTK file.') + default='plot', help='Path to output VTK file.') args = parser.parse_args() # Read data from voxel file From 291aa61f12ba79e9eb3f2a905ee15f5f2041c314 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 8 Oct 2018 15:33:36 -0500 Subject: [PATCH 32/38] Changing to check_call to support Python < 3.5 --- tests/regression_tests/plot_voxel/test.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py index 1e57e557d6..f6af8fa116 100644 --- a/tests/regression_tests/plot_voxel/test.py +++ b/tests/regression_tests/plot_voxel/test.py @@ -1,8 +1,9 @@ import glob import hashlib import os -from subprocess import run +from subprocess import check_call import importlib + import h5py import openmc import pytest @@ -23,8 +24,8 @@ class PlotVoxelTestHarness(TestHarness): def _run_openmc(self): openmc.plot_geometry(openmc_exec=config['exe']) - run(['../../../scripts/openmc-voxel-to-vtk'] + - glob.glob('plot_4.h5'), check=True) + check_call(['../../../scripts/openmc-voxel-to-vtk'] + + glob.glob('plot_4.h5')) def _test_output_created(self): """Make sure *.ppm has been created.""" From 66ba55c3f67ccc73dd1c0bee0eea3983d2df9bd6 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 8 Oct 2018 21:52:12 -0400 Subject: [PATCH 33/38] Fix Fortran character allocation errors for HDF5 --- src/hdf5_interface.F90 | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 40f56d88ef..0c7bdb92d5 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -675,6 +675,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -698,6 +699,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -720,6 +722,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -742,6 +745,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -764,6 +768,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_double_c(obj_id, c_loc(name_), buffer, indep_) else @@ -888,6 +893,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -911,6 +917,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -933,6 +940,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -955,6 +963,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -977,6 +986,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_int_c(obj_id, c_loc(name_), buffer, indep_) else @@ -1025,6 +1035,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_llong_c(obj_id, c_loc(name_), buffer_, indep_) else @@ -1425,6 +1436,7 @@ contains ! If 'name' argument is passed, obj_id is interpreted to be a group and ! 'name' is the name of the dataset we should read from if (present(name)) then + allocate(name_(len_trim(name) + 1)) name_ = to_c_string(name) call read_complex_c(obj_id, c_loc(name_), buffer, indep_) else From e2ed9c7553d6dcc1daafa1ec603c5956a566917f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 24 Sep 2018 13:29:11 -0500 Subject: [PATCH 34/38] Fix various bugs in HDF5 interface / MGXS code (found by valgrind) --- include/openmc/hdf5_interface.h | 4 ++-- src/hdf5_interface.cpp | 10 +++++++--- src/mgxs_interface.cpp | 2 +- 3 files changed, 10 insertions(+), 6 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index f0b3a245c7..a39561d23c 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -193,13 +193,13 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) // Allocate a C char array to get strings auto n = attribute_typesize(obj_id, name); - char buffer[m][n+1]; + char buffer[m][n]; // Read char data in attribute read_attr_string(obj_id, name, n, buffer[0]); for (int i = 0; i < m; ++i) { - vec.emplace_back(&buffer[i][0]); + vec.emplace_back(&buffer[i][0], n); } } diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index cecf283ac1..d72f7d424d 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -392,7 +392,7 @@ object_name(hid_t obj_id) // Read and return name H5Iget_name(obj_id, buffer, size); - return {buffer, size}; + return buffer; } @@ -439,7 +439,9 @@ read_attr_string(hid_t obj_id, const char* name, size_t slen, char* buffer) { // Create datatype for a string hid_t datatype = H5Tcopy(H5T_C_S1); - H5Tset_size(datatype, slen + 1); + H5Tset_size(datatype, slen); + // numpy uses null-padding when writing fixed-length strings + H5Tset_strpad(datatype, H5T_STR_NULLPAD); // Read data into buffer read_attr(obj_id, name, datatype, buffer); @@ -503,7 +505,9 @@ read_string(hid_t obj_id, const char* name, size_t slen, char* buffer, bool inde { // Create datatype for a string hid_t datatype = H5Tcopy(H5T_C_S1); - H5Tset_size(datatype, slen + 1); + H5Tset_size(datatype, slen); + // numpy uses null-padding when writing fixed-length strings + H5Tset_strpad(datatype, H5T_STR_NULLPAD); // Read data into buffer read_dataset(obj_id, name, datatype, buffer, indep); diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 56f3392ff2..0bb464da07 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -208,7 +208,7 @@ void get_name_c(int index, int name_len, char* name) { // First blank out our input string - std::string str(name_len, ' '); + std::string str(name_len - 1, ' '); std::strcpy(name, str.c_str()); // Now get the data and copy to the C-string From 2d1793b7f96d32a1ce697218e8fb9528eadb577c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 8 Oct 2018 15:03:20 -0500 Subject: [PATCH 35/38] Make sure buffer for array of strings doesn't contain uninitialized values --- src/hdf5_interface.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 0c7bdb92d5..842ce48f37 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -1131,6 +1131,7 @@ contains do i = 0, dims(1) - 1 n = len_trim(buffer(i+1)) + 1 buffer_(i*m+1 : i*m+n) = to_c_string(buffer(i+1)) + if (n < m) buffer_(i*m+n : i*m+m) = C_NULL_CHAR end do call write_string_c(group_id, 1, dims, m, to_c_string(name), & From a10c357b013067ed56cac6692b008c406dc2fb49 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 9 Oct 2018 08:36:05 -0500 Subject: [PATCH 36/38] Updating test function name. Updating skip condition for vtk import. --- tests/regression_tests/plot_voxel/test.py | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/tests/regression_tests/plot_voxel/test.py b/tests/regression_tests/plot_voxel/test.py index f6af8fa116..d2767e2f18 100644 --- a/tests/regression_tests/plot_voxel/test.py +++ b/tests/regression_tests/plot_voxel/test.py @@ -2,7 +2,6 @@ import glob import hashlib import os from subprocess import check_call -import importlib import h5py import openmc @@ -12,9 +11,7 @@ from tests.testing_harness import TestHarness from tests.regression_tests import config -pytestmark = pytest.mark.skipif( - importlib.util.find_spec('vtk') is None, - reason="vtk module is not installed.") +vtk = pytest.importorskip('vtk') class PlotVoxelTestHarness(TestHarness): """Specialized TestHarness for running OpenMC voxel plot tests.""" def __init__(self, plot_names): @@ -60,6 +57,6 @@ class PlotVoxelTestHarness(TestHarness): return outstr -def test_plot(): +def test_plot_voxel(): harness = PlotVoxelTestHarness(('plot_4.h5', 'plot.vti')) harness.main() From 596e05c272078444b09a145cebe224058f350a16 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 9 Oct 2018 11:01:30 -0500 Subject: [PATCH 37/38] Fix string lengths when reading HFD5 attributes --- include/openmc/hdf5_interface.h | 4 +++- src/hdf5_interface.F90 | 4 ++-- src/initialize.cpp | 8 +++----- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index a39561d23c..fa64605513 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -1,9 +1,11 @@ #ifndef OPENMC_HDF5_INTERFACE_H #define OPENMC_HDF5_INTERFACE_H +#include // for min #include #include #include +#include // for strlen #include #include #include @@ -199,7 +201,7 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) read_attr_string(obj_id, name, n, buffer[0]); for (int i = 0; i < m; ++i) { - vec.emplace_back(&buffer[i][0], n); + vec.emplace_back(&buffer[i][0], std::min(strlen(buffer[i]), n)); } } diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 842ce48f37..8872a84f6a 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -1362,7 +1362,7 @@ contains ! Allocate a C char array to get strings n = attribute_typesize(obj_id, to_c_string(name)) m = size(buffer) - allocate(buffer_((n+1)*m)) + allocate(buffer_(n*m)) ! Read attribute call read_attr_string_c(obj_id, to_c_string(name), n, buffer_) @@ -1371,7 +1371,7 @@ contains do i = 1, m buffer(i) = '' do j = 1, n - k = (i-1)*(n+1) + j + k = (i-1)*n + j if (buffer_(k) == C_NULL_CHAR) exit buffer(i)(j:j) = buffer_(k) end do diff --git a/src/initialize.cpp b/src/initialize.cpp index 19fc82c859..0c26135dc9 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -115,10 +115,9 @@ parse_command_line(int argc, char* argv[]) // Check what type of file this is hid_t file_id = file_open(argv[i], 'r', true); - size_t len = attribute_typesize(file_id, "filetype"); - read_attr_string(file_id, "filetype", len, buffer); + std::string filetype; + read_attribute(file_id, "filetype", filetype); file_close(file_id); - std::string filetype {buffer}; // Set path and flag for type of run if (filetype == "statepoint") { @@ -141,8 +140,7 @@ parse_command_line(int argc, char* argv[]) // Check file type is a source file file_id = file_open(argv[i+1], 'r', true); - len = attribute_typesize(file_id, "filetype"); - read_attr_string(file_id, "filetype", len, buffer); + read_attribute(file_id, "filetype", filetype); file_close(file_id); if (filetype != "source") { std::string msg {"Second file after restart flag must be a source file"}; From 8be83ff9df8b9c209db5733a91698af9baf91dd6 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 10 Oct 2018 10:24:19 -0500 Subject: [PATCH 38/38] Correcting vtk install condition. --- tools/ci/travis-install.sh | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index cc707d861b..670f28f33b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -23,12 +23,12 @@ fi # Build and install OpenMC executable python tools/ci/travis-install.py -# Install Python API in editable mode -pip install -e .[test] - -# conditionally install vtk -if [ $TRAVIS_PYTHON_VERSION -ne 3.7 ]; then - pip install -e .[vtk] +if [[ $TRAVIS_PYTHON_VERSION == "3.7" ]]; then + # Install Python API in editable mode + pip install -e .[test] +else + # Conditionally install vtk + pip install -e .[test,vtk] fi # For uploading to coveralls