From 89f6931b267988e8969715ded808cf8477f013a9 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 25 Jul 2018 21:05:35 -0400 Subject: [PATCH 01/36] Fix k_eff calc for runs with reset in the middle --- src/api.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/api.F90 b/src/api.F90 index 8afc6061cd..e38b2d0b28 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -276,6 +276,10 @@ contains k_abs_tra = ZERO k_sum(:) = ZERO + ! Set the number of inactive batches (used to compute k_eff for the fission + ! bank). + n_inactive = current_batch + ! Clear active tally lists call active_analog_tallies % clear() call active_tracklength_tallies % clear() From c0d3e37e26761fc1c8c515d6743567cb7631dc71 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 30 Jul 2018 12:54:12 -0400 Subject: [PATCH 02/36] Add n_realizations_gen for robust keff calculation --- openmc/capi/tally.py | 2 +- src/api.F90 | 7 ++----- src/eigenvalue.F90 | 6 +++--- src/output.F90 | 4 ++-- src/simulation.F90 | 8 +++++++- src/state_point.F90 | 7 ++++--- src/tallies/tally.F90 | 4 ++-- src/tallies/tally_header.F90 | 5 +++-- 8 files changed, 24 insertions(+), 19 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index d4d70af718..cb647b8cee 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -116,7 +116,7 @@ def global_tallies(): def num_realizations(): """Number of realizations of global tallies.""" - return c_int32.in_dll(_dll, 'n_realizations').value + return c_int32.in_dll(_dll, 'n_realizations_batch').value class Tally(_FortranObjectWithID): diff --git a/src/api.F90 b/src/api.F90 index e38b2d0b28..cf298512c5 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -267,7 +267,8 @@ contains end if ! Reset global tallies - n_realizations = 0 + n_realizations_batch = 0 + n_realizations_gen = 0 if (allocated(global_tallies)) then global_tallies(:, :) = ZERO end if @@ -276,10 +277,6 @@ contains k_abs_tra = ZERO k_sum(:) = ZERO - ! Set the number of inactive batches (used to compute k_eff for the fission - ! bank). - n_inactive = current_batch - ! Clear active tally lists call active_analog_tallies % clear() call active_tracklength_tallies % clear() diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 7a11d1fb48..9b2cc8ff4b 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -380,7 +380,7 @@ contains ! Determine overall generation and number of active generations i = overall_generation() - n = i - n_inactive*gen_per_batch + n = n_realizations_gen if (n <= 0) then ! For inactive generations, use current generation k as estimate for next @@ -438,13 +438,13 @@ contains ! Make sure we have at least four realizations. Notice that at the end, ! there is a N-3 term in a denominator. - if (n_realizations <= 3) then + if (n_realizations_batch <= 3) then err = -1 return end if ! Initialize variables - n = real(n_realizations, 8) + n = real(n_realizations_batch, 8) ! Copy estimates of k-effective and its variance (not variance of the mean) kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n diff --git a/src/output.F90 b/src/output.F90 index 6197cdc11c..c8fa67f03c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -569,7 +569,7 @@ contains ! display header block for results call header("Results", 4) - n = n_realizations + n = n_realizations_batch if (confidence_intervals) then ! Calculate t-value for confidence intervals @@ -714,7 +714,7 @@ contains ! Calculate t-value for confidence intervals if (confidence_intervals) then alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) + t_value = t_percentile(ONE - alpha/TWO, n_realizations_batch - 1) else t_value = ONE end if diff --git a/src/simulation.F90 b/src/simulation.F90 index ffcfb8bab8..14b6ba688b 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -281,6 +281,9 @@ contains if (entropy_on) call shannon_entropy() ! Collect results and statistics + if (current_batch > n_inactive) then + n_realizations_gen = n_realizations_gen + 1 + end if call calculate_generation_keff() call calculate_average_keff() @@ -325,7 +328,7 @@ contains ! Reset global tally results if (current_batch <= n_inactive) then global_tallies(:,:) = ZERO - n_realizations = 0 + n_realizations_batch = 0 end if if (run_mode == MODE_EIGENVALUE) then @@ -373,6 +376,9 @@ contains if (run_mode == MODE_EIGENVALUE) then do current_gen = 1, gen_per_batch + if (current_batch > n_inactive) then + n_realizations_gen = n_realizations_gen + 1 + end if call calculate_average_keff() ! print out batch keff diff --git a/src/state_point.F90 b/src/state_point.F90 index ebb473d496..eef0cf370b 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -359,7 +359,7 @@ contains elseif (master) then ! Write number of global realizations - call write_dataset(file_id, "n_realizations", n_realizations) + call write_dataset(file_id, "n_realizations", n_realizations_batch) ! Write global tallies call write_dataset(file_id, "global_tallies", global_tallies) @@ -518,7 +518,7 @@ contains if (master) then ! Write number of realizations - call write_dataset(file_id, "n_realizations", n_realizations) + call write_dataset(file_id, "n_realizations", n_realizations_batch) ! Write number of global tallies call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) @@ -756,7 +756,8 @@ contains #endif ! Read number of realizations for global tallies - call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) + call read_dataset(n_realizations_batch, file_id, "n_realizations", & + indep=.true.) ! Read global tally data call read_dataset(global_tallies, file_id, "global_tallies") diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index ed0963a883..9761ef288a 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3784,9 +3784,9 @@ contains ! Increase number of realizations (only used for global tallies) if (reduce_tallies) then - n_realizations = n_realizations + 1 + n_realizations_batch = n_realizations_batch + 1 else - n_realizations = n_realizations + n_procs + n_realizations_batch = n_realizations_batch + n_procs end if ! Accumulate on master only unless run is not reduced then do it on all diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index ee44685924..3418df12a6 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -145,8 +145,9 @@ module tally_header type(VectorInt), public :: active_tallies type(VectorInt), public :: active_surface_tallies - ! Normalization for statistics - integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations + ! Normalization for batch-based and generation-based global tally statistics + integer(C_INT32_T), public, bind(C) :: n_realizations_batch = 0 + integer(C_INT32_T), public, bind(C) :: n_realizations_gen = 0 real(8), public :: total_weight ! total starting particle weight in realization contains From 5af37d4d5badf28d8864588cd66882b90e7585b5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 31 Jul 2018 07:39:53 -0500 Subject: [PATCH 03/36] Use n_realizations for keff calculation --- openmc/capi/tally.py | 2 +- src/api.F90 | 3 +-- src/eigenvalue.F90 | 10 +++++++--- src/output.F90 | 12 ++++++++---- src/simulation.F90 | 8 +------- src/state_point.F90 | 7 +++---- src/tallies/tally.F90 | 4 ++-- src/tallies/tally_header.F90 | 5 ++--- 8 files changed, 25 insertions(+), 26 deletions(-) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index cb647b8cee..d4d70af718 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -116,7 +116,7 @@ def global_tallies(): def num_realizations(): """Number of realizations of global tallies.""" - return c_int32.in_dll(_dll, 'n_realizations_batch').value + return c_int32.in_dll(_dll, 'n_realizations').value class Tally(_FortranObjectWithID): diff --git a/src/api.F90 b/src/api.F90 index cf298512c5..8afc6061cd 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -267,8 +267,7 @@ contains end if ! Reset global tallies - n_realizations_batch = 0 - n_realizations_gen = 0 + n_realizations = 0 if (allocated(global_tallies)) then global_tallies(:, :) = ZERO end if diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 9b2cc8ff4b..d47bc4dab8 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -380,7 +380,11 @@ contains ! Determine overall generation and number of active generations i = overall_generation() - n = n_realizations_gen + if (current_batch > n_inactive) then + n = gen_per_batch*n_realizations + current_gen + else + n = 0 + end if if (n <= 0) then ! For inactive generations, use current generation k as estimate for next @@ -438,13 +442,13 @@ contains ! Make sure we have at least four realizations. Notice that at the end, ! there is a N-3 term in a denominator. - if (n_realizations_batch <= 3) then + if (n_realizations <= 3) then err = -1 return end if ! Initialize variables - n = real(n_realizations_batch, 8) + n = real(n_realizations, 8) ! Copy estimates of k-effective and its variance (not variance of the mean) kv(1) = global_tallies(RESULT_SUM, K_COLLISION) / n diff --git a/src/output.F90 b/src/output.F90 index c8fa67f03c..128080e61c 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -324,7 +324,11 @@ contains ! Determine overall generation and number of active generations i = overall_generation() - n = i - n_inactive*gen_per_batch + if (current_batch > n_inactive) then + n = gen_per_batch*n_realizations + current_gen + else + n = 0 + end if ! write out information about batch and generation write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & @@ -357,7 +361,7 @@ contains ! Determine overall generation and number of active generations i = current_batch*gen_per_batch - n = i - n_inactive*gen_per_batch + n = n_realizations*gen_per_batch ! write out information batch and option independent output write(UNIT=OUTPUT_UNIT, FMT='(2X,A9)', ADVANCE='NO') & @@ -569,7 +573,7 @@ contains ! display header block for results call header("Results", 4) - n = n_realizations_batch + n = n_realizations if (confidence_intervals) then ! Calculate t-value for confidence intervals @@ -714,7 +718,7 @@ contains ! Calculate t-value for confidence intervals if (confidence_intervals) then alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations_batch - 1) + t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) else t_value = ONE end if diff --git a/src/simulation.F90 b/src/simulation.F90 index 14b6ba688b..ffcfb8bab8 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -281,9 +281,6 @@ contains if (entropy_on) call shannon_entropy() ! Collect results and statistics - if (current_batch > n_inactive) then - n_realizations_gen = n_realizations_gen + 1 - end if call calculate_generation_keff() call calculate_average_keff() @@ -328,7 +325,7 @@ contains ! Reset global tally results if (current_batch <= n_inactive) then global_tallies(:,:) = ZERO - n_realizations_batch = 0 + n_realizations = 0 end if if (run_mode == MODE_EIGENVALUE) then @@ -376,9 +373,6 @@ contains if (run_mode == MODE_EIGENVALUE) then do current_gen = 1, gen_per_batch - if (current_batch > n_inactive) then - n_realizations_gen = n_realizations_gen + 1 - end if call calculate_average_keff() ! print out batch keff diff --git a/src/state_point.F90 b/src/state_point.F90 index eef0cf370b..ebb473d496 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -359,7 +359,7 @@ contains elseif (master) then ! Write number of global realizations - call write_dataset(file_id, "n_realizations", n_realizations_batch) + call write_dataset(file_id, "n_realizations", n_realizations) ! Write global tallies call write_dataset(file_id, "global_tallies", global_tallies) @@ -518,7 +518,7 @@ contains if (master) then ! Write number of realizations - call write_dataset(file_id, "n_realizations", n_realizations_batch) + call write_dataset(file_id, "n_realizations", n_realizations) ! Write number of global tallies call write_dataset(file_id, "n_global_tallies", N_GLOBAL_TALLIES) @@ -756,8 +756,7 @@ contains #endif ! Read number of realizations for global tallies - call read_dataset(n_realizations_batch, file_id, "n_realizations", & - indep=.true.) + call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) ! Read global tally data call read_dataset(global_tallies, file_id, "global_tallies") diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 9761ef288a..ed0963a883 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -3784,9 +3784,9 @@ contains ! Increase number of realizations (only used for global tallies) if (reduce_tallies) then - n_realizations_batch = n_realizations_batch + 1 + n_realizations = n_realizations + 1 else - n_realizations_batch = n_realizations_batch + n_procs + n_realizations = n_realizations + n_procs end if ! Accumulate on master only unless run is not reduced then do it on all diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 3418df12a6..ee44685924 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -145,9 +145,8 @@ module tally_header type(VectorInt), public :: active_tallies type(VectorInt), public :: active_surface_tallies - ! Normalization for batch-based and generation-based global tally statistics - integer(C_INT32_T), public, bind(C) :: n_realizations_batch = 0 - integer(C_INT32_T), public, bind(C) :: n_realizations_gen = 0 + ! Normalization for statistics + integer(C_INT32_T), public, bind(C) :: n_realizations = 0 ! # of independent realizations real(8), public :: total_weight ! total starting particle weight in realization contains From b88b9ba83126510a80b0499ace1ae57e507c2696 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 2 Aug 2018 15:52:18 -0500 Subject: [PATCH 04/36] Fix #1028 for restart runs --- src/simulation.F90 | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/src/simulation.F90 b/src/simulation.F90 index ffcfb8bab8..c11b5bc2fb 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -386,6 +386,13 @@ contains end do 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 + end if + ! Write message at end if (current_batch == restart_batch) then call write_message("Resuming simulation...", 6) @@ -448,6 +455,19 @@ contains 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 + n_realizations = 0 + end if + ! Display header if (master) then if (run_mode == MODE_FIXEDSOURCE) then From 6f99f1ae18265f37924afbb84c7ded30946534c0 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 09:38:39 -0500 Subject: [PATCH 05/36] Tally the particle only when it is within the range. --- src/tallies/tally_filter_zernike.F90 | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index d005416138..ad4f592770 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -67,15 +67,17 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - theta = atan2(y, x) + if (r <= 1) then + theta = atan2(y, x) + + ! Get moments for Zernike polynomial orders 0..n + call calc_zn(this % order, r, theta, zn) - ! Get moments for Zernike polynomial orders 0..n - call calc_zn(this % order, r, theta, zn) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn(i)) - end do + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn(i)) + end do + endif end subroutine get_all_bins subroutine to_statepoint(this, filter_group) From cdfa3565c9f95f793206a194563a1948a48e1572 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 12:42:18 -0500 Subject: [PATCH 06/36] Add math function for radial Zernike --- src/math.F90 | 10 ++++++++++ src/math_functions.cpp | 29 +++++++++++++++++++++++++++++ src/math_functions.h | 23 +++++++++++++++++++++++ 3 files changed, 62 insertions(+) diff --git a/src/math.F90 b/src/math.F90 index 33f1ab4d30..87c747f579 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -11,6 +11,7 @@ module math public :: calc_pn public :: calc_rn public :: calc_zn + public :: calc_zn_rad public :: evaluate_legendre public :: rotate_angle public :: maxwell_spectrum @@ -70,6 +71,15 @@ module math real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) end subroutine calc_zn + pure subroutine calc_zn_rad(n, rho, phi, zn_rad) bind(C, name='calc_zn_rad_c') + use ISO_C_BINDING + implicit none + integer(C_INT), value, intent(in) :: n + real(C_DOUBLE), value, intent(in) :: rho + real(C_DOUBLE), value, intent(in) :: phi + real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1) + end subroutine calc_zn + subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING implicit none diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 31f86ea5f6..2a79bb4b24 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -587,6 +587,35 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) { } +void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { + // Calculate R_p0(rho) as Zn_p0(rho) + // Set up the array of the coefficients + int length = int(n/2) + 1; + double zn_rad[length]; + double q = 0; + + // R_00 is always 1 + zn_rad[0] = 1; + + // Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong) + for (int p = 2; p <= n; p += 2) { + int index = int(p/2); + if (p == 2) { + // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) + R_22 = std::pow(rho, 2); + zn_rad[index] = 2 * R_22 - zn_rad[0] + } + else { + double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; + double k2 = 2 * p * (p - 1) * (p - 2); + double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2); + double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.; + zn_rad[index] = + ((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1; + } + } +} + void rotate_angle_c(double uvw[3], double mu, double* phi) { // Copy original directional cosines diff --git a/src/math_functions.h b/src/math_functions.h index 2ceea98d87..98882463df 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -91,6 +91,29 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); +//============================================================================== +//! Calculate only the even radial components of n-th order modified Zernike +//! polynomial moment with azimuthal dependency m = 0 for a given angle +//! (rho, theta) location on the unit disk. +//! +//! Since m = 0, n could only be even orders. Z_q0 = R_q0 +//! +//! This procedure uses the modified Kintner's method for calculating Zernike +//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan, +//! R. (2003). A comparative analysis of algorithms for fast computation of +//! Zernike moments. Pattern Recognition, 36(3), 731-742. +//! The normalization of the polynomials is such that the integral of Z_pq^2 +//! over the unit disk is exactly pi. +//! +//! @param n The maximum order requested +//! @param rho The radial parameter to specify location on the unit disk +//! @param phi The angle parameter to specify location on the unit disk +//! @param zn_rad The requested moments of order 0 to n (inclusive) +//! evaluated at rho and phi when m = 0. +//============================================================================== + +extern "C" void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]); + //============================================================================== //! Rotate the direction cosines through a polar angle whose cosine is mu and //! through an azimuthal angle sampled uniformly. From 1013782e125669bffcebb289707ab01a7dd8369a Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:40:42 -0500 Subject: [PATCH 07/36] Add new filter file in CMakeLists --- CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 15c89142e0..878a41b935 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -377,6 +377,7 @@ add_library(libopenmc SHARED src/tallies/tally_filter_surface.F90 src/tallies/tally_filter_universe.F90 src/tallies/tally_filter_zernike.F90 + src/tallies/tally_filter_zernike_radial.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 From dfa52096c51917741a1d1324f92e9aee326362b5 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:42:14 -0500 Subject: [PATCH 08/36] Add new filter for radial zernike --- src/tallies/tally_filter_zernike_radial.F90 | 203 ++++++++++++++++++++ 1 file changed, 203 insertions(+) create mode 100644 src/tallies/tally_filter_zernike_radial.F90 diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 new file mode 100644 index 0000000000..6eaf97f536 --- /dev/null +++ b/src/tallies/tally_filter_zernike_radial.F90 @@ -0,0 +1,203 @@ +module tally_filter_zernike_radial + + use, intrinsic :: ISO_C_BINDING + + use constants + use error + use hdf5_interface + use math, only: calc_zn_rad + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + use xml_interface + + implicit none + private + +!=============================================================================== +! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a +! particle's position +!=============================================================================== + + type, public, extends(TallyFilter) :: ZernikeRadialFilter + integer :: order + real(8) :: x + real(8) :: y + real(8) :: r + contains + procedure :: from_xml + procedure :: get_all_bins + procedure :: to_statepoint + procedure :: text_label + end type ZernikeRadialFilter + +contains + +!=============================================================================== +! ZernikeRadialFilter methods +!=============================================================================== + + subroutine from_xml(this, node) + class(ZernikeRadialFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + integer :: n + + ! Get center of cylinder and radius + call get_node_value(node, "x", this % x) + call get_node_value(node, "y", this % y) + call get_node_value(node, "r", this % r) + + ! Get specified order + call get_node_value(node, "order", n) + this % order = n + this % n_bins = n/2 + 1 + end subroutine from_xml + + subroutine get_all_bins(this, p, estimator, match) + class(ZernikeRadialFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + real(8) :: x, y, r + real(C_DOUBLE) :: zn_rad(this % n_bins) + + ! Determine normalized (r,theta) positions + x = p % coord(1) % xyz(1) - this % x + y = p % coord(1) % xyz(2) - this % y + r = sqrt(x*x + y*y)/this % r + if (r <= 1) then + + ! Get moments for even order Zernike polynomial orders 0..n + call calc_zn_rad(this % order, r, zn_rad) + + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn_rad(i)) + end do + endif + end subroutine get_all_bins + + subroutine to_statepoint(this, filter_group) + class(ZernikeRadialFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "zernikeradial") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "order", this % order) + call write_dataset(filter_group, "x", this % x) + call write_dataset(filter_group, "y", this % y) + call write_dataset(filter_group, "r", this % r) + end subroutine to_statepoint + + function text_label(this, bin) result(label) + class(ZernikeRadialFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" + ! integer :: n + ! integer :: first, last + + ! do n = 0, this % order, -2 + ! last = (n + 1)*(n + 2)/2 + ! if (bin <= last) then + ! first = last - n + ! m = -n + (bin - first)*2 + ! label = "Zernike expansion, Z" // trim(to_str(n)) // ",0" + ! exit + ! end if + ! end do + end function text_label + +!=============================================================================== +! C API FUNCTIONS +!=============================================================================== + + function openmc_zernike_radial_filter_get_order(index, order) result(err) bind(C) + ! Get the order of an expansion filter + integer(C_INT32_T), value :: index + integer(C_INT), intent(out) :: order + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + order = f % order + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_get_order + + + function openmc_zernike_radial_filter_get_params(index, x, y, r) result(err) bind(C) + ! Get the Zernike filter parameters + integer(C_INT32_T), value :: index + real(C_DOUBLE), intent(out) :: x + real(C_DOUBLE), intent(out) :: y + real(C_DOUBLE), intent(out) :: r + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + x = f % x + y = f % y + r = f % r + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_get_params + + + function openmc_zernike_radial_filter_set_order(index, order) result(err) bind(C) + ! Set the order of an expansion filter + integer(C_INT32_T), value :: index + integer(C_INT), value :: order + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + f % order = order + f % n_bins = order/2 + 1 + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_set_order + + + function openmc_zernike_radial_filter_set_params(index, x, y, r) result(err) bind(C) + ! Set the Zernike filter parameters + integer(C_INT32_T), value :: index + real(C_DOUBLE), intent(in), optional :: x + real(C_DOUBLE), intent(in), optional :: y + real(C_DOUBLE), intent(in), optional :: r + integer(C_INT) :: err + + err = verify_filter(index) + if (err == 0) then + select type (f => filters(index) % obj) + type is (ZernikeRadialFilter) + if (present(x)) f % x = x + if (present(y)) f % y = y + if (present(r)) f % r = r + class default + err = E_INVALID_TYPE + call set_errmsg("Not a Zernike Radial filter.") + end select + end if + end function openmc_zernike_radial_filter_set_params + +end module tally_filter_zernike_radial From f8cac6004980586262732837997d8f8393e20ce9 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:44:21 -0500 Subject: [PATCH 09/36] Add C API for the new radial zernike filters --- openmc/capi/filter.py | 30 ++++++++++++++++++++++++++++-- openmc/capi/math.py | 24 ++++++++++++++++++++++++ 2 files changed, 52 insertions(+), 2 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 391e63f073..621c2a4f6c 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -20,7 +20,7 @@ __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter', 'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', - 'UniverseFilter', 'ZernikeFilter', 'filters'] + 'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters'] # Tally functions _dll.openmc_cell_filter_get_bins.argtypes = [ @@ -95,6 +95,12 @@ _dll.openmc_zernike_filter_get_order.errcheck = _error_handler _dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int] _dll.openmc_zernike_filter_set_order.restype = c_int _dll.openmc_zernike_filter_set_order.errcheck = _error_handler +_dll.openmc_zernike_radial_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_zernike_radial_filter_get_order.restype = c_int +_dll.openmc_zernike_radial_filter_get_order.errcheck = _error_handler +_dll.openmc_zernike_radial_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_zernike_radial_filter_set_order.restype = c_int +_dll.openmc_zernike_radial_filter_set_order.errcheck = _error_handler class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -360,6 +366,25 @@ class ZernikeFilter(Filter): _dll.openmc_zernike_filter_set_order(self._index, order) +class ZernikeRadialFilter(Filter): + filter_type = 'zernikeradial' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_zernike_radial_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_zernike_radial_filter_set_order(self._index, order) + + _FILTER_TYPE_MAP = { 'azimuthal': AzimuthalFilter, 'cell': CellFilter, @@ -380,7 +405,8 @@ _FILTER_TYPE_MAP = { 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, 'universe': UniverseFilter, - 'zernike': ZernikeFilter + 'zernike': ZernikeFilter, + 'zernikeradial': ZernikeRadialFilter } diff --git a/openmc/capi/math.py b/openmc/capi/math.py index d1d7abdf64..4f8bb5cdbe 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -21,6 +21,9 @@ _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)] _dll.calc_zn_c.restype = None _dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)] +_dll.calc_zn_rad_c.restype = None +_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)] + _dll.rotate_angle_c.restype = None _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double, POINTER(c_double)] @@ -154,6 +157,27 @@ def calc_zn(n, rho, phi): _dll.calc_zn_c(n, rho, phi, zn) return zn +def calc_zn_rad(n, rho): + """ Calculate the even orders in n-th order modified Zernike polynomial moment with no azimuthal dependency (m=0) for a given radial location in the unit disk. The normalization of the polynomials is such that the integral of Z_pq*Z_pq over the unit disk is exactly pi + + Parameters + ---------- + n : int + Maximum order + rho : float + Radial location in the unit disk + + Returns + ------- + numpy.ndarray + Corresponding resulting list of coefficients + + """ + + num_bins = n // 2 + 1 + zn_rad = np.zeros(num_bins, dtype=np.float64) + _dll.calc_zn_rad_c(n, rho, zn_rad) + return zn_rad def rotate_angle(uvw0, mu, phi=None): """ Rotates direction cosines through a polar angle whose cosine is From aa28f4a1452089e0e83ff2b82543d31f5f67ccb8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:46:00 -0500 Subject: [PATCH 10/36] Add new python class for radial zernike filter --- openmc/filter.py | 2 +- openmc/filter_expansion.py | 133 +++++++++++++++++++++++++++++++++++++ 2 files changed, 134 insertions(+), 1 deletion(-) diff --git a/openmc/filter.py b/openmc/filter.py index 3180989eb7..ae44024941 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,7 +22,7 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike' + 'sphericalharmonics', 'zernike', 'zernikeradial' ) _CURRENT_NAMES = ( diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index 1970c07710..842d03477f 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -463,3 +463,136 @@ class ZernikeFilter(ExpansionFilter): subelement.text = str(self.r) return element + +class ZernikeRadialFilter(ExpansionFilter): + r"""Score radial Zernike expansion moments in space up to specified order. + + The Zernike polynomials are defined the same as in ZernikeFilter. + + .. math:: + Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0 + + Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0 + + Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0 + + where the radial polynomials are + + .. math:: + R_n^m(\rho) = \sum\limits_{k=0}^{(n-m)/2} \frac{(-1)^k (n-k)!}{k! ( + \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. + + With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk + is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is + 2 if :math:`m` equals 0 and 1 otherwise. + + If there is only radial dependency, the polynomials are integrated over the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. Therefore, for a radial Zernike polynomials up to order of :math:`n`, there are :math:`\frac{n}{2} + 1` terms in total. + + Parameters + ---------- + order : int + Maximum radial Zernike polynomial order + x : float + x-coordinate of center of circle for normalization + y : float + y-coordinate of center of circle for normalization + r : int or None + Radius of circle for normalization + + Attributes + ---------- + order : int + Maximum radial Zernike polynomial order + x : float + x-coordinate of center of circle for normalization + y : float + y-coordinate of center of circle for normalization + r : int or None + Radius of circle for normalization + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins + + """ + + def __init__(self, order, x=0.0, y=0.0, r=1.0, filter_id=None): + super().__init__(order, filter_id) + self.x = x + self.y = y + self.r = r + + def __hash__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + return hash(string) + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tOrder', self.order) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string + + @ExpansionFilter.order.setter + def order(self, order): + ExpansionFilter.order.__set__(self, order) + self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)] + + @property + def x(self): + return self._x + + @x.setter + def x(self, x): + cv.check_type('x', x, Real) + self._x = x + + @property + def y(self): + return self._y + + @y.setter + def y(self, y): + cv.check_type('y', y, Real) + self._y = y + + @property + def r(self): + return self._r + + @r.setter + def r(self, r): + cv.check_type('r', r, Real) + self._r = r + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'].value.decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'].value.decode() + " instead") + + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + order = group['order'].value + x, y, r = group['x'].value, group['y'].value, group['r'].value + + return cls(order, x, y, r, filter_id) + + def to_xml_element(self): + """Return XML Element representing the filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing radial Zernike filter data + + """ + element = super().to_xml_element() + subelement = ET.SubElement(element, 'x') + subelement.text = str(self.x) + subelement = ET.SubElement(element, 'y') + subelement.text = str(self.y) + subelement = ET.SubElement(element, 'r') + subelement.text = str(self.r) + + return element \ No newline at end of file From 76c4ba0c7823899e5c929c3caa7d15179de21d8b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:47:32 -0500 Subject: [PATCH 11/36] Math for radial components of zernike polynomials --- src/math.F90 | 5 ++--- src/math_functions.cpp | 9 ++++----- src/math_functions.h | 2 +- 3 files changed, 7 insertions(+), 9 deletions(-) diff --git a/src/math.F90 b/src/math.F90 index 87c747f579..ab5eff14c1 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -71,14 +71,13 @@ module math real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2) end subroutine calc_zn - pure subroutine calc_zn_rad(n, rho, phi, zn_rad) bind(C, name='calc_zn_rad_c') + pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c') use ISO_C_BINDING implicit none integer(C_INT), value, intent(in) :: n real(C_DOUBLE), value, intent(in) :: rho - real(C_DOUBLE), value, intent(in) :: phi real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1) - end subroutine calc_zn + end subroutine calc_zn_rad subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c') use ISO_C_BINDING diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 2a79bb4b24..53bf296dfe 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -587,11 +587,10 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) { } -void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { +void calc_zn_rad_c(int n, double rho, double zn_rad[]) { // Calculate R_p0(rho) as Zn_p0(rho) // Set up the array of the coefficients - int length = int(n/2) + 1; - double zn_rad[length]; + double q = 0; // R_00 is always 1 @@ -602,8 +601,8 @@ void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]) { int index = int(p/2); if (p == 2) { // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) - R_22 = std::pow(rho, 2); - zn_rad[index] = 2 * R_22 - zn_rad[0] + double R_22 = std::pow(rho, 2); + zn_rad[index] = 2 * R_22 - zn_rad[0]; } else { double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; diff --git a/src/math_functions.h b/src/math_functions.h index 98882463df..0993b1a01b 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -112,7 +112,7 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); //! evaluated at rho and phi when m = 0. //============================================================================== -extern "C" void calc_zn_rad_c(int n, double rho, double phi, double zn_rad[]); +extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]); //============================================================================== //! Rotate the direction cosines through a polar angle whose cosine is mu and From 0f9920719f1f1278ba340eb36839e8289c5c02bf Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:48:10 -0500 Subject: [PATCH 12/36] Update the constant for radial zernike filter --- src/constants.F90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index ab6887610e..f377f55f33 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -343,7 +343,7 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 20 + integer, parameter :: N_FILTER_TYPES = 21 integer, parameter :: & FILTER_UNIVERSE = 1, & FILTER_MATERIAL = 2, & @@ -364,7 +364,8 @@ module constants FILTER_LEGENDRE = 17, & FILTER_SPH_HARMONICS = 18, & FILTER_SPTL_LEGENDRE = 19, & - FILTER_ZERNIKE = 20 + FILTER_ZERNIKE = 20, & + FILTER_ZERNIKE_RADIAL = 21 ! Mesh types integer, parameter :: & From 8e9ed84e111fa01e827d8fc38fc8a96b0e4749a8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 3 Aug 2018 16:49:51 -0500 Subject: [PATCH 13/36] Add the radial zernike filter in all the F90 files --- src/tallies/tally_filter.F90 | 5 +++++ src/tallies/tally_header.F90 | 3 +++ 2 files changed, 8 insertions(+) diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 4b589cedbf..8d8bfadd57 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -26,6 +26,7 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike + use tally_filter_zernike_radial implicit none @@ -88,6 +89,8 @@ contains type_ = 'universe' type is (ZernikeFilter) type_ = 'zernike' + type is (ZernikeRadialFilter) + type_ = 'zernikeradial' end select ! Convert Fortran string to null-terminated C string. We assume the @@ -159,6 +162,8 @@ contains allocate(UniverseFilter :: filters(index) % obj) case ('zernike') allocate(ZernikeFilter :: filters(index) % obj) + case ('zernikeradial') + allocate(ZernikeRadialFilter :: filters(index) % obj) case default err = E_UNASSIGNED call set_errmsg("Unknown filter type: " // trim(type_)) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index ee44685924..224fafc3b8 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -339,6 +339,9 @@ contains type is (ZernikeFilter) j = FILTER_ZERNIKE this % estimator = ESTIMATOR_COLLISION + type is (ZernikeRadialFilter) + j = FILTER_ZERNIKE_RADIAL + this % estimator = ESTIMATOR_COLLISION end select this % find_filter(j) = i end do From 4c1c1c2cf00581e20a677df3a1e825d029d83690 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Mon, 6 Aug 2018 08:37:23 -0500 Subject: [PATCH 14/36] Add unit test for Radial Zernike filter --- tests/unit_tests/test_filters.py | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 488badbfaf..5817c4c3d2 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -100,6 +100,23 @@ def test_zernike(): assert elem.attrib['type'] == 'zernike' assert elem.find('order').text == str(n) +def test_zernike_radial(): + n = 4 + f = openmc.ZernikeRadialFilter(n, 0., 0., 1.) + assert f.order == n + assert f.bins[0] == 'Z0,0' + assert f.bins[-1] == 'Z{},0'.format(n) + assert len(f.bins) == n//2 + 1 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'zernikeradial' + assert elem.find('order').text == str(n) + def test_first_moment(run_in_tmpdir, box_model): plain_tally = openmc.Tally() From f0429f53647bd48212b688e35631e83db55f8cf4 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 6 Aug 2018 19:42:56 -0400 Subject: [PATCH 15/36] Read n_realizations from all processes --- src/state_point.F90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/state_point.F90 b/src/state_point.F90 index ebb473d496..b8d808f2ae 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -741,6 +741,9 @@ contains end if end if + ! Read number of realizations for global tallies + call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) + ! 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 & @@ -754,10 +757,6 @@ contains #else if (master) then #endif - - ! Read number of realizations for global tallies - call read_dataset(n_realizations, file_id, "n_realizations", indep=.true.) - ! Read global tally data call read_dataset(global_tallies, file_id, "global_tallies") From d65b7a11bd73426ef142056ada287841f8db0c4c Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 10:19:29 -0500 Subject: [PATCH 16/36] Some travial modfication to meet the coding standard. --- openmc/capi/math.py | 7 ++++++- openmc/filter_expansion.py | 23 +++++++++++++-------- src/math_functions.cpp | 5 ++--- src/tallies/tally_filter_zernike_radial.F90 | 12 ----------- 4 files changed, 22 insertions(+), 25 deletions(-) diff --git a/openmc/capi/math.py b/openmc/capi/math.py index 4f8bb5cdbe..390ee466ff 100644 --- a/openmc/capi/math.py +++ b/openmc/capi/math.py @@ -157,8 +157,12 @@ def calc_zn(n, rho, phi): _dll.calc_zn_c(n, rho, phi, zn) return zn + def calc_zn_rad(n, rho): - """ Calculate the even orders in n-th order modified Zernike polynomial moment with no azimuthal dependency (m=0) for a given radial location in the unit disk. The normalization of the polynomials is such that the integral of Z_pq*Z_pq over the unit disk is exactly pi + """ Calculate the even orders in n-th order modified Zernike polynomial + moment with no azimuthal dependency (m=0) for a given radial location in + the unit disk. The normalization of the polynomials is such that the + integral of Z_pq*Z_pq over the unit disk is exactly pi. Parameters ---------- @@ -178,6 +182,7 @@ def calc_zn_rad(n, rho): zn_rad = np.zeros(num_bins, dtype=np.float64) _dll.calc_zn_rad_c(n, rho, zn_rad) return zn_rad + def rotate_angle(uvw0, mu, phi=None): """ Rotates direction cosines through a polar angle whose cosine is diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index 842d03477f..a756f7ba4a 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -234,8 +234,8 @@ class SphericalHarmonicsFilter(ExpansionFilter): r"""Score spherical harmonic expansion moments up to specified order. This filter allows you to obtain real spherical harmonic moments of either - the particle's direction or the cosine of the scattering angle. Specifying a - filter with order :math:`\ell` tallies moments for all orders from 0 to + the particle's direction or the cosine of the scattering angle. Specifying + a filter with order :math:`\ell` tallies moments for all orders from 0 to :math:`\ell`. Parameters @@ -342,11 +342,11 @@ class ZernikeFilter(ExpansionFilter): \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk - is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is - 2 if :math:`m` equals 0 and 1 otherwise. + is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where + :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise. - Specifying a filter with order N tallies moments for all :math:`n` from 0 to - N and each value of :math:`m`. The ordering of the Zernike polynomial + Specifying a filter with order N tallies moments for all :math:`n` from 0 + to N and each value of :math:`m`. The ordering of the Zernike polynomial moments follows the ANSI Z80.28 standard, where the one-dimensional index :math:`j` corresponds to the :math:`n` and :math:`m` by @@ -464,6 +464,7 @@ class ZernikeFilter(ExpansionFilter): return element + class ZernikeRadialFilter(ExpansionFilter): r"""Score radial Zernike expansion moments in space up to specified order. @@ -483,10 +484,14 @@ class ZernikeRadialFilter(ExpansionFilter): \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk - is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is - 2 if :math:`m` equals 0 and 1 otherwise. + is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where + :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise. - If there is only radial dependency, the polynomials are integrated over the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. Therefore, for a radial Zernike polynomials up to order of :math:`n`, there are :math:`\frac{n}{2} + 1` terms in total. + If there is only radial dependency, the polynomials are integrated over + the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) + = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. + Therefore, for a radial Zernike polynomials up to order of :math:`n`, + there are :math:`\frac{n}{2} + 1` terms in total. Parameters ---------- diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 3d60ce4d0b..8a36f7b922 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -601,10 +601,9 @@ void calc_zn_rad_c(int n, double rho, double zn_rad[]) { int index = int(p/2); if (p == 2) { // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong) - double R_22 = std::pow(rho, 2); + double R_22 = rho * rho; zn_rad[index] = 2 * R_22 - zn_rad[0]; - } - else { + } else { double k1 = ((p + q) * (p - q) * (p - 2)) / 2.; double k2 = 2 * p * (p - 1) * (p - 2); double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2); diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 index 6eaf97f536..1a4ed5ed55 100644 --- a/src/tallies/tally_filter_zernike_radial.F90 +++ b/src/tallies/tally_filter_zernike_radial.F90 @@ -98,18 +98,6 @@ contains character(MAX_LINE_LEN) :: label label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" - ! integer :: n - ! integer :: first, last - - ! do n = 0, this % order, -2 - ! last = (n + 1)*(n + 2)/2 - ! if (bin <= last) then - ! first = last - n - ! m = -n + (bin - first)*2 - ! label = "Zernike expansion, Z" // trim(to_str(n)) // ",0" - ! exit - ! end if - ! end do end function text_label !=============================================================================== From 87209a8af217f0cf9035c35cac014428e5b2cf73 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:22:11 -0500 Subject: [PATCH 17/36] ZernikeRadial as an extension of Zernike --- src/tallies/tally_filter_zernike.F90 | 119 ++++++++++-- src/tallies/tally_filter_zernike_radial.F90 | 191 -------------------- 2 files changed, 106 insertions(+), 204 deletions(-) delete mode 100644 src/tallies/tally_filter_zernike_radial.F90 diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index ad4f592770..d574c542d5 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -5,7 +5,7 @@ module tally_filter_zernike use constants use error use hdf5_interface - use math, only: calc_zn + use math, only: calc_zn, calc_zn_rad use particle_header, only: Particle use string, only: to_str use tally_filter_header @@ -24,19 +24,32 @@ module tally_filter_zernike real(8) :: y real(8) :: r contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label + procedure :: from_xml => from_xml_zn + procedure :: get_all_bins => get_all_bins_zn + procedure :: to_statepoint => to_statepoint_zn + procedure :: text_label => text_label_zn end type ZernikeFilter +!=============================================================================== +! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a +! particle's position +!=============================================================================== + + type, public, extends(ZernikeFilter) :: ZernikeRadialFilter + contains + procedure :: from_xml => from_xml_zn_rad + procedure :: get_all_bins => get_all_bins_zn_rad + procedure :: to_statepoint => to_statepoint_zn_rad + procedure :: text_label => text_label_zn_rad + end type ZernikeRadialFilter + contains !=============================================================================== ! ZernikeFilter methods !=============================================================================== - subroutine from_xml(this, node) + subroutine from_xml_zn(this, node) class(ZernikeFilter), intent(inout) :: this type(XMLNode), intent(in) :: node @@ -51,9 +64,9 @@ contains call get_node_value(node, "order", n) this % order = n this % n_bins = ((n + 1)*(n + 2))/2 - end subroutine from_xml + end subroutine from_xml_zn - subroutine get_all_bins(this, p, estimator, match) + subroutine get_all_bins_zn(this, p, estimator, match) class(ZernikeFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator @@ -78,9 +91,9 @@ contains call match % weights % push_back(zn(i)) end do endif - end subroutine get_all_bins + end subroutine get_all_bins_zn - subroutine to_statepoint(this, filter_group) + subroutine to_statepoint_zn(this, filter_group) class(ZernikeFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group @@ -90,9 +103,9 @@ contains call write_dataset(filter_group, "x", this % x) call write_dataset(filter_group, "y", this % y) call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint + end subroutine to_statepoint_zn - function text_label(this, bin) result(label) + function text_label_zn(this, bin) result(label) class(ZernikeFilter), intent(in) :: this integer, intent(in) :: bin character(MAX_LINE_LEN) :: label @@ -110,7 +123,74 @@ contains exit end if end do - end function text_label + end function text_label_zn + +!=============================================================================== +! ZernikeRadialFilter methods +!=============================================================================== + + subroutine from_xml_zn_rad(this, node) + class(ZernikeRadialFilter), intent(inout) :: this + type(XMLNode), intent(in) :: node + + integer :: n + + ! Get center of cylinder and radius + call get_node_value(node, "x", this % x) + call get_node_value(node, "y", this % y) + call get_node_value(node, "r", this % r) + + ! Get specified order + call get_node_value(node, "order", n) + this % order = n + this % n_bins = n/2 + 1 + end subroutine from_xml_zn_rad + + subroutine get_all_bins_zn_rad(this, p, estimator, match) + class(ZernikeRadialFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: i + real(8) :: x, y, r + real(C_DOUBLE) :: zn_rad(this % n_bins) + + ! Determine normalized (r,theta) positions + x = p % coord(1) % xyz(1) - this % x + y = p % coord(1) % xyz(2) - this % y + r = sqrt(x*x + y*y)/this % r + if (r <= 1) then + + ! Get moments for even order Zernike polynomial orders 0..n + call calc_zn_rad(this % order, r, zn_rad) + + do i = 1, this % n_bins + call match % bins % push_back(i) + call match % weights % push_back(zn_rad(i)) + end do + endif + end subroutine get_all_bins_zn_rad + + subroutine to_statepoint_zn_rad(this, filter_group) + class(ZernikeRadialFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "zernikeradial") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "order", this % order) + call write_dataset(filter_group, "x", this % x) + call write_dataset(filter_group, "y", this % y) + call write_dataset(filter_group, "r", this % r) + end subroutine to_statepoint_zn_rad + + function text_label_zn_rad(this, bin) result(label) + class(ZernikeRadialFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" + end function text_label_zn_rad !=============================================================================== ! C API FUNCTIONS @@ -127,6 +207,8 @@ contains select type (f => filters(index) % obj) type is (ZernikeFilter) order = f % order + type is (ZernikeRadialFilter) + order = f % order class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -150,6 +232,10 @@ contains x = f % x y = f % y r = f % r + type is (ZernikeRadialFilter) + x = f % x + y = f % y + r = f % r class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -170,6 +256,9 @@ contains type is (ZernikeFilter) f % order = order f % n_bins = ((order + 1)*(order + 2))/2 + type is (ZernikeRadialFilter) + f % order = order + f % n_bins = order/2 + 1 class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") @@ -193,6 +282,10 @@ contains if (present(x)) f % x = x if (present(y)) f % y = y if (present(r)) f % r = r + type is (ZernikeRadialFilter) + if (present(x)) f % x = x + if (present(y)) f % y = y + if (present(r)) f % r = r class default err = E_INVALID_TYPE call set_errmsg("Not a Zernike filter.") diff --git a/src/tallies/tally_filter_zernike_radial.F90 b/src/tallies/tally_filter_zernike_radial.F90 deleted file mode 100644 index 1a4ed5ed55..0000000000 --- a/src/tallies/tally_filter_zernike_radial.F90 +++ /dev/null @@ -1,191 +0,0 @@ -module tally_filter_zernike_radial - - use, intrinsic :: ISO_C_BINDING - - use constants - use error - use hdf5_interface - use math, only: calc_zn_rad - use particle_header, only: Particle - use string, only: to_str - use tally_filter_header - use xml_interface - - implicit none - private - -!=============================================================================== -! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a -! particle's position -!=============================================================================== - - type, public, extends(TallyFilter) :: ZernikeRadialFilter - integer :: order - real(8) :: x - real(8) :: y - real(8) :: r - contains - procedure :: from_xml - procedure :: get_all_bins - procedure :: to_statepoint - procedure :: text_label - end type ZernikeRadialFilter - -contains - -!=============================================================================== -! ZernikeRadialFilter methods -!=============================================================================== - - subroutine from_xml(this, node) - class(ZernikeRadialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - ! Get center of cylinder and radius - call get_node_value(node, "x", this % x) - call get_node_value(node, "y", this % y) - call get_node_value(node, "r", this % r) - - ! Get specified order - call get_node_value(node, "order", n) - this % order = n - this % n_bins = n/2 + 1 - end subroutine from_xml - - subroutine get_all_bins(this, p, estimator, match) - class(ZernikeRadialFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: i - real(8) :: x, y, r - real(C_DOUBLE) :: zn_rad(this % n_bins) - - ! Determine normalized (r,theta) positions - x = p % coord(1) % xyz(1) - this % x - y = p % coord(1) % xyz(2) - this % y - r = sqrt(x*x + y*y)/this % r - if (r <= 1) then - - ! Get moments for even order Zernike polynomial orders 0..n - call calc_zn_rad(this % order, r, zn_rad) - - do i = 1, this % n_bins - call match % bins % push_back(i) - call match % weights % push_back(zn_rad(i)) - end do - endif - end subroutine get_all_bins - - subroutine to_statepoint(this, filter_group) - class(ZernikeRadialFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "zernikeradial") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "order", this % order) - call write_dataset(filter_group, "x", this % x) - call write_dataset(filter_group, "y", this % y) - call write_dataset(filter_group, "r", this % r) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(ZernikeRadialFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" - end function text_label - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_zernike_radial_filter_get_order(index, order) result(err) bind(C) - ! Get the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), intent(out) :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - order = f % order - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_get_order - - - function openmc_zernike_radial_filter_get_params(index, x, y, r) result(err) bind(C) - ! Get the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(out) :: x - real(C_DOUBLE), intent(out) :: y - real(C_DOUBLE), intent(out) :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - x = f % x - y = f % y - r = f % r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_get_params - - - function openmc_zernike_radial_filter_set_order(index, order) result(err) bind(C) - ! Set the order of an expansion filter - integer(C_INT32_T), value :: index - integer(C_INT), value :: order - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - f % order = order - f % n_bins = order/2 + 1 - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_set_order - - - function openmc_zernike_radial_filter_set_params(index, x, y, r) result(err) bind(C) - ! Set the Zernike filter parameters - integer(C_INT32_T), value :: index - real(C_DOUBLE), intent(in), optional :: x - real(C_DOUBLE), intent(in), optional :: y - real(C_DOUBLE), intent(in), optional :: r - integer(C_INT) :: err - - err = verify_filter(index) - if (err == 0) then - select type (f => filters(index) % obj) - type is (ZernikeRadialFilter) - if (present(x)) f % x = x - if (present(y)) f % y = y - if (present(r)) f % r = r - class default - err = E_INVALID_TYPE - call set_errmsg("Not a Zernike Radial filter.") - end select - end if - end function openmc_zernike_radial_filter_set_params - -end module tally_filter_zernike_radial From abdd753d5926c92ef1ea15226084a800bb56d1b7 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:23:45 -0500 Subject: [PATCH 18/36] Delete unecessary names in files --- CMakeLists.txt | 1 - src/tallies/tally_filter.F90 | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5f5706124f..73b3a3851b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -380,7 +380,6 @@ add_library(libopenmc SHARED src/tallies/tally_filter_surface.F90 src/tallies/tally_filter_universe.F90 src/tallies/tally_filter_zernike.F90 - src/tallies/tally_filter_zernike_radial.F90 src/tallies/tally_header.F90 src/tallies/trigger.F90 src/tallies/trigger_header.F90 diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index e2dcdd87d4..efeac180e2 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -27,7 +27,7 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike - use tally_filter_zernike_radial + ! use tally_filter_zernike_radial implicit none From 0b45dc66c4874b3d55f2a19b756ecfec21d87f8f Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:34:29 -0500 Subject: [PATCH 19/36] Change the ZernikeRadial filter on the Python API side as a subclass of Zernike filter. --- openmc/capi/filter.py | 21 ---------- openmc/filter_expansion.py | 78 +------------------------------------- 2 files changed, 1 insertion(+), 98 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 621c2a4f6c..5bd451b300 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -95,12 +95,6 @@ _dll.openmc_zernike_filter_get_order.errcheck = _error_handler _dll.openmc_zernike_filter_set_order.argtypes = [c_int32, c_int] _dll.openmc_zernike_filter_set_order.restype = c_int _dll.openmc_zernike_filter_set_order.errcheck = _error_handler -_dll.openmc_zernike_radial_filter_get_order.argtypes = [c_int32, POINTER(c_int)] -_dll.openmc_zernike_radial_filter_get_order.restype = c_int -_dll.openmc_zernike_radial_filter_get_order.errcheck = _error_handler -_dll.openmc_zernike_radial_filter_set_order.argtypes = [c_int32, c_int] -_dll.openmc_zernike_radial_filter_set_order.restype = c_int -_dll.openmc_zernike_radial_filter_set_order.errcheck = _error_handler class Filter(_FortranObjectWithID): __instances = WeakValueDictionary() @@ -369,21 +363,6 @@ class ZernikeFilter(Filter): class ZernikeRadialFilter(Filter): filter_type = 'zernikeradial' - def __init__(self, order=None, uid=None, new=True, index=None): - super().__init__(uid, new, index) - if order is not None: - self.order = order - - @property - def order(self): - temp_order = c_int() - _dll.openmc_zernike_radial_filter_get_order(self._index, temp_order) - return temp_order.value - - @order.setter - def order(self, order): - _dll.openmc_zernike_radial_filter_set_order(self._index, order) - _FILTER_TYPE_MAP = { 'azimuthal': AzimuthalFilter, diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index a756f7ba4a..ad761f2ff7 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -465,7 +465,7 @@ class ZernikeFilter(ExpansionFilter): return element -class ZernikeRadialFilter(ExpansionFilter): +class ZernikeRadialFilter(ZernikeFilter): r"""Score radial Zernike expansion moments in space up to specified order. The Zernike polynomials are defined the same as in ZernikeFilter. @@ -521,83 +521,7 @@ class ZernikeRadialFilter(ExpansionFilter): """ - def __init__(self, order, x=0.0, y=0.0, r=1.0, filter_id=None): - super().__init__(order, filter_id) - self.x = x - self.y = y - self.r = r - - def __hash__(self): - string = type(self).__name__ + '\n' - string += '{: <16}=\t{}\n'.format('\tOrder', self.order) - return hash(string) - - def __repr__(self): - string = type(self).__name__ + '\n' - string += '{: <16}=\t{}\n'.format('\tOrder', self.order) - string += '{: <16}=\t{}\n'.format('\tID', self.id) - return string - @ExpansionFilter.order.setter def order(self, order): ExpansionFilter.order.__set__(self, order) self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)] - - @property - def x(self): - return self._x - - @x.setter - def x(self, x): - cv.check_type('x', x, Real) - self._x = x - - @property - def y(self): - return self._y - - @y.setter - def y(self, y): - cv.check_type('y', y, Real) - self._y = y - - @property - def r(self): - return self._r - - @r.setter - def r(self, r): - cv.check_type('r', r, Real) - self._r = r - - @classmethod - def from_hdf5(cls, group, **kwargs): - if group['type'].value.decode() != cls.short_name.lower(): - raise ValueError("Expected HDF5 data for filter type '" - + cls.short_name.lower() + "' but got '" - + group['type'].value.decode() + " instead") - - filter_id = int(group.name.split('/')[-1].lstrip('filter ')) - order = group['order'].value - x, y, r = group['x'].value, group['y'].value, group['r'].value - - return cls(order, x, y, r, filter_id) - - def to_xml_element(self): - """Return XML Element representing the filter. - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing radial Zernike filter data - - """ - element = super().to_xml_element() - subelement = ET.SubElement(element, 'x') - subelement.text = str(self.x) - subelement = ET.SubElement(element, 'y') - subelement.text = str(self.y) - subelement = ET.SubElement(element, 'r') - subelement.text = str(self.r) - - return element \ No newline at end of file From a4cdef712d75ce2934b34bc120d99653c37dd32b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Tue, 7 Aug 2018 11:39:44 -0500 Subject: [PATCH 20/36] Remove comment --- src/tallies/tally_filter.F90 | 1 - 1 file changed, 1 deletion(-) diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index efeac180e2..81b1958169 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -27,7 +27,6 @@ module tally_filter use tally_filter_surface use tally_filter_universe use tally_filter_zernike - ! use tally_filter_zernike_radial implicit none From 5a0fc9eebac7e37fdb34d8cd9b98019336c3732e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Aug 2018 10:06:00 -0500 Subject: [PATCH 21/36] Convert whitespace to blanks when reading arrays from XML --- src/xml_interface.F90 | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/src/xml_interface.F90 b/src/xml_interface.F90 index a9dc5224bc..0e08f75d1e 100644 --- a/src/xml_interface.F90 +++ b/src/xml_interface.F90 @@ -261,6 +261,7 @@ contains ! Get value of text node/attribute str = node_value_string(node, name) + call whitespace_to_blanks(str) ! Read numbers into array read(UNIT=str, FMT=*, IOSTAT=stat) array @@ -283,6 +284,7 @@ contains ! Get value of text node/attribute str = node_value_string(node, name) + call whitespace_to_blanks(str) ! Read numbers into array read(UNIT=str, FMT=*, IOSTAT=stat) array @@ -333,4 +335,21 @@ contains end if end subroutine get_node_array_string +!=============================================================================== +! WHITESPACE_TO_BLANKS converts all whitespace to blanks +!=============================================================================== + + subroutine whitespace_to_blanks(str) + character(len=*, kind=C_CHAR), intent(inout) :: str + + integer :: i + + do i = 1, len(str) + select case (str(i:i)) + case (C_NEW_LINE, C_HORIZONTAL_TAB, C_CARRIAGE_RETURN) + str(i:i) = ' ' + end select + end do + end subroutine + end module xml_interface From 1074da87f8a6277276c0cd5abc34afed5b566cf5 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 13:27:41 -0500 Subject: [PATCH 22/36] Made some change accordingly to the suggestions. --- openmc/filter.py | 2 +- src/constants.F90 | 5 +++-- src/math_functions.h | 13 ++++--------- 3 files changed, 8 insertions(+), 12 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 36053ec46c..e63d72412f 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,7 +22,7 @@ _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'particle', 'zernikeradial' + 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle' ) _CURRENT_NAMES = ( diff --git a/src/constants.F90 b/src/constants.F90 index 1bdcec74fb..4de52322e0 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -379,8 +379,9 @@ module constants FILTER_SPH_HARMONICS = 18, & FILTER_SPTL_LEGENDRE = 19, & FILTER_ZERNIKE = 20, & - FILTER_PARTICLE = 21, & - FILTER_ZERNIKE_RADIAL = 22 + FILTER_ZERNIKE_RADIAL = 21, & + FILTER_PARTICLE = 22 + ! Mesh types integer, parameter :: & diff --git a/src/math_functions.h b/src/math_functions.h index e0ecad60ec..7ef39c4461 100644 --- a/src/math_functions.h +++ b/src/math_functions.h @@ -92,18 +92,13 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]); extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]); //============================================================================== -//! Calculate only the even radial components of n-th order modified Zernike -//! polynomial moment with azimuthal dependency m = 0 for a given angle -//! (rho, theta) location on the unit disk. +//! Calculate only the even order components of n-th order modified Zernike +//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho) +//! location on the unit disk. //! //! Since m = 0, n could only be even orders. Z_q0 = R_q0 //! -//! This procedure uses the modified Kintner's method for calculating Zernike -//! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan, -//! R. (2003). A comparative analysis of algorithms for fast computation of -//! Zernike moments. Pattern Recognition, 36(3), 731-742. -//! The normalization of the polynomials is such that the integral of Z_pq^2 -//! over the unit disk is exactly pi. +//! See calc_zn_c for methodology. //! //! @param n The maximum order requested //! @param rho The radial parameter to specify location on the unit disk From a7ce79c62f0cb9cd9f327126435b68c570547fb2 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 13:32:08 -0500 Subject: [PATCH 23/36] ONE instead of 1 --- src/tallies/tally_filter_zernike.F90 | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index d574c542d5..1ece62dbb2 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -80,7 +80,7 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - if (r <= 1) then + if (r <= ONE) then theta = atan2(y, x) ! Get moments for Zernike polynomial orders 0..n @@ -160,7 +160,7 @@ contains x = p % coord(1) % xyz(1) - this % x y = p % coord(1) % xyz(2) - this % y r = sqrt(x*x + y*y)/this % r - if (r <= 1) then + if (r <= ONE) then ! Get moments for even order Zernike polynomial orders 0..n call calc_zn_rad(this % order, r, zn_rad) @@ -228,11 +228,7 @@ contains err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) - type is (ZernikeFilter) - x = f % x - y = f % y - r = f % r - type is (ZernikeRadialFilter) + class is (ZernikeFilter) x = f % x y = f % y r = f % r @@ -278,11 +274,7 @@ contains err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) - type is (ZernikeFilter) - if (present(x)) f % x = x - if (present(y)) f % y = y - if (present(r)) f % r = r - type is (ZernikeRadialFilter) + class is (ZernikeFilter) if (present(x)) f % x = x if (present(y)) f % y = y if (present(r)) f % r = r From d9a7db16777eb6239bcefd1fccf01161eb3fda05 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 25 Jun 2018 21:00:59 -0500 Subject: [PATCH 24/36] Introduce Position/Angle types --- CMakeLists.txt | 1 + src/cell.cpp | 34 ++-- src/cell.h | 11 +- src/geometry.cpp | 59 +++++++ src/geometry.h | 62 +++++++ src/surface.cpp | 449 +++++++++++++++++++++++------------------------ src/surface.h | 130 +++++++------- 7 files changed, 421 insertions(+), 325 deletions(-) create mode 100644 src/geometry.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index ea9bd5969c..0be9ccb325 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -386,6 +386,7 @@ add_library(libopenmc SHARED src/cell.cpp src/initialize.cpp src/finalize.cpp + src/geometry.cpp src/geometry_aux.cpp src/hdf5_interface.cpp src/lattice.cpp diff --git a/src/cell.cpp b/src/cell.cpp index 263191b316..972540fbd2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -279,33 +279,31 @@ Cell::Cell(pugi::xml_node cell_node) //============================================================================== bool -Cell::contains(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains(Position r, Angle a, int32_t on_surface) const { if (simple) { - return contains_simple(xyz, uvw, on_surface); + return contains_simple(r, a, on_surface); } else { - return contains_complex(xyz, uvw, on_surface); + return contains_complex(r, a, on_surface); } } //============================================================================== std::pair -Cell::distance(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::distance(Position r, Angle a, int32_t on_surface) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; for (int32_t token : rpn) { // Ignore this token if it corresponds to an operator rather than a region. - if (token >= OP_UNION) {continue;} + if (token >= OP_UNION) continue; // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {token == on_surface}; - double d {surfaces_c[abs(token)-1]->distance(xyz, uvw, coincident)}; + double d {surfaces_c[abs(token)-1]->distance(r, a, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { @@ -356,8 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const //============================================================================== bool -Cell::contains_simple(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains_simple(Position r, Angle a, int32_t on_surface) const { for (int32_t token : rpn) { if (token < OP_UNION) { @@ -370,7 +367,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3], return false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw); + bool sense = surfaces_c[abs(token)-1]->sense(r, a); if (sense != (token > 0)) {return false;} } } @@ -381,8 +378,7 @@ Cell::contains_simple(const double xyz[3], const double uvw[3], //============================================================================== bool -Cell::contains_complex(const double xyz[3], const double uvw[3], - int32_t on_surface) const +Cell::contains_complex(Position r, Angle a, int32_t on_surface) const { // Make a stack of booleans. We don't know how big it needs to be, but we do // know that rpn.size() is an upper-bound. @@ -413,7 +409,7 @@ Cell::contains_complex(const double xyz[3], const double uvw[3], stack[i_stack] = false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(xyz, uvw);; + bool sense = surfaces_c[abs(token)-1]->sense(r, a);; stack[i_stack] = (sense == (token > 0)); } } @@ -494,12 +490,18 @@ extern "C" { bool cell_simple(Cell *c) {return c->simple;} bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface) - {return c->contains(xyz, uvw, on_surface);} + { + Position r {xyz}; + Angle a {uvw}; + return c->contains(r, a, on_surface); + } void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, double *min_dist, int32_t *i_surf) { - std::pair out = c->distance(xyz, uvw, on_surface); + Position r {xyz}; + Angle a {uvw}; + std::pair out = c->distance(r, a, on_surface); *min_dist = out.first; *i_surf = out.second; } diff --git a/src/cell.h b/src/cell.h index 24413b4fd4..08a6abf3f6 100644 --- a/src/cell.h +++ b/src/cell.h @@ -8,6 +8,7 @@ #include "hdf5.h" #include "pugixml.hpp" +#include "geometry.h" namespace openmc { @@ -97,21 +98,19 @@ public: //! known to be on. This index takes precedence over surface sense //! calculations. bool - contains(const double xyz[3], const double uvw[3], int32_t on_surface) const; + contains(Position r, Angle a, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair - distance(const double xyz[3], const double uvw[3], int32_t on_surface) const; + distance(Position r, Angle a, int32_t on_surface) const; //! \brief Write cell information to an HDF5 group. //! @param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; protected: - bool contains_simple(const double xyz[3], const double uvw[3], - int32_t on_surface) const; - bool contains_complex(const double xyz[3], const double uvw[3], - int32_t on_surface) const; + bool contains_simple(Position r, Angle a, int32_t on_surface) const; + bool contains_complex(Position r, Angle a, int32_t on_surface) const; }; } // namespace openmc diff --git a/src/geometry.cpp b/src/geometry.cpp new file mode 100644 index 0000000000..53bf32b5c1 --- /dev/null +++ b/src/geometry.cpp @@ -0,0 +1,59 @@ +#include "geometry.h" + +namespace openmc { + +Position& +Position::operator+=(Position other) +{ + x += other.x; + y += other.y; + z += other.z; + return *this; +} + +Position& +Position::operator+=(double v) +{ + x += v; + y += v; + z += v; + return *this; +} + +Position& +Position::operator-=(Position other) +{ + x -= other.x; + y -= other.y; + z -= other.z; + return *this; +} + +Position& +Position::operator-=(double v) +{ + x -= v; + y -= v; + z -= v; + return *this; +} + +Position& +Position::operator*=(Position other) +{ + x *= other.x; + y *= other.y; + z *= other.z; + return *this; +} + +Position& +Position::operator*=(double v) +{ + x *= v; + y *= v; + z *= v; + return *this; +} + +} // namespace openmc diff --git a/src/geometry.h b/src/geometry.h index 8837c961a7..792d653aee 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -1,6 +1,68 @@ #ifndef GEOMETRY_H #define GEOMETRY_H +namespace openmc { + extern "C" int openmc_root_universe; +struct Position { + double x = 0.; + double y = 0.; + double z = 0.; + + Position() = default; + Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; + Position(double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + + Position& operator+=(Position); + Position& operator+=(double); + Position& operator-=(Position); + Position& operator-=(double); + Position& operator*=(Position); + Position& operator*=(double); + const double& operator[](int i) const { + switch (i) { + case 0: return x; + case 1: return y; + case 2: return z; + } + } + double& operator[](int i) { + switch (i) { + case 0: return x; + case 1: return y; + case 2: return z; + } + } + + inline double dot(Position other) { + return x*other.x + y*other.y + z*other.z; + } +}; + +inline Position operator+(Position a, Position b) { return a += b; } +inline Position operator+(Position a, double b) { return a += b; } +inline Position operator+(double a, Position b) { return b += a; } + +inline Position operator-(Position a, Position b) { return a -= b; } +inline Position operator-(Position a, double b) { return a -= b; } +inline Position operator-(double a, Position b) { return b -= a; } + +inline Position operator*(Position a, Position b) { return a *= b; } +inline Position operator*(Position a, double b) { return a *= b; } +inline Position operator*(double a, Position b) { return b *= a; } + + +struct Angle : Position { + double& u() { return x; } + double& v() { return y; } + double& w() { return z; } + Angle() = default; + Angle(double u, double v, double w) : Position{u, v, w} { }; + Angle(double uvw[]) : Position{uvw} { }; + Angle(Position r) : Position{r} { }; +}; + +} // namespace openmc + #endif // GEOMETRY_H diff --git a/src/surface.cpp b/src/surface.cpp index 6400da2f77..aec0972b8f 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include "error.h" #include "hdf5_interface.h" @@ -177,38 +178,33 @@ Surface::Surface(pugi::xml_node surf_node) } bool -Surface::sense(const double xyz[3], const double uvw[3]) const +Surface::sense(Position r, Angle a) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. - const double f = evaluate(xyz); + const double f = evaluate(r); // Check which side of surface the point is on. if (std::abs(f) < FP_COINCIDENT) { // Particle may be coincident with this surface. To determine the sense, we // look at the direction of the particle relative to the surface normal (by // default in the positive direction) via their dot product. - double norm[3]; - normal(xyz, norm); - return uvw[0] * norm[0] + uvw[1] * norm[1] + uvw[2] * norm[2] > 0.0; + return a.dot(normal(r)) > 0.0; } return f > 0.0; } -void -Surface::reflect(const double xyz[3], double uvw[3]) const +Angle +Surface::reflect(Position r, Angle a) const { // Determine projection of direction onto normal and squared magnitude of // normal. - double norm[3]; - normal(xyz, norm); - const double projection = norm[0]*uvw[0] + norm[1]*uvw[1] + norm[2]*uvw[2]; - const double magnitude = norm[0]*norm[0] + norm[1]*norm[1] + norm[2]*norm[2]; + Angle n = normal(r); + const double projection = n.dot(a); + const double magnitude = n.dot(n); // Reflect direction according to normal. - uvw[0] -= 2.0 * projection / magnitude * norm[0]; - uvw[1] -= 2.0 * projection / magnitude * norm[1]; - uvw[2] -= 2.0 * projection / magnitude * norm[2]; + return a -= (2.0 * projection / magnitude) * n; } void @@ -261,33 +257,15 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) // The template parameter indicates the axis normal to the plane. template double -axis_aligned_plane_evaluate(const double xyz[3], double offset) +axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset) { - return xyz[i] - offset; -} - -// The template parameter indicates the axis normal to the plane. -template double -axis_aligned_plane_distance(const double xyz[3], const double uvw[3], - bool coincident, double offset) -{ - const double f = offset - xyz[i]; - if (coincident or std::abs(f) < FP_COINCIDENT or uvw[i] == 0.0) return INFTY; - const double d = f / uvw[i]; + const double f = offset - r[i]; + if (coincident or std::abs(f) < FP_COINCIDENT or a[i] == 0.0) return INFTY; + const double d = f / a[i]; if (d < 0.0) return INFTY; return d; } -// The first template parameter indicates the axis normal to the plane. The -// other two parameters indicate the other two axes. -template void -axis_aligned_plane_normal(const double xyz[3], double uvw[3]) -{ - uvw[i1] = 1.0; - uvw[i2] = 0.0; - uvw[i3] = 0.0; -} - //============================================================================== // SurfaceXPlane implementation //============================================================================== @@ -298,20 +276,19 @@ SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0); } -inline double SurfaceXPlane::evaluate(const double xyz[3]) const +double SurfaceXPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<0>(xyz, x0); + return r.x - x0; } -inline double SurfaceXPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<0>(xyz, uvw, coincident, x0); + return axis_aligned_plane_distance<0>(r, a, coincident, x0); } -inline void SurfaceXPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXPlane::normal(Position r) const { - axis_aligned_plane_normal<0, 1, 2>(xyz, uvw); + return {1., 0., 0.}; } void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const @@ -321,26 +298,24 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { - double other_norm[3]; - other->normal(xyz, other_norm); - if (other_norm[0] == 1 and other_norm[1] == 0 and other_norm[2] == 0) { - xyz[0] = x0; + Angle other_n = other->normal(r); + if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 0) { + r.x = x0; return false; } else { // Assume the partner is an YPlane (the only supported partner). Use the - // evaluate function to find y0, then adjust xyz and uvw for rotational + // evaluate function to find y0, then adjust position/angle for rotational // symmetry. - double xyz_test[3] {0, 0, 0}; - double y0 = -other->evaluate(xyz_test); - xyz[1] = xyz[0] - x0 + y0; - xyz[0] = x0; + double y0 = -other->evaluate({0., 0., 0.}); + r.y = r.x - x0 + y0; + r.x = x0; - double u = uvw[0]; - uvw[0] = -uvw[1]; - uvw[1] = u; + double u = a.u(); + a.u() = -a.v(); + a.v() = u; return true; } @@ -349,8 +324,7 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfaceXPlane::bounding_box() const { - BoundingBox out {x0, x0, -INFTY, INFTY, -INFTY, INFTY}; - return out; + return {x0, x0, -INFTY, INFTY, -INFTY, INFTY}; } //============================================================================== @@ -363,20 +337,19 @@ SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, y0); } -inline double SurfaceYPlane::evaluate(const double xyz[3]) const +double SurfaceYPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<1>(xyz, y0); + return r.y - y0; } -inline double SurfaceYPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<1>(xyz, uvw, coincident, y0); + return axis_aligned_plane_distance<1>(r, a, coincident, y0); } -inline void SurfaceYPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYPlane::normal(Position r) const { - axis_aligned_plane_normal<1, 0, 2>(xyz, uvw); + return {0., 1., 0.}; } void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const @@ -386,27 +359,25 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { - double other_norm[3]; - other->normal(xyz, other_norm); - if (other_norm[0] == 0 and other_norm[1] == 1 and other_norm[2] == 0) { + Angle other_n = other->normal(r); + if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) { // The periodic partner is also aligned along y. Just change the y coord. - xyz[1] = y0; + r.y = y0; return false; } else { // Assume the partner is an XPlane (the only supported partner). Use the - // evaluate function to find x0, then adjust xyz and uvw for rotational + // evaluate function to find x0, then adjust position/angle for rotational // symmetry. - double xyz_test[3] {0, 0, 0}; - double x0 = -other->evaluate(xyz_test); - xyz[0] = xyz[1] - y0 + x0; - xyz[1] = y0; + double x0 = -other->evaluate({0., 0., 0.}); + r.x = r.y - y0 + x0; + r.y = y0; - double u = uvw[0]; - uvw[0] = uvw[1]; - uvw[1] = -u; + double u = a.u(); + a.u() = a.v(); + a.v() = -u; return true; } @@ -415,8 +386,7 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfaceYPlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, y0, y0, -INFTY, INFTY}; - return out; + return {-INFTY, INFTY, y0, y0, -INFTY, INFTY}; } //============================================================================== @@ -429,20 +399,19 @@ SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) read_coeffs(surf_node, id, z0); } -inline double SurfaceZPlane::evaluate(const double xyz[3]) const +double SurfaceZPlane::evaluate(Position r) const { - return axis_aligned_plane_evaluate<2>(xyz, z0); + return r.z - z0; } -inline double SurfaceZPlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_plane_distance<2>(xyz, uvw, coincident, z0); + return axis_aligned_plane_distance<2>(r, a, coincident, z0); } -inline void SurfaceZPlane::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZPlane::normal(Position r) const { - axis_aligned_plane_normal<2, 0, 1>(xyz, uvw); + return {0., 0., 1.}; } void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const @@ -452,19 +421,18 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { // Assume the other plane is aligned along z. Just change the z coord. - xyz[2] = z0; + r.z = z0; return false; } BoundingBox SurfaceZPlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, z0, z0}; - return out; + return {-INFTY, INFTY, -INFTY, INFTY, z0, z0}; } //============================================================================== @@ -478,17 +446,16 @@ SurfacePlane::SurfacePlane(pugi::xml_node surf_node) } double -SurfacePlane::evaluate(const double xyz[3]) const +SurfacePlane::evaluate(Position r) const { - return A*xyz[0] + B*xyz[1] + C*xyz[2] - D; + return A*r.x + B*r.y + C*r.z - D; } double -SurfacePlane::distance(const double xyz[3], const double uvw[3], - bool coincident) const +SurfacePlane::distance(Position r, Angle a, bool coincident) const { - const double f = A*xyz[0] + B*xyz[1] + C*xyz[2] - D; - const double projection = A*uvw[0] + B*uvw[1] + C*uvw[2]; + const double f = A*r.x + B*r.y + C*r.z - D; + const double projection = A*a.u() + B*a.v() + C*a.w(); if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) { return INFTY; } else { @@ -498,12 +465,10 @@ SurfacePlane::distance(const double xyz[3], const double uvw[3], } } -void -SurfacePlane::normal(const double xyz[3], double uvw[3]) const +Angle +SurfacePlane::normal(Position r) const { - uvw[0] = A; - uvw[1] = B; - uvw[2] = C; + return {A, B, C}; } void SurfacePlane::to_hdf5_inner(hid_t group_id) const @@ -513,18 +478,18 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const +bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const { // This function assumes the other plane shares this plane's normal direction. // Determine the distance to intersection. - double d = evaluate(xyz) / (A*A + B*B + C*C); + double d = evaluate(r) / (A*A + B*B + C*C); // Move the particle that distance along the normal vector. - xyz[0] -= d * A; - xyz[1] -= d * B; - xyz[2] -= d * C; + r.x -= d * A; + r.y -= d * B; + r.z -= d * C; return false; } @@ -532,8 +497,7 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], BoundingBox SurfacePlane::bounding_box() const { - BoundingBox out {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; - return out; + return {-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; } //============================================================================== @@ -544,28 +508,28 @@ SurfacePlane::bounding_box() const // cylinder. offset1 and offset2 should correspond with i1 and i2, // respectively. template double -axis_aligned_cylinder_evaluate(const double xyz[3], double offset1, +axis_aligned_cylinder_evaluate(Position r, double offset1, double offset2, double radius) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - return xyz1*xyz1 + xyz2*xyz2 - radius*radius; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + return r1*r1 + r2*r2 - radius*radius; } // The first template parameter indicates which axis the cylinder is aligned to. // The other two parameters indicate the other two axes. offset1 and offset2 // should correspond with i2 and i3, respectively. template double -axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], +axis_aligned_cylinder_distance(Position r, Angle u, bool coincident, double offset1, double offset2, double radius) { - const double a = 1.0 - uvw[i1]*uvw[i1]; // u^2 + v^2 + const double a = 1.0 - u[i1]*u[i1]; // u^2 + v^2 if (a == 0.0) return INFTY; - const double xyz2 = xyz[i2] - offset1; - const double xyz3 = xyz[i3] - offset2; - const double k = xyz2 * uvw[i2] + xyz3 * uvw[i3]; - const double c = xyz2*xyz2 + xyz3*xyz3 - radius*radius; + const double r2 = r[i2] - offset1; + const double r3 = r[i3] - offset2; + const double k = r2 * u[i2] + r3 * u[i3]; + const double c = r2*r2 + r3*r3 - radius*radius; const double quad = k*k - a*c; if (quad < 0.0) { @@ -601,13 +565,14 @@ axis_aligned_cylinder_distance(const double xyz[3], const double uvw[3], // The first template parameter indicates which axis the cylinder is aligned to. // The other two parameters indicate the other two axes. offset1 and offset2 // should correspond with i2 and i3, respectively. -template void -axis_aligned_cylinder_normal(const double xyz[3], double uvw[3], double offset1, - double offset2) +template Angle +axis_aligned_cylinder_normal(Position r, double offset1, double offset2) { - uvw[i2] = 2.0 * (xyz[i2] - offset1); - uvw[i3] = 2.0 * (xyz[i3] - offset2); - uvw[i1] = 0.0; + Angle a; + a[i2] = 2.0 * (r[i2] - offset1); + a[i3] = 2.0 * (r[i3] - offset2); + a[i1] = 0.0; + return a; } //============================================================================== @@ -620,21 +585,20 @@ SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, y0, z0, r); } -inline double SurfaceXCylinder::evaluate(const double xyz[3]) const +double SurfaceXCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<1, 2>(xyz, y0, z0, r); + return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r); } -inline double SurfaceXCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<0, 1, 2>(xyz, uvw, coincident, y0, z0, - r); + return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0, + this->r); } -inline void SurfaceXCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<0, 1, 2>(xyz, uvw, y0, z0); + return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0); } @@ -655,21 +619,20 @@ SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, z0, r); } -inline double SurfaceYCylinder::evaluate(const double xyz[3]) const +double SurfaceYCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 2>(xyz, x0, z0, r); + return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r); } -inline double SurfaceYCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<1, 0, 2>(xyz, uvw, coincident, x0, z0, - r); + return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0, + this->r); } -inline void SurfaceYCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<1, 0, 2>(xyz, uvw, x0, z0); + return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0); } void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const @@ -689,21 +652,20 @@ SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, r); } -inline double SurfaceZCylinder::evaluate(const double xyz[3]) const +double SurfaceZCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 1>(xyz, x0, y0, r); + return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r); } -inline double SurfaceZCylinder::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cylinder_distance<2, 0, 1>(xyz, uvw, coincident, x0, y0, - r); + return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0, + this->r); } -inline void SurfaceZCylinder::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZCylinder::normal(Position r) const { - axis_aligned_cylinder_normal<2, 0, 1>(xyz, uvw, x0, y0); + return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0); } void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const @@ -723,22 +685,21 @@ SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r); } -double SurfaceSphere::evaluate(const double xyz[3]) const +double SurfaceSphere::evaluate(Position r) const { - const double x = xyz[0] - x0; - const double y = xyz[1] - y0; - const double z = xyz[2] - z0; - return x*x + y*y + z*z - r*r; + const double x = r.x - x0; + const double y = r.y - y0; + const double z = r.z - z0; + return x*x + y*y + z*z - this->r*this->r; } -double SurfaceSphere::distance(const double xyz[3], const double uvw[3], - bool coincident) const +double SurfaceSphere::distance(Position r, Angle a, bool coincident) const { - const double x = xyz[0] - x0; - const double y = xyz[1] - y0; - const double z = xyz[2] - z0; - const double k = x*uvw[0] + y*uvw[1] + z*uvw[2]; - const double c = x*x + y*y + z*z - r*r; + const double x = r.x - x0; + const double y = r.y - y0; + const double z = r.z - z0; + const double k = x*a.u() + y*a.v() + z*a.w(); + const double c = x*x + y*y + z*z - this->r*this->r; const double quad = k*k - c; if (quad < 0.0) { @@ -770,11 +731,9 @@ double SurfaceSphere::distance(const double xyz[3], const double uvw[3], } } -inline void SurfaceSphere::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceSphere::normal(Position r) const { - uvw[0] = 2.0 * (xyz[0] - x0); - uvw[1] = 2.0 * (xyz[1] - y0); - uvw[2] = 2.0 * (xyz[2] - z0); + return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)}; } void SurfaceSphere::to_hdf5_inner(hid_t group_id) const @@ -792,30 +751,30 @@ void SurfaceSphere::to_hdf5_inner(hid_t group_id) const // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. template double -axis_aligned_cone_evaluate(const double xyz[3], double offset1, +axis_aligned_cone_evaluate(Position r, double offset1, double offset2, double offset3, double radius_sq) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - const double xyz3 = xyz[i3] - offset3; - return xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + const double r3 = r[i3] - offset3; + return r2*r2 + r3*r3 - radius_sq*r1*r1; } // The first template parameter indicates which axis the cone is aligned to. // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. template double -axis_aligned_cone_distance(const double xyz[3], const double uvw[3], +axis_aligned_cone_distance(Position r, Angle u, bool coincident, double offset1, double offset2, double offset3, double radius_sq) { - const double xyz1 = xyz[i1] - offset1; - const double xyz2 = xyz[i2] - offset2; - const double xyz3 = xyz[i3] - offset3; - const double a = uvw[i2]*uvw[i2] + uvw[i3]*uvw[i3] - - radius_sq*uvw[i1]*uvw[i1]; - const double k = xyz2*uvw[i2] + xyz3*uvw[i3] - radius_sq*xyz1*uvw[i1]; - const double c = xyz2*xyz2 + xyz3*xyz3 - radius_sq*xyz1*xyz1; + const double r1 = r[i1] - offset1; + const double r2 = r[i2] - offset2; + const double r3 = r[i3] - offset3; + const double a = u[i2]*u[i2] + u[i3]*u[i3] + - radius_sq*u[i1]*u[i1]; + const double k = r2*u[i2] + r3*u[i3] - radius_sq*r1*u[i1]; + const double c = r2*r2 + r3*r3 - radius_sq*r1*r1; double quad = k*k - a*c; double d; @@ -858,13 +817,15 @@ axis_aligned_cone_distance(const double xyz[3], const double uvw[3], // The first template parameter indicates which axis the cone is aligned to. // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. -template void -axis_aligned_cone_normal(const double xyz[3], double uvw[3], double offset1, - double offset2, double offset3, double radius_sq) +template Angle +axis_aligned_cone_normal(Position r, double offset1, double offset2, + double offset3, double radius_sq) { - uvw[i1] = -2.0 * radius_sq * (xyz[i1] - offset1); - uvw[i2] = 2.0 * (xyz[i2] - offset2); - uvw[i3] = 2.0 * (xyz[i3] - offset3); + Angle a; + a[i1] = -2.0 * radius_sq * (r[i1] - offset1); + a[i2] = 2.0 * (r[i2] - offset2); + a[i3] = 2.0 * (r[i3] - offset3); + return a; } //============================================================================== @@ -877,21 +838,20 @@ SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceXCone::evaluate(const double xyz[3]) const +double SurfaceXCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<0, 1, 2>(xyz, x0, y0, z0, r_sq); + return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq); } -inline double SurfaceXCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceXCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<0, 1, 2>(xyz, uvw, coincident, x0, y0, z0, + return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0, r_sq); } -inline void SurfaceXCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceXCone::normal(Position r) const { - axis_aligned_cone_normal<0, 1, 2>(xyz, uvw, x0, y0, z0, r_sq); + return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq); } void SurfaceXCone::to_hdf5_inner(hid_t group_id) const @@ -911,21 +871,20 @@ SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceYCone::evaluate(const double xyz[3]) const +double SurfaceYCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<1, 0, 2>(xyz, y0, x0, z0, r_sq); + return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq); } -inline double SurfaceYCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceYCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<1, 0, 2>(xyz, uvw, coincident, y0, x0, z0, + return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0, r_sq); } -inline void SurfaceYCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceYCone::normal(Position r) const { - axis_aligned_cone_normal<1, 0, 2>(xyz, uvw, y0, x0, z0, r_sq); + return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq); } void SurfaceYCone::to_hdf5_inner(hid_t group_id) const @@ -945,21 +904,20 @@ SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) read_coeffs(surf_node, id, x0, y0, z0, r_sq); } -inline double SurfaceZCone::evaluate(const double xyz[3]) const +double SurfaceZCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<2, 0, 1>(xyz, z0, x0, y0, r_sq); + return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq); } -inline double SurfaceZCone::distance(const double xyz[3], - const double uvw[3], bool coincident) const +double SurfaceZCone::distance(Position r, Angle a, bool coincident) const { - return axis_aligned_cone_distance<2, 0, 1>(xyz, uvw, coincident, z0, x0, y0, + return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0, r_sq); } -inline void SurfaceZCone::normal(const double xyz[3], double uvw[3]) const +Angle SurfaceZCone::normal(Position r) const { - axis_aligned_cone_normal<2, 0, 1>(xyz, uvw, z0, x0, y0, r_sq); + return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq); } void SurfaceZCone::to_hdf5_inner(hid_t group_id) const @@ -980,26 +938,25 @@ SurfaceQuadric::SurfaceQuadric(pugi::xml_node surf_node) } double -SurfaceQuadric::evaluate(const double xyz[3]) const +SurfaceQuadric::evaluate(Position r) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; + const double x = r.x; + const double y = r.y; + const double z = r.z; return x*(A*x + D*y + G) + y*(B*y + E*z + H) + z*(C*z + F*x + J) + K; } double -SurfaceQuadric::distance(const double xyz[3], - const double uvw[3], bool coincident) const +SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; - const double &u = uvw[0]; - const double &v = uvw[1]; - const double &w = uvw[2]; + const double &x = r.x; + const double &y = r.y; + const double &z = r.z; + const double &u = ang.u(); + const double &v = ang.v(); + const double &w = ang.w(); const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w; const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) + @@ -1045,15 +1002,15 @@ SurfaceQuadric::distance(const double xyz[3], return d; } -void -SurfaceQuadric::normal(const double xyz[3], double uvw[3]) const +Angle +SurfaceQuadric::normal(Position r) const { - const double &x = xyz[0]; - const double &y = xyz[1]; - const double &z = xyz[2]; - uvw[0] = 2.0*A*x + D*y + F*z + G; - uvw[1] = 2.0*B*y + D*x + E*z + H; - uvw[2] = 2.0*C*z + E*y + F*x + J; + const double &x = r.x; + const double &y = r.y; + const double &z = r.z; + return {2.0*A*x + D*y + F*z + G, + 2.0*B*y + D*x + E*z + H, + 2.0*C*z + E*y + F*x + J}; } void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const @@ -1263,10 +1220,24 @@ extern "C" { int surface_bc(Surface *surf) {return surf->bc;} void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) - {surf->reflect(xyz, uvw);} + { + Position r {xyz}; + Angle a {uvw}; + a = surf->reflect(r, a); + + uvw[0] = a.u(); + uvw[1] = a.v(); + uvw[2] = a.w(); + } void surface_normal(Surface *surf, double xyz[3], double uvw[3]) - {return surf->normal(xyz, uvw);} + { + Position r {xyz}; + Angle a = surf->normal(r); + uvw[0] = a.u(); + uvw[1] = a.v(); + uvw[2] = a.w(); + } void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} @@ -1275,7 +1246,21 @@ extern "C" { bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3], double uvw[3]) - {return surf->periodic_translate(other, xyz, uvw);} + { + Position r {xyz}; + Angle a {uvw}; + bool rotational = surf->periodic_translate(other, r, a); + + // Copy back to arrays + xyz[0] = r.x; + xyz[1] = r.y; + xyz[2] = r.z; + uvw[0] = a.u(); + uvw[1] = a.v(); + uvw[2] = a.w(); + + return rotational; + } void free_memory_surfaces_c() { diff --git a/src/surface.h b/src/surface.h index 666fc0dbae..be3ad52ab9 100644 --- a/src/surface.h +++ b/src/surface.h @@ -9,6 +9,7 @@ #include "pugixml.hpp" #include "constants.h" +#include "geometry.h" namespace openmc { @@ -65,39 +66,38 @@ public: virtual ~Surface() {} //! Determine which side of a surface a point lies on. - //! @param xyz[3] The 3D Cartesian coordinate of a point. - //! @param uvw[3] A direction used to "break ties" and pick a sense when the + //! @param r The 3D Cartesian coordinate of a point. + //! @param o A direction used to "break ties" and pick a sense when the //! point is very close to the surface. //! @return true if the point is on the "positive" side of the surface and //! false otherwise. - bool sense(const double xyz[3], const double uvw[3]) const; + bool sense(Position r, Angle a) const; //! Determine the direction of a ray reflected from the surface. - //! @param xyz[3] The point at which the ray is incident. - //! @param uvw[3] A direction. This is both an input and an output parameter. + //! @param r The point at which the ray is incident. + //! @param o A direction. This is both an input and an output parameter. //! It specifies the icident direction on input and the reflected direction //! on output. - void reflect(const double xyz[3], double uvw[3]) const; + Angle reflect(Position r, Angle a) const; //! Evaluate the equation describing the surface. //! //! Surfaces can be described by some function f(x, y, z) = 0. This member //! function evaluates that mathematical function. - //! @param xyz[3] A 3D Cartesian coordinate. - virtual double evaluate(const double xyz[3]) const = 0; + //! @param r A 3D Cartesian coordinate. + virtual double evaluate(Position r) const = 0; //! Compute the distance between a point and the surface along a ray. - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] The direction of the ray. + //! @param r A 3D Cartesian coordinate. + //! @param o The direction of the ray. //! @param coincident A hint to the code that the given point should lie //! exactly on the surface. - virtual double distance(const double xyz[3], const double uvw[3], - bool coincident) const = 0; + virtual double distance(Position r, Angle a, bool coincident) const = 0; //! Compute the local outward normal direction of the surface. - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] This output argument provides the normal. - virtual void normal(const double xyz[3], double uvw[3]) const = 0; + //! @param r A 3D Cartesian coordinate. + //! @return Normal direction + virtual Angle normal(Position r) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. //! @param group_id An HDF5 group id. @@ -125,14 +125,14 @@ public: //! Translate a particle onto this surface from a periodic partner surface. //! @param other A pointer to the partner surface in this periodic BC. - //! @param xyz[3] A point on the partner surface that will be translated onto + //! @param r A point on the partner surface that will be translated onto //! this surface. - //! @param uvw[3] A direction that will be rotated for systems with rotational + //! @param a A direction that will be rotated for systems with rotational //! periodicity. //! @return true if this surface and its partner make a rotationally-periodic //! boundary condition. - virtual bool periodic_translate(PeriodicSurface *other, double xyz[3], - double uvw[3]) const = 0; + virtual bool periodic_translate(const PeriodicSurface *other, Position& r, + Angle& a) const = 0; //! Get the bounding box for this surface. virtual BoundingBox bounding_box() const = 0; @@ -149,12 +149,11 @@ class SurfaceXPlane : public PeriodicSurface double x0; public: explicit SurfaceXPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], bool coincident) - const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -170,12 +169,11 @@ class SurfaceYPlane : public PeriodicSurface double y0; public: explicit SurfaceYPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -191,12 +189,11 @@ class SurfaceZPlane : public PeriodicSurface double z0; public: explicit SurfaceZPlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -212,12 +209,11 @@ class SurfacePlane : public PeriodicSurface double A, B, C, D; public: explicit SurfacePlane(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], bool coincident) - const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) + bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) const; BoundingBox bounding_box() const; }; @@ -234,10 +230,9 @@ class SurfaceXCylinder : public Surface double y0, z0, r; public: explicit SurfaceXCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -253,10 +248,9 @@ class SurfaceYCylinder : public Surface double x0, z0, r; public: explicit SurfaceYCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -272,10 +266,9 @@ class SurfaceZCylinder : public Surface double x0, y0, r; public: explicit SurfaceZCylinder(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -291,10 +284,9 @@ class SurfaceSphere : public Surface double x0, y0, z0, r; public: explicit SurfaceSphere(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -310,10 +302,9 @@ class SurfaceXCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceXCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -329,10 +320,9 @@ class SurfaceYCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceYCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -348,10 +338,9 @@ class SurfaceZCone : public Surface double x0, y0, z0, r_sq; public: explicit SurfaceZCone(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -367,10 +356,9 @@ class SurfaceQuadric : public Surface double A, B, C, D, E, F, G, H, J, K; public: explicit SurfaceQuadric(pugi::xml_node surf_node); - double evaluate(const double xyz[3]) const; - double distance(const double xyz[3], const double uvw[3], - bool coincident) const; - void normal(const double xyz[3], double uvw[3]) const; + double evaluate(Position r) const; + double distance(Position r, Angle a, bool coincident) const; + Angle normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; From f7649a61dc0a67d4d4d48b4e544cd0f0c6aca1a7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 26 Jun 2018 14:03:54 -0500 Subject: [PATCH 25/36] Use Position/Angle in lattice classes --- src/geometry.h | 4 +- src/hdf5_interface.h | 9 +++ src/lattice.cpp | 130 ++++++++++++++++++++----------------------- src/lattice.h | 41 +++++++------- 4 files changed, 93 insertions(+), 91 deletions(-) diff --git a/src/geometry.h b/src/geometry.h index 792d653aee..e069b9315b 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -12,7 +12,7 @@ struct Position { Position() = default; Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; - Position(double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; Position& operator+=(Position); Position& operator+=(double); @@ -59,7 +59,7 @@ struct Angle : Position { double& w() { return z; } Angle() = default; Angle(double u, double v, double w) : Position{u, v, w} { }; - Angle(double uvw[]) : Position{uvw} { }; + Angle(const double uvw[]) : Position{uvw} { }; Angle(Position r) : Position{r} { }; }; diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index b0d2cc2c19..8b1f7a434d 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -11,6 +11,8 @@ #include #include +#include "geometry.h" + namespace openmc { @@ -181,5 +183,12 @@ write_dataset(hid_t obj_id, const char* name, const std::array& buffer) write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } +template<> inline void +write_dataset(hid_t obj_id, const char* name, Position r) +{ + std::array buffer {r.x, r.y, r.z}; + write_dataset(obj_id, name, buffer); +} + } // namespace openmc #endif //HDF5_INTERFACE_H diff --git a/src/lattice.cpp b/src/lattice.cpp index c95762a63d..c67fc419f0 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -223,15 +223,14 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -RectLattice::distance(const double xyz[3], const double uvw[3], - const int i_xyz[3]) const +RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const { // Get short aliases to the coordinates. - double x {xyz[0]}; - double y {xyz[1]}; - double z {xyz[2]}; - double u {uvw[0]}; - double v {uvw[1]}; + double x = r.x; + double y = r.y; + double z = r.z; + double u = a.u(); + double v = a.v(); // Determine the oncoming edge. double x0 {copysign(0.5 * pitch[0], u)}; @@ -264,7 +263,7 @@ RectLattice::distance(const double xyz[3], const double uvw[3], // Top and bottom sides if (is_3d) { - double w {uvw[2]}; + double w {a.w()}; double z0 {copysign(0.5 * pitch[2], w)}; if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { double this_d = (z0 - z) / w; @@ -285,13 +284,13 @@ RectLattice::distance(const double xyz[3], const double uvw[3], //============================================================================== std::array -RectLattice::get_indices(const double xyz[3]) const +RectLattice::get_indices(Position r) const { - int ix {static_cast(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1}; - int iy {static_cast(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1}; + int ix {static_cast(std::ceil((r.x - lower_left.x) / pitch.x))-1}; + int iy {static_cast(std::ceil((r.y - lower_left.y) / pitch.y))-1}; int iz; if (is_3d) { - iz = static_cast(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1; + iz = static_cast(std::ceil((r.z - lower_left.z) / pitch.z))-1; } else { iz = 0; } @@ -300,18 +299,15 @@ RectLattice::get_indices(const double xyz[3]) const //============================================================================== -std::array -RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +Position +RectLattice::get_local_position(Position r, const int i_xyz[3]) const { - std::array local_xyz; - local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]); - local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]); + r.x -= (lower_left.x + (i_xyz[0] + 0.5)*pitch.x); + r.y -= (lower_left.y + (i_xyz[1] + 0.5)*pitch.y); if (is_3d) { - local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]); - } else { - local_xyz[2] = global_xyz[2]; + r.z -= (lower_left.z + (i_xyz[2] + 0.5)*pitch.z); } - return local_xyz; + return r; } //============================================================================== @@ -583,12 +579,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -HexLattice::distance(const double xyz[3], const double uvw[3], - const int i_xyz[3]) const +HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const { // Compute the direction on the hexagonal basis. - double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0; - double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0; + double beta_dir = a.u() * std::sqrt(3.0) / 2.0 + a.v() / 2.0; + double gamma_dir = a.u() * std::sqrt(3.0) / 2.0 - a.v() / 2.0; // Note that hexagonal lattice distance calculations are performed // using the particle's coordinates relative to the neighbor lattice @@ -600,15 +595,15 @@ HexLattice::distance(const double xyz[3], const double uvw[3], double d {INFTY}; std::array lattice_trans; double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge - std::array xyz_t; + Position r_t; if (beta_dir > 0) { const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0; + double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { d = (edge - beta) / beta_dir; if (beta_dir > 0) { @@ -622,12 +617,12 @@ HexLattice::distance(const double xyz[3], const double uvw[3], edge = -copysign(0.5*pitch[0], gamma_dir); if (gamma_dir > 0) { const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0; + double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { double this_d = (edge - gamma) / gamma_dir; if (this_d < d) { @@ -641,18 +636,18 @@ HexLattice::distance(const double xyz[3], const double uvw[3], } // Upper and lower sides. - edge = -copysign(0.5*pitch[0], uvw[1]); - if (uvw[1] > 0) { + edge = -copysign(0.5*pitch[0], a.v()); + if (a.v() > 0) { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) { - double this_d = (edge - xyz_t[1]) / uvw[1]; + if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) { + double this_d = (edge - r_t.y) / a.v(); if (this_d < d) { - if (uvw[1] > 0) { + if (a.v() > 0) { lattice_trans = {0, 1, 0}; } else { lattice_trans = {0, -1, 0}; @@ -663,8 +658,8 @@ HexLattice::distance(const double xyz[3], const double uvw[3], // Top and bottom sides if (is_3d) { - double z {xyz[2]}; - double w {uvw[2]}; + double z = r.z; + double w = a.w(); double z0 {copysign(0.5 * pitch[1], w)}; if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { double this_d = (z0 - z) / w; @@ -686,24 +681,24 @@ HexLattice::distance(const double xyz[3], const double uvw[3], //============================================================================== std::array -HexLattice::get_indices(const double xyz[3]) const +HexLattice::get_indices(Position r) const { // Offset the xyz by the lattice center. - double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]}; - if (is_3d) {xyz_o[2] -= center[2];} + Position r_o {r.x - center.x, r.y - center.y, r.z}; + if (is_3d) {r_o.z -= center.z;} // Index the z direction. std::array out; if (is_3d) { - out[2] = static_cast(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial))-1; + out[2] = static_cast(std::ceil(r_o.z / pitch[1] + 0.5 * n_axial))-1; } else { out[2] = 0; } // Convert coordinates into skewed bases. The (x, alpha) basis is used to // find the index of the global coordinates to within 4 cells. - double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0); - out[0] = static_cast(std::floor(xyz_o[0] + double alpha = r_o.y - r_o.x / std::sqrt(3.0); + out[0] = static_cast(std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch[0]))); out[1] = static_cast(std::floor(alpha / pitch[0])); @@ -725,8 +720,8 @@ HexLattice::get_indices(const double xyz[3]) const for (int i = 0; i < 2; i++) { for (int j = 0; j < 2; j++) { int i_xyz[3] {out[0] + j, out[1] + i, 0}; - std::array xyz_t = get_local_xyz(xyz, i_xyz); - double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1]; + Position r_t = get_local_position(r, i_xyz); + double d = r_t.x*r_t.x + r_t.y*r_t.y; if (d < d_min) { d_min = d; k_min = k; @@ -751,26 +746,19 @@ HexLattice::get_indices(const double xyz[3]) const //============================================================================== -std::array -HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +Position +HexLattice::get_local_position(Position r, const int i_xyz[3]) const { - std::array local_xyz; - // x_l = x_g - (center + pitch_x*cos(30)*index_x) - local_xyz[0] = global_xyz[0] - (center[0] - + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); + r.x -= (center.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - local_xyz[1] = global_xyz[1] - (center[1] - + (i_xyz[1] - n_rings + 1) * pitch[0] - + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); + r.y -= (center.y + (i_xyz[1] - n_rings + 1) * pitch[0] + + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); if (is_3d) { - local_xyz[2] = global_xyz[2] - center[2] - + (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; - } else { - local_xyz[2] = global_xyz[2]; + r.z -= center.z - (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; } - return local_xyz; + return r; } //============================================================================== @@ -908,7 +896,9 @@ extern "C" { void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], const int i_xyz[3], double *d, int lattice_trans[3]) { - std::pair> ld {lat->distance(xyz, uvw, i_xyz)}; + Position r {xyz}; + Angle a {uvw}; + std::pair> ld {lat->distance(r, a, i_xyz)}; *d = ld.first; lattice_trans[0] = ld.second[0]; lattice_trans[1] = ld.second[1]; @@ -917,7 +907,8 @@ extern "C" { void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) { - std::array inds = lat->get_indices(xyz); + Position r {xyz}; + std::array inds = lat->get_indices(r); i_xyz[0] = inds[0]; i_xyz[1] = inds[1]; i_xyz[2] = inds[2]; @@ -926,10 +917,11 @@ extern "C" { void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3], const int i_xyz[3], double local_xyz[3]) { - std::array xyz = lat->get_local_xyz(global_xyz, i_xyz); - local_xyz[0] = xyz[0]; - local_xyz[1] = xyz[1]; - local_xyz[2] = xyz[2]; + Position global {global_xyz}; + Position local = lat->get_local_position(global, i_xyz); + local_xyz[0] = local.x; + local_xyz[1] = local.y; + local_xyz[2] = local.z; } int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) diff --git a/src/lattice.h b/src/lattice.h index bc706a5a3f..a6c4496621 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -8,6 +8,7 @@ #include #include "constants.h" +#include "geometry.h" #include "hdf5.h" #include "pugixml.hpp" @@ -75,26 +76,26 @@ public: virtual bool are_valid_indices(const int i_xyz[3]) const = 0; //! \brief Find the next lattice surface crossing - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] A 3D Cartesian direction. + //! @param r A 3D Cartesian coordinate. + //! @param a A 3D Cartesian direction. //! @param i_xyz[3] The indices for a lattice tile. //! @return The distance to the next crossing and an array indicating how the //! lattice indices would change after crossing that boundary. virtual std::pair> - distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const + distance(Position r, Angle a, const int i_xyz[3]) const = 0; //! \brief Find the lattice tile indices for a given point. - //! @param xyz[3] A 3D Cartesian coordinate. + //! @param r A 3D Cartesian coordinate. //! @return An array containing the indices of a lattice tile. - virtual std::array get_indices(const double xyz[3]) const = 0; + virtual std::array get_indices(Position r) const = 0; //! \brief Get coordinates local to a lattice tile. - //! @param global_xyz[3] A 3D Cartesian coordinate. + //! @param r A 3D Cartesian coordinate. //! @param i_xyz[3] The indices for a lattice tile. //! @return Local 3D Cartesian coordinates. - virtual std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0; + virtual Position + get_local_position(Position r, const int i_xyz[3]) const = 0; //! \brief Check flattened lattice index. //! @param indx The index for a lattice tile. @@ -193,12 +194,12 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + distance(Position r, Angle a, const int i_xyz[3]) const; - std::array get_indices(const double xyz[3]) const; + std::array get_indices(Position r) const; - std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + Position + get_local_position(Position r, const int i_xyz[3]) const; int32_t& offset(int map, const int i_xyz[3]); @@ -207,9 +208,9 @@ public: void to_hdf5_inner(hid_t group_id) const; private: - std::array n_cells; //!< Number of cells along each axis - std::array lower_left; //!< Global lower-left corner of the lattice - std::array pitch; //!< Lattice tile width along each axis + std::array n_cells; //!< Number of cells along each axis + Position lower_left; //!< Global lower-left corner of the lattice + Position pitch; //!< Lattice tile width along each axis // Convenience aliases int &nx {n_cells[0]}; @@ -233,12 +234,12 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + distance(Position r, Angle a, const int i_xyz[3]) const; - std::array get_indices(const double xyz[3]) const; + std::array get_indices(Position r) const; - std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + Position + get_local_position(Position r, const int i_xyz[3]) const; bool is_valid_index(int indx) const; @@ -251,7 +252,7 @@ public: private: int n_rings; //!< Number of radial tile positions int n_axial; //!< Number of axial tile positions - std::array center; //!< Global center of lattice + Position center; //!< Global center of lattice std::array pitch; //!< Lattice tile width and height }; From 0fba0e287c3c75d4c49bc93f9ec2d1b89764cee8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Aug 2018 13:36:42 -0500 Subject: [PATCH 26/36] Settle on Position/Direction, declared in position.h --- CMakeLists.txt | 2 +- src/cell.cpp | 26 ++-- src/cell.h | 11 +- src/geometry.h | 64 +--------- src/hdf5_interface.h | 2 +- src/lattice.cpp | 58 +++++---- src/lattice.h | 11 +- src/{geometry.cpp => position.cpp} | 6 +- src/position.h | 74 ++++++++++++ src/surface.cpp | 184 ++++++++++++++--------------- src/surface.h | 94 +++++++-------- 11 files changed, 275 insertions(+), 257 deletions(-) rename src/{geometry.cpp => position.cpp} (75%) create mode 100644 src/position.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 0be9ccb325..c94f4df5fd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -386,7 +386,6 @@ add_library(libopenmc SHARED src/cell.cpp src/initialize.cpp src/finalize.cpp - src/geometry.cpp src/geometry_aux.cpp src/hdf5_interface.cpp src/lattice.cpp @@ -396,6 +395,7 @@ add_library(libopenmc SHARED src/mgxs_interface.cpp src/particle.cpp src/plot.cpp + src/position.cpp src/pugixml/pugixml_c.cpp src/random_lcg.cpp src/scattdata.cpp diff --git a/src/cell.cpp b/src/cell.cpp index 972540fbd2..acf1473bb8 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -279,19 +279,19 @@ Cell::Cell(pugi::xml_node cell_node) //============================================================================== bool -Cell::contains(Position r, Angle a, int32_t on_surface) const +Cell::contains(Position r, Direction u, int32_t on_surface) const { if (simple) { - return contains_simple(r, a, on_surface); + return contains_simple(r, u, on_surface); } else { - return contains_complex(r, a, on_surface); + return contains_complex(r, u, on_surface); } } //============================================================================== std::pair -Cell::distance(Position r, Angle a, int32_t on_surface) const +Cell::distance(Position r, Direction u, int32_t on_surface) const { double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; @@ -303,7 +303,7 @@ Cell::distance(Position r, Angle a, int32_t on_surface) const // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {token == on_surface}; - double d {surfaces_c[abs(token)-1]->distance(r, a, coincident)}; + double d {surfaces_c[abs(token)-1]->distance(r, u, coincident)}; // Check if this distance is the new minimum. if (d < min_dist) { @@ -354,7 +354,7 @@ Cell::to_hdf5(hid_t cell_group) const //============================================================================== bool -Cell::contains_simple(Position r, Angle a, int32_t on_surface) const +Cell::contains_simple(Position r, Direction u, int32_t on_surface) const { for (int32_t token : rpn) { if (token < OP_UNION) { @@ -367,7 +367,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const return false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(r, a); + bool sense = surfaces_c[abs(token)-1]->sense(r, u); if (sense != (token > 0)) {return false;} } } @@ -378,7 +378,7 @@ Cell::contains_simple(Position r, Angle a, int32_t on_surface) const //============================================================================== bool -Cell::contains_complex(Position r, Angle a, int32_t on_surface) const +Cell::contains_complex(Position r, Direction u, int32_t on_surface) const { // Make a stack of booleans. We don't know how big it needs to be, but we do // know that rpn.size() is an upper-bound. @@ -409,7 +409,7 @@ Cell::contains_complex(Position r, Angle a, int32_t on_surface) const stack[i_stack] = false; } else { // Note the off-by-one indexing - bool sense = surfaces_c[abs(token)-1]->sense(r, a);; + bool sense = surfaces_c[abs(token)-1]->sense(r, u);; stack[i_stack] = (sense == (token > 0)); } } @@ -492,16 +492,16 @@ extern "C" { bool cell_contains(Cell *c, double xyz[3], double uvw[3], int32_t on_surface) { Position r {xyz}; - Angle a {uvw}; - return c->contains(r, a, on_surface); + Direction u {uvw}; + return c->contains(r, u, on_surface); } void cell_distance(Cell *c, double xyz[3], double uvw[3], int32_t on_surface, double *min_dist, int32_t *i_surf) { Position r {xyz}; - Angle a {uvw}; - std::pair out = c->distance(r, a, on_surface); + Direction u {uvw}; + std::pair out = c->distance(r, u, on_surface); *min_dist = out.first; *i_surf = out.second; } diff --git a/src/cell.h b/src/cell.h index 08a6abf3f6..10265beb1b 100644 --- a/src/cell.h +++ b/src/cell.h @@ -8,7 +8,8 @@ #include "hdf5.h" #include "pugixml.hpp" -#include "geometry.h" + +#include "position.h" namespace openmc { @@ -98,19 +99,19 @@ public: //! known to be on. This index takes precedence over surface sense //! calculations. bool - contains(Position r, Angle a, int32_t on_surface) const; + contains(Position r, Direction u, int32_t on_surface) const; //! Find the oncoming boundary of this cell. std::pair - distance(Position r, Angle a, int32_t on_surface) const; + distance(Position r, Direction u, int32_t on_surface) const; //! \brief Write cell information to an HDF5 group. //! @param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; protected: - bool contains_simple(Position r, Angle a, int32_t on_surface) const; - bool contains_complex(Position r, Angle a, int32_t on_surface) const; + bool contains_simple(Position r, Direction u, int32_t on_surface) const; + bool contains_complex(Position r, Direction u, int32_t on_surface) const; }; } // namespace openmc diff --git a/src/geometry.h b/src/geometry.h index e069b9315b..09e7eda8eb 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -1,68 +1,10 @@ -#ifndef GEOMETRY_H -#define GEOMETRY_H +#ifndef OPENMC_GEOMETRY_H +#define OPENMC_GEOMETRY_H namespace openmc { extern "C" int openmc_root_universe; -struct Position { - double x = 0.; - double y = 0.; - double z = 0.; - - Position() = default; - Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; - Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; - - Position& operator+=(Position); - Position& operator+=(double); - Position& operator-=(Position); - Position& operator-=(double); - Position& operator*=(Position); - Position& operator*=(double); - const double& operator[](int i) const { - switch (i) { - case 0: return x; - case 1: return y; - case 2: return z; - } - } - double& operator[](int i) { - switch (i) { - case 0: return x; - case 1: return y; - case 2: return z; - } - } - - inline double dot(Position other) { - return x*other.x + y*other.y + z*other.z; - } -}; - -inline Position operator+(Position a, Position b) { return a += b; } -inline Position operator+(Position a, double b) { return a += b; } -inline Position operator+(double a, Position b) { return b += a; } - -inline Position operator-(Position a, Position b) { return a -= b; } -inline Position operator-(Position a, double b) { return a -= b; } -inline Position operator-(double a, Position b) { return b -= a; } - -inline Position operator*(Position a, Position b) { return a *= b; } -inline Position operator*(Position a, double b) { return a *= b; } -inline Position operator*(double a, Position b) { return b *= a; } - - -struct Angle : Position { - double& u() { return x; } - double& v() { return y; } - double& w() { return z; } - Angle() = default; - Angle(double u, double v, double w) : Position{u, v, w} { }; - Angle(const double uvw[]) : Position{uvw} { }; - Angle(Position r) : Position{r} { }; -}; - } // namespace openmc -#endif // GEOMETRY_H +#endif // OPENMC_GEOMETRY_H diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index 8b1f7a434d..b0466e613d 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -11,7 +11,7 @@ #include #include -#include "geometry.h" +#include "position.h" namespace openmc { diff --git a/src/lattice.cpp b/src/lattice.cpp index c67fc419f0..73dcb98667 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -223,25 +223,23 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const +RectLattice::distance(Position r, Direction u, const int i_xyz[3]) const { // Get short aliases to the coordinates. double x = r.x; double y = r.y; double z = r.z; - double u = a.u(); - double v = a.v(); // Determine the oncoming edge. - double x0 {copysign(0.5 * pitch[0], u)}; - double y0 {copysign(0.5 * pitch[1], v)}; + double x0 {copysign(0.5 * pitch[0], u.x)}; + double y0 {copysign(0.5 * pitch[1], u.y)}; // Left and right sides double d {INFTY}; std::array lattice_trans; - if ((std::abs(x - x0) > FP_PRECISION) && u != 0) { - d = (x0 - x) / u; - if (u > 0) { + if ((std::abs(x - x0) > FP_PRECISION) && u.x != 0) { + d = (x0 - x) / u.x; + if (u.x > 0) { lattice_trans = {1, 0, 0}; } else { lattice_trans = {-1, 0, 0}; @@ -249,11 +247,11 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const } // Front and back sides - if ((std::abs(y - y0) > FP_PRECISION) && v != 0) { - double this_d = (y0 - y) / v; + if ((std::abs(y - y0) > FP_PRECISION) && u.y != 0) { + double this_d = (y0 - y) / u.y; if (this_d < d) { d = this_d; - if (v > 0) { + if (u.y > 0) { lattice_trans = {0, 1, 0}; } else { lattice_trans = {0, -1, 0}; @@ -263,13 +261,12 @@ RectLattice::distance(Position r, Angle a, const int i_xyz[3]) const // Top and bottom sides if (is_3d) { - double w {a.w()}; - double z0 {copysign(0.5 * pitch[2], w)}; - if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { - double this_d = (z0 - z) / w; + double z0 {copysign(0.5 * pitch[2], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; if (this_d < d) { d = this_d; - if (w > 0) { + if (u.z > 0) { lattice_trans = {0, 0, 1}; } else { lattice_trans = {0, 0, -1}; @@ -579,11 +576,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const +HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const { // Compute the direction on the hexagonal basis. - double beta_dir = a.u() * std::sqrt(3.0) / 2.0 + a.v() / 2.0; - double gamma_dir = a.u() * std::sqrt(3.0) / 2.0 - a.v() / 2.0; + double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; + double gamma_dir = u.x * std::sqrt(3.0) / 2.0 - u.y / 2.0; // Note that hexagonal lattice distance calculations are performed // using the particle's coordinates relative to the neighbor lattice @@ -636,18 +633,18 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const } // Upper and lower sides. - edge = -copysign(0.5*pitch[0], a.v()); - if (a.v() > 0) { + edge = -copysign(0.5*pitch[0], u.y); + if (u.y > 0) { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } - if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) { - double this_d = (edge - r_t.y) / a.v(); + if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) { + double this_d = (edge - r_t.y) / u.y; if (this_d < d) { - if (a.v() > 0) { + if (u.y > 0) { lattice_trans = {0, 1, 0}; } else { lattice_trans = {0, -1, 0}; @@ -659,13 +656,12 @@ HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const // Top and bottom sides if (is_3d) { double z = r.z; - double w = a.w(); - double z0 {copysign(0.5 * pitch[1], w)}; - if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { - double this_d = (z0 - z) / w; + double z0 {copysign(0.5 * pitch[1], u.z)}; + if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { + double this_d = (z0 - z) / u.z; if (this_d < d) { d = this_d; - if (w > 0) { + if (u.z > 0) { lattice_trans = {0, 0, 1}; } else { lattice_trans = {0, 0, -1}; @@ -897,8 +893,8 @@ extern "C" { const int i_xyz[3], double *d, int lattice_trans[3]) { Position r {xyz}; - Angle a {uvw}; - std::pair> ld {lat->distance(r, a, i_xyz)}; + Direction u {uvw}; + std::pair> ld {lat->distance(r, u, i_xyz)}; *d = ld.first; lattice_trans[0] = ld.second[0]; lattice_trans[1] = ld.second[1]; diff --git a/src/lattice.h b/src/lattice.h index a6c4496621..5bd15a6126 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -7,11 +7,12 @@ #include #include -#include "constants.h" -#include "geometry.h" #include "hdf5.h" #include "pugixml.hpp" +#include "constants.h" +#include "position.h" + namespace openmc { @@ -82,7 +83,7 @@ public: //! @return The distance to the next crossing and an array indicating how the //! lattice indices would change after crossing that boundary. virtual std::pair> - distance(Position r, Angle a, const int i_xyz[3]) const + distance(Position r, Direction u, const int i_xyz[3]) const = 0; //! \brief Find the lattice tile indices for a given point. @@ -194,7 +195,7 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(Position r, Angle a, const int i_xyz[3]) const; + distance(Position r, Direction u, const int i_xyz[3]) const; std::array get_indices(Position r) const; @@ -234,7 +235,7 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(Position r, Angle a, const int i_xyz[3]) const; + distance(Position r, Direction u, const int i_xyz[3]) const; std::array get_indices(Position r) const; diff --git a/src/geometry.cpp b/src/position.cpp similarity index 75% rename from src/geometry.cpp rename to src/position.cpp index 53bf32b5c1..85e511041e 100644 --- a/src/geometry.cpp +++ b/src/position.cpp @@ -1,7 +1,11 @@ -#include "geometry.h" +#include "position.h" namespace openmc { +//============================================================================== +// Position implementation +//============================================================================== + Position& Position::operator+=(Position other) { diff --git a/src/position.h b/src/position.h new file mode 100644 index 0000000000..808e0bcf67 --- /dev/null +++ b/src/position.h @@ -0,0 +1,74 @@ +#ifndef OPENMC_POSITION_H +#define OPENMC_POSITION_H + +namespace openmc { + +//============================================================================== +//! Type representing a position in Cartesian coordinates +//============================================================================== + +struct Position { + // Constructors + Position() = default; + Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; + Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + + // Unary operators + Position& operator+=(Position); + Position& operator+=(double); + Position& operator-=(Position); + Position& operator-=(double); + Position& operator*=(Position); + Position& operator*=(double); + const double& operator[](int i) const { + switch (i) { + case 0: return x; + case 1: return y; + case 2: return z; + } + } + double& operator[](int i) { + switch (i) { + case 0: return x; + case 1: return y; + case 2: return z; + } + } + + // Other member functions + + //! Dot product of two vectors + //! \param[in] other Vector to take dot product with + //! \result Resulting dot product + inline double dot(Position other) { + return x*other.x + y*other.y + z*other.z; + } + + // Data members + double x = 0.; + double y = 0.; + double z = 0.; +}; + +// Binary operators +inline Position operator+(Position a, Position b) { return a += b; } +inline Position operator+(Position a, double b) { return a += b; } +inline Position operator+(double a, Position b) { return b += a; } + +inline Position operator-(Position a, Position b) { return a -= b; } +inline Position operator-(Position a, double b) { return a -= b; } +inline Position operator-(double a, Position b) { return b -= a; } + +inline Position operator*(Position a, Position b) { return a *= b; } +inline Position operator*(Position a, double b) { return a *= b; } +inline Position operator*(double a, Position b) { return b *= a; } + +//============================================================================== +//! Type representing a vector direction in Cartesian coordinates +//============================================================================== + +using Direction = Position; + +} // namespace openmc + +#endif // OPENMC_POSITION_H \ No newline at end of file diff --git a/src/surface.cpp b/src/surface.cpp index aec0972b8f..8258861526 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -178,7 +178,7 @@ Surface::Surface(pugi::xml_node surf_node) } bool -Surface::sense(Position r, Angle a) const +Surface::sense(Position r, Direction u) const { // Evaluate the surface equation at the particle's coordinates to determine // which side the particle is on. @@ -189,22 +189,22 @@ Surface::sense(Position r, Angle a) const // Particle may be coincident with this surface. To determine the sense, we // look at the direction of the particle relative to the surface normal (by // default in the positive direction) via their dot product. - return a.dot(normal(r)) > 0.0; + return u.dot(normal(r)) > 0.0; } return f > 0.0; } -Angle -Surface::reflect(Position r, Angle a) const +Direction +Surface::reflect(Position r, Direction u) const { // Determine projection of direction onto normal and squared magnitude of // normal. - Angle n = normal(r); - const double projection = n.dot(a); + Direction n = normal(r); + const double projection = n.dot(u); const double magnitude = n.dot(n); // Reflect direction according to normal. - return a -= (2.0 * projection / magnitude) * n; + return u -= (2.0 * projection / magnitude) * n; } void @@ -257,11 +257,11 @@ PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) // The template parameter indicates the axis normal to the plane. template double -axis_aligned_plane_distance(Position r, Angle a, bool coincident, double offset) +axis_aligned_plane_distance(Position r, Direction u, bool coincident, double offset) { const double f = offset - r[i]; - if (coincident or std::abs(f) < FP_COINCIDENT or a[i] == 0.0) return INFTY; - const double d = f / a[i]; + if (coincident or std::abs(f) < FP_COINCIDENT or u[i] == 0.0) return INFTY; + const double d = f / u[i]; if (d < 0.0) return INFTY; return d; } @@ -281,12 +281,12 @@ double SurfaceXPlane::evaluate(Position r) const return r.x - x0; } -double SurfaceXPlane::distance(Position r, Angle a, bool coincident) const +double SurfaceXPlane::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_plane_distance<0>(r, a, coincident, x0); + return axis_aligned_plane_distance<0>(r, u, coincident, x0); } -Angle SurfaceXPlane::normal(Position r) const +Direction SurfaceXPlane::normal(Position r) const { return {1., 0., 0.}; } @@ -299,23 +299,23 @@ void SurfaceXPlane::to_hdf5_inner(hid_t group_id) const } bool SurfaceXPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Angle& a) const + Direction& u) const { - Angle other_n = other->normal(r); - if (other_n.u() == 1 and other_n.v() == 0 and other_n.w() == 0) { + Direction other_n = other->normal(r); + if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) { r.x = x0; return false; } else { // Assume the partner is an YPlane (the only supported partner). Use the - // evaluate function to find y0, then adjust position/angle for rotational + // evaluate function to find y0, then adjust position/Direction for rotational // symmetry. double y0 = -other->evaluate({0., 0., 0.}); r.y = r.x - x0 + y0; r.x = x0; - double u = a.u(); - a.u() = -a.v(); - a.v() = u; + double ux = u.x; + u.x = -u.y; + u.y = ux; return true; } @@ -342,12 +342,12 @@ double SurfaceYPlane::evaluate(Position r) const return r.y - y0; } -double SurfaceYPlane::distance(Position r, Angle a, bool coincident) const +double SurfaceYPlane::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_plane_distance<1>(r, a, coincident, y0); + return axis_aligned_plane_distance<1>(r, u, coincident, y0); } -Angle SurfaceYPlane::normal(Position r) const +Direction SurfaceYPlane::normal(Position r) const { return {0., 1., 0.}; } @@ -360,24 +360,24 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const } bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Angle& a) const + Direction& u) const { - Angle other_n = other->normal(r); - if (other_n.u() == 0 and other_n.v() == 1 and other_n.w() == 0) { + Direction other_n = other->normal(r); + if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) { // The periodic partner is also aligned along y. Just change the y coord. r.y = y0; return false; } else { // Assume the partner is an XPlane (the only supported partner). Use the - // evaluate function to find x0, then adjust position/angle for rotational + // evaluate function to find x0, then adjust position/Direction for rotational // symmetry. double x0 = -other->evaluate({0., 0., 0.}); r.x = r.y - y0 + x0; r.y = y0; - double u = a.u(); - a.u() = a.v(); - a.v() = -u; + double ux = u.x; + u.x = u.y; + u.y = -ux; return true; } @@ -404,12 +404,12 @@ double SurfaceZPlane::evaluate(Position r) const return r.z - z0; } -double SurfaceZPlane::distance(Position r, Angle a, bool coincident) const +double SurfaceZPlane::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_plane_distance<2>(r, a, coincident, z0); + return axis_aligned_plane_distance<2>(r, u, coincident, z0); } -Angle SurfaceZPlane::normal(Position r) const +Direction SurfaceZPlane::normal(Position r) const { return {0., 0., 1.}; } @@ -422,7 +422,7 @@ void SurfaceZPlane::to_hdf5_inner(hid_t group_id) const } bool SurfaceZPlane::periodic_translate(const PeriodicSurface *other, Position& r, - Angle& a) const + Direction& u) const { // Assume the other plane is aligned along z. Just change the z coord. r.z = z0; @@ -452,10 +452,10 @@ SurfacePlane::evaluate(Position r) const } double -SurfacePlane::distance(Position r, Angle a, bool coincident) const +SurfacePlane::distance(Position r, Direction u, bool coincident) const { const double f = A*r.x + B*r.y + C*r.z - D; - const double projection = A*a.u() + B*a.v() + C*a.w(); + const double projection = A*u.x + B*u.y + C*u.z; if (coincident or std::abs(f) < FP_COINCIDENT or projection == 0.0) { return INFTY; } else { @@ -465,7 +465,7 @@ SurfacePlane::distance(Position r, Angle a, bool coincident) const } } -Angle +Direction SurfacePlane::normal(Position r) const { return {A, B, C}; @@ -479,7 +479,7 @@ void SurfacePlane::to_hdf5_inner(hid_t group_id) const } bool SurfacePlane::periodic_translate(const PeriodicSurface *other, Position& r, - Angle& a) const + Direction& u) const { // This function assumes the other plane shares this plane's normal direction. @@ -520,7 +520,7 @@ axis_aligned_cylinder_evaluate(Position r, double offset1, // The other two parameters indicate the other two axes. offset1 and offset2 // should correspond with i2 and i3, respectively. template double -axis_aligned_cylinder_distance(Position r, Angle u, +axis_aligned_cylinder_distance(Position r, Direction u, bool coincident, double offset1, double offset2, double radius) { const double a = 1.0 - u[i1]*u[i1]; // u^2 + v^2 @@ -565,14 +565,14 @@ axis_aligned_cylinder_distance(Position r, Angle u, // The first template parameter indicates which axis the cylinder is aligned to. // The other two parameters indicate the other two axes. offset1 and offset2 // should correspond with i2 and i3, respectively. -template Angle +template Direction axis_aligned_cylinder_normal(Position r, double offset1, double offset2) { - Angle a; - a[i2] = 2.0 * (r[i2] - offset1); - a[i3] = 2.0 * (r[i3] - offset2); - a[i1] = 0.0; - return a; + Direction u; + u[i2] = 2.0 * (r[i2] - offset1); + u[i3] = 2.0 * (r[i3] - offset2); + u[i1] = 0.0; + return u; } //============================================================================== @@ -590,13 +590,13 @@ double SurfaceXCylinder::evaluate(Position r) const return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r); } -double SurfaceXCylinder::distance(Position r, Angle a, bool coincident) const +double SurfaceXCylinder::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cylinder_distance<0, 1, 2>(r, a, coincident, y0, z0, + return axis_aligned_cylinder_distance<0, 1, 2>(r, u, coincident, y0, z0, this->r); } -Angle SurfaceXCylinder::normal(Position r) const +Direction SurfaceXCylinder::normal(Position r) const { return axis_aligned_cylinder_normal<0, 1, 2>(r, y0, z0); } @@ -624,13 +624,13 @@ double SurfaceYCylinder::evaluate(Position r) const return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r); } -double SurfaceYCylinder::distance(Position r, Angle a, bool coincident) const +double SurfaceYCylinder::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cylinder_distance<1, 0, 2>(r, a, coincident, x0, z0, + return axis_aligned_cylinder_distance<1, 0, 2>(r, u, coincident, x0, z0, this->r); } -Angle SurfaceYCylinder::normal(Position r) const +Direction SurfaceYCylinder::normal(Position r) const { return axis_aligned_cylinder_normal<1, 0, 2>(r, x0, z0); } @@ -657,13 +657,13 @@ double SurfaceZCylinder::evaluate(Position r) const return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r); } -double SurfaceZCylinder::distance(Position r, Angle a, bool coincident) const +double SurfaceZCylinder::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cylinder_distance<2, 0, 1>(r, a, coincident, x0, y0, + return axis_aligned_cylinder_distance<2, 0, 1>(r, u, coincident, x0, y0, this->r); } -Angle SurfaceZCylinder::normal(Position r) const +Direction SurfaceZCylinder::normal(Position r) const { return axis_aligned_cylinder_normal<2, 0, 1>(r, x0, y0); } @@ -693,12 +693,12 @@ double SurfaceSphere::evaluate(Position r) const return x*x + y*y + z*z - this->r*this->r; } -double SurfaceSphere::distance(Position r, Angle a, bool coincident) const +double SurfaceSphere::distance(Position r, Direction u, bool coincident) const { const double x = r.x - x0; const double y = r.y - y0; const double z = r.z - z0; - const double k = x*a.u() + y*a.v() + z*a.w(); + const double k = x*u.x + y*u.y + z*u.z; const double c = x*x + y*y + z*z - this->r*this->r; const double quad = k*k - c; @@ -731,7 +731,7 @@ double SurfaceSphere::distance(Position r, Angle a, bool coincident) const } } -Angle SurfaceSphere::normal(Position r) const +Direction SurfaceSphere::normal(Position r) const { return {2.0*(r.x - x0), 2.0*(r.y - y0), 2.0*(r.z - z0)}; } @@ -764,7 +764,7 @@ axis_aligned_cone_evaluate(Position r, double offset1, // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. template double -axis_aligned_cone_distance(Position r, Angle u, +axis_aligned_cone_distance(Position r, Direction u, bool coincident, double offset1, double offset2, double offset3, double radius_sq) { @@ -817,15 +817,15 @@ axis_aligned_cone_distance(Position r, Angle u, // The first template parameter indicates which axis the cone is aligned to. // The other two parameters indicate the other two axes. offset1, offset2, // and offset3 should correspond with i1, i2, and i3, respectively. -template Angle +template Direction axis_aligned_cone_normal(Position r, double offset1, double offset2, double offset3, double radius_sq) { - Angle a; - a[i1] = -2.0 * radius_sq * (r[i1] - offset1); - a[i2] = 2.0 * (r[i2] - offset2); - a[i3] = 2.0 * (r[i3] - offset3); - return a; + Direction u; + u[i1] = -2.0 * radius_sq * (r[i1] - offset1); + u[i2] = 2.0 * (r[i2] - offset2); + u[i3] = 2.0 * (r[i3] - offset3); + return u; } //============================================================================== @@ -843,13 +843,13 @@ double SurfaceXCone::evaluate(Position r) const return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq); } -double SurfaceXCone::distance(Position r, Angle a, bool coincident) const +double SurfaceXCone::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cone_distance<0, 1, 2>(r, a, coincident, x0, y0, z0, + return axis_aligned_cone_distance<0, 1, 2>(r, u, coincident, x0, y0, z0, r_sq); } -Angle SurfaceXCone::normal(Position r) const +Direction SurfaceXCone::normal(Position r) const { return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq); } @@ -876,13 +876,13 @@ double SurfaceYCone::evaluate(Position r) const return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq); } -double SurfaceYCone::distance(Position r, Angle a, bool coincident) const +double SurfaceYCone::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cone_distance<1, 0, 2>(r, a, coincident, y0, x0, z0, + return axis_aligned_cone_distance<1, 0, 2>(r, u, coincident, y0, x0, z0, r_sq); } -Angle SurfaceYCone::normal(Position r) const +Direction SurfaceYCone::normal(Position r) const { return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq); } @@ -909,13 +909,13 @@ double SurfaceZCone::evaluate(Position r) const return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq); } -double SurfaceZCone::distance(Position r, Angle a, bool coincident) const +double SurfaceZCone::distance(Position r, Direction u, bool coincident) const { - return axis_aligned_cone_distance<2, 0, 1>(r, a, coincident, z0, x0, y0, + return axis_aligned_cone_distance<2, 0, 1>(r, u, coincident, z0, x0, y0, r_sq); } -Angle SurfaceZCone::normal(Position r) const +Direction SurfaceZCone::normal(Position r) const { return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq); } @@ -949,14 +949,14 @@ SurfaceQuadric::evaluate(Position r) const } double -SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const +SurfaceQuadric::distance(Position r, Direction ang, bool coincident) const { const double &x = r.x; const double &y = r.y; const double &z = r.z; - const double &u = ang.u(); - const double &v = ang.v(); - const double &w = ang.w(); + const double &u = ang.x; + const double &v = ang.y; + const double &w = ang.z; const double a = A*u*u + B*v*v + C*w*w + D*u*v + E*v*w + F*u*w; const double k = (A*u*x + B*v*y + C*w*z + 0.5*(D*(u*y + v*x) + @@ -1002,7 +1002,7 @@ SurfaceQuadric::distance(Position r, Angle ang, bool coincident) const return d; } -Angle +Direction SurfaceQuadric::normal(Position r) const { const double &x = r.x; @@ -1222,21 +1222,21 @@ extern "C" { void surface_reflect(Surface *surf, double xyz[3], double uvw[3]) { Position r {xyz}; - Angle a {uvw}; - a = surf->reflect(r, a); + Direction u {uvw}; + u = surf->reflect(r, u); - uvw[0] = a.u(); - uvw[1] = a.v(); - uvw[2] = a.w(); + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; } void surface_normal(Surface *surf, double xyz[3], double uvw[3]) { Position r {xyz}; - Angle a = surf->normal(r); - uvw[0] = a.u(); - uvw[1] = a.v(); - uvw[2] = a.w(); + Direction u = surf->normal(r); + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; } void surface_to_hdf5(Surface *surf, hid_t group) {surf->to_hdf5(group);} @@ -1248,16 +1248,16 @@ extern "C" { double uvw[3]) { Position r {xyz}; - Angle a {uvw}; - bool rotational = surf->periodic_translate(other, r, a); + Direction u {uvw}; + bool rotational = surf->periodic_translate(other, r, u); // Copy back to arrays xyz[0] = r.x; xyz[1] = r.y; xyz[2] = r.z; - uvw[0] = a.u(); - uvw[1] = a.v(); - uvw[2] = a.w(); + uvw[0] = u.x; + uvw[1] = u.y; + uvw[2] = u.z; return rotational; } diff --git a/src/surface.h b/src/surface.h index be3ad52ab9..7d639c57a8 100644 --- a/src/surface.h +++ b/src/surface.h @@ -9,7 +9,7 @@ #include "pugixml.hpp" #include "constants.h" -#include "geometry.h" +#include "position.h" namespace openmc { @@ -71,14 +71,14 @@ public: //! point is very close to the surface. //! @return true if the point is on the "positive" side of the surface and //! false otherwise. - bool sense(Position r, Angle a) const; + bool sense(Position r, Direction u) const; //! Determine the direction of a ray reflected from the surface. //! @param r The point at which the ray is incident. //! @param o A direction. This is both an input and an output parameter. //! It specifies the icident direction on input and the reflected direction //! on output. - Angle reflect(Position r, Angle a) const; + Direction reflect(Position r, Direction u) const; //! Evaluate the equation describing the surface. //! @@ -92,12 +92,12 @@ public: //! @param o The direction of the ray. //! @param coincident A hint to the code that the given point should lie //! exactly on the surface. - virtual double distance(Position r, Angle a, bool coincident) const = 0; + virtual double distance(Position r, Direction u, bool coincident) const = 0; //! Compute the local outward normal direction of the surface. //! @param r A 3D Cartesian coordinate. //! @return Normal direction - virtual Angle normal(Position r) const = 0; + virtual Direction normal(Position r) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. //! @param group_id An HDF5 group id. @@ -132,7 +132,7 @@ public: //! @return true if this surface and its partner make a rotationally-periodic //! boundary condition. virtual bool periodic_translate(const PeriodicSurface *other, Position& r, - Angle& a) const = 0; + Direction& u) const = 0; //! Get the bounding box for this surface. virtual BoundingBox bounding_box() const = 0; @@ -150,10 +150,10 @@ class SurfaceXPlane : public PeriodicSurface public: explicit SurfaceXPlane(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) + bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) const; BoundingBox bounding_box() const; }; @@ -170,10 +170,10 @@ class SurfaceYPlane : public PeriodicSurface public: explicit SurfaceYPlane(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) + bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) const; BoundingBox bounding_box() const; }; @@ -190,10 +190,10 @@ class SurfaceZPlane : public PeriodicSurface public: explicit SurfaceZPlane(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) + bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) const; BoundingBox bounding_box() const; }; @@ -210,10 +210,10 @@ class SurfacePlane : public PeriodicSurface public: explicit SurfacePlane(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; - bool periodic_translate(const PeriodicSurface *other, Position& r, Angle& a) + bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) const; BoundingBox bounding_box() const; }; @@ -231,8 +231,8 @@ class SurfaceXCylinder : public Surface public: explicit SurfaceXCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -249,8 +249,8 @@ class SurfaceYCylinder : public Surface public: explicit SurfaceYCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -267,8 +267,8 @@ class SurfaceZCylinder : public Surface public: explicit SurfaceZCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -285,8 +285,8 @@ class SurfaceSphere : public Surface public: explicit SurfaceSphere(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -303,8 +303,8 @@ class SurfaceXCone : public Surface public: explicit SurfaceXCone(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -321,8 +321,8 @@ class SurfaceYCone : public Surface public: explicit SurfaceYCone(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -339,8 +339,8 @@ class SurfaceZCone : public Surface public: explicit SurfaceZCone(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -357,8 +357,8 @@ class SurfaceQuadric : public Surface public: explicit SurfaceQuadric(pugi::xml_node surf_node); double evaluate(Position r) const; - double distance(Position r, Angle a, bool coincident) const; - Angle normal(Position r) const; + double distance(Position r, Direction u, bool coincident) const; + Direction normal(Position r) const; void to_hdf5_inner(hid_t group_id) const; }; @@ -367,19 +367,19 @@ public: //============================================================================== extern "C" { - Surface* surface_pointer(int surf_ind); - int surface_id(Surface *surf); - int surface_bc(Surface *surf); - bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); - void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); - double surface_distance(Surface *surf, double xyz[3], double uvw[3], - bool coincident); - void surface_normal(Surface *surf, double xyz[3], double uvw[3]); - void surface_to_hdf5(Surface *surf, hid_t group); - int surface_i_periodic(PeriodicSurface *surf); - bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, - double xyz[3], double uvw[3]); - void free_memory_surfaces_c(); + Surface* surface_pointer(int surf_ind); + int surface_id(Surface *surf); + int surface_bc(Surface *surf); + bool surface_sense(Surface *surf, double xyz[3], double uvw[3]); + void surface_reflect(Surface *surf, double xyz[3], double uvw[3]); + double surface_distance(Surface *surf, double xyz[3], double uvw[3], + bool coincident); + void surface_normal(Surface *surf, double xyz[3], double uvw[3]); + void surface_to_hdf5(Surface *surf, hid_t group); + int surface_i_periodic(PeriodicSurface *surf); + bool surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, + double xyz[3], double uvw[3]); + void free_memory_surfaces_c(); } } // namespace openmc From c72ff04afb79219f172c082409e936511eb7d4d0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Aug 2018 13:46:53 -0500 Subject: [PATCH 27/36] Make write_dataset for Position an overload (not specialization) --- src/hdf5_interface.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index b0466e613d..6faba4078f 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -183,8 +183,8 @@ write_dataset(hid_t obj_id, const char* name, const std::array& buffer) write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } -template<> inline void -write_dataset(hid_t obj_id, const char* name, Position r) +inline void +write_dataset(hid_t obj_id, const char* name, Position r) { std::array buffer {r.x, r.y, r.z}; write_dataset(obj_id, name, buffer); From 134ad969d832525beff391b4f54214b6d3c0d0f1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Aug 2018 13:54:33 -0500 Subject: [PATCH 28/36] Fix some comments, include guards --- src/cell.h | 14 +++++++------- src/hdf5_interface.h | 6 +++--- src/lattice.h | 44 ++++++++++++++++++++++---------------------- src/surface.h | 38 +++++++++++++++++++------------------- 4 files changed, 51 insertions(+), 51 deletions(-) diff --git a/src/cell.h b/src/cell.h index 10265beb1b..062251b48c 100644 --- a/src/cell.h +++ b/src/cell.h @@ -1,5 +1,5 @@ -#ifndef CELL_H -#define CELL_H +#ifndef OPENMC_CELL_H +#define OPENMC_CELL_H #include #include @@ -92,10 +92,10 @@ public: //! provides a performance benefit for the common case. In //! contains_complex, we evaluate the RPN expression using a stack, similar to //! how a RPN calculator would work. - //! @param xyz[3] The 3D Cartesian coordinate to check. - //! @param uvw[3] A direction used to "break ties" the coordinates are very + //! \param r The 3D Cartesian coordinate to check. + //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. - //! @param on_surface The signed index of a surface that the coordinate is + //! \param on_surface The signed index of a surface that the coordinate is //! known to be on. This index takes precedence over surface sense //! calculations. bool @@ -106,7 +106,7 @@ public: distance(Position r, Direction u, int32_t on_surface) const; //! \brief Write cell information to an HDF5 group. - //! @param group_id An HDF5 group id. + //! \param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; protected: @@ -115,4 +115,4 @@ protected: }; } // namespace openmc -#endif // CELL_H +#endif // OPENMC_CELL_H diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index 6faba4078f..152c52c1cd 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -1,5 +1,5 @@ -#ifndef HDF5_INTERFACE_H -#define HDF5_INTERFACE_H +#ifndef OPENMC_HDF5_INTERFACE_H +#define OPENMC_HDF5_INTERFACE_H #include "hdf5.h" #include "hdf5_hl.h" @@ -191,4 +191,4 @@ write_dataset(hid_t obj_id, const char* name, Position r) } } // namespace openmc -#endif //HDF5_INTERFACE_H +#endif // OPENMC_HDF5_INTERFACE_H diff --git a/src/lattice.h b/src/lattice.h index 5bd15a6126..94d7b43878 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -1,5 +1,5 @@ -#ifndef LATTICE_H -#define LATTICE_H +#ifndef OPENMC_LATTICE_H +#define OPENMC_LATTICE_H #include #include @@ -71,54 +71,54 @@ public: int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map); //! \brief Check lattice indices. - //! @param i_xyz[3] The indices for a lattice tile. - //! @return true if the given indices fit within the lattice bounds. False + //! \param i_xyz[3] The indices for a lattice tile. + //! \return true if the given indices fit within the lattice bounds. False //! otherwise. virtual bool are_valid_indices(const int i_xyz[3]) const = 0; //! \brief Find the next lattice surface crossing - //! @param r A 3D Cartesian coordinate. - //! @param a A 3D Cartesian direction. - //! @param i_xyz[3] The indices for a lattice tile. - //! @return The distance to the next crossing and an array indicating how the + //! \param r A 3D Cartesian coordinate. + //! \param u A 3D Cartesian direction. + //! \param i_xyz[3] The indices for a lattice tile. + //! \return The distance to the next crossing and an array indicating how the //! lattice indices would change after crossing that boundary. virtual std::pair> distance(Position r, Direction u, const int i_xyz[3]) const = 0; //! \brief Find the lattice tile indices for a given point. - //! @param r A 3D Cartesian coordinate. - //! @return An array containing the indices of a lattice tile. + //! \param r A 3D Cartesian coordinate. + //! \return An array containing the indices of a lattice tile. virtual std::array get_indices(Position r) const = 0; //! \brief Get coordinates local to a lattice tile. - //! @param r A 3D Cartesian coordinate. - //! @param i_xyz[3] The indices for a lattice tile. - //! @return Local 3D Cartesian coordinates. + //! \param r A 3D Cartesian coordinate. + //! \param i_xyz[3] The indices for a lattice tile. + //! \return Local 3D Cartesian coordinates. virtual Position get_local_position(Position r, const int i_xyz[3]) const = 0; //! \brief Check flattened lattice index. - //! @param indx The index for a lattice tile. - //! @return true if the given index fit within the lattice bounds. False + //! \param indx The index for a lattice tile. + //! \return true if the given index fit within the lattice bounds. False //! otherwise. virtual bool is_valid_index(int indx) const {return (indx >= 0) && (indx < universes.size());} //! \brief Get the distribcell offset for a lattice tile. - //! @param The map index for the target cell. - //! @param i_xyz[3] The indices for a lattice tile. - //! @return Distribcell offset i.e. the largest instance number for the target + //! \param The map index for the target cell. + //! \param i_xyz[3] The indices for a lattice tile. + //! \return Distribcell offset i.e. the largest instance number for the target //! cell found in the geometry tree under this lattice tile. virtual int32_t& offset(int map, const int i_xyz[3]) = 0; //! \brief Convert an array index to a useful human-readable string. - //! @param indx The index for a lattice tile. - //! @return A string representing the lattice tile. + //! \param indx The index for a lattice tile. + //! \return A string representing the lattice tile. virtual std::string index_to_string(int indx) const = 0; //! \brief Write lattice information to an HDF5 group. - //! @param group_id An HDF5 group id. + //! \param group_id An HDF5 group id. void to_hdf5(hid_t group_id) const; protected: @@ -258,4 +258,4 @@ private: }; } // namespace openmc -#endif // LATTICE_H +#endif // OPENMC_LATTICE_H diff --git a/src/surface.h b/src/surface.h index 7d639c57a8..b267af05ca 100644 --- a/src/surface.h +++ b/src/surface.h @@ -1,5 +1,5 @@ -#ifndef SURFACE_H -#define SURFACE_H +#ifndef OPENMC_SURFACE_H +#define OPENMC_SURFACE_H #include #include // For numeric_limits @@ -66,16 +66,16 @@ public: virtual ~Surface() {} //! Determine which side of a surface a point lies on. - //! @param r The 3D Cartesian coordinate of a point. - //! @param o A direction used to "break ties" and pick a sense when the + //! \param r The 3D Cartesian coordinate of a point. + //! \param u A direction used to "break ties" and pick a sense when the //! point is very close to the surface. - //! @return true if the point is on the "positive" side of the surface and + //! \return true if the point is on the "positive" side of the surface and //! false otherwise. bool sense(Position r, Direction u) const; //! Determine the direction of a ray reflected from the surface. - //! @param r The point at which the ray is incident. - //! @param o A direction. This is both an input and an output parameter. + //! \param r The point at which the ray is incident. + //! \param u A direction. This is both an input and an output parameter. //! It specifies the icident direction on input and the reflected direction //! on output. Direction reflect(Position r, Direction u) const; @@ -84,23 +84,23 @@ public: //! //! Surfaces can be described by some function f(x, y, z) = 0. This member //! function evaluates that mathematical function. - //! @param r A 3D Cartesian coordinate. + //! \param r A 3D Cartesian coordinate. virtual double evaluate(Position r) const = 0; //! Compute the distance between a point and the surface along a ray. - //! @param r A 3D Cartesian coordinate. - //! @param o The direction of the ray. - //! @param coincident A hint to the code that the given point should lie + //! \param r A 3D Cartesian coordinate. + //! \param u The direction of the ray. + //! \param coincident A hint to the code that the given point should lie //! exactly on the surface. virtual double distance(Position r, Direction u, bool coincident) const = 0; //! Compute the local outward normal direction of the surface. - //! @param r A 3D Cartesian coordinate. - //! @return Normal direction + //! \param r A 3D Cartesian coordinate. + //! \return Normal direction virtual Direction normal(Position r) const = 0; //! Write all information needed to reconstruct the surface to an HDF5 group. - //! @param group_id An HDF5 group id. + //! \param group_id An HDF5 group id. //TODO: this probably needs to include i_periodic for PeriodicSurface void to_hdf5(hid_t group_id) const; @@ -124,12 +124,12 @@ public: explicit PeriodicSurface(pugi::xml_node surf_node); //! Translate a particle onto this surface from a periodic partner surface. - //! @param other A pointer to the partner surface in this periodic BC. - //! @param r A point on the partner surface that will be translated onto + //! \param other A pointer to the partner surface in this periodic BC. + //! \param r A point on the partner surface that will be translated onto //! this surface. - //! @param a A direction that will be rotated for systems with rotational + //! \param u A direction that will be rotated for systems with rotational //! periodicity. - //! @return true if this surface and its partner make a rotationally-periodic + //! \return true if this surface and its partner make a rotationally-periodic //! boundary condition. virtual bool periodic_translate(const PeriodicSurface *other, Position& r, Direction& u) const = 0; @@ -383,4 +383,4 @@ extern "C" { } } // namespace openmc -#endif // SURFACE_H +#endif // OPENMC_SURFACE_H From 1bd1a1b855fccb26cecf7300032bcc4a5a90a092 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Wed, 8 Aug 2018 14:52:09 -0500 Subject: [PATCH 29/36] Remove some redundancy --- src/tallies/tally_filter_zernike.F90 | 39 ++++++++++++++++------------ 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 1ece62dbb2..96e016d2ee 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -13,6 +13,10 @@ module tally_filter_zernike implicit none private + public :: openmc_zernike_filter_get_order + public :: openmc_zernike_filter_get_params + public :: openmc_zernike_filter_set_order + public :: openmc_zernike_filter_set_params !=============================================================================== ! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position @@ -24,6 +28,7 @@ module tally_filter_zernike real(8) :: y real(8) :: r contains + procedure :: calc_n_bins => calc_n_bins_zn procedure :: from_xml => from_xml_zn procedure :: get_all_bins => get_all_bins_zn procedure :: to_statepoint => to_statepoint_zn @@ -37,7 +42,8 @@ module tally_filter_zernike type, public, extends(ZernikeFilter) :: ZernikeRadialFilter contains - procedure :: from_xml => from_xml_zn_rad + procedure :: calc_n_bins => calc_n_bins_zn_rad + ! Inherit from_xml from ZernikeFilter procedure :: get_all_bins => get_all_bins_zn_rad procedure :: to_statepoint => to_statepoint_zn_rad procedure :: text_label => text_label_zn_rad @@ -49,6 +55,15 @@ contains ! ZernikeFilter methods !=============================================================================== + function calc_n_bins_zn(this) result(num_n_bins) + class(ZernikeFilter), intent(in) :: this + integer :: n + integer :: num_n_bins + + n = this % order + num_n_bins = ((n+1) * (n+2))/2 + end function calc_n_bins_zn + subroutine from_xml_zn(this, node) class(ZernikeFilter), intent(inout) :: this type(XMLNode), intent(in) :: node @@ -63,7 +78,7 @@ contains ! Get specified order call get_node_value(node, "order", n) this % order = n - this % n_bins = ((n + 1)*(n + 2))/2 + this % n_bins = this % calc_n_bins() end subroutine from_xml_zn subroutine get_all_bins_zn(this, p, estimator, match) @@ -129,22 +144,14 @@ contains ! ZernikeRadialFilter methods !=============================================================================== - subroutine from_xml_zn_rad(this, node) - class(ZernikeRadialFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - + function calc_n_bins_zn_rad(this) result(num_n_bins) + class(ZernikeRadialFilter), intent(in) :: this integer :: n + integer :: num_n_bins - ! Get center of cylinder and radius - call get_node_value(node, "x", this % x) - call get_node_value(node, "y", this % y) - call get_node_value(node, "r", this % r) - - ! Get specified order - call get_node_value(node, "order", n) - this % order = n - this % n_bins = n/2 + 1 - end subroutine from_xml_zn_rad + n = this % order + num_n_bins = n/2 + 1 + end function calc_n_bins_zn_rad subroutine get_all_bins_zn_rad(this, p, estimator, match) class(ZernikeRadialFilter), intent(in) :: this From 98916c301725973f7c4f3a06b9cb87d1598572a4 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 8 Aug 2018 17:24:26 -0500 Subject: [PATCH 30/36] Change r -> radius in cylinder/sphere/cone types --- src/surface.cpp | 62 ++++++++++++++++++++++++------------------------- src/surface.h | 21 ++++++++--------- 2 files changed, 41 insertions(+), 42 deletions(-) diff --git a/src/surface.cpp b/src/surface.cpp index 8258861526..94c24c82db 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -582,18 +582,18 @@ axis_aligned_cylinder_normal(Position r, double offset1, double offset2) SurfaceXCylinder::SurfaceXCylinder(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, y0, z0, r); + read_coeffs(surf_node, id, y0, z0, radius); } double SurfaceXCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, this->r); + return axis_aligned_cylinder_evaluate<1, 2>(r, y0, z0, radius); } double SurfaceXCylinder::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<0, 1, 2>(r, u, coincident, y0, z0, - this->r); + radius); } Direction SurfaceXCylinder::normal(Position r) const @@ -605,7 +605,7 @@ Direction SurfaceXCylinder::normal(Position r) const void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cylinder", false); - std::array coeffs {{y0, z0, r}}; + std::array coeffs {{y0, z0, radius}}; write_dataset(group_id, "coefficients", coeffs); } @@ -616,18 +616,18 @@ void SurfaceXCylinder::to_hdf5_inner(hid_t group_id) const SurfaceYCylinder::SurfaceYCylinder(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, z0, r); + read_coeffs(surf_node, id, x0, z0, radius); } double SurfaceYCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, this->r); + return axis_aligned_cylinder_evaluate<0, 2>(r, x0, z0, radius); } double SurfaceYCylinder::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<1, 0, 2>(r, u, coincident, x0, z0, - this->r); + radius); } Direction SurfaceYCylinder::normal(Position r) const @@ -638,7 +638,7 @@ Direction SurfaceYCylinder::normal(Position r) const void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cylinder", false); - std::array coeffs {{x0, z0, r}}; + std::array coeffs {{x0, z0, radius}}; write_dataset(group_id, "coefficients", coeffs); } @@ -649,18 +649,18 @@ void SurfaceYCylinder::to_hdf5_inner(hid_t group_id) const SurfaceZCylinder::SurfaceZCylinder(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, y0, r); + read_coeffs(surf_node, id, x0, y0, radius); } double SurfaceZCylinder::evaluate(Position r) const { - return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, this->r); + return axis_aligned_cylinder_evaluate<0, 1>(r, x0, y0, radius); } double SurfaceZCylinder::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cylinder_distance<2, 0, 1>(r, u, coincident, x0, y0, - this->r); + radius); } Direction SurfaceZCylinder::normal(Position r) const @@ -671,7 +671,7 @@ Direction SurfaceZCylinder::normal(Position r) const void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cylinder", false); - std::array coeffs {{x0, y0, r}}; + std::array coeffs {{x0, y0, radius}}; write_dataset(group_id, "coefficients", coeffs); } @@ -682,7 +682,7 @@ void SurfaceZCylinder::to_hdf5_inner(hid_t group_id) const SurfaceSphere::SurfaceSphere(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, y0, z0, r); + read_coeffs(surf_node, id, x0, y0, z0, radius); } double SurfaceSphere::evaluate(Position r) const @@ -690,7 +690,7 @@ double SurfaceSphere::evaluate(Position r) const const double x = r.x - x0; const double y = r.y - y0; const double z = r.z - z0; - return x*x + y*y + z*z - this->r*this->r; + return x*x + y*y + z*z - radius*radius; } double SurfaceSphere::distance(Position r, Direction u, bool coincident) const @@ -699,7 +699,7 @@ double SurfaceSphere::distance(Position r, Direction u, bool coincident) const const double y = r.y - y0; const double z = r.z - z0; const double k = x*u.x + y*u.y + z*u.z; - const double c = x*x + y*y + z*z - this->r*this->r; + const double c = x*x + y*y + z*z - radius*radius; const double quad = k*k - c; if (quad < 0.0) { @@ -739,7 +739,7 @@ Direction SurfaceSphere::normal(Position r) const void SurfaceSphere::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "sphere", false); - std::array coeffs {{x0, y0, z0, r}}; + std::array coeffs {{x0, y0, z0, radius}}; write_dataset(group_id, "coefficients", coeffs); } @@ -835,29 +835,29 @@ axis_aligned_cone_normal(Position r, double offset1, double offset2, SurfaceXCone::SurfaceXCone(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, y0, z0, r_sq); + read_coeffs(surf_node, id, x0, y0, z0, radius_sq); } double SurfaceXCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, r_sq); + return axis_aligned_cone_evaluate<0, 1, 2>(r, x0, y0, z0, radius_sq); } double SurfaceXCone::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<0, 1, 2>(r, u, coincident, x0, y0, z0, - r_sq); + radius_sq); } Direction SurfaceXCone::normal(Position r) const { - return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, r_sq); + return axis_aligned_cone_normal<0, 1, 2>(r, x0, y0, z0, radius_sq); } void SurfaceXCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "x-cone", false); - std::array coeffs {{x0, y0, z0, r_sq}}; + std::array coeffs {{x0, y0, z0, radius_sq}}; write_dataset(group_id, "coefficients", coeffs); } @@ -868,29 +868,29 @@ void SurfaceXCone::to_hdf5_inner(hid_t group_id) const SurfaceYCone::SurfaceYCone(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, y0, z0, r_sq); + read_coeffs(surf_node, id, x0, y0, z0, radius_sq); } double SurfaceYCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, r_sq); + return axis_aligned_cone_evaluate<1, 0, 2>(r, y0, x0, z0, radius_sq); } double SurfaceYCone::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<1, 0, 2>(r, u, coincident, y0, x0, z0, - r_sq); + radius_sq); } Direction SurfaceYCone::normal(Position r) const { - return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, r_sq); + return axis_aligned_cone_normal<1, 0, 2>(r, y0, x0, z0, radius_sq); } void SurfaceYCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "y-cone", false); - std::array coeffs {{x0, y0, z0, r_sq}}; + std::array coeffs {{x0, y0, z0, radius_sq}}; write_dataset(group_id, "coefficients", coeffs); } @@ -901,29 +901,29 @@ void SurfaceYCone::to_hdf5_inner(hid_t group_id) const SurfaceZCone::SurfaceZCone(pugi::xml_node surf_node) : Surface(surf_node) { - read_coeffs(surf_node, id, x0, y0, z0, r_sq); + read_coeffs(surf_node, id, x0, y0, z0, radius_sq); } double SurfaceZCone::evaluate(Position r) const { - return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, r_sq); + return axis_aligned_cone_evaluate<2, 0, 1>(r, z0, x0, y0, radius_sq); } double SurfaceZCone::distance(Position r, Direction u, bool coincident) const { return axis_aligned_cone_distance<2, 0, 1>(r, u, coincident, z0, x0, y0, - r_sq); + radius_sq); } Direction SurfaceZCone::normal(Position r) const { - return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, r_sq); + return axis_aligned_cone_normal<2, 0, 1>(r, z0, x0, y0, radius_sq); } void SurfaceZCone::to_hdf5_inner(hid_t group_id) const { write_string(group_id, "type", "z-cone", false); - std::array coeffs {{x0, y0, z0, r_sq}}; + std::array coeffs {{x0, y0, z0, radius_sq}}; write_dataset(group_id, "coefficients", coeffs); } diff --git a/src/surface.h b/src/surface.h index b267af05ca..62e0111ca0 100644 --- a/src/surface.h +++ b/src/surface.h @@ -74,10 +74,9 @@ public: bool sense(Position r, Direction u) const; //! Determine the direction of a ray reflected from the surface. - //! \param r The point at which the ray is incident. - //! \param u A direction. This is both an input and an output parameter. - //! It specifies the icident direction on input and the reflected direction - //! on output. + //! \param[in] r The point at which the ray is incident. + //! \param[in] u Incident direction of the ray + //! \return Outgoing direction of the ray Direction reflect(Position r, Direction u) const; //! Evaluate the equation describing the surface. @@ -227,7 +226,7 @@ public: class SurfaceXCylinder : public Surface { - double y0, z0, r; + double y0, z0, radius; public: explicit SurfaceXCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -245,7 +244,7 @@ public: class SurfaceYCylinder : public Surface { - double x0, z0, r; + double x0, z0, radius; public: explicit SurfaceYCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -263,7 +262,7 @@ public: class SurfaceZCylinder : public Surface { - double x0, y0, r; + double x0, y0, radius; public: explicit SurfaceZCylinder(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -281,7 +280,7 @@ public: class SurfaceSphere : public Surface { - double x0, y0, z0, r; + double x0, y0, z0, radius; public: explicit SurfaceSphere(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -299,7 +298,7 @@ public: class SurfaceXCone : public Surface { - double x0, y0, z0, r_sq; + double x0, y0, z0, radius_sq; public: explicit SurfaceXCone(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -317,7 +316,7 @@ public: class SurfaceYCone : public Surface { - double x0, y0, z0, r_sq; + double x0, y0, z0, radius_sq; public: explicit SurfaceYCone(pugi::xml_node surf_node); double evaluate(Position r) const; @@ -335,7 +334,7 @@ public: class SurfaceZCone : public Surface { - double x0, y0, z0, r_sq; + double x0, y0, z0, radius_sq; public: explicit SurfaceZCone(pugi::xml_node surf_node); double evaluate(Position r) const; From fd6c8906aca1c3cc1600030e71277725757f20f1 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 09:42:36 -0500 Subject: [PATCH 31/36] Change order --- src/tallies/tally_header.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index c6f3772e1f..7793739a49 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -339,11 +339,11 @@ contains type is (ZernikeFilter) j = FILTER_ZERNIKE this % estimator = ESTIMATOR_COLLISION - type is (ParticleFilter) - j = FILTER_PARTICLE type is (ZernikeRadialFilter) j = FILTER_ZERNIKE_RADIAL this % estimator = ESTIMATOR_COLLISION + type is (ParticleFilter) + j = FILTER_PARTICLE end select this % find_filter(j) = i end do From 960211de0cee2a04ce538e8dd87dbdc56f4e67a2 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 11:05:16 -0500 Subject: [PATCH 32/36] Add ZernikeRadialFilter in Sphinx --- docs/source/pythonapi/base.rst | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b9eb7033f1..752cf417e8 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -122,6 +122,7 @@ Constructing Tallies openmc.SpatialLegendreFilter openmc.SphericalHarmonicsFilter openmc.ZernikeFilter + openmc.ZernikeRadialFilter openmc.ParticleFilter openmc.Mesh openmc.Trigger From 11723c8263e17161a11e3d88645d2710f4c6c36a Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 11:15:22 -0500 Subject: [PATCH 33/36] Fortran index starts from 1. --- src/tallies/tally_filter_zernike.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 96e016d2ee..b680fc06d0 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -133,7 +133,7 @@ contains if (bin <= last) then first = last - n m = -n + (bin - first)*2 - label = "Zernike expansion, Z" // trim(to_str(n)) // "," & + label = "Zernike expansion, Z" // trim(to_str((n-1))) // "," & // trim(to_str(m)) exit end if From 7f015db3434f909c0d0e62d06d0f543bfc7eb1a8 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Thu, 9 Aug 2018 12:41:30 -0500 Subject: [PATCH 34/36] Change the description of ZernikeRadialFilter --- openmc/filter_expansion.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index ad761f2ff7..7debebd6b9 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -466,32 +466,30 @@ class ZernikeFilter(ExpansionFilter): class ZernikeRadialFilter(ZernikeFilter): - r"""Score radial Zernike expansion moments in space up to specified order. + r"""Score the :math:`m = 0` (radial variation only) Zernike moments up to + specified order. - The Zernike polynomials are defined the same as in ZernikeFilter. + The Zernike polynomials are defined the same as in :class:`ZernikeFilter`. .. math:: - Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0 - Z_n^{m}(\rho, \theta) = R_n^{m}(\rho) \sin (m\theta), \quad m < 0 - - Z_n^{m}(\rho, \theta) = R_n^{m}(\rho), \quad m = 0 + Z_n^{0}(\rho, \theta) = R_n^{0}(\rho) where the radial polynomials are .. math:: - R_n^m(\rho) = \sum\limits_{k=0}^{(n-m)/2} \frac{(-1)^k (n-k)!}{k! ( - \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}. + R_n^{0}(\rho) = \sum\limits_{k=0}^{n/2} \frac{(-1)^k (n-k)!}{k! (( + \frac{n}{2} - k)!)^{2}} \rho^{n-2k}. - With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk - is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where - :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise. + With this definition, the integral of :math:`(Z_n^0)^2` over the unit disk + is :math:`\frac{\pi}{n+1}`. If there is only radial dependency, the polynomials are integrated over the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders. Therefore, for a radial Zernike polynomials up to order of :math:`n`, - there are :math:`\frac{n}{2} + 1` terms in total. + there are :math:`\frac{n}{2} + 1` terms in total. The indexing is from the + lowest even order (0) to highest even order. Parameters ---------- From e2774801a5f509b30e539cb7f257d583aec1480b Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Fri, 10 Aug 2018 09:22:08 -0500 Subject: [PATCH 35/36] Fixed the typo --- src/tallies/tally_filter_zernike.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index b680fc06d0..73fda7fade 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -133,7 +133,7 @@ contains if (bin <= last) then first = last - n m = -n + (bin - first)*2 - label = "Zernike expansion, Z" // trim(to_str((n-1))) // "," & + label = "Zernike expansion, Z" // trim(to_str(n)) // "," & // trim(to_str(m)) exit end if @@ -196,7 +196,7 @@ contains integer, intent(in) :: bin character(MAX_LINE_LEN) :: label - label = "Zernike expansion, Z" // trim(to_str(2*bin)) // ",0" + label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0" end function text_label_zn_rad !=============================================================================== From daf57e663ac10773274cb95b193a058b7fa08548 Mon Sep 17 00:00:00 2001 From: Zhuoran Han Date: Mon, 13 Aug 2018 08:58:39 -0500 Subject: [PATCH 36/36] Fixed some issues --- openmc/capi/filter.py | 2 +- src/tallies/tally_filter_zernike.F90 | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/openmc/capi/filter.py b/openmc/capi/filter.py index 5bd451b300..8d9d794d16 100644 --- a/openmc/capi/filter.py +++ b/openmc/capi/filter.py @@ -360,7 +360,7 @@ class ZernikeFilter(Filter): _dll.openmc_zernike_filter_set_order(self._index, order) -class ZernikeRadialFilter(Filter): +class ZernikeRadialFilter(ZernikeFilter): filter_type = 'zernikeradial' diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 index 73fda7fade..04515cf318 100644 --- a/src/tallies/tally_filter_zernike.F90 +++ b/src/tallies/tally_filter_zernike.F90 @@ -55,13 +55,13 @@ contains ! ZernikeFilter methods !=============================================================================== - function calc_n_bins_zn(this) result(num_n_bins) + function calc_n_bins_zn(this) result(n_bins) class(ZernikeFilter), intent(in) :: this integer :: n - integer :: num_n_bins + integer :: n_bins n = this % order - num_n_bins = ((n+1) * (n+2))/2 + n_bins = ((n+1) * (n+2))/2 end function calc_n_bins_zn subroutine from_xml_zn(this, node)