Merge remote-tracking branch 'upstream/develop' into mgxs

This commit is contained in:
Will Boyd 2015-08-09 21:23:45 -07:00
commit 5ee36fc2a8
48 changed files with 1154 additions and 722 deletions

View file

@ -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
------

View file

@ -1,3 +1,36 @@
from collections import Iterable
from numbers import Integral, Real
import numpy as np
def _isinstance(value, expected_type):
"""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):
"""Ensure that an object is of an expected type. Optionally, if the object is
iterable, check that each element is of a particular type.
@ -16,20 +49,97 @@ 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__)
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 there should be before
reaching the ultimately contained items
"""
# 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 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,
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.

View file

@ -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

View file

@ -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 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.
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

12
openmc/temp.py Normal file
View file

@ -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))

View file

@ -7,7 +7,7 @@ import sys
import numpy as np
import openmc
from openmc.checkvalue import check_type, check_length, check_greater_than
import openmc.checkvalue as cv
if sys.version_info[0] >= 3:
basestring = str
@ -111,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
@ -148,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):
@ -432,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):
@ -671,24 +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)
self._universes = np.asarray(universes, dtype=Universe)
cv.check_iterable_type('lattice universes', universes, Universe,
min_depth=2, max_depth=3)
self._universes = universes
def get_unique_universes(self):
"""Determine all unique universes in the lattice
@ -701,13 +702,21 @@ 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]
univs[u._id] = u
else:
for i in range(len(self._universes[k][j])):
u = self._universes[k][j][i]
assert isinstance(u, Universe)
univs[u._id] = u
for universe in unique_universes:
universes[universe._id] = universe
univs[self._outer._id] = self._outer
return universes
return univs
def get_all_nuclides(self):
"""Return all nuclides contained in the lattice
@ -825,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):
@ -1041,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
@ -1091,14 +1100,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 +1115,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:

View file

@ -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

View file

@ -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)

View file

@ -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

View file

@ -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)

View file

@ -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 <upper_right> and <width> 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 <upper_right> and <width> 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 <width> must be the same as the &
&number of entries on <lower_left>.")
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 <upper_right> must be the same &
&as the number of entries on <lower_left>.")
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

View file

@ -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()

View file

@ -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()

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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)
@ -394,23 +394,23 @@ 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.
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, &
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
@ -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
@ -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
@ -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

View file

@ -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)

View file

@ -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 <material>
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 <name> and <xs> 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 <upper_right> and <width> 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 <upper_right> and <width> 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 <width> must be the same as &
&the number of entries on <lower_left>.")
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 <upper_right> must be the &
&same as the number of entries on <lower_left>.")
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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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)

View file

@ -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

View file

@ -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)

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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
!===============================================================================

View file

@ -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

View file

@ -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), &

View file

@ -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

View file

@ -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

View file

@ -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.

View file

@ -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

View file

@ -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

View file

@ -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

View file

@ -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

232
tests/check_source.py Executable file
View file

@ -0,0 +1,232 @@
#!/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 is 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)

View file

@ -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