From 83a2be6ec51f793a1ac300ecf2c940f7a6327476 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 21 Jul 2015 22:10:15 -0600 Subject: [PATCH 1/8] Make source code style checking script --- tests/check_source.py | 235 ++++++++++++++++++++++++++++++++++++++++++ tests/travis.sh | 1 + 2 files changed, 236 insertions(+) create mode 100755 tests/check_source.py diff --git a/tests/check_source.py b/tests/check_source.py new file mode 100755 index 0000000000..490eb9c7c3 --- /dev/null +++ b/tests/check_source.py @@ -0,0 +1,235 @@ +#!/usr/bin/env python + +from __future__ import print_function + +import glob +from string import whitespace +from sys import exit +from textwrap import TextWrapper + + +################################################################################ +# Fortran reading/parsing backend. +################################################################################ + + +class LineOfCode(object): + """Contains and provides basic info about a line of Fortran 90 code.""" + def __init__(self, number, content, is_continuation=False, + string_terminator=None): + # Initialize member variables. + self.number = number + self.content = content + self.is_continuation = is_continuation + self.is_continued = False + assert (string_terminator == None or + string_terminator == "'" or + string_terminator == '"') + self.initial_string_terminator = string_terminator + self.final_string_terminator = None + self.is_comment_only = False + + # Parse the string. + self.parse() + + def __repr__(self): + rep = 'LineOfCode: line number = {0:d}\n'.format(self.number) + rep += '\tis_continuation = ' + str(self.is_continuation) + '\n' + rep += '\tis_continued = ' + str(self.is_continued) + '\n' + rep += ('\tinitial_string_terminator = ' + + str(self.initial_string_terminator) + '\n') + rep += ('\tfinal_string_terminator = ' + + str(self.final_string_terminator) + '\n') + rep += '\tis_comment_only = ' + str(self.is_comment_only) + '\n' + rep += '\tContent:\n' + self.content + return rep + + def __str__(self): + return repr(self) + + def parse(self): + # Initialize the string variables. + if self.initial_string_terminator == "'": + in_a_single_quote = True + in_a_double_quote = False + elif self.initial_string_terminator == '"': + in_a_single_quote = False + in_a_double_quote = True + else: + in_a_single_quote = False + in_a_double_quote = False + + # Initialize other variables. + in_a_comment = False + ends_with_ampersand = False + has_real_content = False + + # Parse through the line. + for char in self.content: + # Check for the start or end of strings. + if char == "'" and in_a_single_quote: + in_a_single_quote = False + elif char == "'" and not in_a_double_quote: + in_a_single_quote = True + elif char == '"' and in_a_double_quote: + in_a_double_quote = False + elif char == '"' and not in_a_single_quote: + in_a_double_quote = True + + # Check for the start of a comment. + if (char == "!" and not in_a_single_quote + and not in_a_double_quote): + in_a_comment = True + break # Don't care about anything in a comment. + + # Check for a continuation ampersand. + if char == "&": + ends_with_ampersand = True + elif all([char != w for w in whitespace]): + ends_with_ampersand = False + + # Check to see if there is any real content (non-whitespace, not in + # a comment) + if (not in_a_comment) and all([char != w for w in whitespace]): + has_real_content = True + + # Make sure nothing went terribly wrong. + assert not (in_a_single_quote and in_a_double_quote) + assert not (in_a_single_quote and in_a_comment) + assert not (in_a_double_quote and in_a_comment) + + # Is this line continued onto the next? + if ends_with_ampersand: + self.is_continued = True + + # Is a multiline string continued on the next line? + if in_a_single_quote: + assert self.is_continued + self.final_string_terminator = "'" + if in_a_double_quote: + assert self.is_continued + self.final_string_terminator = '"' + + # Is this line a comment-only line? + if (not has_real_content) and in_a_comment: + self.is_comment_only = True + + + def get_indent(self): + """Return the number of indentation spaces prefixing this line.""" + return len(self.content) - len(self.content.lstrip(' ')) + + + def contains_whitespace_only(self): + """Return True/False if all characters in the line are whitespace.""" + if len(self.content.strip('\n')) == 0: return False + is_ws = [ any([char == ws for ws in whitespace]) + for char in self.content.strip('\n') ] + return all(is_ws) + + + def has_trailing_whitespace(self): + """Return True/False if this line ends with a whitespace character.""" + stripped = self.content.strip('\n') + if len(stripped) == 0: return False + return any([stripped[-1] == ws for ws in whitespace]) + + def contains_tab(self): + """Return True/False if a tab character appears in the line.""" + return '\t' in self.content + + +def read_lines_of_code(fname): + line_num = 0 + cont = False + str_term = None + with open(fname) as fin: + for line in fin: + loc = LineOfCode(line_num, line, is_continuation=cont, + string_terminator=str_term) + cont = loc.is_continued + str_term = loc.final_string_terminator + line_num += 1 + yield loc + + +################################################################################ +# Error checking. +################################################################################ + + +def print_error(fname, line_number, err_msg): + header = "Error in file {0}, line {1}:".format(fname, line_number + 1) + + tw = TextWrapper(width=80, initial_indent=' '*4, subsequent_indent=' '*4) + body = '\n'.join(tw.wrap(err_msg)) + + print(header + '\n' + body + '\n') + + +def check_source(fname): + """Make sure the given Fortran source file meets OpenMC standards. + + If errors are found, messages will be printed to the screen describing the + error. The function will return True if no errors were found or False + otherwise. + """ + good_code = True + base_indent = 0 + + for loc in read_lines_of_code(fname): + # Check for extra whitespace errors. + if loc.contains_whitespace_only(): + good_code = False + print_error(fname, loc.number, "Line contains whitespace " \ + "characters but no content. Please remove whitespace.") + elif loc.has_trailing_whitespace(): + good_code = False + print_error(fname, loc.number, "Line has trailing whitespace after"\ + " the content. Please remove trailing whitespace.") + if loc.contains_tab(): + good_code = False + print_error(fname, loc.number, "Line contains a tab character. " \ + "Please replace with single whitespace characters.") + + # Check indentation. + current_indent = loc.get_indent() + if ((not loc.is_continuation) and (not loc.is_comment_only) and + (not loc.contains_whitespace_only()) and current_indent % 2 != 0): + good_code = False + print_error(fname, loc.number, "Line is indented an odd number of "\ + "spaces. All non-continuation lines should be indented an "\ + "even number of spaces.") + if loc.is_continuation and current_indent < base_indent + 5: + good_code = False + print_error(fname, loc.number, "Continuation lines must be "\ + "indented by at least 5 spaces, but this line is indented {0}"\ + " spaces.".format(current_indent - base_indent)) + + # Set base indentation for next lines. + if not loc.is_continuation: + base_indent = current_indent + + return good_code + + +################################################################################ +# Main loop. +################################################################################ + + +if __name__ == '__main__': + # Get a list of the F90 source files. + f90_files = glob.glob('../src/*.F90') + + # Make sure we found the source files. + assert len(f90_files) != 0, 'No .F90 source files found.' + + # Make sure all the F90 source files meet our standards. + good_code = [check_source(fname) for fname in f90_files] + if not all(good_code): + print("ERROR: The Fortran source code does not meet OpenMC's standards") + exit(-1) + else: + print("SUCCESS! Looks like the Fortran source meets our standards") + exit(0) diff --git a/tests/travis.sh b/tests/travis.sh index d591428686..fac8d793a6 100755 --- a/tests/travis.sh +++ b/tests/travis.sh @@ -3,6 +3,7 @@ set -ev # Run all debug tests +./check_source.py if [ "$TRAVIS_PULL_REQUEST" != "false" ]; then ./run_tests.py -C "^basic-debug$|^hdf5-debug$|^mpi-omp-debug$|^phdf5-omp-debug$" -j 2 -s else From 0f273878108e779626a3dc57d0228c6e91c14740 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 21 Jul 2015 23:39:43 -0600 Subject: [PATCH 2/8] Make source code match style guide --- src/ace.F90 | 27 +++-- src/cmfd_data.F90 | 30 ++--- src/cmfd_execute.F90 | 2 +- src/cmfd_header.F90 | 10 +- src/cmfd_input.F90 | 18 +-- src/cmfd_loss_operator.F90 | 20 +-- src/cmfd_prod_operator.F90 | 4 +- src/cmfd_solver.F90 | 6 +- src/cross_section.F90 | 10 +- src/dict_header.F90 | 22 ++-- src/doppler.F90 | 205 +++++++++++++++---------------- src/eigenvalue.F90 | 12 +- src/endf_header.F90 | 16 +-- src/geometry_header.F90 | 69 ++++++----- src/global.F90 | 6 +- src/hdf5_interface.F90 | 2 +- src/hdf5_summary.F90 | 20 +-- src/initialize.F90 | 4 +- src/input_xml.F90 | 238 ++++++++++++++++++------------------ src/list_header.F90 | 18 +-- src/math.F90 | 241 +++++++++++++++++++------------------ src/matrix_header.F90 | 26 ++-- src/mpiio_interface.F90 | 28 ++--- src/output.F90 | 51 ++++---- src/output_interface.F90 | 28 ++--- src/physics.F90 | 34 +++--- src/plot.F90 | 2 +- src/plot_header.F90 | 4 +- src/ppmlib.F90 | 12 +- src/progress_header.F90 | 20 +-- src/random_lcg.F90 | 2 +- src/search.F90 | 6 +- src/state_point.F90 | 28 ++--- src/string.F90 | 4 +- src/tally_header.F90 | 6 +- src/timer_header.F90 | 12 +- src/trigger.F90 | 6 +- src/trigger_header.F90 | 6 +- src/vector_header.F90 | 12 +- src/xml_interface.F90 | 50 ++++---- 40 files changed, 669 insertions(+), 648 deletions(-) diff --git a/src/ace.F90 b/src/ace.F90 index c1845fffb4..0dd92ea2c7 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -88,16 +88,16 @@ contains nuclides(i_nuclide) % resonant = .true. nuclides(i_nuclide) % name_0K = nuclides_0K(n) % name_0K nuclides(i_nuclide) % name_0K = trim(nuclides(i_nuclide) % & - & name_0K) + & name_0K) nuclides(i_nuclide) % scheme = nuclides_0K(n) % scheme nuclides(i_nuclide) % scheme = trim(nuclides(i_nuclide) % & - & scheme) + & scheme) nuclides(i_nuclide) % E_min = nuclides_0K(n) % E_min nuclides(i_nuclide) % E_max = nuclides_0K(n) % E_max if (.not. already_read % contains(nuclides(i_nuclide) % & - & name_0K)) then + & name_0K)) then i_listing = xs_listing_dict % get_key(nuclides(i_nuclide) % & - & name_0K) + & name_0K) call read_ace_table(i_nuclide, i_listing) end if exit @@ -446,9 +446,10 @@ contains if (nuc % elastic_0K(i) < ZERO) nuc % elastic_0K(i) = ZERO ! build xs cdf - xs_cdf_sum = xs_cdf_sum + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) & - & + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO & - & * (nuc % energy_0K(i+1) - nuc % energy_0K(i)) + xs_cdf_sum = xs_cdf_sum & + + (sqrt(nuc % energy_0K(i)) * nuc % elastic_0K(i) & + + sqrt(nuc % energy_0K(i+1)) * nuc % elastic_0K(i+1)) / TWO & + * (nuc % energy_0K(i+1) - nuc % energy_0K(i)) nuc % xs_cdf(i) = xs_cdf_sum end do @@ -752,31 +753,31 @@ contains ! Set flag and allocate space for Tab1 to store yield rxn % multiplicity_with_E = .true. allocate(rxn % multiplicity_E) - + XSS_index = JXS(11) + rxn % multiplicity - 101 NR = nint(XSS(XSS_index)) rxn % multiplicity_E % n_regions = NR - + ! allocate space for ENDF interpolation parameters if (NR > 0) then allocate(rxn % multiplicity_E % nbt(NR)) allocate(rxn % multiplicity_E % int(NR)) end if - + ! read ENDF interpolation parameters XSS_index = XSS_index + 1 if (NR > 0) then rxn % multiplicity_E % nbt = get_int(NR) rxn % multiplicity_E % int = get_int(NR) end if - + ! allocate space for yield data XSS_index = XSS_index + 2*NR NE = nint(XSS(XSS_index)) rxn % multiplicity_E % n_pairs = NE allocate(rxn % multiplicity_E % x(NE)) allocate(rxn % multiplicity_E % y(NE)) - + ! read yield data XSS_index = XSS_index + 1 rxn % multiplicity_E % x = get_real(NE) @@ -1480,7 +1481,7 @@ contains table % inelastic_data(i) % e_out_pdf(j) = XSS(XSS_index + 2) table % inelastic_data(i) % e_out_cdf(j) = XSS(XSS_index + 3) table % inelastic_data(i) % mu(:, j) = & - XSS(XSS_index + 4: XSS_index + 4 + NMU - 1) + XSS(XSS_index + 4: XSS_index + 4 + NMU - 1) XSS_index = XSS_index + 4 + NMU - 1 end do end do diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 70a18bac78..2b17784196 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -104,23 +104,23 @@ contains cmfd % keff_bal = ZERO - ! Begin loop around tallies - TAL: do ital = 1, n_cmfd_tallies + ! Begin loop around tallies + TAL: do ital = 1, n_cmfd_tallies - ! Associate tallies and mesh - t => cmfd_tallies(ital) - i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1) - m => meshes(i_mesh) + ! Associate tallies and mesh + t => cmfd_tallies(ital) + i_mesh = t % filters(t % find_filter(FILTER_MESH)) % int_bins(1) + m => meshes(i_mesh) - i_filter_mesh = t % find_filter(FILTER_MESH) - i_filter_ein = t % find_filter(FILTER_ENERGYIN) - i_filter_eout = t % find_filter(FILTER_ENERGYOUT) - i_filter_surf = t % find_filter(FILTER_SURFACE) + i_filter_mesh = t % find_filter(FILTER_MESH) + i_filter_ein = t % find_filter(FILTER_ENERGYIN) + i_filter_eout = t % find_filter(FILTER_ENERGYOUT) + i_filter_surf = t % find_filter(FILTER_SURFACE) - ! Begin loop around space - ZLOOP: do k = 1,nz + ! Begin loop around space + ZLOOP: do k = 1,nz - YLOOP: do j = 1,ny + YLOOP: do j = 1,ny XLOOP: do i = 1,nx @@ -600,7 +600,7 @@ contains dtilde = (2*cell_dc*neig_dc)/(neig_hxyz(xyz_idx)*cell_dc + & cell_hxyz(xyz_idx)*neig_dc) - end if + end if end if @@ -797,7 +797,7 @@ contains ! Calculate albedo if ((shift_idx == 1 .and. current(2*l ) < 1.0e-10_8) .or. & - (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then + (shift_idx == -1 .and. current(2*l-1) < 1.0e-10_8)) then albedo = ONE else albedo = (current(2*l-1)/current(2*l))**(shift_idx) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index e620c1be38..c55d206cbc 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -121,7 +121,7 @@ contains ! Allocate cmfd source if not already allocated and allocate buffer if (.not. allocated(cmfd % cmfd_src)) & - allocate(cmfd % cmfd_src(ng,nx,ny,nz)) + allocate(cmfd % cmfd_src(ng,nx,ny,nz)) ! Reset cmfd source to 0 cmfd % cmfd_src = ZERO diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index c9bc0f32af..e743cc928e 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -1,4 +1,4 @@ -module cmfd_header +module cmfd_header use constants, only: CMFD_NOACCEL, ZERO, ONE @@ -17,7 +17,7 @@ module cmfd_header ! Core overlay map integer, allocatable :: coremap(:,:,:) integer, allocatable :: indexmap(:,:) - integer :: mat_dim = CMFD_NOACCEL + integer :: mat_dim = CMFD_NOACCEL ! Energy grid real(8), allocatable :: egrid(:) @@ -51,7 +51,7 @@ module cmfd_header ! Source sites in each mesh box real(8), allocatable :: sourcecounts(:,:,:,:) - ! Weight adjustment factors + ! Weight adjustment factors real(8), allocatable :: weightfactors(:,:,:,:) ! Eigenvector/eigenvalue from cmfd run @@ -118,7 +118,7 @@ contains if (.not. allocated(this % nfissxs)) allocate(this % nfissxs(ng,ng,nx,ny,nz)) if (.not. allocated(this % diffcof)) allocate(this % diffcof(ng,nx,ny,nz)) - ! Allocate dtilde and dhat + ! Allocate dtilde and dhat if (.not. allocated(this % dtilde)) allocate(this % dtilde(6,ng,nx,ny,nz)) if (.not. allocated(this % dhat)) allocate(this % dhat(6,ng,nx,ny,nz)) @@ -167,7 +167,7 @@ contains end subroutine allocate_cmfd !=============================================================================== -! DEALLOCATE_CMFD frees all memory of cmfd type +! DEALLOCATE_CMFD frees all memory of cmfd type !=============================================================================== subroutine deallocate_cmfd(this) diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 169c763957..2d44d3e9bc 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -120,7 +120,7 @@ contains allocate(cmfd % coremap(cmfd % indices(1), cmfd % indices(2), & cmfd % indices(3))) if (get_arraysize_integer(node_mesh, "map") /= & - product(cmfd % indices(1:3))) then + product(cmfd % indices(1:3))) then call fatal_error('CMFD coremap not to correct dimensions') end if allocate(iarray(get_arraysize_integer(node_mesh, "map"))) @@ -156,7 +156,7 @@ contains call get_node_value(doc, "dhat_reset", temp_str) temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - dhat_reset = .true. + dhat_reset = .true. end if ! Set monitoring @@ -180,7 +180,7 @@ contains call get_node_value(doc, "run_adjoint", temp_str) temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - cmfd_run_adjoint = .true. + cmfd_run_adjoint = .true. end if ! Batch to begin cmfd @@ -289,7 +289,7 @@ contains ! Determine number of dimensions for mesh n = get_arraysize_integer(node_mesh, "dimension") if (n /= 2 .and. n /= 3) then - call fatal_error("Mesh must be two or three dimensions.") + call fatal_error("Mesh must be two or three dimensions.") end if m % n_dimension = n @@ -318,14 +318,14 @@ contains ! Make sure both upper-right or width were specified if (check_for_node(node_mesh, "upper_right") .and. & - check_for_node(node_mesh, "width")) then + check_for_node(node_mesh, "width")) then call fatal_error("Cannot specify both and on a & &tally mesh.") end if ! Make sure either upper-right or width was specified if (.not.check_for_node(node_mesh, "upper_right") .and. & - .not.check_for_node(node_mesh, "width")) then + .not.check_for_node(node_mesh, "width")) then call fatal_error("Must specify either and on a & &tally mesh.") end if @@ -333,7 +333,7 @@ contains if (check_for_node(node_mesh, "width")) then ! Check to ensure width has same dimensions if (get_arraysize_double(node_mesh, "width") /= & - get_arraysize_double(node_mesh, "lower_left")) then + get_arraysize_double(node_mesh, "lower_left")) then call fatal_error("Number of entries on must be the same as the & &number of entries on .") end if @@ -351,7 +351,7 @@ contains elseif (check_for_node(node_mesh, "upper_right")) then ! Check to ensure width has same dimensions if (get_arraysize_double(node_mesh, "upper_right") /= & - get_arraysize_double(node_mesh, "lower_left")) then + get_arraysize_double(node_mesh, "lower_left")) then call fatal_error("Number of entries on must be the same & &as the number of entries on .") end if @@ -548,7 +548,7 @@ contains ! Deallocate filter bins deallocate(filters(1) % int_bins) if (check_for_node(node_mesh, "energy")) & - deallocate(filters(2) % real_bins) + deallocate(filters(2) % real_bins) end do diff --git a/src/cmfd_loss_operator.F90 b/src/cmfd_loss_operator.F90 index b89f0f5772..7dfeed1dbd 100644 --- a/src/cmfd_loss_operator.F90 +++ b/src/cmfd_loss_operator.F90 @@ -54,7 +54,7 @@ contains n_e = 4*(nx - 2) + 4*(ny - 2) + 4*(nz - 2) ! define # of edges n_s = 2*(nx - 2)*(ny - 2) + 2*(nx - 2)*(nz - 2) & + 2*(ny - 2)*(nz - 2) ! define # of sides - n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors + n_i = nx*ny*nz - (n_c + n_e + n_s) ! define # of interiors nz_c = ng*n_c*(4 + ng - 1) ! define # nonzero corners nz_e = ng*n_e*(5 + ng - 1) ! define # nonzero edges nz_s = ng*n_s*(6 + ng - 1) ! define # nonzero sides @@ -131,7 +131,7 @@ contains ! Check for global boundary if (bound(l) /= nxyz(xyz_idx,dir_idx)) then - ! Check for coremap + ! Check for coremap if (cmfd_coremap) then ! Check for neighbor that is non-acceleartred @@ -139,11 +139,11 @@ contains CMFD_NOACCEL) then ! Get neighbor matrix index - call indices_to_matrix(g,neig_idx(1), neig_idx(2), & + call indices_to_matrix(g,neig_idx(1), neig_idx(2), & neig_idx(3), neig_mat_idx, ng, nx, ny) ! Record nonzero - nnz = nnz + 1 + nnz = nnz + 1 end if @@ -218,11 +218,11 @@ contains real(8) :: jn ! direction dependent leakage coeff to neig real(8) :: jo(6) ! leakage coeff in front of cell flux real(8) :: jnet ! net leakage from jo - real(8) :: val ! temporary variable before saving to + real(8) :: val ! temporary variable before saving to ! Check for adjoint adjoint_calc = .false. - if (present(adjoint)) adjoint_calc = adjoint + if (present(adjoint)) adjoint_calc = adjoint ! Get maximum number of cells in each direction nx = cmfd%indices(1) @@ -257,11 +257,11 @@ contains dhat = ZERO end if - ! Create boundary vector + ! Create boundary vector bound = (/i,i,j,j,k,k/) ! Begin loop over leakages - ! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z + ! 1=-x, 2=+x, 3=-y, 4=+y, 5=-z, 6=+z LEAK: do l = 1,6 ! Define (x,y,z) and (-,+) indices @@ -350,7 +350,7 @@ contains scattxshg = cmfd%scattxs(h, g, i, j, k) end if - ! Negate the scattering xs + ! Negate the scattering xs val = -scattxshg ! Record value in matrix @@ -358,7 +358,7 @@ contains end do SCATTR - end do ROWS + end do ROWS ! CSR requires n+1 row call loss_matrix % new_row() diff --git a/src/cmfd_prod_operator.F90 b/src/cmfd_prod_operator.F90 index e6d055ced8..7779c24551 100644 --- a/src/cmfd_prod_operator.F90 +++ b/src/cmfd_prod_operator.F90 @@ -97,7 +97,7 @@ contains end if - ! loop around all other groups + ! loop around all other groups NFISS: do h = 1, ng @@ -122,7 +122,7 @@ contains end do NFISS - end do ROWS + end do ROWS ! CSR requires n+1 row call prod_matrix % new_row() diff --git a/src/cmfd_solver.F90 b/src/cmfd_solver.F90 index ee5127a374..84f0016d3d 100644 --- a/src/cmfd_solver.F90 +++ b/src/cmfd_solver.F90 @@ -65,7 +65,7 @@ contains ! Check for physical adjoint if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'physical') & - physical_adjoint = .true. + physical_adjoint = .true. ! Start timer for build call time_cmfdbuild % start() @@ -75,7 +75,7 @@ contains ! Check for mathematical adjoint calculation if (adjoint_calc .and. trim(cmfd_adjoint_type) == 'math') & - call compute_adjoint() + call compute_adjoint() ! Stop timer for build call time_cmfdbuild % stop() @@ -286,7 +286,7 @@ contains ROWS: do irow = 1, loss % n COLS: do icol = loss % get_row(irow), loss % get_row(irow + 1) - 1 if (loss % get_col(icol) == prod % get_col(jcol) .and. & - jcol < prod % get_row(irow + 1)) then + jcol < prod % get_row(irow + 1)) then loss % val(icol) = loss % val(icol) - ONE/k_s*prod % val(jcol) jcol = jcol + 1 end if diff --git a/src/cross_section.F90 b/src/cross_section.F90 index c8a567938e..b937b03a15 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -446,7 +446,7 @@ contains ! Calculate elastic cross section/factor elastic = ZERO if (urr % prob(i_energy, URR_ELASTIC, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then + urr % prob(i_energy + 1, URR_ELASTIC, i_up) > ZERO) then elastic = exp((ONE - f) * log(urr % prob(i_energy, URR_ELASTIC, & i_low)) + f * log(urr % prob(i_energy + 1, URR_ELASTIC, & i_up))) @@ -455,7 +455,7 @@ contains ! Calculate fission cross section/factor fission = ZERO if (urr % prob(i_energy, URR_FISSION, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then + urr % prob(i_energy + 1, URR_FISSION, i_up) > ZERO) then fission = exp((ONE - f) * log(urr % prob(i_energy, URR_FISSION, & i_low)) + f * log(urr % prob(i_energy + 1, URR_FISSION, & i_up))) @@ -464,7 +464,7 @@ contains ! Calculate capture cross section/factor capture = ZERO if (urr % prob(i_energy, URR_N_GAMMA, i_low) > ZERO .and. & - urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then + urr % prob(i_energy + 1, URR_N_GAMMA, i_up) > ZERO) then capture = exp((ONE - f) * log(urr % prob(i_energy, URR_N_GAMMA, & i_low)) + f * log(urr % prob(i_energy + 1, URR_N_GAMMA, & i_up))) @@ -571,11 +571,11 @@ contains ! calculate interpolation factor f = (E - nuc % energy_0K(i_grid)) & - & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) + & / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid)) ! Calculate microscopic nuclide elastic cross section xs_out = (ONE - f) * nuc % elastic_0K(i_grid) & - & + f * nuc % elastic_0K(i_grid + 1) + & + f * nuc % elastic_0K(i_grid + 1) end function elastic_xs_0K diff --git a/src/dict_header.F90 b/src/dict_header.F90 index 7da2015fd3..ffe488d525 100644 --- a/src/dict_header.F90 +++ b/src/dict_header.F90 @@ -56,7 +56,7 @@ module dict_header !=============================================================================== ! DICT* is a dictionary of (key,value) pairs with convenience methods as ! type-bound procedures. DictCharInt has character(*) keys and integer values, -! and DictIntInt has integer keys and values. +! and DictIntInt has integer keys and values. !=============================================================================== type, public :: DictCharInt @@ -105,7 +105,7 @@ contains if (associated(elem)) then elem % value = value else - ! Get hash + ! Get hash hash = dict_hash_key_ci(key) ! Create new element @@ -135,7 +135,7 @@ contains if (associated(elem)) then elem % value = value else - ! Get hash + ! Get hash hash = dict_hash_key_ii(key) ! Create new element @@ -156,13 +156,13 @@ contains !=============================================================================== function dict_get_elem_ci(this, key) result(elem) - + class(DictCharInt) :: this character(*), intent(in) :: key type(ElemKeyValueCI), pointer :: elem integer :: hash - + ! Check for dictionary not being allocated if (.not. associated(this % table)) then allocate(this % table(HASH_SIZE)) @@ -178,13 +178,13 @@ contains end function dict_get_elem_ci function dict_get_elem_ii(this, key) result(elem) - + class(DictIntInt) :: this integer, intent(in) :: key type(ElemKeyValueII), pointer :: elem integer :: hash - + ! Check for dictionary not being allocated if (.not. associated(this % table)) then allocate(this % table(HASH_SIZE)) @@ -279,7 +279,7 @@ contains character(*), intent(in) :: key integer :: val - + integer :: i val = 0 @@ -298,7 +298,7 @@ contains integer, intent(in) :: key integer :: val - + val = 0 ! Added the absolute val on val-1 since the sum in the do loop is @@ -420,7 +420,7 @@ contains integer :: i type(ElemKeyValueII), pointer :: current type(ElemKeyValueII), pointer :: next - + if (associated(this % table)) then do i = 1, size(this % table) current => this % table(i) % list @@ -440,5 +440,5 @@ contains end if end subroutine dict_clear_ii - + end module dict_header diff --git a/src/doppler.F90 b/src/doppler.F90 index e888eb2b18..85ab6da3c0 100644 --- a/src/doppler.F90 +++ b/src/doppler.F90 @@ -52,62 +52,115 @@ contains ! Loop over incoming neutron energies ENERGY_NEUTRON: do i = 1, n - sigma = ZERO - y = x(i) - y_sq = y*y - y_inv = ONE / y - y_inv_sq = y_inv / y + sigma = ZERO + y = x(i) + y_sq = y*y + y_inv = ONE / y + y_inv_sq = y_inv / y - ! ======================================================================= - ! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4 + ! ======================================================================= + ! EVALUATE FIRST TERM FROM x(k) - y = 0 to -4 - k = i - a = ZERO - call calculate_F(F_a, a) + k = i + a = ZERO + call calculate_F(F_a, a) - do while (a >= -4.0 .and. k > 1) - ! Move to next point - F_b = F_a - k = k - 1 - a = x(k) - y + do while (a >= -4.0 .and. k > 1) + ! Move to next point + F_b = F_a + k = k - 1 + a = x(k) - y - ! Calculate F and H functions - call calculate_F(F_a, a) - H = F_a - F_b + ! Calculate F and H functions + call calculate_F(F_a, a) + H = F_a - F_b - ! Calculate A(k), B(k), and slope terms - Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0) - Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0) - slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2) + ! Calculate A(k), B(k), and slope terms + Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0) + Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0) + slope = (xs(k+1) - xs(k)) / (x(k+1)**2 - x(k)**2) - ! Add contribution to broadened cross section - sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk - end do + ! Add contribution to broadened cross section + sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk + end do - ! ======================================================================= - ! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE + ! ======================================================================= + ! EXTEND CROSS SECTION TO 0 ASSUMING 1/V SHAPE - if (k == 1 .and. a >= -4.0) then - ! Since x = 0, this implies that a = -y - F_b = F_a - a = -y + if (k == 1 .and. a >= -4.0) then + ! Since x = 0, this implies that a = -y + F_b = F_a + a = -y - ! Calculate F and H functions - call calculate_F(F_a, a) - H = F_a - F_b + ! Calculate F and H functions + call calculate_F(F_a, a) + H = F_a - F_b - ! Add contribution to broadened cross section - sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0)) - end if + ! Add contribution to broadened cross section + sigma = sigma + xs(k)*x(k)*(y_inv_sq*H(1) + y_inv*H(0)) + end if - ! ======================================================================= - ! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4 + ! ======================================================================= + ! EVALUATE FIRST TERM FROM x(k) - y = 0 to 4 - k = i - b = ZERO - call calculate_F(F_b, b) + k = i + b = ZERO + call calculate_F(F_b, b) - do while (b <= 4.0 .and. k < n) + do while (b <= 4.0 .and. k < n) + ! Move to next point + F_a = F_b + k = k + 1 + b = x(k) - y + + ! Calculate F and H functions + call calculate_F(F_b, b) + H = F_a - F_b + + ! Calculate A(k), B(k), and slope terms + Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0) + Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0) + slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2) + + ! Add contribution to broadened cross section + sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk + end do + + ! ======================================================================= + ! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE + + if (k == n .and. b <= 4.0) then + ! Calculate F function at last energy point + a = x(k) - y + call calculate_F(F_a, a) + + ! Add contribution to broadened cross section + sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0)) + end if + + ! ======================================================================= + ! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4 + + if (y <= 4.0) then + ! Swap signs on y + y = -y + y_inv = -y_inv + k = 1 + + ! Calculate a and b based on 0 and x(1) + a = -y + b = x(k) - y + + ! Calculate F and H functions + call calculate_F(F_a, a) + call calculate_F(F_b, b) + H = F_a - F_b + + ! Add contribution to broadened cross section + sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0)) + + ! Now progress forward doing the remainder of the second term + do while (b <= 4.0) ! Move to next point F_a = F_b k = k + 1 @@ -119,69 +172,17 @@ contains ! Calculate A(k), B(k), and slope terms Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0) - Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0) + Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) & + + y_sq*H(0) slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2) ! Add contribution to broadened cross section - sigma = sigma + Ak*(xs(k) - slope*x(k)**2) + slope*Bk - end do + sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk + end do + end if - ! ======================================================================= - ! EXTEND CROSS SECTION TO INFINITY ASSUMING CONSTANT SHAPE - - if (k == n .and. b <= 4.0) then - ! Calculate F function at last energy point - a = x(k) - y - call calculate_F(F_a, a) - - ! Add contribution to broadened cross section - sigma = sigma + xs(k) * (y_inv_sq*F_a(2) + 2.0*y_inv*F_a(1) + F_a(0)) - end if - - ! ======================================================================= - ! EVALUATE SECOND TERM FROM x(k) + y = 0 to +4 - - if (y <= 4.0) then - ! Swap signs on y - y = -y - y_inv = -y_inv - k = 1 - - ! Calculate a and b based on 0 and x(1) - a = -y - b = x(k) - y - - ! Calculate F and H functions - call calculate_F(F_a, a) - call calculate_F(F_b, b) - H = F_a - F_b - - ! Add contribution to broadened cross section - sigma = sigma - xs(k) * x(k) * (y_inv_sq*H(1) + y_inv*H(0)) - - ! Now progress forward doing the remainder of the second term - do while (b <= 4.0) - ! Move to next point - F_a = F_b - k = k + 1 - b = x(k) - y - - ! Calculate F and H functions - call calculate_F(F_b, b) - H = F_a - F_b - - ! Calculate A(k), B(k), and slope terms - Ak = y_inv_sq*H(2) + 2.0*y_inv*H(1) + H(0) - Bk = y_inv_sq*H(4) + 4.0*y_inv*H(3) + 6.0*H(2) + 4.0*y*H(1) + y_sq*H(0) - slope = (xs(k) - xs(k-1)) / (x(k)**2 - x(k-1)**2) - - ! Add contribution to broadened cross section - sigma = sigma - Ak*(xs(k) - slope*x(k)**2) - slope*Bk - end do - end if - - ! Set broadened cross section - sigmaNew(i) = sigma + ! Set broadened cross section + sigmaNew(i) = sigma end do ENERGY_NEUTRON diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index aefb597d0e..1bbcbe850e 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -259,7 +259,7 @@ contains ! Write out source point if it's been specified for this batch if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. & - source_write) then + source_write) then call write_source_point() end if @@ -880,8 +880,8 @@ contains !$omp do ordered schedule(static) do i = 1, n_threads !$omp ordered - master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank) - total = total + n_bank + master_fission_bank(total+1:total+n_bank) = fission_bank(1:n_bank) + total = total + n_bank !$omp end ordered end do !$omp end do @@ -891,10 +891,10 @@ contains ! Now copy the shared fission bank sites back to the master thread's copy. if (thread_id == 0) then - n_bank = total - fission_bank(1:n_bank) = master_fission_bank(1:n_bank) + n_bank = total + fission_bank(1:n_bank) = master_fission_bank(1:n_bank) else - n_bank = 0 + n_bank = 0 end if !$omp end parallel diff --git a/src/endf_header.F90 b/src/endf_header.F90 index a2e8b10de4..54af0f7383 100644 --- a/src/endf_header.F90 +++ b/src/endf_header.F90 @@ -3,7 +3,7 @@ module endf_header implicit none !=============================================================================== -! TAB1 represents a one-dimensional interpolable function +! TAB1 represents a one-dimensional interpolable function !=============================================================================== type Tab1 @@ -13,28 +13,28 @@ module endf_header integer :: n_pairs ! # of pairs of (x,y) values real(8), allocatable :: x(:) ! values of abscissa real(8), allocatable :: y(:) ! values of ordinate - + ! Type-Bound procedures contains procedure :: clear => tab1_clear ! deallocates a Tab1 Object. end type Tab1 - + contains - + !=============================================================================== ! TAB1_CLEAR deallocates the items in Tab1 !=============================================================================== subroutine tab1_clear(this) - + class(Tab1), intent(inout) :: this ! The Tab1 to clear - + if (allocated(this % nbt)) & deallocate(this % nbt, this % int) - + if (allocated(this % x)) & deallocate(this % x, this % y) - + end subroutine tab1_clear end module endf_header diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 9a5778c173..93e5dd1fdb 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -9,13 +9,13 @@ module geometry_header !=============================================================================== type Universe - integer :: id ! Unique ID - integer :: type ! Type - integer :: n_cells ! # of cells within - integer, allocatable :: cells(:) ! List of cells within - real(8) :: x0 ! Translation in x-coordinate - real(8) :: y0 ! Translation in y-coordinate - real(8) :: z0 ! Translation in z-coordinate + integer :: id ! Unique ID + integer :: type ! Type + integer :: n_cells ! # of cells within + integer, allocatable :: cells(:) ! List of cells within + real(8) :: x0 ! Translation in x-coordinate + real(8) :: y0 ! Translation in y-coordinate + real(8) :: z0 ! Translation in z-coordinate end type Universe !=============================================================================== @@ -119,14 +119,14 @@ module geometry_header !=============================================================================== type Surface - integer :: id ! Unique ID - character(len=52) :: name = "" ! User-defined name - integer :: type ! Type of surface - real(8), allocatable :: coeffs(:) ! Definition of surface - integer, allocatable :: & - neighbor_pos(:), & ! List of cells on positive side - neighbor_neg(:) ! List of cells on negative side - integer :: bc ! Boundary condition + integer :: id ! Unique ID + character(len=52) :: name = "" ! User-defined name + integer :: type ! Type of surface + real(8), allocatable :: coeffs(:) ! Definition of surface + integer, allocatable :: & + neighbor_pos(:), & ! List of cells on positive side + neighbor_neg(:) ! List of cells on negative side + integer :: bc ! Boundary condition end type Surface !=============================================================================== @@ -134,24 +134,29 @@ module geometry_header !=============================================================================== type Cell - integer :: id ! Unique ID - character(len=52) :: name = "" ! User-defined name - integer :: type ! Type of cell (normal, universe, lattice) - integer :: universe ! universe # this cell is in - integer :: fill ! universe # filling this cell - integer :: instances ! number of instances of this cell in the geom - integer :: material ! Material within cell (0 for universe) - integer :: n_surfaces ! Number of surfaces within - integer, allocatable :: offset (:) ! Distribcell offset for tally counter - integer, allocatable :: & - & surfaces(:) ! List of surfaces bounding cell -- note that - ! parentheses, union, etc operators will be listed - ! here too + integer :: id ! Unique ID + character(len=52) :: name = "" ! User-defined name + integer :: type ! Type of cell (normal, universe, + ! lattice) + integer :: universe ! universe # this cell is in + integer :: fill ! universe # filling this cell + integer :: instances ! number of instances of this cell in + ! the geom + integer :: material ! Material within cell (0 for + ! universe) + integer :: n_surfaces ! Number of surfaces within + integer, allocatable :: offset (:) ! Distribcell offset for tally + ! counter + integer, allocatable :: & + & surfaces(:) ! List of surfaces bounding cell + ! -- note that parentheses, union, + ! etc operators will be listed here + ! too - ! Rotation matrix and translation vector - real(8), allocatable :: translation(:) - real(8), allocatable :: rotation(:) - real(8), allocatable :: rotation_matrix(:,:) + ! Rotation matrix and translation vector + real(8), allocatable :: translation(:) + real(8), allocatable :: rotation(:) + real(8), allocatable :: rotation_matrix(:,:) end type Cell ! array index of universe 0 diff --git a/src/global.F90 b/src/global.F90 index f7d9f184e1..a5ff7453ea 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -529,9 +529,9 @@ contains ! Deallocate entropy mesh if (associated(entropy_mesh)) then if (allocated(entropy_mesh % lower_left)) & - deallocate(entropy_mesh % lower_left) + deallocate(entropy_mesh % lower_left) if (allocated(entropy_mesh % upper_right)) & - deallocate(entropy_mesh % upper_right) + deallocate(entropy_mesh % upper_right) if (allocated(entropy_mesh % width)) deallocate(entropy_mesh % width) deallocate(entropy_mesh) end if @@ -541,7 +541,7 @@ contains if (associated(ufs_mesh)) then if (allocated(ufs_mesh % lower_left)) deallocate(ufs_mesh % lower_left) if (allocated(ufs_mesh % upper_right)) & - deallocate(ufs_mesh % upper_right) + deallocate(ufs_mesh % upper_right) if (allocated(ufs_mesh % width)) deallocate(ufs_mesh % width) deallocate(ufs_mesh) end if diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index c3e5b44ba8..34b72196ee 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -7,7 +7,7 @@ module hdf5_interface use, intrinsic :: ISO_C_BINDING #ifdef MPI - use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL + use message_passing, only: MPI_COMM_WORLD, MPI_INFO_NULL #endif implicit none diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index 161246f1c9..d4953d076e 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -186,7 +186,7 @@ contains call su % write_data("lattice", "fill_type", & group="geometry/cells/cell " // trim(to_str(c % id))) call su % write_data(lattices(c % fill) % obj % id, "lattice", & - group="geometry/cells/cell " // trim(to_str(c % id))) + group="geometry/cells/cell " // trim(to_str(c % id))) end select ! Write list of bounding surfaces @@ -373,7 +373,7 @@ contains ! Write lattice universes. allocate(lattice_universes(lat % n_cells(1), lat % n_cells(2), & - &lat % n_cells(3))) + &lat % n_cells(3))) do j = 1, lat % n_cells(1) do k = 1, lat % n_cells(2) do m = 1, lat % n_cells(3) @@ -399,13 +399,13 @@ contains ! Write lattice center, pitch and outer universe. if (lat % is_3d) then - call su % write_data(lat % center, "center", length=3, & + call su % write_data(lat % center, "center", length=3, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) else - call su % write_data(lat % center, "center", length=2, & + call su % write_data(lat % center, "center", length=2, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) end if - + if (lat % is_3d) then call su % write_data(lat % pitch, "pitch", length=3, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) @@ -429,7 +429,7 @@ contains ! Write lattice universes. allocate(lattice_universes(2*lat % n_rings - 1, 2*lat % n_rings - 1, & - &lat % n_axial)) + &lat % n_axial)) do m = 1, lat % n_axial do k = 1, 2*lat % n_rings - 1 do j = 1, 2*lat % n_rings - 1 @@ -647,12 +647,12 @@ contains // "/filter " // trim(to_str(j))) case(FILTER_ENERGYIN) call su % write_data("energy", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) case(FILTER_ENERGYOUT) call su % write_data("energyout", "type_name", & - group="tallies/tally " // trim(to_str(t % id)) & - // "/filter " // trim(to_str(j))) + group="tallies/tally " // trim(to_str(t % id)) & + // "/filter " // trim(to_str(j))) end select end do FILTER_LOOP diff --git a/src/initialize.F90 b/src/initialize.F90 index 99a66e2ce8..409d531a38 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -916,9 +916,9 @@ contains thread_id = omp_get_thread_num() if (thread_id == 0) then - allocate(fission_bank(3*work)) + allocate(fission_bank(3*work)) else - allocate(fission_bank(3*work/n_threads)) + allocate(fission_bank(3*work/n_threads)) end if !$omp end parallel allocate(master_fission_bank(3*work), STAT=alloc_err) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9305313511..5c939a394f 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -472,7 +472,7 @@ contains ! Check for type of energy distribution type = '' if (check_for_node(node_dist, "type")) & - call get_node_value(node_dist, "type", type) + call get_node_value(node_dist, "type", type) select case (to_lower(type)) case ('monoenergetic') external_source % type_energy = SRC_ENERGY_MONO @@ -798,7 +798,7 @@ contains if (.not. source_separate) then do i = 1, n_source_points if (.not. statepoint_batch % contains(sourcepoint_batch % & - get_item(i))) then + get_item(i))) then call fatal_error('Sourcepoint batches are not a subset& & of statepoint batches.') end if @@ -811,7 +811,7 @@ contains call get_node_value(doc, "no_reduce", temp_str) temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - reduce_tallies = .false. + reduce_tallies = .false. end if ! Check if the user has specified to use confidence intervals for @@ -884,11 +884,11 @@ contains &// trim(to_str(i)) // " in settings.xml file!") end if call get_node_value(node_scatterer, "nuclide", & - nuclides_0K(i) % nuclide) + nuclides_0K(i) % nuclide) if (check_for_node(node_scatterer, "method")) then call get_node_value(node_scatterer, "method", & - nuclides_0K(i) % scheme) + nuclides_0K(i) % scheme) end if ! check to make sure xs name for which method is applied is given @@ -898,7 +898,7 @@ contains &// " given in cross_sections.xml") end if call get_node_value(node_scatterer, "xs_label", & - nuclides_0K(i) % name) + nuclides_0K(i) % name) ! check to make sure 0K xs name for which method is applied is given if (.not. check_for_node(node_scatterer, "xs_label_0K")) then @@ -906,11 +906,11 @@ contains &// trim(to_str(i)) // " given in cross_sections.xml") end if call get_node_value(node_scatterer, "xs_label_0K", & - nuclides_0K(i) % name_0K) + nuclides_0K(i) % name_0K) if (check_for_node(node_scatterer, "E_min")) then call get_node_value(node_scatterer, "E_min", & - nuclides_0K(i) % E_min) + nuclides_0K(i) % E_min) end if ! check that E_min is non-negative @@ -921,7 +921,7 @@ contains if (check_for_node(node_scatterer, "E_max")) then call get_node_value(node_scatterer, "E_max", & - nuclides_0K(i) % E_max) + nuclides_0K(i) % E_max) end if ! check that E_max is not less than E_min @@ -1081,7 +1081,7 @@ contains ! Read material word = '' if (check_for_node(node_cell, "material")) & - call get_node_value(node_cell, "material", word) + call get_node_value(node_cell, "material", word) select case(to_lower(word)) case ('void') c % material = MATERIAL_VOID @@ -1248,7 +1248,7 @@ contains ! Copy and interpret surface type word = '' if (check_for_node(node_surf, "type")) & - call get_node_value(node_surf, "type", word) + call get_node_value(node_surf, "type", word) select case(to_lower(word)) case ('x-plane') s % type = SURF_PX @@ -1306,7 +1306,7 @@ contains ! Boundary conditions word = '' if (check_for_node(node_surf, "boundary")) & - call get_node_value(node_surf, "boundary", word) + call get_node_value(node_surf, "boundary", word) select case (to_lower(word)) case ('transmission', 'transmit', '') s % bc = BC_TRANSMIT @@ -1447,7 +1447,7 @@ contains do k = 0, n_y - 1 do j = 1, n_x lat % universes(j, n_y - k, m) = & - &temp_int_array(j + n_x*k + n_x*n_y*(m-1)) + &temp_int_array(j + n_x*k + n_x*n_y*(m-1)) end do end do end do @@ -1713,7 +1713,7 @@ contains ! Copy default cross section if present if (check_for_node(doc, "default_xs")) & - call get_node_value(doc, "default_xs", default_xs) + call get_node_value(doc, "default_xs", default_xs) ! Get pointer to list of XML call get_node_list(doc, "material", node_mat_list) @@ -1844,7 +1844,7 @@ contains ! store full name call get_node_value(node_nuc, "name", temp_str) if (check_for_node(node_nuc, "xs")) & - call get_node_value(node_nuc, "xs", name) + call get_node_value(node_nuc, "xs", name) name = trim(temp_str) // "." // trim(name) ! save name and density to list @@ -1853,7 +1853,7 @@ contains ! Check if no atom/weight percents were specified or if both atom and ! weight percents were specified if (.not.check_for_node(node_nuc, "ao") .and. & - .not.check_for_node(node_nuc, "wo")) then + .not.check_for_node(node_nuc, "wo")) then call fatal_error("No atom or weight percent specified for nuclide " & &// trim(name)) elseif (check_for_node(node_nuc, "ao") .and. & @@ -1903,7 +1903,7 @@ contains ! Check if no atom/weight percents were specified or if both atom and ! weight percents were specified if (.not.check_for_node(node_ele, "ao") .and. & - .not.check_for_node(node_ele, "wo")) then + .not.check_for_node(node_ele, "wo")) then call fatal_error("No atom or weight percent specified for element " & &// trim(name)) elseif (check_for_node(node_ele, "ao") .and. & @@ -2010,7 +2010,7 @@ contains ! Determine name of S(a,b) table if (.not.check_for_node(node_sab, "name") .or. & - .not.check_for_node(node_sab, "xs")) then + .not.check_for_node(node_sab, "xs")) then call fatal_error("Need to specify and for S(a,b) & &table.") end if @@ -2157,7 +2157,7 @@ contains call get_node_value(doc, "assume_separate", temp_str) temp_str = to_lower(temp_str) if (trim(temp_str) == 'true' .or. trim(temp_str) == '1') & - assume_separate = .true. + assume_separate = .true. end if ! ========================================================================== @@ -2185,7 +2185,7 @@ contains ! Read mesh type temp_str = '' if (check_for_node(node_mesh, "type")) & - call get_node_value(node_mesh, "type", temp_str) + call get_node_value(node_mesh, "type", temp_str) select case (to_lower(temp_str)) case ('rect', 'rectangle', 'rectangular') m % type = LATTICE_RECT @@ -2227,14 +2227,14 @@ contains ! Make sure both upper-right or width were specified if (check_for_node(node_mesh, "upper_right") .and. & - check_for_node(node_mesh, "width")) then + check_for_node(node_mesh, "width")) then call fatal_error("Cannot specify both and on a & &tally mesh.") end if ! Make sure either upper-right or width was specified if (.not.check_for_node(node_mesh, "upper_right") .and. & - .not.check_for_node(node_mesh, "width")) then + .not.check_for_node(node_mesh, "width")) then call fatal_error("Must specify either and on a & &tally mesh.") end if @@ -2242,7 +2242,7 @@ contains if (check_for_node(node_mesh, "width")) then ! Check to ensure width has same dimensions if (get_arraysize_double(node_mesh, "width") /= & - get_arraysize_double(node_mesh, "lower_left")) then + get_arraysize_double(node_mesh, "lower_left")) then call fatal_error("Number of entries on must be the same as & &the number of entries on .") end if @@ -2260,7 +2260,7 @@ contains elseif (check_for_node(node_mesh, "upper_right")) then ! Check to ensure width has same dimensions if (get_arraysize_double(node_mesh, "upper_right") /= & - get_arraysize_double(node_mesh, "lower_left")) then + get_arraysize_double(node_mesh, "lower_left")) then call fatal_error("Number of entries on must be the & &same as the number of entries on .") end if @@ -2323,7 +2323,7 @@ contains ! Copy tally name if (check_for_node(node_tal, "name")) & - call get_node_value(node_tal, "name", t % name) + call get_node_value(node_tal, "name", t % name) ! ======================================================================= ! READ DATA FOR FILTERS @@ -2345,13 +2345,13 @@ contains ! Convert filter type to lower case temp_str = '' if (check_for_node(node_filt, "type")) & - call get_node_value(node_filt, "type", temp_str) + call get_node_value(node_filt, "type", temp_str) temp_str = to_lower(temp_str) ! Determine number of bins if (check_for_node(node_filt, "bins")) then if (trim(temp_str) == 'energy' .or. & - trim(temp_str) == 'energyout') then + trim(temp_str) == 'energyout') then n_words = get_arraysize_double(node_filt, "bins") else n_words = get_arraysize_integer(node_filt, "bins") @@ -2618,7 +2618,7 @@ contains if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then n_order_pos = scan(score_name,'0123456789') n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) + score_name(n_order_pos:(len_trim(score_name)))),4) if (n_order > MAX_ANG_ORDER) then ! User requested too many orders; throw a warning and set to the ! maximum order. @@ -2661,7 +2661,7 @@ contains if (starts_with(score_name,trim(MOMENT_STRS(imomstr)))) then n_order_pos = scan(score_name,'0123456789') n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) + score_name(n_order_pos:(len_trim(score_name)))),4) if (n_order > MAX_ANG_ORDER) then ! User requested too many orders; throw a warning and set to the ! maximum order. @@ -2684,7 +2684,7 @@ contains if (starts_with(score_name,trim(MOMENT_N_STRS(imomstr)))) then n_order_pos = scan(score_name,'0123456789') n_order = int(str_to_int( & - score_name(n_order_pos:(len_trim(score_name)))),4) + score_name(n_order_pos:(len_trim(score_name)))),4) if (n_order > MAX_ANG_ORDER) then ! User requested too many orders; throw a warning and set to the ! maximum order. @@ -3006,7 +3006,7 @@ contains ! Get the trigger type - "variance", "std_dev" or "rel_err" if (check_for_node(node_trigger, "type")) then call get_node_value(node_trigger, "type", temp_str) - temp_str = to_lower(temp_str) + temp_str = to_lower(temp_str) else call fatal_error("Must specify trigger type for tally " // & trim(to_str(t % id)) // " in tally XML file.") @@ -3105,7 +3105,7 @@ contains ! Increment the overall trigger index trig_ind = trig_ind + 1 - end if + end if end do SCORE_LOOP ! Deallocate the list of tally scores used to create triggers @@ -3223,7 +3223,7 @@ contains ! Copy plot type temp_str = 'slice' if (check_for_node(node_plot, "type")) & - call get_node_value(node_plot, "type", temp_str) + call get_node_value(node_plot, "type", temp_str) temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("slice") @@ -3238,7 +3238,7 @@ contains ! Set output file path filename = trim(to_str(pl % id)) // "_plot" if (check_for_node(node_plot, "filename")) & - call get_node_value(node_plot, "filename", filename) + call get_node_value(node_plot, "filename", filename) select case (pl % type) case (PLOT_TYPE_SLICE) pl % path_plot = trim(path_input) // trim(filename) // ".ppm" @@ -3283,7 +3283,7 @@ contains if (pl % type == PLOT_TYPE_SLICE) then temp_str = 'xy' if (check_for_node(node_plot, "basis")) & - call get_node_value(node_plot, "basis", temp_str) + call get_node_value(node_plot, "basis", temp_str) temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("xy") @@ -3338,7 +3338,7 @@ contains ! Copy plot color type and initialize all colors randomly temp_str = "cell" if (check_for_node(node_plot, "color")) & - call get_node_value(node_plot, "color", temp_str) + call get_node_value(node_plot, "color", temp_str) temp_str = to_lower(temp_str) select case (trim(temp_str)) case ("cell") @@ -3452,7 +3452,7 @@ contains ! Ensure that there is a linewidth for this meshlines specification if (check_for_node(node_meshlines, "linewidth")) then call get_node_value(node_meshlines, "linewidth", & - pl % meshlines_width) + pl % meshlines_width) else call fatal_error("Must specify a linewidth for meshlines & &specification in plot " // trim(to_str(pl % id))) @@ -3468,7 +3468,7 @@ contains end if call get_node_array(node_meshlines, "color", & - pl % meshlines_color % rgb) + pl % meshlines_color % rgb) else pl % meshlines_color % rgb = (/ 0, 0, 0 /) @@ -3494,8 +3494,8 @@ contains end if i_mesh = cmfd_tallies(1) % & - filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % & - int_bins(1) + filters(cmfd_tallies(1) % find_filter(FILTER_MESH)) % & + int_bins(1) pl % meshlines_mesh => meshes(i_mesh) case ('entropy') @@ -3654,9 +3654,9 @@ contains ! Check if cross_sections.xml exists inquire(FILE=path_cross_sections, EXIST=file_exists) if (.not. file_exists) then - ! Could not find cross_sections.xml file - call fatal_error("Cross sections XML file '" & - &// trim(path_cross_sections) // "' does not exist!") + ! Could not find cross_sections.xml file + call fatal_error("Cross sections XML file '" & + &// trim(path_cross_sections) // "' does not exist!") end if call write_message("Reading cross sections XML file...", 5) @@ -3665,28 +3665,28 @@ contains call open_xmldoc(doc, path_cross_sections) if (check_for_node(doc, "directory")) then - ! Copy directory information if present - call get_node_value(doc, "directory", directory) + ! Copy directory information if present + call get_node_value(doc, "directory", directory) else - ! If no directory is listed in cross_sections.xml, by default select the - ! directory in which the cross_sections.xml file resides - i = index(path_cross_sections, "/", BACK=.true.) - directory = path_cross_sections(1:i) + ! If no directory is listed in cross_sections.xml, by default select the + ! directory in which the cross_sections.xml file resides + i = index(path_cross_sections, "/", BACK=.true.) + directory = path_cross_sections(1:i) end if ! determine whether binary/ascii temp_str = '' if (check_for_node(doc, "filetype")) & - call get_node_value(doc, "filetype", temp_str) + call get_node_value(doc, "filetype", temp_str) if (trim(temp_str) == 'ascii') then - filetype = ASCII + filetype = ASCII elseif (trim(temp_str) == 'binary') then - filetype = BINARY + filetype = BINARY elseif (len_trim(temp_str) == 0) then - filetype = ASCII + filetype = ASCII else - call fatal_error("Unknown filetype in cross_sections.xml: " & - &// trim(temp_str)) + call fatal_error("Unknown filetype in cross_sections.xml: " & + &// trim(temp_str)) end if ! copy default record length and entries for binary files @@ -3701,83 +3701,83 @@ contains ! Allocate xs_listings array if (n_listings == 0) then - call fatal_error("No ACE table listings present in cross_sections.xml & - &file!") + call fatal_error("No ACE table listings present in cross_sections.xml & + &file!") else - allocate(xs_listings(n_listings)) + allocate(xs_listings(n_listings)) end if do i = 1, n_listings - listing => xs_listings(i) + listing => xs_listings(i) - ! Get pointer to ace table XML node - call get_list_item(node_ace_list, i, node_ace) + ! Get pointer to ace table XML node + call get_list_item(node_ace_list, i, node_ace) - ! copy a number of attributes - call get_node_value(node_ace, "name", listing % name) - if (check_for_node(node_ace, "alias")) & - call get_node_value(node_ace, "alias", listing % alias) - call get_node_value(node_ace, "zaid", listing % zaid) - call get_node_value(node_ace, "awr", listing % awr) - if (check_for_node(node_ace, "temperature")) & - call get_node_value(node_ace, "temperature", listing % kT) - call get_node_value(node_ace, "location", listing % location) + ! copy a number of attributes + call get_node_value(node_ace, "name", listing % name) + if (check_for_node(node_ace, "alias")) & + call get_node_value(node_ace, "alias", listing % alias) + call get_node_value(node_ace, "zaid", listing % zaid) + call get_node_value(node_ace, "awr", listing % awr) + if (check_for_node(node_ace, "temperature")) & + call get_node_value(node_ace, "temperature", listing % kT) + call get_node_value(node_ace, "location", listing % location) - ! determine type of cross section - if (ends_with(listing % name, 'c')) then - listing % type = ACE_NEUTRON - elseif (ends_with(listing % name, 't')) then - listing % type = ACE_THERMAL - end if + ! determine type of cross section + if (ends_with(listing % name, 'c')) then + listing % type = ACE_NEUTRON + elseif (ends_with(listing % name, 't')) then + listing % type = ACE_THERMAL + end if - ! set filetype, record length, and number of entries - if (check_for_node(node_ace, "filetype")) then - temp_str = '' - call get_node_value(node_ace, "filetype", temp_str) - if (temp_str == 'ascii') then - listing % filetype = ASCII - else if (temp_str == 'binary') then - listing % filetype = BINARY - end if - else - listing % filetype = filetype - end if + ! set filetype, record length, and number of entries + if (check_for_node(node_ace, "filetype")) then + temp_str = '' + call get_node_value(node_ace, "filetype", temp_str) + if (temp_str == 'ascii') then + listing % filetype = ASCII + else if (temp_str == 'binary') then + listing % filetype = BINARY + end if + else + listing % filetype = filetype + end if - ! Set record length and entries for binary files - if (filetype == BINARY) then - listing % recl = recl - listing % entries = entries - end if + ! Set record length and entries for binary files + if (filetype == BINARY) then + listing % recl = recl + listing % entries = entries + end if - ! determine metastable state - if (.not.check_for_node(node_ace, "metastable")) then - listing % metastable = .false. - else - listing % metastable = .true. - end if + ! determine metastable state + if (.not.check_for_node(node_ace, "metastable")) then + listing % metastable = .false. + else + listing % metastable = .true. + end if - ! determine path of cross section table - if (check_for_node(node_ace, "path")) then - call get_node_value(node_ace, "path", temp_str) - else - call fatal_error("Path missing for isotope " // listing % name) - end if + ! determine path of cross section table + if (check_for_node(node_ace, "path")) then + call get_node_value(node_ace, "path", temp_str) + else + call fatal_error("Path missing for isotope " // listing % name) + end if - if (starts_with(temp_str, '/')) then - listing % path = trim(temp_str) - else - if (ends_with(directory,'/')) then - listing % path = trim(directory) // trim(temp_str) - else - listing % path = trim(directory) // '/' // trim(temp_str) - end if - end if + if (starts_with(temp_str, '/')) then + listing % path = trim(temp_str) + else + if (ends_with(directory,'/')) then + listing % path = trim(directory) // trim(temp_str) + else + listing % path = trim(directory) // '/' // trim(temp_str) + end if + end if - ! create dictionary entry for both name and alias - call xs_listing_dict % add_key(to_lower(listing % name), i) - if (check_for_node(node_ace, "alias")) then - call xs_listing_dict % add_key(to_lower(listing % alias), i) - end if + ! create dictionary entry for both name and alias + call xs_listing_dict % add_key(to_lower(listing % name), i) + if (check_for_node(node_ace, "alias")) then + call xs_listing_dict % add_key(to_lower(listing % alias), i) + end if end do ! Check that 0K nuclides are listed in the cross_sections.xml file diff --git a/src/list_header.F90 b/src/list_header.F90 index a9311ec64a..8020b96b17 100644 --- a/src/list_header.F90 +++ b/src/list_header.F90 @@ -45,7 +45,7 @@ module list_header private integer :: count = 0 ! Number of elements in list - ! Used in get_item for fast sequential lookups + ! Used in get_item for fast sequential lookups integer :: last_index = huge(0) type(ListElemInt), pointer :: last_elem => null() @@ -67,7 +67,7 @@ module list_header private integer :: count = 0 ! Number of elements in list - ! Used in get_item for fast sequential lookups + ! Used in get_item for fast sequential lookups integer :: last_index = huge(0) type(ListElemReal), pointer :: last_elem => null() @@ -89,7 +89,7 @@ module list_header private integer :: count = 0 ! Number of elements in list - ! Used in get_item for fast sequential lookups + ! Used in get_item for fast sequential lookups integer :: last_index = huge(0) type(ListElemChar), pointer :: last_elem => null() @@ -198,7 +198,7 @@ contains type(ListElemInt), pointer :: current => null() type(ListElemInt), pointer :: next => null() - + if (this % count > 0) then current => this % head do while (associated(current)) @@ -224,7 +224,7 @@ contains type(ListElemReal), pointer :: current => null() type(ListElemReal), pointer :: next => null() - + if (this % count > 0) then current => this % head do while (associated(current)) @@ -250,7 +250,7 @@ contains type(ListElemChar), pointer :: current => null() type(ListElemChar), pointer :: next => null() - + if (this % count > 0) then current => this % head do while (associated(current)) @@ -518,7 +518,7 @@ contains ! Allocate new element allocate(new_elem) new_elem % data = data - + ! Put it before the i-th element new_elem % prev => elem % prev new_elem % next => elem @@ -573,7 +573,7 @@ contains ! Allocate new element allocate(new_elem) new_elem % data = data - + ! Put it before the i-th element new_elem % prev => elem % prev new_elem % next => elem @@ -628,7 +628,7 @@ contains ! Allocate new element allocate(new_elem) new_elem % data = data - + ! Put it before the i-th element new_elem % prev => elem % prev new_elem % next => elem diff --git a/src/math.F90 b/src/math.F90 index 58406bda18..826172fa24 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -143,20 +143,20 @@ contains pnx = 7.875_8 * (x ** 5) - 8.75_8 * x * x * x + 1.875 * x case(6) pnx = 14.4375_8 * (x ** 6) - 19.6875_8 * (x ** 4) + & - 6.5625_8 * x * x - 0.3125_8 + 6.5625_8 * x * x - 0.3125_8 case(7) pnx = 26.8125_8 * (x ** 7) - 43.3125_8 * (x ** 5) + & - 19.6875_8 * x * x * x - 2.1875_8 * x + 19.6875_8 * x * x * x - 2.1875_8 * x case(8) pnx = 50.2734375_8 * (x ** 8) - 93.84375_8 * (x ** 6) + & - 54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8 + 54.140625 * (x ** 4) - 9.84375_8 * x * x + 0.2734375_8 case(9) pnx = 94.9609375_8 * (x ** 9) - 201.09375_8 * (x ** 7) + & - 140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x + 140.765625_8 * (x ** 5) - 36.09375_8 * x * x * x + 2.4609375_8 * x case(10) pnx = 180.42578125_8 * (x ** 10) - 427.32421875_8 * (x ** 8) + & - 351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + & - 13.53515625_8 * x * x - 0.24609375_8 + 351.9140625_8 * (x ** 6) - 117.3046875_8 * (x ** 4) + & + 13.53515625_8 * x * x - 0.24609375_8 case default pnx = ONE ! correct for case(0), incorrect for the rest end select @@ -215,12 +215,12 @@ contains rn(2) = 1.93649167310371_8 * w*(w2m1) * sin(TWO*phi) ! l = 3, m = -1 rn(3) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * & - sin(phi) + sin(phi) ! l = 3, m = 0 rn(4) = 2.5_8 * w**3 - 1.5_8 * w ! l = 3, m = 1 rn(5) = 0.408248290463863_8*sqrt(w2m1)*((15.0_8/TWO)*w**2 - THREE/TWO) * & - cos(phi) + cos(phi) ! l = 3, m = 2 rn(6) = 1.93649167310371_8 * w*(w2m1) * cos(TWO*phi) ! l = 3, m = 3 @@ -232,18 +232,18 @@ contains rn(2) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * sin(THREE* phi) ! l = 4, m = -2 rn(3) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & - sin(TWO*phi) + sin(TWO*phi) ! l = 4, m = -1 - rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * & - sin(phi) + rn(4) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)& + * sin(phi) ! l = 4, m = 0 rn(5) = 4.375_8 * w**4 - 3.75_8 * w**2 + 0.375_8 ! l = 4, m = 1 - rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w) * & - cos(phi) + rn(6) = 0.316227766016838_8*sqrt(w2m1)*((35.0_8/TWO)*w**3 - 15.0_8/TWO*w)& + * cos(phi) ! l = 4, m = 2 rn(7) = 0.074535599249993_8 * (w2m1)*((105.0_8/TWO)*w**2 - 15.0_8/TWO) * & - cos(TWO*phi) + cos(TWO*phi) ! l = 4, m = 3 rn(8) = 2.09165006633519_8 * w*(w2m1)**(THREE/TWO) * cos(THREE* phi) ! l = 4, m = 4 @@ -255,24 +255,26 @@ contains rn(2) = 2.21852991866236_8 * w*(w2m1)**2 * sin(4.0_8*phi) ! l = 5, m = -3 rn(3) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* & - ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi) + ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * sin(THREE*phi) ! l = 5, m = -2 - rn(4) = 0.0487950036474267_8 * (w2m1)*((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * & - sin(TWO*phi) + rn(4) = 0.0487950036474267_8 * (w2m1) & + * ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * sin(TWO*phi) ! l = 5, m = -1 rn(5) = 0.258198889747161_8*sqrt(w2m1)* & - ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * sin(phi) + ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) & + * sin(phi) ! l = 5, m = 0 rn(6) = 7.875_8 * w**5 - 8.75_8 * w**3 + 1.875_8 * w ! l = 5, m = 1 rn(7) = 0.258198889747161_8*sqrt(w2m1)* & - ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) * cos(phi) + ((315.0_8/8.0_8)*w**4 - 105.0_8/4.0_8 * w**2 + 15.0_8/8.0_8) & + * cos(phi) ! l = 5, m = 2 rn(8) = 0.0487950036474267_8 * (w2m1)* & - ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi) + ((315.0_8/TWO)*w**3 - 105.0_8/TWO*w) * cos(TWO*phi) ! l = 5, m = 3 rn(9) = 0.00996023841111995_8 * (w2m1)**(THREE/TWO)* & - ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi) + ((945.0_8 /TWO)*w**2 - 105.0_8/TWO) * cos(THREE*phi) ! l = 5, m = 4 rn(10) = 2.21852991866236_8 * w*(w2m1)**2 * cos(4.0_8*phi) ! l = 5, m = 5 @@ -284,30 +286,34 @@ contains rn(2) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * sin(5.0_8*phi) ! l = 6, m = -4 rn(3) = 0.00104990131391452_8 * (w2m1)**2 * & - ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi) + ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * sin(4.0_8*phi) ! l = 6, m = -3 rn(4) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * & - ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi) + ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * sin(THREE*phi) ! l = 6, m = -2 rn(5) = 0.0345032779671177_8 * (w2m1) * & - ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * sin(TWO*phi) + ((3465.0_8/8.0_8)*w**4 - 945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) & + * sin(TWO*phi) ! l = 6, m = -1 rn(6) = 0.218217890235992_8*sqrt(w2m1) * & - ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * sin(phi) + ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) & + * sin(phi) ! l = 6, m = 0 rn(7) = 14.4375_8 * w**6 - 19.6875_8 * w**4 + 6.5625_8 * w**2 - 0.3125_8 ! l = 6, m = 1 rn(8) = 0.218217890235992_8*sqrt(w2m1) * & - ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) * cos(phi) + ((693.0_8/8.0_8)*w**5- 315.0_8/4.0_8 * w**3 + (105.0_8/8.0_8)*w) & + * cos(phi) ! l = 6, m = 2 rn(9) = 0.0345032779671177_8 * (w2m1) * & - ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) * cos(TWO*phi) + ((3465.0_8/8.0_8)*w**4 -945.0_8/4.0_8 * w**2 + 105.0_8/8.0_8) & + * cos(TWO*phi) ! l = 6, m = 3 rn(10) = 0.00575054632785295_8 * (w2m1)**(THREE/TWO) * & - ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi) + ((3465.0_8/TWO)*w**3 - 945.0_8/TWO*w) * cos(THREE*phi) ! l = 6, m = 4 rn(11) = 0.00104990131391452_8 * (w2m1)**2 * & - ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi) + ((10395.0_8/TWO)*w**2 - 945.0_8/TWO) * cos(4.0_8*phi) ! l = 6, m = 5 rn(12) = 2.32681380862329_8 * w*(w2m1)**(5.0_8/TWO) * cos(5.0_8*phi) ! l = 6, m = 6 @@ -319,42 +325,43 @@ contains rn(2) = 2.42182459624969_8 * w*(w2m1)**3 * sin(6.0_8*phi) ! l = 7, m = -5 rn(3) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & - ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi) + ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * sin(5.0_8*phi) ! l = 7, m = -4 rn(4) = 0.000548293079133141_8 * (w2m1)**2* & - ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi) + ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * sin(4.0_8*phi) ! l = 7, m = -3 rn(5) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* & - ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & - sin(THREE*phi) + ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & + sin(THREE*phi) ! l = 7, m = -2 rn(6) = 0.025717224993682_8 * (w2m1)* & - ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & - sin(TWO*phi) + ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & + sin(TWO*phi) ! l = 7, m = -1 rn(7) = 0.188982236504614_8*sqrt(w2m1)* & - ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & - (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi) + ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & + (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * sin(phi) ! l = 7, m = 0 - rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 * w + rn(8) = 26.8125_8 * w**7 - 43.3125_8 * w**5 + 19.6875_8 * w**3 -2.1875_8 & + * w ! l = 7, m = 1 rn(9) = 0.188982236504614_8*sqrt(w2m1)* & - ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & - (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi) + ((3003.0_8/16.0_8)*w**6 - 3465.0_8/16.0_8 * w**4 + & + (945.0_8/16.0_8)*w**2 - 35.0_8/16.0_8) * cos(phi) ! l = 7, m = 2 rn(10) = 0.025717224993682_8 * (w2m1)* & - ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & - cos(TWO*phi) + ((9009.0_8/8.0_8)*w**5 -3465.0_8/4.0_8 * w**3 + (945.0_8/8.0_8)*w)* & + cos(TWO*phi) ! l = 7, m = 3 rn(11) = 0.00363696483726654_8 * (w2m1)**(THREE/TWO)* & - ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & - cos(THREE*phi) + ((45045.0_8/8.0_8)*w**4 - 10395.0_8/4.0_8 * w**2 + 945.0_8/8.0_8)* & + cos(THREE*phi) ! l = 7, m = 4 rn(12) = 0.000548293079133141_8 * (w2m1)**2 * & - ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi) + ((45045.0_8/TWO)*w**3 - 10395.0_8/TWO*w) * cos(4.0_8*phi) ! l = 7, m = 5 rn(13) = 9.13821798555235d-5*(w2m1)**(5.0_8/TWO)* & - ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi) + ((135135.0_8/TWO)*w**2 - 10395.0_8/TWO) * cos(5.0_8*phi) ! l = 7, m = 6 rn(14) = 2.42182459624969_8 * w*(w2m1)**3 * cos(6.0_8*phi) ! l = 7, m = 7 @@ -366,48 +373,50 @@ contains rn(2) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * sin(7.0_8*phi) ! l = 8, m = -6 rn(3) = 6.77369783729086d-6*(w2m1)**3* & - ((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi) + ((2027025.0_8/TWO)*w**2 - 135135.0_8/TWO) * sin(6.0_8*phi) ! l = 8, m = -5 rn(4) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)* & - ((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi) + ((675675.0_8/TWO)*w**3 - 135135.0_8/TWO*w) * sin(5.0_8*phi) ! l = 8, m = -4 rn(5) = 0.000316557156832328_8 * (w2m1)**2* & - ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * sin(4.0_8*phi) + ((675675.0_8/8.0_8)*w**4 - 135135.0_8/4.0_8 * w**2 & + + 10395.0_8/8.0_8) * sin(4.0_8*phi) ! l = 8, m = -3 rn(6) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* & - ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + (10395.0_8/8.0_8)*w) * sin(THREE*phi) + ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 & + + (10395.0_8/8.0_8)*w) * sin(THREE*phi) ! l = 8, m = -2 rn(7) = 0.0199204768222399_8 * (w2m1)* & - ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + & - (10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi) + ((45045.0_8/16.0_8)*w**6- 45045.0_8/16.0_8 * w**4 + & + (10395.0_8/16.0_8)*w**2 - 315.0_8/16.0_8) * sin(TWO*phi) ! l = 8, m = -1 rn(8) = 0.166666666666667_8*sqrt(w2m1)* & - ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & - (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi) + ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & + (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * sin(phi) ! l = 8, m = 0 rn(9) = 50.2734375_8 * w**8 - 93.84375_8 * w**6 + 54.140625_8 * w**4 -& - 9.84375_8 * w**2 + 0.2734375_8 + 9.84375_8 * w**2 + 0.2734375_8 ! l = 8, m = 1 rn(10) = 0.166666666666667_8*sqrt(w2m1)* & - ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & - (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi) + ((6435.0_8/16.0_8)*w**7 - 9009.0_8/16.0_8 * w**5 + & + (3465.0_8/16.0_8)*w**3 - 315.0_8/16.0_8 * w) * cos(phi) ! l = 8, m = 2 rn(11) = 0.0199204768222399_8 * (w2m1)*((45045.0_8/16.0_8)*w**6- & - 45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - & - 315.0_8/16.0_8) * cos(TWO*phi) + 45045.0_8/16.0_8 * w**4 + (10395.0_8/16.0_8)*w**2 - & + 315.0_8/16.0_8) * cos(TWO*phi) ! l = 8, m = 3 rn(12) = 0.00245204119306875_8 * (w2m1)**(THREE/TWO)* & - ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + & - (10395.0_8/8.0_8)*w) * cos(THREE*phi) + ((135135.0_8/8.0_8)*w**5 - 45045.0_8/4.0_8 * w**3 + & + (10395.0_8/8.0_8)*w) * cos(THREE*phi) ! l = 8, m = 4 rn(13) = 0.000316557156832328_8 * (w2m1)**2*((675675.0_8/8.0_8)*w**4 - & - 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi) + 135135.0_8/4.0_8 * w**2 + 10395.0_8/8.0_8) * cos(4.0_8*phi) ! l = 8, m = 5 - rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 - & - 135135.0_8/TWO*w) * cos(5.0_8*phi) + rn(14) = 4.38985792528482d-5*(w2m1)**(5.0_8/TWO)*((675675.0_8/TWO)*w**3 -& + 135135.0_8/TWO*w) * cos(5.0_8*phi) ! l = 8, m = 6 rn(15) = 6.77369783729086d-6*(w2m1)**3*((2027025.0_8/TWO)*w**2 - & - 135135.0_8/TWO) * cos(6.0_8*phi) + 135135.0_8/TWO) * cos(6.0_8*phi) ! l = 8, m = 7 rn(16) = 2.50682661696018_8 * w*(w2m1)**(7.0_8/TWO) * cos(7.0_8*phi) ! l = 8, m = 8 @@ -419,54 +428,56 @@ contains rn(2) = 2.58397773170915_8 * w*(w2m1)**4 * sin(8.0_8*phi) ! l = 9, m = -7 rn(3) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & - ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi) + ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * sin(7.0_8*phi) ! l = 9, m = -6 rn(4) = 3.02928976464514d-6*(w2m1)**3* & - ((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi) + ((11486475.0_8/TWO)*w**3 - 2027025.0_8/TWO*w) * sin(6.0_8*phi) ! l = 9, m = -5 rn(5) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)* & - ((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + & - 135135.0_8/8.0_8) * sin(5.0_8*phi) + ((11486475.0_8/8.0_8)*w**4 - 2027025.0_8/4.0_8 * w**2 + & + 135135.0_8/8.0_8) * sin(5.0_8*phi) ! l = 9, m = -4 rn(6) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & - 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi) + 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * sin(4.0_8*phi) ! l = 9, m = -3 rn(7) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)* & - ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + & - (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi) + ((765765.0_8/16.0_8)*w**6 - 675675.0_8/16.0_8 * w**4 + & + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8) * sin(THREE*phi) ! l = 9, m = -2 rn(8) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7- & - 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/16.0_8 * w)* & - sin(TWO*phi) + 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 & + - 3465.0_8/16.0_8 * w) * sin(TWO*phi) ! l = 9, m = -1 rn(9) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & - 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * sin(phi) + 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 - 3465.0_8/32.0_8 & + * w**2 + 315.0_8/128.0_8) * sin(phi) ! l = 9, m = 0 rn(10) = 94.9609375_8 * w**9 - 201.09375_8 * w**7 + 140.765625_8 * w**5- & - 36.09375_8 * w**3 + 2.4609375_8 * w + 36.09375_8 * w**3 + 2.4609375_8 * w ! l = 9, m = 1 rn(11) = 0.149071198499986_8*sqrt(w2m1)*((109395.0_8/128.0_8)*w**8 - & - 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 * w**2 + 315.0_8/128.0_8) * cos(phi) + 45045.0_8/32.0_8 * w**6 + (45045.0_8/64.0_8)*w**4 -3465.0_8/32.0_8 & + * w**2 + 315.0_8/128.0_8) * cos(phi) ! l = 9, m = 2 rn(12) = 0.0158910431540932_8 * (w2m1)*((109395.0_8/16.0_8)*w**7 - & - 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 - 3465.0_8/ 16.0_8 * w) * & - cos(TWO*phi) + 135135.0_8/16.0_8 * w**5 + (45045.0_8/16.0_8)*w**3 & + - 3465.0_8/ 16.0_8 * w) * cos(TWO*phi) ! l = 9, m = 3 - rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)*w**6 - & - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 - 3465.0_8/16.0_8)* & - cos(THREE*phi) + rn(13) = 0.00173385495536766_8 * (w2m1)**(THREE/TWO)*((765765.0_8/16.0_8)& + *w**6 - 675675.0_8/16.0_8 * w**4 + (135135.0_8/16.0_8)*w**2 & + - 3465.0_8/16.0_8)* cos(THREE*phi) ! l = 9, m = 4 rn(14) = 0.000196320414650061_8 * (w2m1)**2*((2297295.0_8/8.0_8)*w**5 - & - 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi) + 675675.0_8/4.0_8 * w**3 + (135135.0_8/8.0_8)*w) * cos(4.0_8*phi) ! l = 9, m = 5 rn(15) = 2.34647776186144d-5*(w2m1)**(5.0_8/TWO)*((11486475.0_8/8.0_8)* & - w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi) + w**4 - 2027025.0_8/4.0_8 * w**2 + 135135.0_8/8.0_8) * cos(5.0_8*phi) ! l = 9, m = 6 rn(16) = 3.02928976464514d-6*(w2m1)**3*((11486475.0_8/TWO)*w**3 - & - 2027025.0_8/TWO*w) * cos(6.0_8*phi) + 2027025.0_8/TWO*w) * cos(6.0_8*phi) ! l = 9, m = 7 rn(17) = 4.37240315267812d-7*(w2m1)**(7.0_8/TWO)* & - ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi) + ((34459425.0_8/TWO)*w**2 - 2027025.0_8/TWO) * cos(7.0_8*phi) ! l = 9, m = 8 rn(18) = 2.58397773170915_8 * w*(w2m1)**4 * cos(8.0_8*phi) ! l = 9, m = 9 @@ -478,65 +489,65 @@ contains rn(2) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * sin(9.0_8*phi) ! l = 10, m = -8 rn(3) = 2.49953651452314d-8*(w2m1)**4*((654729075.0_8/TWO)*w**2 - & - 34459425.0_8/TWO) * sin(8.0_8*phi) + 34459425.0_8/TWO) * sin(8.0_8*phi) ! l = 10, m = -7 rn(4) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & - ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi) + ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * sin(7.0_8*phi) ! l = 10, m = -6 rn(5) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & - 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi) + 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * sin(6.0_8*phi) ! l = 10, m = -5 rn(6) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & - ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & - (2027025.0_8/8.0_8)*w) * sin(5.0_8*phi) + ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & + (2027025.0_8/8.0_8)*w) * sin(5.0_8*phi) ! l = 10, m = -4 rn(7) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - & - 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & - 45045.0_8/16.0_8) * sin(4.0_8*phi) + 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & + 45045.0_8/16.0_8) * sin(4.0_8*phi) ! l = 10, m = -3 rn(8) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* & - ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & - (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi) + ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & + (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * sin(THREE*phi) ! l = 10, m = -2 rn(9) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & - 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - & - 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi) + 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 - & + 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * sin(TWO*phi) ! l = 10, m = -1 rn(10) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & - 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - & - 15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi) + 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 - & + 15015.0_8/32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * sin(phi) ! l = 10, m = 0 - rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 * w**6 - & - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8 + rn(11) = 180.42578125_8 * w**10 - 427.32421875_8 * w**8 +351.9140625_8 & + * w**6 - 117.3046875_8 * w**4 + 13.53515625_8 * w**2 -0.24609375_8 ! l = 10, m = 1 rn(12) = 0.134839972492648_8*sqrt(w2m1)*((230945.0_8/128.0_8)*w**9 - & - 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ & - 32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi) + 109395.0_8/32.0_8 * w**7 + (135135.0_8/64.0_8)*w**5 -15015.0_8/ & + 32.0_8 * w**3 + (3465.0_8/128.0_8)*w) * cos(phi) ! l = 10, m = 2 rn(13) = 0.012974982402692_8 * (w2m1)*((2078505.0_8/128.0_8)*w**8 - & - 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -& - 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi) + 765765.0_8/32.0_8 * w**6 + (675675.0_8/64.0_8)*w**4 -& + 45045.0_8/32.0_8 * w**2 + 3465.0_8/128.0_8) * cos(TWO*phi) ! l = 10, m = 3 rn(14) = 0.00127230170115096_8 * (w2m1)**(THREE/TWO)* & - ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & - (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi) + ((2078505.0_8/16.0_8)*w**7 - 2297295.0_8/16.0_8 * w**5 + & + (675675.0_8/16.0_8)*w**3 - 45045.0_8/16.0_8 * w) * cos(THREE*phi) ! l = 10, m = 4 - rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 - & - 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & - 45045.0_8/16.0_8) * cos(4.0_8*phi) + rn(15) = 0.000128521880085575_8 * (w2m1)**2*((14549535.0_8/16.0_8)*w**6 -& + 11486475.0_8/16.0_8 * w**4 + (2027025.0_8/16.0_8)*w**2 - & + 45045.0_8/16.0_8) * cos(4.0_8*phi) ! l = 10, m = 5 rn(16) = 1.35473956745817d-5*(w2m1)**(5.0_8/TWO)* & - ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & - (2027025.0_8/8.0_8)*w) * cos(5.0_8*phi) + ((43648605.0_8/8.0_8)*w**5 - 11486475.0_8/4.0_8 * w**3 + & + (2027025.0_8/8.0_8)*w) * cos(5.0_8*phi) ! l = 10, m = 6 rn(17) = 1.51464488232257d-6*(w2m1)**3*((218243025.0_8/8.0_8)*w**4 - & - 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi) + 34459425.0_8/4.0_8 * w**2 + 2027025.0_8/8.0_8) * cos(6.0_8*phi) ! l = 10, m = 7 rn(18) = 1.83677671621093d-7*(w2m1)**(7.0_8/TWO)* & - ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi) + ((218243025.0_8/TWO)*w**3 - 34459425.0_8/TWO*w) * cos(7.0_8*phi) ! l = 10, m = 8 rn(19) = 2.49953651452314d-8*(w2m1)**4* & - ((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi) + ((654729075.0_8/TWO)*w**2 - 34459425.0_8/TWO) * cos(8.0_8*phi) ! l = 10, m = 9 rn(20) = 2.65478475211798_8 * w*(w2m1)**(9.0_8/TWO) * cos(9.0_8*phi) ! l = 10, m = 10 diff --git a/src/matrix_header.F90 b/src/matrix_header.F90 index 116b6a962e..b2cd3b1291 100644 --- a/src/matrix_header.F90 +++ b/src/matrix_header.F90 @@ -11,19 +11,19 @@ module matrix_header integer, allocatable :: row(:) ! csr row vector integer, allocatable :: col(:) ! column vector real(8), allocatable :: val(:) ! matrix value vector - contains - procedure :: create => matrix_create - procedure :: destroy => matrix_destroy - procedure :: add_value => matrix_add_value - procedure :: new_row => matrix_new_row - procedure :: assemble => matrix_assemble - procedure :: get_row => matrix_get_row - procedure :: get_col => matrix_get_col - procedure :: vector_multiply => matrix_vector_multiply - procedure :: search_indices => matrix_search_indices - procedure :: write => matrix_write - procedure :: copy => matrix_copy - procedure :: transpose => matrix_transpose + contains + procedure :: create => matrix_create + procedure :: destroy => matrix_destroy + procedure :: add_value => matrix_add_value + procedure :: new_row => matrix_new_row + procedure :: assemble => matrix_assemble + procedure :: get_row => matrix_get_row + procedure :: get_col => matrix_get_col + procedure :: vector_multiply => matrix_vector_multiply + procedure :: search_indices => matrix_search_indices + procedure :: write => matrix_write + procedure :: copy => matrix_copy + procedure :: transpose => matrix_transpose end type matrix contains diff --git a/src/mpiio_interface.F90 b/src/mpiio_interface.F90 index cae01a68e9..b9652c8ae4 100644 --- a/src/mpiio_interface.F90 +++ b/src/mpiio_interface.F90 @@ -129,13 +129,13 @@ contains integer, intent(inout) :: buffer ! read data to here logical, intent(in) :: collect ! collective I/O - if (collect) then - call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpiio_err) - else - call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, & - MPI_STATUS_IGNORE, mpiio_err) - end if + if (collect) then + call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + else + call MPI_FILE_READ(fh, buffer, 1, MPI_INTEGER, & + MPI_STATUS_IGNORE, mpiio_err) + end if end subroutine mpi_read_integer @@ -341,13 +341,13 @@ contains real(8), intent(inout) :: buffer ! read data to here logical, intent(in) :: collect ! collective I/O - if (collect) then - call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpiio_err) - else - call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, & - MPI_STATUS_IGNORE, mpiio_err) - end if + if (collect) then + call MPI_FILE_READ_ALL(fh, buffer, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + else + call MPI_FILE_READ(fh, buffer, 1, MPI_REAL8, & + MPI_STATUS_IGNORE, mpiio_err) + end if end subroutine mpi_read_double diff --git a/src/output.F90 b/src/output.F90 index 49908fed6c..cd08e46110 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1372,7 +1372,7 @@ contains if (cmfd_run) then write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" if (cmfd_display /= '') & - write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" + write(UNIT=ou, FMT='(A8,3X)', ADVANCE='NO') "========" end if write(UNIT=ou, FMT=*) @@ -1432,20 +1432,20 @@ contains ! write out cmfd keff if it is active and other display info if (cmfd_on) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % k_cmfd(current_batch) + cmfd % k_cmfd(current_batch) select case(trim(cmfd_display)) case('entropy') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % entropy(current_batch) + cmfd % entropy(current_batch) case('balance') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % balance(current_batch) + cmfd % balance(current_batch) case('source') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % src_cmp(current_batch) + cmfd % src_cmp(current_batch) case('dominance') write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - cmfd % dom(current_batch) + cmfd % dom(current_batch) end select end if @@ -1567,7 +1567,7 @@ contains speed_inactive = real(n_particles * (n_inactive - restart_batch) * & gen_per_batch) / time_inactive % elapsed speed_active = real(n_particles * n_active * gen_per_batch) / & - time_active % elapsed + time_active % elapsed else speed_inactive = ZERO speed_active = real(n_particles * (n_batches - restart_batch) * & @@ -1884,21 +1884,22 @@ contains select case(t % score_bins(k)) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) score_name = 'P' // trim(to_str(t % moment_order(k))) // " " // & - score_names(abs(t % score_bins(k))) + score_names(abs(t % score_bins(k))) write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) case (SCORE_SCATTER_PN, SCORE_NU_SCATTER_PN) score_index = score_index - 1 do n_order = 0, t % moment_order(k) score_index = score_index + 1 score_name = 'P' // trim(to_str(n_order)) // " " //& - score_names(abs(t % score_bins(k))) + score_names(abs(t % score_bins(k))) write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) & + % sum_sq)) end do k = k + t % moment_order(k) case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & @@ -1908,11 +1909,13 @@ contains do nm_order = -n_order, n_order score_index = score_index + 1 score_name = 'Y' // trim(to_str(n_order)) // ',' // & - trim(to_str(nm_order)) // " " // score_names(abs(t % score_bins(k))) + trim(to_str(nm_order)) // " " & + // score_names(abs(t % score_bins(k))) write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index)& + % sum_sq)) end do end do k = k + (t % moment_order(k) + 1)**2 - 1 @@ -1923,9 +1926,9 @@ contains score_name = score_names(abs(t % score_bins(k))) end if write(UNIT=UNIT_TALLY, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(t % results(score_index,filter_index) % sum), & - trim(to_str(t % results(score_index,filter_index) % sum_sq)) + repeat(" ", indent), score_name, & + to_str(t % results(score_index,filter_index) % sum), & + trim(to_str(t % results(score_index,filter_index) % sum_sq)) end select end do indent = indent - 2 @@ -2334,8 +2337,8 @@ contains ! Loop over lattice coordinates do k = 1, n_x - do l = 1, n_y - do m = 1, n_z + do l = 1, n_y + do m = 1, n_z if (final >= lat % offset(map, k, l, m) + offset) then if (k == n_x .and. l == n_y .and. m == n_z) then diff --git a/src/output_interface.F90 b/src/output_interface.F90 index 9a4341135b..05fcd0010b 100644 --- a/src/output_interface.F90 +++ b/src/output_interface.F90 @@ -31,7 +31,7 @@ module output_interface #endif #endif logical :: serial ! Serial I/O when using MPI/PHDF5 - contains + contains generic, public :: write_data => write_double, & write_double_1Darray, & write_double_2Darray, & @@ -111,7 +111,7 @@ contains self % serial = serial else self % serial = .true. - end if + end if #ifdef HDF5 # ifdef MPI @@ -153,7 +153,7 @@ contains self % serial = serial else self % serial = .true. - end if + end if #ifdef HDF5 # ifdef MPI @@ -201,18 +201,18 @@ contains #ifdef HDF5 # ifdef MPI - call hdf5_file_close(self % hdf5_fh) + call hdf5_file_close(self % hdf5_fh) # else - call hdf5_file_close(self % hdf5_fh) + call hdf5_file_close(self % hdf5_fh) # endif #elif MPI - if (self % serial) then - close(UNIT=self % unit_fh) - else - call mpi_close_file(self % mpi_fh) - end if + if (self % serial) then + close(UNIT=self % unit_fh) + else + call mpi_close_file(self % mpi_fh) + end if #else - close(UNIT=self % unit_fh) + close(UNIT=self % unit_fh) #endif end subroutine file_close @@ -475,7 +475,7 @@ contains call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length) else call hdf5_read_double_1Darray_parallel(self % hdf5_grp, name_, buffer, & - length, collect_) + length, collect_) end if # else call hdf5_read_double_1Darray(self % hdf5_grp, name_, buffer, length) @@ -663,8 +663,8 @@ contains if (self % serial) then call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length) else - call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, length, & - collect_) + call hdf5_write_double_3Darray_parallel(self % hdf5_grp, name_, buffer, & + length, collect_) end if # else call hdf5_write_double_3Darray(self % hdf5_grp, name_, buffer, length) diff --git a/src/physics.F90 b/src/physics.F90 index 6f9457ff68..9ca59a9872 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -428,7 +428,7 @@ contains ! Sample velocity of target nucleus if (.not. micro_xs(i_nuclide) % use_ptable) then call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, & - & micro_xs(i_nuclide) % elastic) + & micro_xs(i_nuclide) % elastic) else v_t = ZERO end if @@ -582,7 +582,7 @@ contains ! accompanying PDF and CDF is utilized) if ((sab % secondary_mode == SAB_SECONDARY_EQUAL) .or. & - (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then + (sab % secondary_mode == SAB_SECONDARY_SKEWED)) then if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then ! All bins equally likely @@ -843,16 +843,16 @@ contains ! interpolate xs since we're not exactly at the energy indices xs_low = nuc % elastic_0K(i_E_low) m = (nuc % elastic_0K(i_E_low + 1) - xs_low) & - & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low)) xs_up = nuc % elastic_0K(i_E_up) m = (nuc % elastic_0K(i_E_up + 1) - xs_up) & - & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up)) ! get max 0K xs value over range of practical relative energies xs_max = max(xs_low, & - & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) + & maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up - 1)), xs_up) reject = .true. @@ -898,26 +898,26 @@ contains ! cdf value at lower bound attainable energy if (i_E_low > 1) then m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) & - & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) cdf_low = nuc % xs_cdf(i_E_low - 1) & - & + m * (E_low - nuc % energy_0K(i_E_low)) + & + m * (E_low - nuc % energy_0K(i_E_low)) else m = nuc % xs_cdf(i_E_low) & - & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) + & / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low)) cdf_low = m * (E_low - nuc % energy_0K(i_E_low)) if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO end if ! cdf value at upper bound attainable energy m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) & - & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) + & / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up)) cdf_up = nuc % xs_cdf(i_E_up - 1) & - & + m * (E_up - nuc % energy_0K(i_E_up)) + & + m * (E_up - nuc % energy_0K(i_E_up)) ! values used to sample the Maxwellian E_mode = kT p_mode = TWO * sqrt(E_mode / pi) * sqrt((ONE / kT)**3) & - & * exp(-E_mode / kT) + & * exp(-E_mode / kT) E_t_max = 16.0_8 * E_mode reject = .true. @@ -927,7 +927,7 @@ contains ! perform Maxwellian rejection sampling E_t = E_t_max * prn()**2 p_t = TWO * sqrt(E_t / pi) * sqrt((ONE / kT)**3) & - & * exp(-E_t / kT) + & * exp(-E_t / kT) R_speed = p_t / p_mode if (prn() < R_speed) then @@ -935,12 +935,12 @@ contains ! sample a relative energy using the xs cdf cdf_rel = cdf_low + prn() * (cdf_up - cdf_low) i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), & - & i_E_up - i_E_low + 2, cdf_rel) + & i_E_up - i_E_low + 2, cdf_rel) E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1) m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) & - & - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & - & / (nuc % energy_0K(i_E_low + i_E_rel) & - & - nuc % energy_0K(i_E_low + i_E_rel - 1)) + & - nuc % xs_cdf(i_E_low + i_E_rel - 2)) & + & / (nuc % energy_0K(i_E_low + i_E_rel) & + & - nuc % energy_0K(i_E_low + i_E_rel - 1)) E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m ! perform rejection sampling on cosine between @@ -1026,7 +1026,7 @@ contains ! Determine rejection probability accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) & - /(beta_vn + beta_vt) + /(beta_vn + beta_vt) ! Perform rejection sampling on vt and mu if (prn() < accept_prob) exit diff --git a/src/plot.F90 b/src/plot.F90 index c6341e3a9b..e507b6093b 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -250,7 +250,7 @@ contains do j = ijk_ll(inner), ijk_ur(inner) ! check if we're in the mesh for this ijk if (i > 0 .and. i <= m % dimension(outer) .and. & - j > 0 .and. j <= m % dimension(inner)) then + j > 0 .and. j <= m % dimension(inner)) then ! get xyz's of lower left and upper right of this mesh cell xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 7ce0aa4500..8dc725d9d4 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -24,7 +24,7 @@ module plot_header integer :: color_by ! quantity to color regions by real(8) :: origin(3) ! xyz center of plot location real(8) :: width(3) ! xyz widths of plot - integer :: basis ! direction of plot slice + integer :: basis ! direction of plot slice integer :: pixels(3) ! pixel width/height of plot slice integer :: meshlines_width ! pixel width of meshlines integer :: level ! universe depth to plot the cells of @@ -37,7 +37,7 @@ module plot_header ! Plot type integer, parameter :: PLOT_TYPE_SLICE = 1 integer, parameter :: PLOT_TYPE_VOXEL = 2 - + ! Plot level integer, parameter :: PLOT_LEVEL_LOWEST = -1 diff --git a/src/ppmlib.F90 b/src/ppmlib.F90 index b9df00e1ca..ebd8105bef 100644 --- a/src/ppmlib.F90 +++ b/src/ppmlib.F90 @@ -1,7 +1,7 @@ module ppmlib implicit none - + !=============================================================================== ! Image holds RGB information for output PPM image !=============================================================================== @@ -12,7 +12,7 @@ module ppmlib end type Image contains - + !=============================================================================== ! INIT_IMAGE initializes the Image derived type !=============================================================================== @@ -55,7 +55,7 @@ contains !=============================================================================== ! DEALLOCATE_IMAGE !=============================================================================== - + subroutine deallocate_image(img) type(Image) :: img @@ -70,7 +70,7 @@ contains ! INSIDE_IMAGE determines whether a point (x,y) is inside the image !=============================================================================== - + function inside_image(img, x, y) result(inside) type(Image), intent(in) :: img @@ -83,7 +83,7 @@ contains (x >= 0) .and. (y >= 0)) inside = .true. end function inside_image - + !=============================================================================== ! VALID_IMAGE checks whether the image has a width and height and if its color ! arrays are allocated @@ -123,5 +123,5 @@ contains end if end subroutine set_pixel - + end module ppmlib diff --git a/src/progress_header.F90 b/src/progress_header.F90 index d0ee7563a5..5f462e0fac 100644 --- a/src/progress_header.F90 +++ b/src/progress_header.F90 @@ -22,12 +22,12 @@ module progress_header private character(len=72) :: bar="???% | " // & " |" - contains - procedure :: set_value => bar_set_value + contains + procedure :: set_value => bar_set_value end type ProgressBar contains - + !=============================================================================== ! IS_TERMINAL checks if output is currently being output to a terminal. Relies ! on a POSIX implementation of isatty, and defaults to false if that is not @@ -46,7 +46,7 @@ contains #endif end function is_terminal - + !=============================================================================== ! BAR_SET_VALUE prints the progress bar without advancing. The value is ! specified as percent completion, from 0 to 100. If the value is ever set to @@ -57,7 +57,7 @@ contains class(ProgressBar), intent(inout) :: self real(8), intent(in) :: val - + integer :: i if (.not. is_terminal()) return @@ -84,19 +84,19 @@ contains write(OUTPUT_UNIT, '(A1,A1,A72)', ADVANCE='no') '+', char(13), self % bar flush(OUTPUT_UNIT) - + if (val >= 100.) then - + ! make new line write(OUTPUT_UNIT, "(A)") "" flush(OUTPUT_UNIT) - + ! reset the bar in case we want to use this instance again self % bar = "???% | " // & " |" - + end if - + end subroutine bar_set_value end module progress_header diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index e6a28587df..1ebe651d3c 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -93,7 +93,7 @@ contains end do end subroutine set_particle_seed - + !=============================================================================== ! PRN_SKIP advances the random number seed 'n' times from the current seed !=============================================================================== diff --git a/src/search.F90 b/src/search.F90 index 8de05a955b..dab7fa67ca 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -62,7 +62,7 @@ contains n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then call fatal_error("Reached maximum number of iterations on binary & - &search.") + &search.") end if end do @@ -114,7 +114,7 @@ contains n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then call fatal_error("Reached maximum number of iterations on binary & - &search.") + &search.") end if end do @@ -166,7 +166,7 @@ contains n_iteration = n_iteration + 1 if (n_iteration == MAX_ITERATION) then call fatal_error("Reached maximum number of iterations on binary & - &search.") + &search.") end if end do diff --git a/src/state_point.F90 b/src/state_point.F90 index 0ba0931767..6b983231f6 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -52,7 +52,7 @@ contains ! Set filename for state point filename = trim(path_output) // 'statepoint.' // & - & zero_padded(current_batch, count_digits(n_max_batches)) + & zero_padded(current_batch, count_digits(n_max_batches)) ! Append appropriate extension #ifdef HDF5 @@ -256,7 +256,7 @@ contains group="tallies/tally " // trim(to_str(tally % id)) // & "/filter " // to_str(j)) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then + tally % filters(j) % type == FILTER_ENERGYOUT) then call sp % write_data(tally % filters(j) % real_bins, "bins", & group="tallies/tally " // trim(to_str(tally % id)) // & "/filter " // to_str(j), & @@ -309,8 +309,8 @@ contains do n_order = 0, tally % moment_order(k) moment_name = 'P' // trim(to_str(n_order)) call sp % write_data(moment_name, "order" // trim(to_str(k)), & - group="tallies/tally " // trim(to_str(tally % id)) // & - "/moments") + group="tallies/tally " // trim(to_str(tally % id)) // & + "/moments") k = k + 1 end do case (SCORE_SCATTER_YN, SCORE_NU_SCATTER_YN, SCORE_FLUX_YN, & @@ -318,7 +318,7 @@ contains do n_order = 0, tally % moment_order(k) do nm_order = -n_order, n_order moment_name = 'Y' // trim(to_str(n_order)) // ',' // & - trim(to_str(nm_order)) + trim(to_str(nm_order)) call sp % write_data(moment_name, "order" // & trim(to_str(k)), & group="tallies/tally " // trim(to_str(tally % id)) // & @@ -412,7 +412,7 @@ contains ! Set filename filename = trim(path_output) // 'source.' // & - & zero_padded(current_batch, count_digits(n_max_batches)) + & zero_padded(current_batch, count_digits(n_max_batches)) #ifdef HDF5 filename = trim(filename) // '.h5' @@ -434,7 +434,7 @@ contains ! Set filename for state point filename = trim(path_output) // 'statepoint.' // & - & zero_padded(current_batch, count_digits(n_max_batches)) + & zero_padded(current_batch, count_digits(n_max_batches)) #ifdef HDF5 filename = trim(filename) // '.h5' #else @@ -607,12 +607,12 @@ contains tally % results(:,:) % sum_sq = tally_temp(2,:,:) end if - ! Put in temporary tally result - allocate(tallyresult_temp(m,n)) - tallyresult_temp(:,:) % sum = tally_temp(1,:,:) - tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:) + ! Put in temporary tally result + allocate(tallyresult_temp(m,n)) + tallyresult_temp(:,:) % sum = tally_temp(1,:,:) + tallyresult_temp(:,:) % sum_sq = tally_temp(2,:,:) - ! Write reduced tally results to file + ! Write reduced tally results to file call sp % write_tally_result(tally % results, "results", & group="tallies/tally " // trim(to_str(tally % id)), n1=m, n2=n) @@ -687,7 +687,7 @@ contains call sp % read_data(int_array(2), "version_minor") call sp % read_data(int_array(3), "version_release") if (int_array(1) /= VERSION_MAJOR .or. int_array(2) /= VERSION_MINOR & - .or. int_array(3) /= VERSION_RELEASE) then + .or. int_array(3) /= VERSION_RELEASE) then if (master) call warning("State point file was created with a different & &version of OpenMC.") end if @@ -843,7 +843,7 @@ contains group="tallies/tally " // trim(to_str(curr_key)) // & "/filter " // to_str(j)) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & - tally % filters(j) % type == FILTER_ENERGYOUT) then + tally % filters(j) % type == FILTER_ENERGYOUT) then call sp % read_data(tally % filters(j) % real_bins, "bins", & group="tallies/tally " // trim(to_str(curr_key)) // & "/filter " // to_str(j), & diff --git a/src/string.F90 b/src/string.F90 index ebe2310f24..5f57bb3855 100644 --- a/src/string.F90 +++ b/src/string.F90 @@ -7,7 +7,7 @@ module string implicit none interface to_str - module procedure int4_to_str, int8_to_str, real_to_str + module procedure int4_to_str, int8_to_str, real_to_str end interface contains @@ -50,7 +50,7 @@ contains n = n + 1 if (i_end - i_start + 1 > len(words(n))) then if (master) call warning("The word '" // string(i_start:i_end) & - &// "' is longer than the space allocated for it.") + &// "' is longer than the space allocated for it.") end if words(n) = string(i_start:i_end) ! reset indices diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 43ec868c54..d20c3fea20 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -124,7 +124,7 @@ module tally_header ! Number of realizations of tally random variables integer :: n_realizations = 0 - + ! Tally precision triggers integer :: n_triggers = 0 ! # of triggers type(TriggerObject), allocatable :: triggers(:) ! Array of triggers @@ -197,10 +197,10 @@ module tally_header this % reset = .false. this % n_realizations = 0 - + if (allocated(this % triggers)) & deallocate (this % triggers) - + this % n_triggers = 0 end subroutine tallyobject_clear diff --git a/src/timer_header.F90 b/src/timer_header.F90 index 0dde85d0c2..0bf1b7aef5 100644 --- a/src/timer_header.F90 +++ b/src/timer_header.F90 @@ -15,11 +15,11 @@ module timer_header logical :: running = .false. ! is timer running? integer :: start_counts = 0 ! counts when started real(8), public :: elapsed = ZERO ! total time elapsed in seconds - contains - procedure :: start => timer_start - procedure :: get_value => timer_get_value - procedure :: stop => timer_stop - procedure :: reset => timer_reset + contains + procedure :: start => timer_start + procedure :: get_value => timer_get_value + procedure :: stop => timer_stop + procedure :: reset => timer_reset end type Timer contains @@ -83,7 +83,7 @@ contains !=============================================================================== subroutine timer_reset(self) - + class(Timer), intent(inout) :: self self % running = .false. diff --git a/src/trigger.F90 b/src/trigger.F90 index 6909598db2..4a16cd5cab 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -43,7 +43,7 @@ contains ! When trigger threshold is reached, write information if (satisfy_triggers) then call write_message("Triggers satisfied for batch " // & - trim(to_str(current_batch))) + trim(to_str(current_batch))) ! When trigger is not reached write convergence info for user elseif (name == "eigenvalue") then @@ -83,7 +83,7 @@ contains ! and finds the maximum uncertainty/threshold ratio for all triggers !=============================================================================== - subroutine check_tally_triggers(max_ratio, tally_id, name) + subroutine check_tally_triggers(max_ratio, tally_id, name) ! Variables to reflect distance to trigger convergence criteria real(8), intent(inout) :: max_ratio ! max uncertainty/thresh ratio @@ -299,7 +299,7 @@ contains ! precision trigger(s). !=============================================================================== - subroutine compute_tally_current(t, trigger) + subroutine compute_tally_current(t, trigger) integer :: i ! mesh index for x integer :: j ! mesh index for y diff --git a/src/trigger_header.F90 b/src/trigger_header.F90 index 8289ca6154..e137829bd0 100644 --- a/src/trigger_header.F90 +++ b/src/trigger_header.F90 @@ -14,16 +14,16 @@ module trigger_header character(len=52) :: score_name ! the name of the score integer :: score_index ! the index of the score real(8) :: variance=0.0 ! temp variance container - real(8) :: std_dev =0.0 ! temp std. dev. container + real(8) :: std_dev =0.0 ! temp std. dev. container real(8) :: rel_err =0.0 ! temp rel. err. container end type TriggerObject - + !=============================================================================== ! KTRIGGER describes a user-specified precision trigger for k-effective !=============================================================================== type KTrigger integer :: trigger_type = 0 - real(8) :: threshold = 0 + real(8) :: threshold = 0 end type KTrigger end module trigger_header diff --git a/src/vector_header.F90 b/src/vector_header.F90 index 424aedd686..1ada0fba74 100644 --- a/src/vector_header.F90 +++ b/src/vector_header.F90 @@ -9,12 +9,12 @@ module vector_header integer :: n ! number of rows/cols in matrix real(8), allocatable :: data(:) ! where vector data is stored real(8), pointer :: val(:) ! pointer to vector data - contains - procedure :: create => vector_create - procedure :: destroy => vector_destroy - procedure :: add_value => vector_add_value - procedure :: copy => vector_copy - ! TODO: procedure :: write => vector_write + contains + procedure :: create => vector_create + procedure :: destroy => vector_destroy + procedure :: add_value => vector_add_value + procedure :: copy => vector_copy + ! TODO: procedure :: write => vector_write end type Vector contains diff --git a/src/xml_interface.F90 b/src/xml_interface.F90 index 1a24bc14b3..4ae05ee281 100644 --- a/src/xml_interface.F90 +++ b/src/xml_interface.F90 @@ -92,7 +92,7 @@ contains ! Check if node exists if (associated(temp_ptr)) return - ! Check for a sub-element + ! Check for a sub-element elem_list => getChildrenByTagName(ptr, trim(node_name)) ! Get the length of the list @@ -122,7 +122,7 @@ contains ! Set found to false found_ = .false. - ! Check for a sub-element + ! Check for a sub-element elem_list => getChildrenByTagName(in_ptr, trim(node_name)) ! Get the length of the list @@ -147,7 +147,7 @@ contains type(Node), pointer, intent(in) :: in_ptr type(NodeList), pointer, intent(out) :: out_ptr - ! Check for a sub-element + ! Check for a sub-element out_ptr => getChildrenByTagName(in_ptr, trim(node_name)) end subroutine get_node_list @@ -176,7 +176,7 @@ contains type(NodeList), pointer, intent(in) :: in_ptr type(Node), pointer, intent(out) :: out_ptr - ! Check for a sub-element + ! Check for a sub-element out_ptr => item(in_ptr, idx - 1) end subroutine get_list_item @@ -199,7 +199,7 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if @@ -210,7 +210,7 @@ contains else call extractDataContent(temp_ptr, result_int) end if - + end subroutine get_node_value_integer !=============================================================================== @@ -231,18 +231,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_long) else call extractDataContent(temp_ptr, result_long) end if - + end subroutine get_node_value_long !=============================================================================== @@ -263,18 +263,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_double) else call extractDataContent(temp_ptr, result_double) end if - + end subroutine get_node_value_double !=============================================================================== @@ -295,18 +295,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_int) else call extractDataContent(temp_ptr, result_int) end if - + end subroutine get_node_array_integer !=============================================================================== @@ -327,18 +327,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_double) else call extractDataContent(temp_ptr, result_double) end if - + end subroutine get_node_array_double !=============================================================================== @@ -359,18 +359,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " & &// getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_string) else call extractDataContent(temp_ptr, result_string) end if - + end subroutine get_node_array_string !=============================================================================== @@ -391,18 +391,18 @@ contains call get_node(ptr, node_name, temp_ptr, node_type, found) ! Leave if it was not found - if (.not. found) then + if (.not. found) then call fatal_error("Node " // node_name // " not part of Node " // & getNodeName(ptr) // ".") end if - + ! Extract value if (node_type == ATTR_NODE) then call extractDataAttribute(ptr, node_name, result_str) else call extractDataContent(temp_ptr, result_str) end if - + end subroutine get_node_value_string !=============================================================================== @@ -513,7 +513,7 @@ contains ! Check if node exists if (associated(out_ptr)) return - ! Check for a sub-element + ! Check for a sub-element elem_list => getChildrenByTagName(in_ptr, trim(node_name)) ! Get the length of the list From a8b3a0b2cd5da23633d7ee6c69879bfaf6c99378 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 26 Jul 2015 12:23:21 -0600 Subject: [PATCH 3/8] Fix some PEP-8 for #423 --- tests/check_source.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/tests/check_source.py b/tests/check_source.py index 490eb9c7c3..3a79d44c28 100755 --- a/tests/check_source.py +++ b/tests/check_source.py @@ -22,7 +22,7 @@ class LineOfCode(object): self.content = content self.is_continuation = is_continuation self.is_continued = False - assert (string_terminator == None or + assert (string_terminator is None or string_terminator == "'" or string_terminator == '"') self.initial_string_terminator = string_terminator @@ -114,12 +114,10 @@ class LineOfCode(object): if (not has_real_content) and in_a_comment: self.is_comment_only = True - def get_indent(self): """Return the number of indentation spaces prefixing this line.""" return len(self.content) - len(self.content.lstrip(' ')) - def contains_whitespace_only(self): """Return True/False if all characters in the line are whitespace.""" if len(self.content.strip('\n')) == 0: return False @@ -127,7 +125,6 @@ class LineOfCode(object): for char in self.content.strip('\n') ] return all(is_ws) - def has_trailing_whitespace(self): """Return True/False if this line ends with a whitespace character.""" stripped = self.content.strip('\n') @@ -181,15 +178,15 @@ def check_source(fname): # Check for extra whitespace errors. if loc.contains_whitespace_only(): good_code = False - print_error(fname, loc.number, "Line contains whitespace " \ + print_error(fname, loc.number, "Line contains whitespace " "characters but no content. Please remove whitespace.") elif loc.has_trailing_whitespace(): good_code = False - print_error(fname, loc.number, "Line has trailing whitespace after"\ + print_error(fname, loc.number, "Line has trailing whitespace after" " the content. Please remove trailing whitespace.") if loc.contains_tab(): good_code = False - print_error(fname, loc.number, "Line contains a tab character. " \ + print_error(fname, loc.number, "Line contains a tab character. " "Please replace with single whitespace characters.") # Check indentation. @@ -197,13 +194,13 @@ def check_source(fname): if ((not loc.is_continuation) and (not loc.is_comment_only) and (not loc.contains_whitespace_only()) and current_indent % 2 != 0): good_code = False - print_error(fname, loc.number, "Line is indented an odd number of "\ - "spaces. All non-continuation lines should be indented an "\ + print_error(fname, loc.number, "Line is indented an odd number of " + "spaces. All non-continuation lines should be indented an " "even number of spaces.") if loc.is_continuation and current_indent < base_indent + 5: good_code = False - print_error(fname, loc.number, "Continuation lines must be "\ - "indented by at least 5 spaces, but this line is indented {0}"\ + print_error(fname, loc.number, "Continuation lines must be " + "indented by at least 5 spaces, but this line is indented {0}" " spaces.".format(current_indent - base_indent)) # Set base indentation for next lines. From 315236ea58b3fe1d7f016f66a6799ffdad8017a9 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 27 Jul 2015 20:24:10 -0600 Subject: [PATCH 4/8] Fix PyAPI HexLattice.get_unique_universes() closes #413 --- openmc/checkvalue.py | 78 ++++++++++++++++++++++++++++++++++++++++++++ openmc/universe.py | 37 ++++++++++++++------- 2 files changed, 104 insertions(+), 11 deletions(-) diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index 01613908ba..f9e05e3f8c 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -1,3 +1,5 @@ +from collections import Iterable + def check_type(name, value, expected_type, expected_iter_type=None): """Ensure that an object is of an expected type. Optionally, if the object is iterable, check that each element is of a particular type. @@ -30,6 +32,82 @@ def check_type(name, value, expected_type, expected_iter_type=None): raise ValueError(msg) +def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1): + """Ensure that an object is an iterable containing an expected type. + + Parameters + ---------- + name : str + Description of value being checked + value : Iterable + Iterable, possibly of other iterables, that should ultimately contain + the expected type + expected_type : type + type that the iterable should contain + min_depth : int + The minimum number of layers of nested iterables there should be before + reaching the ultimately contained items + max_depth : int + The maximum number of layers of nested iterables ... + """ + # Initialize the tree at the very first item. + tree = [value] + index = [0] + + # Traverse the tree. + while index[0] != len(tree[0]): + # If we are done with this level of the tree, go to the next branch on + # the level above this one. + if index[-1] == len(tree[-1]): + del index[-1] + del tree[-1] + index[-1] += 1 + continue + + # Get a string representation of the current index in case we raise an + # exception. + form = '[' + '{:d}, '*(len(index)-1) + '{:d}]' + ind_str = form.format(*index) + + # What is the current item we are looking at? + current_item = tree[-1][index[-1]] + + # If this item is of the expected type, then we've reached the bottom + # level of this branch. + if isinstance(current_item, expected_type): + # Is this deep enough? + if len(tree) < min_depth: + msg = 'Error setting {0}: The item at {1} does not meet the ' \ + 'minimum depth of {2}'.format(name, ind_str, min_depth) + raise ValueError(msg) + + # This item is okay. Move on to the next item. + index[-1] += 1 + + # If this item is not of the expected type, then it's either an error or + # another level of the tree that we need to pursue deeper. + else: + if isinstance(current_item, Iterable): + # The tree goes deeper here, let's explore it. + tree.append(current_item) + index.append(0) + + # But first, have we exceeded the max depth? + if len(tree) > max_depth: + msg = 'Error setting {0}: Found an iterable at {1}, items '\ + 'in that iterable excceed the maximum depth of {2}' \ + .format(name, ind_str, max_depth) + raise ValueError(msg) + + else: + # This item is completely unexected. + msg = "Error setting {0}: Items must be of type '{1}', but " \ + "item at {2} is of type '{3}'"\ + .format(name, expected_type.__name__, ind_str, + type(current_item).__name__) + raise ValueError(msg) + + def check_length(name, value, length_min, length_max=None): """Ensure that a sized object has length within a given range. diff --git a/openmc/universe.py b/openmc/universe.py index 84ce3696ee..f424b1a944 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -7,7 +7,8 @@ import sys import numpy as np import openmc -from openmc.checkvalue import check_type, check_length, check_greater_than +from openmc.checkvalue import check_type, check_iterable_type, check_length, \ + check_greater_than if sys.version_info[0] >= 3: basestring = str @@ -687,8 +688,11 @@ class Lattice(object): @universes.setter def universes(self, universes): - check_type('lattice universes', universes, Iterable) - self._universes = np.asarray(universes, dtype=Universe) + #check_type('lattice universes', universes, Iterable) + check_iterable_type('lattice universes', universes, Universe, + min_depth=2, max_depth=3) + self._universes = universes + #self._universes = np.asarray(universes, dtype=Universe) def get_unique_universes(self): """Determine all unique universes in the lattice @@ -701,13 +705,24 @@ class Lattice(object): """ - unique_universes = np.unique(self._universes.ravel()) - universes = {} + univs = dict() + for k in range(len(self._universes)): + for j in range(len(self._universes[k])): + if isinstance(self._universes[k][j], Universe): + u = self._universes[k][j] + if u._id not in univs: + univs[u._id] = u + else: + for i in range(len(self._universes[k][j])): + u = self._universes[k][j][i] + assert isinstance(u, Universe) + if u._id not in univs: + univs[u._id] = u - for universe in unique_universes: - universes[universe._id] = universe + if self._outer._id not in univs: + univs[self._outer._id] = self._outer - return universes + return univs def get_all_nuclides(self): """Return all nuclides contained in the lattice @@ -1091,14 +1106,14 @@ class HexLattice(Lattice): # Set the number of axial positions. if n_dims == 3: - self.num_axial = self._universes.shape[0] + self.num_axial = len(self._universes) else: self._num_axial = None # Set the number of rings and make sure this number is consistent for # all axial positions. if n_dims == 3: - self.num_rings = len(self._universes[0]) + self.num_rings = len(self._universes) for rings in self._universes: if len(rings) != self._num_rings: msg = 'HexLattice ID={0:d} has an inconsistent number of ' \ @@ -1106,7 +1121,7 @@ class HexLattice(Lattice): raise ValueError(msg) else: - self.num_rings = self._universes.shape[0] + self.num_rings = len(self._universes) # Make sure there are the correct number of elements in each ring. if n_dims == 3: From 2cda46e38467f15078e9b1cee5dd788fbae411e1 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 1 Aug 2015 22:00:38 -0600 Subject: [PATCH 5/8] Fix Numpy type issue Numpy integer types aren't recognized as an instance of Integral, at least not with numpy 1.8.2 --- openmc/checkvalue.py | 18 +++++++++++++++--- openmc/filter.py | 21 ++++++--------------- 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index f9e05e3f8c..5b5f302e5a 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -1,4 +1,16 @@ from collections import Iterable +from numbers import Integral + +import numpy as np + +def _isinstance(value, expected_type): + """A Numpy-aware replacement for isinstance""" + np_ints = [np.intc, np.intp, np.int8, np.int16, np.int32, np.int64, + np.uint8, np.uint16, np.uint32, np.uint64] + if expected_type is Integral: + types = np_ints + [Integral] + return any(isinstance(value, t) for t in types) + return isinstance(value, expected_type) def check_type(name, value, expected_type, expected_iter_type=None): """Ensure that an object is of an expected type. Optionally, if the object is @@ -18,14 +30,14 @@ def check_type(name, value, expected_type, expected_iter_type=None): """ - if not isinstance(value, expected_type): + if not _isinstance(value, expected_type): msg = 'Unable to set {0} to {1} which is not of type {2}'.format( name, value, expected_type.__name__) raise ValueError(msg) if expected_iter_type: for item in value: - if not isinstance(item, expected_iter_type): + if not _isinstance(item, expected_iter_type): msg = 'Unable to set {0} to {1} since each item must be ' \ 'of type {2}'.format(name, value, expected_iter_type.__name__) @@ -74,7 +86,7 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1): # If this item is of the expected type, then we've reached the bottom # level of this branch. - if isinstance(current_item, expected_type): + if _isinstance(current_item, expected_type): # Is this deep enough? if len(tree) < min_depth: msg = 'Error setting {0}: The item at {1} does not meet the ' \ diff --git a/openmc/filter.py b/openmc/filter.py index 09928e235a..ca5d430a95 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -6,7 +6,8 @@ import numpy as np from openmc import Mesh from openmc.constants import * -from openmc.checkvalue import check_type +from openmc.checkvalue import check_type, check_iterable_type, \ + check_greater_than class Filter(object): """A filter used to constrain a tally to a specific criterion, e.g. only tally @@ -134,15 +135,9 @@ class Filter(object): if self._type in ['cell', 'cellborn', 'surface', 'material', 'universe', 'distribcell']: + check_iterable_type('filter bins', bins, Integral) for edge in bins: - if not isinstance(edge, Integral): - msg = 'Unable to add bin "{0}" to a {1} Filter since ' \ - 'it is not an integer'.format(edge, self._type) - raise ValueError(msg) - elif edge < 0: - msg = 'Unable to add bin "{0}" to a {1} Filter since ' \ - 'it is negative'.format(edge, self._type) - raise ValueError(msg) + check_greater_than('filter bin', edge, 0, equality=True) elif self._type in ['energy', 'energyout']: for edge in bins: @@ -185,12 +180,8 @@ class Filter(object): # FIXME @num_bins.setter def num_bins(self, num_bins): - if not isinstance(num_bins, Integral) or num_bins < 0: - msg = 'Unable to set the number of bins "{0}" for a {1} Filter ' \ - 'since it is not a positive ' \ - 'integer'.format(num_bins, self._type) - raise ValueError(msg) - + check_type('filter num_bins', num_bins, Integral) + check_greater_than('filter num_bins', num_bins, 0, equality=True) self._num_bins = num_bins @mesh.setter From 4ded3a39e17b2c646f07902a8672f70b6ba605b8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 2 Aug 2015 12:24:06 -0600 Subject: [PATCH 6/8] Forbid trailing whitespace in the style guide --- docs/source/devguide/styleguide.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/source/devguide/styleguide.rst b/docs/source/devguide/styleguide.rst index 25c58ae773..c8644e45cc 100644 --- a/docs/source/devguide/styleguide.rst +++ b/docs/source/devguide/styleguide.rst @@ -153,6 +153,8 @@ Avoid extraneous whitespace in the following situations: Yes: if (variable == 2) then No: if ( variable==2 ) then +Do not leave trailing whitespace at the end of a line. + ------ Python ------ From 1dca4d842bfff20f6ea3cb8a9a5faab2734f8792 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 5 Aug 2015 19:10:03 -0600 Subject: [PATCH 7/8] Fix hex lattices in HDF5 summary --- openmc/summary.py | 72 ++++++++++++++++++++++++++++++++++++++------ src/hdf5_summary.F90 | 8 ++--- 2 files changed, 67 insertions(+), 13 deletions(-) diff --git a/openmc/summary.py b/openmc/summary.py index c0673ef233..2165d114ad 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -411,16 +411,70 @@ class Summary(object): if outer != -22: lattice.outer = self.universes[outer] - # Build array of Universe pointers for the Lattice - universes = \ - np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe) + # Build array of Universe pointers for the Lattice. Note that + # we need to convert between the HDF5's square array of + # (x, alpha, z) to the PythonAPI's format of a ragged nested + # nested list of (z, ring, theta). + universes = [] + for z in range(lattice.num_axial): + # Add a list for this axial level. + universes.append([]) + x = lattice.num_rings - 1 + a = 2*lattice.num_rings - 2 + for r in range(lattice.num_rings - 1, 0, -1): + # Add a list for this ring. + universes[-1].append([]) - for i in range(universe_ids.shape[0]): - for j in range(universe_ids.shape[1]): - for k in range(universe_ids.shape[2]): - if universe_ids[i, j, k] != -1: - universes[i, j, k] = self.get_universe_by_id( - universe_ids[i, j, k]) + # Climb down the top-right. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + x += 1 + a -= 1 + + # Climb down the right. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + a -= 1 + + # Climb down the bottom-right. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + x -= 1 + + # Climb up the bottom-left. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + x -= 1 + a += 1 + + # Climb up the left. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + a += 1 + + # Climb up the top-left. + for i in range(r): + universes[-1][-1].append(universe_ids[z, a, x]) + x += 1 + + # Move down to the next ring. + a -= 1 + + # Convert the ids into Universe objects. + universes[-1][-1] = [self.get_universe_by_id(u_id) + for u_id in universes[-1][-1]] + + # Handle the degenerate center ring separately. + u_id = universe_ids[z, a, x] + universes[-1].append([self.get_universe_by_id(u_id)]) + + # Add the universes to the lattice. + if len(pitch) == 2: + # Lattice is 3D + lattice.universes = universes + else: + # Lattice is 2D; extract the only axial level + lattice.universes = universes[0] if offset_size > 0: lattice.offsets = offsets diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index 161246f1c9..15caf2edfa 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -394,7 +394,7 @@ contains ! Write number of lattice cells. call su % write_data(lat % n_rings, "n_rings", & group="geometry/lattices/lattice " // trim(to_str(lat % id))) - call su % write_data(lat % n_rings, "n_axial", & + call su % write_data(lat % n_axial, "n_axial", & group="geometry/lattices/lattice " // trim(to_str(lat % id))) ! Write lattice center, pitch and outer universe. @@ -407,10 +407,10 @@ contains end if if (lat % is_3d) then - call su % write_data(lat % pitch, "pitch", length=3, & + call su % write_data(lat % pitch, "pitch", length=2, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) else - call su % write_data(lat % pitch, "pitch", length=2, & + call su % write_data(lat % pitch, "pitch", length=1, & group="geometry/lattices/lattice " // trim(to_str(lat % id))) end if @@ -447,7 +447,7 @@ contains end do end do call su % write_data(lattice_universes, "universes", & - &length=(/lat % n_axial, 2*lat % n_rings-1, 2*lat % n_rings-1/), & + &length=(/2*lat % n_rings-1, 2*lat % n_rings-1, lat % n_axial/), & &group="geometry/lattices/lattice " // trim(to_str(lat % id))) deallocate(lattice_universes) end select From 65679c22207e0a1c09435fe3ea7229863406d5c3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 7 Aug 2015 20:36:21 -0600 Subject: [PATCH 8/8] Small fixes for #425 --- openmc/checkvalue.py | 38 +++++++++++++++----- openmc/summary.py | 4 +-- openmc/temp.py | 12 +++++++ openmc/universe.py | 84 ++++++++++++++++++++------------------------ 4 files changed, 82 insertions(+), 56 deletions(-) create mode 100644 openmc/temp.py diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index d3f782ac6f..787052f5d6 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -1,15 +1,34 @@ from collections import Iterable -from numbers import Integral +from numbers import Integral, Real import numpy as np def _isinstance(value, expected_type): - """A Numpy-aware replacement for isinstance""" - np_ints = [np.intc, np.intp, np.int8, np.int16, np.int32, np.int64, - np.uint8, np.uint16, np.uint32, np.uint64] - if expected_type is Integral: - types = np_ints + [Integral] - return any(isinstance(value, t) for t in types) + """A Numpy-aware replacement for isinstance + + This function will be obsolete when Numpy v. >= 1.9 is established. + """ + + # Declare numpy numeric types. + np_ints = (np.int_, np.intc, np.intp, np.int8, np.int16, np.int32, np.int64, + np.uint8, np.uint16, np.uint32, np.uint64) + np_floats = (np.float_, np.float16, np.float32, np.float64) + + # Include numpy integers, if necessary. + if type(expected_type) is tuple: + if Integral in expected_type: + expected_type = expected_type + np_ints + elif expected_type is Integral: + expected_type = (Integral, ) + np_ints + + # Include numpy floats, if necessary. + if type(expected_type) is tuple: + if Real in expected_type: + expected_type = expected_type + np_floats + elif expected_type is Real: + expected_type = (Real, ) + np_floats + + # Now, make the instance check. return isinstance(value, expected_type) def check_type(name, value, expected_type, expected_iter_type=None): @@ -60,7 +79,8 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1): The minimum number of layers of nested iterables there should be before reaching the ultimately contained items max_depth : int - The maximum number of layers of nested iterables ... + The maximum number of layers of nested iterables there should be before + reaching the ultimately contained items """ # Initialize the tree at the very first item. tree = [value] @@ -112,7 +132,7 @@ def check_iterable_type(name, value, expected_type, min_depth=1, max_depth=1): raise ValueError(msg) else: - # This item is completely unexected. + # This item is completely unexpected. msg = "Error setting {0}: Items must be of type '{1}', but " \ "item at {2} is of type '{3}'"\ .format(name, expected_type.__name__, ind_str, diff --git a/openmc/summary.py b/openmc/summary.py index 2165d114ad..7f9d5387e9 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -413,8 +413,8 @@ class Summary(object): # Build array of Universe pointers for the Lattice. Note that # we need to convert between the HDF5's square array of - # (x, alpha, z) to the PythonAPI's format of a ragged nested - # nested list of (z, ring, theta). + # (x, alpha, z) to the Python API's format of a ragged nested + # list of (z, ring, theta). universes = [] for z in range(lattice.num_axial): # Add a list for this axial level. diff --git a/openmc/temp.py b/openmc/temp.py new file mode 100644 index 0000000000..91f6082996 --- /dev/null +++ b/openmc/temp.py @@ -0,0 +1,12 @@ +from checkvalue import * +from checkvalue import _isinstance + +import numpy as np + +zs = np.zeros((2,)) + +print _isinstance(zs[0], Integral) +print _isinstance(zs[0], Real) +print _isinstance(zs[0], (Integral, Real)) + +print check_iterable_type('thing', zs, (Real, Integral)) diff --git a/openmc/universe.py b/openmc/universe.py index 2785574008..75b6737465 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -7,8 +7,7 @@ import sys import numpy as np import openmc -from openmc.checkvalue import check_type, check_iterable_type, check_length, \ - check_greater_than +import openmc.checkvalue as cv if sys.version_info[0] >= 3: basestring = str @@ -112,13 +111,13 @@ class Cell(object): self._id = AUTO_CELL_ID AUTO_CELL_ID += 1 else: - check_type('cell ID', cell_id, Integral) - check_greater_than('cell ID', cell_id, 0) + cv.check_type('cell ID', cell_id, Integral) + cv.check_greater_than('cell ID', cell_id, 0) self._id = cell_id @name.setter def name(self, name): - check_type('cell name', name, basestring) + cv.check_type('cell name', name, basestring) self._name = name @fill.setter @@ -149,19 +148,19 @@ class Cell(object): @rotation.setter def rotation(self, rotation): - check_type('cell rotation', rotation, Iterable, Real) - check_length('cell rotation', rotation, 3) + cv.check_type('cell rotation', rotation, Iterable, Real) + cv.check_length('cell rotation', rotation, 3) self._rotation = rotation @translation.setter def translation(self, translation): - check_type('cell translation', translation, Iterable, Real) - check_length('cell translation', translation, 3) + cv.check_type('cell translation', translation, Iterable, Real) + cv.check_length('cell translation', translation, 3) self._translation = translation @offsets.setter def offsets(self, offsets): - check_type('cell offsets', offsets, Iterable) + cv.check_type('cell offsets', offsets, Iterable) self._offsets = offsets def add_surface(self, surface, halfspace): @@ -433,13 +432,13 @@ class Universe(object): self._id = AUTO_UNIVERSE_ID AUTO_UNIVERSE_ID += 1 else: - check_type('universe ID', universe_id, Integral) - check_greater_than('universe ID', universe_id, 0, True) + cv.check_type('universe ID', universe_id, Integral) + cv.check_greater_than('universe ID', universe_id, 0, True) self._id = universe_id @name.setter def name(self, name): - check_type('universe name', name, basestring) + cv.check_type('universe name', name, basestring) self._name = name def add_cell(self, cell): @@ -672,27 +671,25 @@ class Lattice(object): self._id = AUTO_UNIVERSE_ID AUTO_UNIVERSE_ID += 1 else: - check_type('lattice ID', lattice_id, Integral) - check_greater_than('lattice ID', lattice_id, 0) + cv.check_type('lattice ID', lattice_id, Integral) + cv.check_greater_than('lattice ID', lattice_id, 0) self._id = lattice_id @name.setter def name(self, name): - check_type('lattice name', name, basestring) + cv.check_type('lattice name', name, basestring) self._name = name @outer.setter def outer(self, outer): - check_type('outer universe', outer, Universe) + cv.check_type('outer universe', outer, Universe) self._outer = outer @universes.setter def universes(self, universes): - #check_type('lattice universes', universes, Iterable) - check_iterable_type('lattice universes', universes, Universe, - min_depth=2, max_depth=3) + cv.check_iterable_type('lattice universes', universes, Universe, + min_depth=2, max_depth=3) self._universes = universes - #self._universes = np.asarray(universes, dtype=Universe) def get_unique_universes(self): """Determine all unique universes in the lattice @@ -710,17 +707,14 @@ class Lattice(object): for j in range(len(self._universes[k])): if isinstance(self._universes[k][j], Universe): u = self._universes[k][j] - if u._id not in univs: - univs[u._id] = u + univs[u._id] = u else: for i in range(len(self._universes[k][j])): u = self._universes[k][j][i] assert isinstance(u, Universe) - if u._id not in univs: - univs[u._id] = u + univs[u._id] = u - if self._outer._id not in univs: - univs[self._outer._id] = self._outer + univs[self._outer._id] = self._outer return univs @@ -840,29 +834,29 @@ class RectLattice(Lattice): @dimension.setter def dimension(self, dimension): - check_type('lattice dimension', dimension, Iterable, Integral) - check_length('lattice dimension', dimension, 2, 3) + cv.check_type('lattice dimension', dimension, Iterable, Integral) + cv.check_length('lattice dimension', dimension, 2, 3) for dim in dimension: - check_greater_than('lattice dimension', dim, 0) + cv.check_greater_than('lattice dimension', dim, 0) self._dimension = dimension @lower_left.setter def lower_left(self, lower_left): - check_type('lattice lower left corner', lower_left, Iterable, Real) - check_length('lattice lower left corner', lower_left, 2, 3) + cv.check_type('lattice lower left corner', lower_left, Iterable, Real) + cv.check_length('lattice lower left corner', lower_left, 2, 3) self._lower_left = lower_left @offsets.setter def offsets(self, offsets): - check_type('lattice offsets', offsets, Iterable) + cv.check_type('lattice offsets', offsets, Iterable) self._offsets = offsets @Lattice.pitch.setter def pitch(self, pitch): - check_type('lattice pitch', pitch, Iterable, Real) - check_length('lattice pitch', pitch, 2, 3) + cv.check_type('lattice pitch', pitch, Iterable, Real) + cv.check_length('lattice pitch', pitch, 2, 3) for dim in pitch: - check_greater_than('lattice pitch', dim, 0.0) + cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch def get_offset(self, path, filter_offset): @@ -1056,28 +1050,28 @@ class HexLattice(Lattice): @num_rings.setter def num_rings(self, num_rings): - check_type('number of rings', num_rings, Integral) - check_greater_than('number of rings', num_rings, 0) + cv.check_type('number of rings', num_rings, Integral) + cv.check_greater_than('number of rings', num_rings, 0) self._num_rings = num_rings @num_axial.setter def num_axial(self, num_axial): - check_type('number of axial', num_axial, Integral) - check_greater_than('number of axial', num_axial, 0) + cv.check_type('number of axial', num_axial, Integral) + cv.check_greater_than('number of axial', num_axial, 0) self._num_axial = num_axial @center.setter def center(self, center): - check_type('lattice center', center, Iterable, Real) - check_length('lattice center', center, 2, 3) + cv.check_type('lattice center', center, Iterable, Real) + cv.check_length('lattice center', center, 2, 3) self._center = center @Lattice.pitch.setter def pitch(self, pitch): - check_type('lattice pitch', pitch, Iterable, Real) - check_length('lattice pitch', pitch, 1, 2) + cv.check_type('lattice pitch', pitch, Iterable, Real) + cv.check_length('lattice pitch', pitch, 1, 2) for dim in pitch: - check_greater_than('lattice pitch', dim, 0) + cv.check_greater_than('lattice pitch', dim, 0) self._pitch = pitch @Lattice.universes.setter