crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

merged Paul's commits to CMFD branch

Browse source
This commit is contained in:
Bryan Herman 2011-11-06 21:01:58 -05:00
commit 81d91ea28f
55 changed files with 4102 additions and 1330 deletions
Download patch file Download diff file Expand all files Collapse all files

25
benchmarks/pu-met-fast-001/geometry.xml
View file

@ -1,25 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!--
==========================================
Description: Bare Pu-239 Jezebel Benchmark
Case: PU-MET-FAST-001
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<cell uid="1">
<universe>0</universe>
<material>1</material>
<surfaces>-1</surfaces>
</cell>
<surface uid="1">
<type>sphere</type>
<coeffs>0. 0. 0. 6.3849</coeffs>
<boundary>vacuum</boundary>
</surface>
</geometry>

21
benchmarks/pu-met-fast-001/materials.xml
View file

@ -1,21 +0,0 @@
<?xml version="1.0"?>
<materials>
<!--
==========================================
Description: Bare Pu-239 Jezebel Benchmark
Case: PU-MET-FAST-001
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<material uid="1">
<density value="0.04029014" units="atom/b-cm" />
<nuclide name="Pu-239" xs="03c" ao="0.037047" />
<nuclide name="Pu-240" xs="03c" ao="0.0017512" />
<nuclide name="Pu-241" xs="03c" ao="0.00011674" />
<nuclide name="Ga-nat" xs="03c" ao="0.0013752" />
</material>
</materials>

28
benchmarks/pu-met-fast-001/settings.xml
View file

@ -1,28 +0,0 @@
<?xml version="1.0"?>
<settings>
<!--
==========================================
Description: Bare Pu-239 Jezebel Benchmark
Case: PU-MET-FAST-001
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<xslibrary path="/opt/serpent/xsdata/endfb7/xsdata" />
<criticality>
<cycles>3000</cycles>
<inactive>20</inactive>
<particles>10000</particles>
</criticality>
<verbosity value="7" />
<source>
<type>box</type>
<coeffs>-1 -1 -1 1 1 1</coeffs>
</source>
</settings>

13
benchmarks/pu-met-fast-001/tallies.xml
View file

@ -1,13 +0,0 @@
<?xml version="1.0"?>
<tallies>
<tally id="1">
<filters>
<cell>1</cell>
</filters>
<macros>total scatter absorption fission nu-fission</macros>
</tally>
</tallies>

25
benchmarks/pu-met-fast-002/geometry.xml
View file

@ -1,25 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!--
==========================================
Description: Bare Pu-240 Jezebel Benchmark
Case: PU-MET-FAST-002
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<cell uid="1">
<universe>0</universe>
<material>1</material>
<surfaces>-1</surfaces>
</cell>
<surface uid="1">
<type>sphere</type>
<coeffs>0. 0. 0. 6.6595</coeffs>
<boundary>vacuum</boundary>
</surface>
</geometry>

22
benchmarks/pu-met-fast-002/materials.xml
View file

@ -1,22 +0,0 @@
<?xml version="1.0"?>
<materials>
<!--
==========================================
Description: Bare Pu-240 Jezebel Benchmark
Case: PU-MET-FAST-002
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<material uid="1">
<density value="0.04055292" units="atom/b-cm" />
<nuclide name="Pu-239" xs="03c" ao="0.029934" />
<nuclide name="Pu-240" xs="03c" ao="0.0078754" />
<nuclide name="Pu-241" xs="03c" ao="0.0012146" />
<nuclide name="Pu-242" xs="03c" ao="0.00015672" />
<nuclide name="Ga-nat" xs="03c" ao="0.0013722" />
</material>
</materials>

28
benchmarks/pu-met-fast-002/settings.xml
View file

@ -1,28 +0,0 @@
<?xml version="1.0"?>
<settings>
<!--
==========================================
Description: Bare Pu-240 Jezebel Benchmark
Case: PU-MET-FAST-002
Written By: Paul Romano
Date: 9/1/2011
==========================================
-->
<xslibrary path="/opt/serpent/xsdata/endfb7/xsdata" />
<criticality>
<cycles>3000</cycles>
<inactive>20</inactive>
<particles>10000</particles>
</criticality>
<verbosity value="7" />
<source>
<type>box</type>
<coeffs>-1 -1 -1 1 1 1</coeffs>
</source>
</settings>

34
benchmarks/pu-met-fast-005/geometry.xml
View file

@ -1,34 +0,0 @@
<?xml version="1.0"?>
<geometry>
<!--
===============================================================
Description: Unmoderated Plutonium Metal Button Array Benchmark
Case: PU-MET-FAST-005
Written By: Paul Romano
Date: 9/1/2011
===============================================================
-->
<cell uid="1">
<universe>0</universe>
<material>1</material>
<surfaces>-1</surfaces>
</cell>
<cell uid="2">
<universe>0</universe>
<material>2</material>
<surfaces>1 -2</surfaces>
</cell>
<surface uid="1">
<type>sphere</type>
<coeffs>0. 0. 0. 5.0419</coeffs>
</surface>
<surface uid="2">
<type>sphere</type>
<coeffs>0. 0. 0. 9.7409</coeffs>
<boundary>vacuum</boundary>
</surface>
</geometry>

28
benchmarks/pu-met-fast-005/materials.xml
View file

@ -1,28 +0,0 @@
<?xml version="1.0"?>
<materials>
<!--
===============================================================
Description: Unmoderated Plutonium Metal Button Array Benchmark
Case: PU-MET-FAST-005
Written By: Paul Romano
Date: 9/1/2011
===============================================================
-->
<material uid="1">
<density value="0.04070346" units="atom/b-cm" />
<nuclide name="Pu-239" xs="03c" ao="0.037291" />
<nuclide name="Pu-240" xs="03c" ao="0.0019277" />
<nuclide name="Pu-241" xs="03c" ao="0.00012196" />
<nuclide name="Ga-nat" xs="03c" ao="0.0013628" />
</material>
<material uid="2">
<density value="0.06605308" units="atom/b-cm" />
<nuclide name="W-nat" xs="03c" ao="0.051468" />
<nuclide name="Ni-nat" xs="03c" ao="0.0097124" />
<nuclide name="Cu-nat" xs="03c" ao="0.0040774" />
<nuclide name="Zr-nat" xs="03c" ao="0.00079528" />
</material>
</materials>

28
benchmarks/pu-met-fast-005/settings.xml
View file

@ -1,28 +0,0 @@
<?xml version="1.0"?>
<settings>
<!--
===============================================================
Description: Unmoderated Plutonium Metal Button Array Benchmark
Case: PU-MET-FAST-005
Written By: Paul Romano
Date: 9/1/2011
===============================================================
-->
<xslibrary path="/opt/serpent/xsdata/endfb7/xsdata" />
<criticality>
<cycles>3000</cycles>
<inactive>20</inactive>
<particles>10000</particles>
</criticality>
<verbosity value="7" />
<source>
<type>box</type>
<coeffs>-1 -1 -1 1 1 1</coeffs>
</source>
</settings>

2705
cross_sections_serpent.xml Normal file
View file

File diff suppressed because it is too large Load diff

16
examples/basic/geometry.xml
View file

@ -2,16 +2,14 @@
<geometry>
<!-- Definition of Cells -->
<cell uid="1" universe="0" fill="37" surfaces="-2 4 -5" />
<cell uid="100" universe="37" material="40" surfaces="-1" />
<cell uid="101" universe="37" material="41" surfaces="1" />
<cell uid="2" universe="0" material="41" surfaces = "2 -3 4 -5" />
<cell id="1" universe="0" fill="37" surfaces="-2" />
<cell id="100" universe="37" material="40" surfaces="-1" />
<cell id="101" universe="37" material="41" surfaces="1" />
<cell id="2" universe="0" material="41" surfaces = "2 -3" />
<!-- Defition of Surfaces -->
<surface uid="1" type="z-cylinder" coeffs="0 0 7" />
<surface uid="2" type="z-cylinder" coeffs="0 0 9" />
<surface uid="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
<surface uid="4" type="z-plane" coeffs="-50" boundary="vacuum" />
<surface uid="5" type="z-plane" coeffs="50" boundary="vacuum" />
<surface id="1" type="z-cylinder" coeffs="0 0 7" />
<surface id="2" type="z-cylinder" coeffs="0 0 9" />
<surface id="3" type="z-cylinder" coeffs="0 0 11" boundary="vacuum" />
</geometry>

4
examples/basic/materials.xml
View file

@ -2,12 +2,12 @@
<materials>
<!-- Definition of materials -->
<material uid="40">
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" ao="1.0" />
</material>
<material uid="41">
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H-1" xs="03c" ao="2.0" />
<nuclide name="O-16" xs="03c" ao="1.0" />

7
examples/basic/settings.xml
View file

@ -2,7 +2,9 @@
<settings>
<!-- Path for cross section library -->
<xslibrary path="/opt/serpent/xsdata/endfb7/xsdata" />
<cross_sections>
/home/paulromano/openmc/cross_sections_serpent.xml
</cross_sections>
<!-- Parameters for criticality calculation -->
<criticality>
@ -11,9 +13,6 @@
<particles>10000</particles>
</criticality>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<type>box</type>

18
examples/lattice/geometry.xml
View file

@ -1,23 +1,23 @@
<?xml version="1.0"?>
<geometry>
<cell uid="1">
<cell id="1">
<universe>0</universe>
<fill>111</fill>
<surfaces>1 -2 3 -4</surfaces>
</cell>
<cell uid="100">
<cell id="100">
<universe>3</universe>
<material>40</material>
<surfaces>-5</surfaces>
</cell>
<cell uid="101">
<cell id="101">
<universe>3</universe>
<material>41</material>
<surfaces>5</surfaces>
</cell>
<lattice uid="111">
<lattice id="111">
<type>rectangular</type>
<dimension>20 20</dimension>
<origin>-10.0 -10.0</origin>
@ -46,11 +46,11 @@
</universes>
</lattice>
<surface uid="1" type="x-plane" coeffs="-9.9" boundary="vacuum" />
<surface uid="2" type="x-plane" coeffs="9.9" boundary="vacuum" />
<surface uid="3" type="y-plane" coeffs="-9.9" boundary="vacuum" />
<surface uid="4" type="y-plane" coeffs="9.9" boundary="vacuum" />
<surface uid="5">
<surface id="1" type="x-plane" coeffs="-9.9" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="9.9" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-9.9" boundary="vacuum" />
<surface id="4" type="y-plane" coeffs="9.9" boundary="vacuum" />
<surface id="5">
<type>z-cylinder</type>
<coeffs>0.0 0.0 0.4</coeffs>
</surface>

6
examples/lattice/materials.xml
View file

@ -2,16 +2,16 @@
<materials>
<!-- Definition of materials -->
<material uid="40">
<material id="40">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" ao="1.0" />
</material>
<material uid="41">
<material id="41">
<density value="1.0" units="g/cc" />
<nuclide name="H-1" xs="03c" ao="2.0" />
<nuclide name="O-16" xs="03c" ao="1.0" />
<sab uid="H-H2O" xs="10t" />
<sab name="lwe7" xs="00t" />
</material>
</materials>

7
examples/lattice/settings.xml
View file

@ -2,7 +2,9 @@
<settings>
<!-- Path for cross section library -->
<xslibrary path="/opt/serpent/xsdata/endfb7/xsdata" />
<cross_sections>
/home/paulromano/openmc/cross_sections_serpent.xml
</cross_sections>
<!-- Parameters for criticality calculation -->
<criticality>
@ -11,9 +13,6 @@
<particles>10000</particles>
</criticality>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<type>box</type>

14
examples/reflective/geometry.xml
View file

@ -2,18 +2,18 @@
<geometry>
<!-- Definition of Cells -->
<cell uid="1">
<cell id="1">
<universe>0</universe>
<material>1</material>
<surfaces>1 -2 3 -4 5 -6</surfaces>
</cell>
<!-- Defition of Surfaces -->
<surface uid="1" type="x-plane" coeffs="-1" boundary="reflective" />
<surface uid="2" type="x-plane" coeffs="1" boundary="reflective" />
<surface uid="3" type="y-plane" coeffs="-1" boundary="reflective" />
<surface uid="4" type="y-plane" coeffs="1" boundary="reflective" />
<surface uid="5" type="z-plane" coeffs="-1" boundary="reflective" />
<surface uid="6" type="z-plane" coeffs="1" boundary="reflective" />
<surface id="1" type="x-plane" coeffs="-1" boundary="vacuum" />
<surface id="2" type="x-plane" coeffs="1" boundary="vacuum" />
<surface id="3" type="y-plane" coeffs="-1" boundary="reflective" />
<surface id="4" type="y-plane" coeffs="1" boundary="reflective" />
<surface id="5" type="z-plane" coeffs="-1" boundary="reflective" />
<surface id="6" type="z-plane" coeffs="1" boundary="reflective" />
</geometry>

2
examples/reflective/materials.xml
View file

@ -2,7 +2,7 @@
<materials>
<!-- Definition of materials -->
<material uid="1">
<material id="1">
<density value="4.5" units="g/cc" />
<nuclide name="U-235" xs="03c" ao="1.0" />
</material>

7
examples/reflective/settings.xml
View file

@ -2,7 +2,9 @@
<settings>
<!-- Path for cross section library -->
<xslibrary path="/opt/serpent/xsdata/sss_endfb7u.xsdata" />
<cross_sections>
/home/paulromano/openmc/cross_sections_serpent.xml
</cross_sections>
<!-- Parameters for criticality calculation -->
<criticality>
@ -11,9 +13,6 @@
<particles>10000</particles>
</criticality>
<!-- How verbose output should be -->
<verbosity value="7" />
<!-- Starting source -->
<source>
<type>box</type>

9
src/DEPENDENCIES
View file

@ -13,6 +13,8 @@ cross_section.o: material_header.o
cross_section.o: output.o
cross_section.o: string.o
cross_section_header.o: constants.o
datatypes.o: datatypes_header.o
doppler.o: constants.o
@ -35,6 +37,7 @@ fileio.o: string.o
fission.o: constants.o
fission.o: cross_section_header.o
fission.o: error.o
fission.o: global.o
fission.o: interpolation.o
fission.o: search.o
@ -46,6 +49,7 @@ geometry.o: global.o
geometry.o: output.o
geometry.o: particle_header.o
geometry.o: string.o
geometry.o: tally.o
global.o: bank_header.o
global.o: constants.o
@ -87,6 +91,7 @@ input_xml.o: mesh_header.o
input_xml.o: output.o
input_xml.o: string.o
input_xml.o: tally_header.o
input_xml.o: xml-fortran/templates/cross_sections_t.o
input_xml.o: xml-fortran/templates/geometry_t.o
input_xml.o: xml-fortran/templates/materials_t.o
input_xml.o: xml-fortran/templates/settings_t.o
@ -94,6 +99,7 @@ input_xml.o: xml-fortran/templates/tallies_t.o
interpolation.o: constants.o
interpolation.o: error.o
interpolation.o: global.o
interpolation.o: search.o
logging.o: constants.o
@ -160,6 +166,7 @@ plot.o: particle_header.o
search.o: constants.o
search.o: error.o
search.o: global.o
source.o: bank_header.o
source.o: constants.o
@ -172,9 +179,9 @@ source.o: physics.o
string.o: constants.o
string.o: error.o
string.o: global.o
tally.o: constants.o
tally.o: cross_section.o
tally.o: error.o
tally.o: global.o
tally.o: mesh.o

2
src/bank_header.f90
View file

@ -9,7 +9,7 @@ module bank_header
!===============================================================================
type Bank
integer(8) :: uid ! Unique ID
integer(8) :: id ! Unique ID
real(8) :: xyz(3) ! location of bank particle
real(8) :: uvw(3) ! diretional cosines
real(8) :: E ! energy

7
src/constants.f90
View file

@ -166,6 +166,12 @@ module constants
N_PT = 116, &
N_DA = 117
! ACE table types
integer, parameter :: &
ACE_NEUTRON = 1, & ! continuous-energy neutron
ACE_THERMAL = 2, & ! thermal S(a,b) scattering data
ACE_DOSIMETRY = 3 ! dosimetry cross sections
! Tally macro reactions
integer, parameter :: N_MACRO_TYPES = 15
integer, parameter :: &
@ -223,6 +229,7 @@ module constants
integer, parameter :: MAX_WORDS = 500
integer, parameter :: MAX_LINE_LEN = 250
integer, parameter :: MAX_WORD_LEN = 150
integer, parameter :: MAX_FILE_LEN = 255
! Unit numbers
integer, parameter :: UNIT_LOG = 11 ! unit # for writing log file

265
src/cross_section.f90
View file

@ -1,7 +1,7 @@
module cross_section
use constants
use cross_section_header, only: Nuclide, Reaction, SAB_Table, xsData
use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing
use datatypes, only: dict_create, dict_add_key, dict_get_key, &
dict_has_key, dict_delete
use datatypes_header, only: DictionaryCI
@ -10,7 +10,7 @@ module cross_section
use fileio, only: read_line, read_data, skip_lines
use global
use material_header, only: Material
use output, only: message
use output, only: write_message
use string, only: split_string, str_to_int, str_to_real, &
lower_case, int_to_str
@ -36,11 +36,10 @@ contains
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
integer :: index ! index in xsdatas array
integer :: index ! index in xs_listings array
integer :: index_nuclides ! index in nuclides
integer :: index_sab ! index in sab_tables
character(10) :: key ! name of isotope, e.g. 92235.03c
character(MAX_LINE_LEN) :: msg ! output/error message
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(SAB_Table), pointer :: sab => null()
@ -106,8 +105,8 @@ contains
do i = 1, n_materials
mat => materials(i)
do j = 1, mat % n_nuclides
! Get index in xsdatas array for this nuclide
index = mat % xsdata(j)
! Get index in xs_listings array for this nuclide
index = mat % xs_listing(j)
! Get name of nuclide
key = mat % names(j)
@ -127,13 +126,13 @@ contains
key = mat % sab_name
if (.not. dict_has_key(temp_dict, key)) then
! Find the entry in xsdatas for this table
if (dict_has_key(xsdata_dict, key)) then
index = dict_get_key(xsdata_dict, key)
! Find the entry in xs_listings for this table
if (dict_has_key(xs_listing_dict, key)) then
index = dict_get_key(xs_listing_dict, key)
else
msg = "Cannot find cross-section " // trim(key) // " in specified &
&xsdata file."
call fatal_error(msg)
message = "Cannot find cross-section " // trim(key) // &
" in specified cross_sections.xml file."
call fatal_error()
end if
! Read the table and add entry to dictionary
@ -156,9 +155,9 @@ contains
! Check to make sure S(a,b) table matched a nuclide
if (mat % sab_nuclide == 0) then
msg = "S(a,b) table " // trim(mat % sab_name) // " did not match " &
// "any nuclide on material " // trim(int_to_str(mat % uid))
call fatal_error(msg)
message = "S(a,b) table " // trim(mat % sab_name) // " did not match " &
// "any nuclide on material " // trim(int_to_str(mat % id))
call fatal_error()
end if
end if
end do
@ -177,7 +176,7 @@ contains
subroutine read_ACE_continuous(index_table, index)
integer, intent(in) :: index_table ! index in nuclides array
integer, intent(in) :: index ! index in xsdatas array
integer, intent(in) :: index ! index in xs_listings array
integer :: in = 7 ! unit to read from
integer :: ioError ! error status for file access
@ -188,56 +187,57 @@ contains
logical :: file_exists ! does ACE library exist?
logical :: found_xs ! did we find table in library?
character(7) :: readable ! is ACE library readable?
character(MAX_LINE_LEN) :: msg ! output/error message
character(MAX_LINE_LEN) :: line ! single line to read
character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line
character(MAX_WORD_LEN) :: filename ! name of ACE library file
character(MAX_FILE_LEN) :: filename ! name of ACE library file
character(10) :: tablename ! name of cross section table
type(Nuclide), pointer :: nuc => null()
! Check to make sure index in nuclides array and xsdata arrays are valid
! Check to make sure index in nuclides array and xs_listings arrays are
! valid
if (index_table > size(nuclides)) then
msg = "Index of table to read is greater than length of nuclides."
call fatal_error(msg)
elseif (index > size(xsdatas)) then
msg = "Index of xsdata entry is greater than length of xsdatas."
call fatal_error(msg)
message = "Index of table to read is greater than length of nuclides."
call fatal_error()
elseif (index > size(xs_listings)) then
message = "Index of xs_listing entry is greater than length of " // &
"xs_listings."
call fatal_error()
end if
filename = xsdatas(index)%path
tablename = xsdatas(index)%id
filename = xs_listings(index) % path
tablename = xs_listings(index) % name
nuc => nuclides(index_table)
! Check if input file exists and is readable
inquire(FILE=filename, EXIST=file_exists, READ=readable)
if (.not. file_exists) then
msg = "ACE library '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "ACE library '" // trim(filename) // "' does not exist!"
call fatal_error()
elseif (readable(1:3) == 'NO') then
msg = "ACE library '" // trim(filename) // "' is not readable! &
message = "ACE library '" // trim(filename) // "' is not readable! &
&Change file permissions with chmod command."
call fatal_error(msg)
call fatal_error()
end if
! display message
msg = "Loading ACE cross section table: " // tablename
call message(msg, 6)
message = "Loading ACE cross section table: " // tablename
call write_message(6)
! open file
open(file=filename, unit=in, status='old', &
& action='read', iostat=ioError)
if (ioError /= 0) then
msg = "Error while opening file: " // filename
call fatal_error(msg)
message = "Error while opening file: " // filename
call fatal_error()
end if
found_xs = .false.
do while (.not. found_xs)
call read_line(in, line, ioError)
if (ioError < 0) then
msg = "Could not find ACE table " // tablename // "."
call fatal_error(msg)
message = "Could not find ACE table " // tablename // "."
call fatal_error()
end if
call split_string(line, words, n)
if (trim(words(1)) == trim(tablename)) then
@ -1073,7 +1073,7 @@ contains
subroutine read_ACE_thermal(index_table, index)
integer, intent(in) :: index_table ! index in sab_tables array
integer, intent(in) :: index ! index in xsdatas array
integer, intent(in) :: index ! index in xs_listings array
integer :: in = 7 ! unit to read from
integer :: ioError ! error status for file access
@ -1084,47 +1084,46 @@ contains
logical :: file_exists ! does ACE library exist?
logical :: found_xs ! did we find table in library?
character(7) :: readable ! is ACE library readable?
character(MAX_LINE_LEN) :: msg ! output/error message
character(MAX_LINE_LEN) :: line ! single line to read
character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line
character(MAX_WORD_LEN) :: filename ! name of ACE library file
character(10) :: tablename ! name of cross section table
type(SAB_Table), pointer :: table => null()
filename = xsdatas(index)%path
tablename = xsdatas(index)%id
filename = xs_listings(index) % path
tablename = xs_listings(index) % name
table => sab_tables(index_table)
! Check if input file exists and is readable
inquire(FILE=filename, EXIST=file_exists, READ=readable)
if (.not. file_exists) then
msg = "ACE library '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "ACE library '" // trim(filename) // "' does not exist!"
call fatal_error()
elseif (readable(1:3) == 'NO') then
msg = "ACE library '" // trim(filename) // "' is not readable! &
message = "ACE library '" // trim(filename) // "' is not readable! &
&Change file permissions with chmod command."
call fatal_error(msg)
call fatal_error()
end if
! display message
msg = "Loading ACE cross section table: " // tablename
call message(msg, 6)
message = "Loading ACE cross section table: " // tablename
call write_message(6)
! open file
open(file=filename, unit=in, status='old', &
& action='read', iostat=ioError)
if (ioError /= 0) then
msg = "Error while opening file: " // filename
call fatal_error(msg)
message = "Error while opening file: " // filename
call fatal_error()
end if
found_xs = .false.
do while (.not. found_xs)
call read_line(in, line, ioError)
if (ioError < 0) then
msg = "Could not find ACE table " // tablename // "."
call fatal_error(msg)
message = "Could not find ACE table " // tablename // "."
call fatal_error()
end if
call split_string(line, words, n)
if (trim(words(1)) == trim(tablename)) then
@ -1226,6 +1225,8 @@ contains
! scattering
NE_out = NXS(4)
NMU = NXS(3) + 1
table % n_inelastic_e_out = NE_out
table % n_inelastic_mu = NMU
allocate(table % inelastic_e_out(NE_out, NE_in))
allocate(table % inelastic_mu(NMU, NE_out, NE_in))
@ -1320,168 +1321,4 @@ contains
end function get_real
!===============================================================================
! GET_MACRO_XS
!===============================================================================
function get_macro_xs(p, mat, MT) result(xs)
type(Particle), pointer :: p
type(Material), pointer :: mat
integer, intent(in) :: MT
real(8) :: xs
integer :: i, j
integer :: n_nuclides
integer :: IE
real(8) :: density_i
real(8) :: sigma_i
real(8) :: f
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
! initialize xs
xs = ZERO
! loop over all nuclides in material
n_nuclides = mat % n_nuclides
do i = 1, n_nuclides
nuc => nuclides(mat % nuclide(i))
! determine nuclide atom density
density_i = mat % atom_density(i)
! search nuclide energy grid
IE = nuc%grid_index(p % IE)
f = (p%E - nuc%energy(IE))/(nuc%energy(IE+1) - nuc%energy(IE))
! handle special case of total cross section
if (MT == 1) then
xs = xs + mat % density * (ONE-f) * nuc%total(IE) + &
& f * (nuc%total(IE+1))
cycle
end if
! loop over reactions in isotope
do j = 1, nuc % n_reaction
rxn => nuc % reactions(i)
! check for matching MT
if (MT /= rxn % MT) cycle
! if energy is below threshold for this reaction, skip it
if (IE < rxn % IE) cycle
! add to cumulative probability
sigma_i = (ONE-f) * rxn%sigma(IE-rxn%IE+1) + &
& f * (rxn%sigma(IE-rxn%IE+2))
end do
! calculate nuclide macroscopic cross-section
xs = xs + density_i * sigma_i
end do
end function get_macro_xs
!===============================================================================
! READ_XSDATA reads the data in a SERPENT xsdata file and builds a dictionary to
! find cross-section information later on.
!===============================================================================
subroutine read_xsdata(path)
character(*), intent(in) :: path
type(xsData), pointer :: iso => null()
character(MAX_LINE_LEN) :: line
character(MAX_LINE_LEN) :: msg
character(MAX_WORD_LEN) :: words(MAX_WORDS)
character(MAX_WORD_LEN) :: filename
integer :: n
integer :: in = 7
logical :: file_exists
character(7) :: readable
integer :: count
integer :: index
integer :: ioError
msg = "Reading cross-section summary file..."
call message(msg, 5)
! Construct filename
filename = trim(path)
! Check if xsdata exists and is readable
inquire(FILE=filename, EXIST=file_exists, READ=readable)
if (.not. file_exists) then
msg = "Cross section summary '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
elseif (readable(1:3) == 'NO') then
msg = "Cross section summary '" // trim(filename) // "' is not readable!" &
& // "Change file permissions with chmod command."
call fatal_error(msg)
end if
! open xsdata file
open(FILE=filename, UNIT=in, STATUS='old', &
& ACTION='read', IOSTAT=ioError)
if (ioError /= 0) then
msg = "Error while opening file: " // filename
call fatal_error(msg)
end if
! determine how many lines
count = 0
do
read(UNIT=in, FMT='(A)', IOSTAT=ioError) line
if (ioError < 0) then
! reached end of file
exit
elseif (ioError > 0) then
msg = "Unknown error while reading file: " // filename
close(UNIT=in)
call fatal_error(msg)
end if
count = count + 1
end do
allocate(xsdatas(count))
! read actual lines
index = 0
rewind(in)
do
read(UNIT=in, FMT='(A)', IOSTAT=ioError) line
if (ioError < 0) exit
index = index + 1
call split_string(line, words, n)
if (n == 0) cycle ! skip blank line
! Check to make sure there are enough arguments
if (n < 9) then
msg = "Not enough arguments on xsdata line: " // line
close(UNIT=in)
call fatal_error(msg)
end if
iso => xsdatas(index)
! store data
iso%alias = words(1)
iso%id = words(2)
iso%type = str_to_int(words(3))
iso%zaid = str_to_int(words(4))
iso%isomeric = str_to_int(words(5))
iso%awr = str_to_real(words(6))
iso%temp = str_to_real(words(7))
iso%binary = str_to_int(words(8))
iso%path = words(9)
! create dictionary entry
call dict_add_key(xsdata_dict, iso%alias, index)
end do
close(UNIT=in)
end subroutine read_xsdata
end module cross_section

16
src/cross_section_header.f90
View file

@ -1,5 +1,7 @@
module cross_section_header
use constants, only: MAX_FILE_LEN
implicit none
!===============================================================================
@ -148,20 +150,20 @@ module cross_section_header
end type SAB_Table
!===============================================================================
! XSDATA contains data read in from a SERPENT xsdata file
! XSLISTING contains data read from a cross_sections.xml file
!===============================================================================
type xsData
type XsListing
character(10) :: name
character(10) :: alias
character(10) :: id
integer :: type
integer :: zaid
integer :: isomeric
real(8) :: awr
real(8) :: temp
integer :: binary
character(150) :: path
end type xsData
logical :: metastable
logical :: binary
character(MAX_FILE_LEN) :: path
end type XsListing
!===============================================================================
! NUCLIDEMICROXS contains cached microscopic cross sections for a

6
src/datatypes.f90
View file

@ -16,9 +16,9 @@ module datatypes
implicit none
integer, parameter :: HASH_SIZE = 4993
integer, parameter :: HASH_MULTIPLIER = 31
integer, parameter :: DICT_NULL = 0
integer, parameter, private :: HASH_SIZE = 4993
integer, parameter, private :: HASH_MULTIPLIER = 31
integer, parameter, private :: DICT_NULL = 0
!===============================================================================
! LIST Interfaces -- these allow one to use a single subroutine or function for

7
src/energy_grid.f90
View file

@ -4,7 +4,7 @@ module energy_grid
use datatypes, only: list_insert, list_size, list_delete
use datatypes_header, only: ListReal
use global
use output, only: message
use output, only: write_message
contains
@ -23,10 +23,9 @@ contains
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
integer :: i, j
character(MAX_LINE_LEN) :: msg
msg = "Creating unionized energy grid..."
call message(msg, 5)
message = "Creating unionized energy grid..."
call write_message(5)
! loop over all materials
do i = 1, n_materials

30
src/error.f90
View file

@ -2,7 +2,7 @@ module error
use ISO_FORTRAN_ENV
use global, only: master, free_memory
use global, only: master, free_memory, message
implicit none
@ -17,24 +17,22 @@ contains
! stream.
!===============================================================================
subroutine warning(msg)
subroutine warning()
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
integer :: n_lines ! number of lines
integer :: i ! loop index for lines
! Only allow master to print to screen
if (.not. master) return
write(ou, fmt='(1X,A9)', advance='no') 'WARNING: '
n_lines = (len_trim(msg)-1)/70 + 1
n_lines = (len_trim(message)-1)/70 + 1
do i = 1, n_lines
if (i == 1) then
write(ou, fmt='(A70)') msg(70*(i-1)+1:70*i)
write(ou, fmt='(A70)') message(70*(i-1)+1:70*i)
else
write(ou, fmt='(10X,A70)') msg(70*(i-1)+1:70*i)
write(ou, fmt='(10X,A70)') message(70*(i-1)+1:70*i)
end if
end do
@ -46,23 +44,21 @@ contains
! the program is aborted.
!===============================================================================
subroutine fatal_error(msg)
subroutine fatal_error()
character(*), intent(in) :: msg
integer :: n_lines
integer :: i
integer :: n_lines ! number of lines
integer :: i ! loop index over lines
! Only allow master to print to screen
if (master) then
write(eu, fmt='(1X,A7)', advance='no') 'ERROR: '
n_lines = (len_trim(msg)-1)/72 + 1
n_lines = (len_trim(message)-1)/72 + 1
do i = 1, n_lines
if (i == 1) then
write(eu, fmt='(A72)') msg(72*(i-1)+1:72*i)
write(eu, fmt='(A72)') message(72*(i-1)+1:72*i)
else
write(eu, fmt='(7X,A72)') msg(72*(i-1)+1:72*i)
write(eu, fmt='(7X,A72)') message(72*(i-1)+1:72*i)
end if
end do
write(eu,*)

12
src/fileio.f90
View file

@ -89,8 +89,8 @@ contains
integer, intent(in) :: n_lines ! number of lines to skip
integer, intent(out) :: ioError ! error status
integer :: i ! index for number of lines
character(MAX_LINE_LEN) :: tmp ! single line
integer :: i ! index for number of lines
character(MAX_LINE_LEN) :: tmp ! single line
do i = 1, n_lines
read(UNIT=unit, FMT='(A)', IOSTAT=ioError) tmp
@ -110,8 +110,8 @@ contains
integer, intent(in) :: lines ! total number of lines
integer, intent(in) :: words_per_line ! number of words per line
integer :: i ! line index
integer :: loc ! locator for array
integer :: i ! line index
integer :: loc ! locator for array
loc = 0
do i = 1, lines
@ -137,8 +137,8 @@ contains
integer, intent(in) :: lines ! total number of lines
integer, intent(in) :: words_per_line ! number of words per line
integer :: i ! line index
integer :: loc ! locator for array
integer :: i ! line index
integer :: loc ! locator for array
loc = 0
do i = 1, lines

18
src/fission.f90
View file

@ -3,6 +3,7 @@ module fission
use constants
use cross_section_header, only: Nuclide
use error, only: fatal_error
use global, only: message
use interpolation, only: interpolate_tab1
use search, only: binary_search
@ -21,14 +22,13 @@ contains
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of total neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
character(MAX_LINE_LEN) :: msg ! error message
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_t_type == NU_NONE) then
msg = "No neutron emission data for table: " // nuc % name
call fatal_error(msg)
message = "No neutron emission data for table: " // nuc % name
call fatal_error()
elseif (nuc % nu_t_type == NU_POLYNOMIAL) then
! determine number of coefficients
NC = int(nuc % nu_t_data(1))
@ -57,9 +57,9 @@ contains
real(8), intent(in) :: E ! energy of incoming neutron
real(8) :: nu ! number of prompt neutrons emitted per fission
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
integer :: i ! loop index
integer :: NC ! number of polynomial coefficients
real(8) :: c ! polynomial coefficient
if (nuc % nu_p_type == NU_NONE) then
! since no prompt or delayed data is present, this means all neutron

125
src/geometry.f90
View file

@ -5,9 +5,10 @@ module geometry
use error, only: fatal_error
use geometry_header, only: Cell, Surface, Universe, Lattice
use global
use output, only: message
use output, only: write_message
use particle_header, only: Particle
use string, only: int_to_str
use tally, only: score_surface_current
implicit none
@ -95,14 +96,12 @@ contains
type(Particle), pointer :: p ! pointer to particle
logical, intent(inout) :: found ! particle found?
character(MAX_LINE_LEN) :: msg ! error message
integer :: i ! index over cells
integer :: x, y
type(Cell), pointer :: c ! pointer to cell
type(Lattice), pointer :: lat ! pointer to lattice
type(Universe), pointer :: lower_univ ! if particle is in lower
! universe, use this pointer
! to call recursively
type(Universe), pointer :: lower_univ ! if particle is in lower universe,
! use this pointer to call recursively
found = .false.
@ -118,7 +117,7 @@ contains
! set particle attributes
p % cell = univ % cells(i)
p % universe = dict_get_key(universe_dict, univ % uid)
p % universe = dict_get_key(universe_dict, univ % id)
p % material = c % material
exit
elseif (c % type == CELL_FILL) then
@ -127,8 +126,8 @@ contains
if (found) then
exit
else
msg = "Could not locate particle in universe: "
call fatal_error(msg)
message = "Could not locate particle in universe: "
call fatal_error()
end if
elseif (c % type == CELL_LATTICE) then
! Set current lattice
@ -153,9 +152,9 @@ contains
if (found) then
exit
else
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call fatal_error(msg)
message = "Could not locate particle in lattice: " &
& // int_to_str(lat % id)
call fatal_error()
end if
end if
end if
@ -191,7 +190,6 @@ contains
real(8) :: dot_prod ! dot product of direction and normal
real(8) :: norm ! "norm" of surface normal
logical :: found ! particle found in universe?
character(MAX_LINE_LEN) :: msg ! output/error message?
type(Surface), pointer :: surf => null()
type(Cell), pointer :: c => null()
type(Lattice), pointer :: lat => null()
@ -199,18 +197,31 @@ contains
surf => surfaces(abs(p % surface))
if (verbosity >= 10) then
msg = " Crossing surface " // trim(int_to_str(surf % uid))
call message(msg)
message = " Crossing surface " // trim(int_to_str(surf % id))
call write_message()
end if
if (surf % bc == BC_VACUUM .and. (.not. plotting)) then
! =======================================================================
! PARTICLE LEAKS OUT OF PROBLEM
! Kill particle
p % alive = .false.
! Score any surface current tallies -- note that the particle is moved
! forward slightly so that if the mesh boundary is on the surface, it is
! still processed
! TODO: Find a better solution to score surface currents than physically
! moving the particle forward slightly
p % xyz = p % xyz + 1e-6 * p % uvw
call score_surface_current(p)
! Display message
if (verbosity >= 10) then
msg = " Leaked out of surface " // trim(int_to_str(surf % uid))
call message(msg)
message = " Leaked out of surface " // trim(int_to_str(surf % id))
call write_message()
end if
return
@ -247,9 +258,9 @@ contains
p % uvw = (/ u, v, w /)
case (SURF_CYL_X)
! Find y-y0, z-z0 and dot product of direction and surface normal
y = p % xyz(2) - surf % coeffs(2)
z = p % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
y = p % xyz(2) - surf % coeffs(1)
z = p % xyz(3) - surf % coeffs(2)
R = surf % coeffs(3)
dot_prod = v*y + w*z
! Reflect direction according to normal
@ -261,8 +272,8 @@ contains
case (SURF_CYL_Y)
! Find x-x0, z-z0 and dot product of direction and surface normal
x = p % xyz(1) - surf % coeffs(1)
z = p % xyz(3) - surf % coeffs(3)
R = surf % coeffs(4)
z = p % xyz(3) - surf % coeffs(2)
R = surf % coeffs(3)
dot_prod = u*x + w*z
! Reflect direction according to normal
@ -275,7 +286,7 @@ contains
! Find x-x0, y-y0 and dot product of direction and surface normal
x = p % xyz(1) - surf % coeffs(1)
y = p % xyz(2) - surf % coeffs(2)
R = surf % coeffs(4)
R = surf % coeffs(3)
dot_prod = u*x + v*y
! Reflect direction according to normal
@ -301,19 +312,26 @@ contains
! Set vector
p % uvw = (/ u, v, w /)
case default
msg = "Reflection not supported for surface " // &
trim(int_to_str(surf % uid))
call fatal_error(msg)
message = "Reflection not supported for surface " // &
trim(int_to_str(surf % id))
call fatal_error()
end select
! Reassign particle's cell and surface
p % cell = last_cell
p % surface = -p % surface
! Score surface currents since reflection causes the direction of the
! particle to change
call score_surface_current(p)
! Set previous coordinate going slightly past surface crossing
p % last_xyz = p % xyz + 1e-6 * p % uvw
! Diagnostic message
if (verbosity >= 10) then
msg = " Reflected from surface " // trim(int_to_str(surf%uid))
call message(msg)
message = " Reflected from surface " // trim(int_to_str(surf%id))
call write_message()
end if
return
end if
@ -332,8 +350,8 @@ contains
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in universe: "
call fatal_error(msg)
message = "Could not locate particle in universe: "
call fatal_error()
end if
elseif (c % type == CELL_LATTICE) then
! Set current lattice
@ -356,9 +374,9 @@ contains
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call fatal_error(msg)
message = "Could not locate particle in lattice: " // &
trim(int_to_str(lat % id))
call fatal_error()
end if
else
! set current pointers
@ -379,8 +397,8 @@ contains
lower_univ => universes(c % fill)
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in universe: "
call fatal_error(msg)
message = "Could not locate particle in universe: "
call fatal_error()
end if
elseif (c % type == CELL_LATTICE) then
! Set current lattice
@ -403,9 +421,9 @@ contains
call find_cell(lower_univ, p, found)
if (.not. found) then
msg = "Could not locate particle in lattice: " &
& // int_to_str(lat % uid)
call fatal_error(msg)
message = "Could not locate particle in lattice: " // &
trim(int_to_str(lat % id))
call fatal_error()
end if
else
! set current pointers
@ -430,9 +448,9 @@ contains
end do
! Couldn't find next cell anywhere!
msg = "After particle crossed surface " // trim(int_to_str(p%surface)) // &
& ", it could not be located in any cell and it did not leak."
call fatal_error(msg)
message = "After particle crossed surface " // trim(int_to_str(p%surface)) &
// ", it could not be located in any cell and it did not leak."
call fatal_error()
end subroutine cross_surface
@ -459,13 +477,12 @@ contains
real(8) :: x0 ! half the width of lattice element
real(8) :: y0 ! half the height of lattice element
logical :: found ! particle found in cell?
character(MAX_LINE_LEN) :: msg ! output/error message
type(Lattice), pointer :: lat
type(Universe), pointer :: univ
if (verbosity >= 10) then
msg = " Crossing lattice"
call message(msg)
message = " Crossing lattice"
call write_message()
end if
lat => lattices(p % lattice)
@ -530,11 +547,11 @@ contains
i_x = p % index_x
i_y = p % index_y
if (i_x < 1 .or. i_x > lat % n_x) then
msg = "Reached edge of lattice."
call fatal_error(msg)
message = "Reached edge of lattice."
call fatal_error()
elseif (i_y < 1 .or. i_y > lat % n_y) then
msg = "Reached edge of lattice."
call fatal_error(msg)
message = "Reached edge of lattice."
call fatal_error()
end if
! Find universe for next lattice element
@ -543,8 +560,8 @@ contains
! Find cell in next lattice element
call find_cell(univ, p, found)
if (.not. found) then
msg = "Could not locate particle in universe: "
call fatal_error(msg)
message = "Could not locate particle in universe: "
call fatal_error()
end if
end subroutine cross_lattice
@ -578,7 +595,6 @@ contains
real(8) :: a,b,c,k ! quadratic equation coefficients
real(8) :: quad ! discriminant of quadratic equation
logical :: on_surface ! is particle on surface?
character(MAX_LINE_LEN) :: msg ! output/error message
type(Cell), pointer :: cell_p => null()
type(Cell), pointer :: parent_p => null()
type(Surface), pointer :: surf_p => null()
@ -882,8 +898,8 @@ contains
end if
case (SURF_GQ)
msg = "Surface distance not yet implement for general quadratic."
call fatal_error(msg)
message = "Surface distance not yet implement for general quadratic."
call fatal_error()
end select
@ -1114,12 +1130,11 @@ contains
integer, allocatable :: count_positive(:) ! # of cells on positive side
integer, allocatable :: count_negative(:) ! # of cells on negative side
logical :: positive ! positive side specified in surface list
character(MAX_LINE_LEN) :: msg ! output/error message
type(Cell), pointer :: c
type(Surface), pointer :: surf
msg = "Building neighboring cells lists for each surface..."
call message(msg, 4)
message = "Building neighboring cells lists for each surface..."
call write_message(4)
allocate(count_positive(n_surfaces))
allocate(count_negative(n_surfaces))

11
src/geometry_header.f90
View file

@ -7,7 +7,7 @@ module geometry_header
!===============================================================================
type Universe
integer :: uid ! Unique ID
integer :: id ! Unique ID
integer :: type ! Type
integer :: level ! Level of universe (0=base)
integer :: n_cells ! # of cells within
@ -15,7 +15,6 @@ module geometry_header
real(8) :: x0 ! Translation in x-coordinate
real(8) :: y0 ! Translation in y-coordinate
real(8) :: z0 ! Translation in z-coordinate
integer, allocatable :: tallies(:)
end type Universe
!===============================================================================
@ -24,7 +23,7 @@ module geometry_header
!===============================================================================
type Lattice
integer :: uid ! Universe number for lattice
integer :: id ! Universe number for lattice
integer :: type ! Type of lattice (rectangular, hex, etc)
integer :: level ! Level of lattice
integer :: n_x ! number of lattice cells in x-direction
@ -34,7 +33,6 @@ module geometry_header
real(8) :: width_x ! width of lattice cell
real(8) :: width_y ! width of lattice cell
integer, allocatable :: element(:,:) ! specified universes
integer, allocatable :: tallies(:)
end type Lattice
!===============================================================================
@ -43,7 +41,7 @@ module geometry_header
!===============================================================================
type Surface
integer :: uid ! Unique ID
integer :: id ! Unique ID
integer :: type ! Type of surface
real(8), allocatable :: coeffs(:) ! Definition of surface
integer, allocatable :: &
@ -57,7 +55,7 @@ module geometry_header
!===============================================================================
type Cell
integer :: uid ! Unique ID
integer :: id ! Unique ID
integer :: type ! Type of cell (normal, universe, lattice)
integer :: universe ! universe # this cell is in
integer :: fill ! universe # filling this cell
@ -68,7 +66,6 @@ module geometry_header
& surfaces(:) ! List of surfaces bounding cell -- note that
! parentheses, union, etc operators will be listed
! here too
integer, allocatable :: tallies(:)
end type Cell
! array index of universe 0

166
src/global.f90
View file

@ -2,8 +2,8 @@ module global
use bank_header, only: Bank
use constants
use cross_section_header, only: Nuclide, SAB_Table, xsData, NuclideMicroXS, &
MaterialMacroXS
use cross_section_header, only: Nuclide, SAB_Table, xsListing, &
NuclideMicroXS, MaterialMacroXS
use datatypes_header, only: DictionaryII, DictionaryCI
use geometry_header, only: Cell, Universe, Lattice, Surface
use material_header, only: Material
@ -20,60 +20,79 @@ module global
implicit none
save
! ============================================================================
! GEOMETRY-RELATED VARIABLES
! Main arrays
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
type(xsData), allocatable, target :: xsdatas(:)
type(StructuredMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
type(Cell), allocatable, target :: cells(:)
type(Universe), allocatable, target :: universes(:)
type(Lattice), allocatable, target :: lattices(:)
type(Surface), allocatable, target :: surfaces(:)
type(Material), allocatable, target :: materials(:)
! Size of main arrays
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
integer :: n_meshes ! # of structured meshes
integer :: n_tallies ! # of tallies
integer :: n_cells ! # of cells
integer :: n_universes ! # of universes
integer :: n_lattices ! # of lattices
integer :: n_surfaces ! # of surfaces
integer :: n_materials ! # of materials
! These dictionaries provide a fast lookup mechanism -- the key is the
! user-specified identifier and the value is the index in the corresponding
! array
type(DictionaryII), pointer :: cell_dict
type(DictionaryII), pointer :: universe_dict
type(DictionaryII), pointer :: lattice_dict
type(DictionaryII), pointer :: surface_dict
type(DictionaryII), pointer :: material_dict
type(DictionaryII), pointer :: mesh_dict
type(DictionaryII), pointer :: tally_dict
type(DictionaryCI), pointer :: xsdata_dict
type(DictionaryCI), pointer :: nuclide_dict
type(DictionaryCI), pointer :: sab_dict
type(DictionaryII), pointer :: cell_dict => null()
type(DictionaryII), pointer :: universe_dict => null()
type(DictionaryII), pointer :: lattice_dict => null()
type(DictionaryII), pointer :: surface_dict => null()
type(DictionaryII), pointer :: material_dict => null()
type(DictionaryII), pointer :: mesh_dict => null()
type(DictionaryII), pointer :: tally_dict => null()
! ============================================================================
! CROSS SECTION RELATED VARIABLES
! Cross section arrays
type(Nuclide), allocatable, target :: nuclides(:)
type(SAB_Table), allocatable, target :: sab_tables(:)
integer :: n_nuclides_total
integer :: n_sab_tables
type(Nuclide), allocatable, target :: nuclides(:) ! Nuclide cross-sections
type(SAB_Table), allocatable, target :: sab_tables(:) ! S(a,b) tables
type(XsListing), allocatable, target :: xs_listings(:) ! cross_sections.xml listings
! Cross section caches
type(NuclideMicroXS), allocatable :: micro_xs(:)
type(MaterialMacroXS) :: material_xs
type(NuclideMicroXS), allocatable :: micro_xs(:) ! Cache for each nuclide
type(MaterialMacroXS) :: material_xs ! Cache for current material
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
integer :: n_nuclides_total ! Number of nuclide cross section tables
integer :: n_sab_tables ! Number of S(a,b) thermal scattering tables
! Dictionaries to look up cross sections and listings
type(DictionaryCI), pointer :: nuclide_dict => null()
type(DictionaryCI), pointer :: sab_dict => null()
type(DictionaryCI), pointer :: xs_listing_dict => null()
! Unionized energy grid
integer :: n_grid ! number of points on unionized grid
real(8), allocatable :: e_grid(:) ! energies on unionized grid
! Histories/cycles/etc for both external source and criticality
integer(8) :: n_particles ! # of particles (per cycle for criticality)
integer :: n_cycles ! # of cycles
integer :: n_inactive ! # of inactive cycles
! ============================================================================
! TALLY-RELATED VARIABLES
type(StructuredMesh), allocatable, target :: meshes(:)
type(TallyObject), allocatable, target :: tallies(:)
! Tally map structure
type(TallyMap), allocatable :: tally_maps(:)
integer :: n_meshes ! # of structured meshes
integer :: n_tallies ! # of tallies
! Flag for turning tallies on
logical :: tallies_on
! ============================================================================
! CRITICALITY SIMULATION VARIABLES
integer(8) :: n_particles = 10000 ! # of particles per cycle
integer :: n_cycles = 500 ! # of cycles
integer :: n_inactive = 50 ! # of inactive cycles
! External source
type(ExtSource), target :: external_source
@ -88,71 +107,60 @@ module global
integer(8) :: source_index ! index for source particles
! cycle keff
real(8) :: keff
real(8) :: keff = ONE
real(8) :: keff_std
! Flag for turning tallies on
logical :: tallies_on
! ============================================================================
! PARALLEL PROCESSING VARIABLES
! Parallel processing variables
integer :: n_procs ! number of processes
integer :: rank ! rank of process
logical :: master ! master process?
logical :: mpi_enabled ! is MPI in use and initialized?
! Timing variables
! ============================================================================
! TIMING VARIABLES
type(Timer) :: time_total ! timer for total run
type(Timer) :: time_init ! timer for initialization
type(Timer) :: time_intercycle ! timer for intercycle synchronization
type(Timer) :: time_inactive ! timer for inactive cycles
type(Timer) :: time_compute ! timer for computation
! Paths to input file, cross section data, etc
character(MAX_WORD_LEN) :: &
& path_input, &
& path_xsdata
! ===========================================================================
! VARIANCE REDUCTION VARIABLES
! Problem type
integer :: problem_type
! The verbosity controls how much information will be printed to the
! screen and in logs
integer :: verbosity
! Variance reduction options
logical :: survival_biasing = .false.
real(8) :: weight_cutoff = 0.25
real(8) :: weight_survive = 1.0
! Plotting options
! ============================================================================
! PLOTTING VARIABLES
logical :: plotting = .false.
real(8) :: plot_origin(3)
real(8) :: plot_width(2)
real(8) :: plot_basis(6)
real(8) :: pixel
! ============================================================================
! MISCELLANEOUS VARIABLES
character(MAX_WORD_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections ! Path to cross_sections.xml
! Message used in message/warning/fatal_error
character(MAX_LINE_LEN) :: message
! Problem type
integer :: problem_type = PROB_CRITICALITY
! The verbosity controls how much information will be printed to the
! screen and in logs
integer :: verbosity = 7
contains
!===============================================================================
! SET_DEFAULTS gives default values for many global parameters
!===============================================================================
subroutine set_defaults()
! Default problem type is external source
problem_type = PROB_SOURCE
! Default number of particles
n_particles = 10000
! Default verbosity
verbosity = 5
! Defualt multiplication factor
keff = ONE
end subroutine set_defaults
!===============================================================================
! FREE_MEMORY deallocates all allocatable arrays in the program, namely the
! cells, surfaces, materials, and sources
@ -170,10 +178,10 @@ contains
if (allocated(materials)) deallocate(materials)
if (allocated(lattices)) deallocate(lattices)
! Deallocate cross section data
if (allocated(xsdatas)) deallocate(xsdatas)
! Deallocate cross section data and listings
if (allocated(nuclides)) deallocate(nuclides)
if (allocated(sab_tables)) deallocate(sab_tables)
if (allocated(xs_listings)) deallocate(xs_listings)
! Deallocate energy grid
if (allocated(e_grid)) deallocate(e_grid)

205
src/initialize.f90
View file

@ -1,16 +1,17 @@
module initialize
use constants
use cross_section, only: read_xs, read_xsdata
use cross_section, only: read_xs
use datatypes, only: dict_create, dict_add_key, dict_get_key, &
dict_has_key, DICT_NULL, dict_keys
dict_has_key, dict_keys
use datatypes_header, only: ListKeyValueII, DictionaryII
use energy_grid, only: unionized_grid, original_indices
use error, only: fatal_error
use geometry, only: neighbor_lists
use geometry_header, only: Cell, Surface, Universe, Lattice, BASE_UNIVERSE
use global
use input_xml, only: read_input_xml, cells_in_univ_dict
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
use logging, only: create_log
use mcnp_random, only: RN_init_problem
use mpi_routines, only: setup_mpi
@ -55,9 +56,6 @@ contains
! 63-bit RNGs.
call RN_init_problem(3, 0_8, 0_8, 0_8, 0)
! Set default values for settings
call set_defaults()
! set up dictionaries
call create_dictionaries()
@ -81,10 +79,10 @@ contains
if (.not. plotting) then
! Read cross section summary file to determine what files contain
! cross-sections
call read_xsdata(path_xsdata)
call read_cross_sections_xml()
! With the AWRs from the xsdata, change all material specifications so that
! they contain atom percents summing to 1
! With the AWRs from the xs_listings, change all material specifications
! so that they contain atom percents summing to 1
call normalize_ao()
! Read ACE-format cross sections
@ -185,7 +183,7 @@ contains
call dict_create(material_dict)
call dict_create(mesh_dict)
call dict_create(tally_dict)
call dict_create(xsdata_dict)
call dict_create(xs_listing_dict)
call dict_create(nuclide_dict)
call dict_create(sab_dict)
@ -214,18 +212,18 @@ contains
! We also need to allocate the cell count lists for each universe. The logic
! for this is a little more convoluted. In universe_dict, the (key,value)
! pairs are the uid of the universe and the index in the array. In
! ucount_dict, it's the uid of the universe and the number of cells.
! pairs are the id of the universe and the index in the array. In
! ucount_dict, it's the id of the universe and the number of cells.
key_list => dict_keys(universe_dict)
do while (associated(key_list))
! find index of universe in universes array
index = key_list%data%value
univ => universes(index)
univ % uid = key_list%data%key
univ % id = key_list%data%key
! check for lowest level universe
if (univ % uid == 0) BASE_UNIVERSE = index
if (univ % id == 0) BASE_UNIVERSE = index
! find cell count for this universe
count = dict_get_key(cells_in_univ_dict, key_list%data%key)
@ -273,8 +271,7 @@ contains
integer :: j ! index over surface list
integer :: k
integer :: index ! index in surfaces/materials array
integer :: uid ! user-specified uid
character(MAX_LINE_LEN) :: msg ! output/error message
integer :: id ! user-specified id
type(Cell), pointer :: c => null()
type(Lattice), pointer :: l => null()
type(TallyObject), pointer :: t => null()
@ -285,15 +282,15 @@ contains
c => cells(i)
do j = 1, c % n_surfaces
uid = c % surfaces(j)
if (uid < OP_DIFFERENCE) then
if (dict_has_key(surface_dict, abs(uid))) then
index = dict_get_key(surface_dict, abs(uid))
c % surfaces(j) = sign(index, uid)
id = c % surfaces(j)
if (id < OP_DIFFERENCE) then
if (dict_has_key(surface_dict, abs(id))) then
index = dict_get_key(surface_dict, abs(id))
c % surfaces(j) = sign(index, id)
else
msg = "Could not find surface " // trim(int_to_str(abs(uid))) // &
& " specified on cell " // trim(int_to_str(c % uid))
call fatal_error(msg)
message = "Could not find surface " // trim(int_to_str(abs(id))) // &
& " specified on cell " // trim(int_to_str(c % id))
call fatal_error()
end if
end if
end do
@ -301,40 +298,40 @@ contains
! =======================================================================
! ADJUST UNIVERSE INDEX FOR EACH CELL
uid = c % universe
if (dict_has_key(universe_dict, uid)) then
c % universe = dict_get_key(universe_dict, uid)
id = c % universe
if (dict_has_key(universe_dict, id)) then
c % universe = dict_get_key(universe_dict, id)
else
msg = "Could not find universe " // trim(int_to_str(uid)) // &
" specified on cell " // trim(int_to_str(c % uid))
call fatal_error(msg)
message = "Could not find universe " // trim(int_to_str(id)) // &
" specified on cell " // trim(int_to_str(c % id))
call fatal_error()
end if
! =======================================================================
! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL
uid = c % material
if (uid /= 0) then
if (dict_has_key(material_dict, uid)) then
id = c % material
if (id /= 0) then
if (dict_has_key(material_dict, id)) then
c % type = CELL_NORMAL
c % material = dict_get_key(material_dict, uid)
c % material = dict_get_key(material_dict, id)
else
msg = "Could not find material " // trim(int_to_str(uid)) // &
" specified on cell " // trim(int_to_str(c % uid))
call fatal_error(msg)
message = "Could not find material " // trim(int_to_str(id)) // &
" specified on cell " // trim(int_to_str(c % id))
call fatal_error()
end if
else
uid = c % fill
if (dict_has_key(universe_dict, uid)) then
id = c % fill
if (dict_has_key(universe_dict, id)) then
c % type = CELL_FILL
c % fill = dict_get_key(universe_dict, uid)
elseif (dict_has_key(lattice_dict, uid)) then
c % fill = dict_get_key(universe_dict, id)
elseif (dict_has_key(lattice_dict, id)) then
c % type = CELL_LATTICE
c % fill = dict_get_key(lattice_dict, uid)
c % fill = dict_get_key(lattice_dict, id)
else
msg = "Specified fill " // trim(int_to_str(uid)) // " on cell " // &
trim(int_to_str(c % uid)) // " is neither a universe nor a lattice."
call fatal_error(msg)
message = "Specified fill " // trim(int_to_str(id)) // " on cell " // &
trim(int_to_str(c % id)) // " is neither a universe nor a lattice."
call fatal_error()
end if
end if
end do
@ -346,13 +343,13 @@ contains
l => lattices(i)
do j = 1, l % n_x
do k = 1, l % n_y
uid = l % element(j,k)
if (dict_has_key(universe_dict, uid)) then
l % element(j,k) = dict_get_key(universe_dict, uid)
id = l % element(j,k)
if (dict_has_key(universe_dict, id)) then
l % element(j,k) = dict_get_key(universe_dict, id)
else
msg = "Invalid universe number " // trim(int_to_str(uid)) &
// " specified on lattice " // trim(int_to_str(l % uid))
call fatal_error(msg)
message = "Invalid universe number " // trim(int_to_str(id)) &
// " specified on lattice " // trim(int_to_str(l % id))
call fatal_error()
end if
end do
end do
@ -366,13 +363,13 @@ contains
if (t % n_bins(T_CELL) > 0) then
do j = 1, t % n_bins(T_CELL)
uid = t % cell_bins(j) % scalar
if (dict_has_key(cell_dict, uid)) then
t % cell_bins(j) % scalar = dict_get_key(cell_dict, uid)
id = t % cell_bins(j) % scalar
if (dict_has_key(cell_dict, id)) then
t % cell_bins(j) % scalar = dict_get_key(cell_dict, id)
else
msg = "Could not find cell " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find cell " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end do
end if
@ -382,13 +379,13 @@ contains
if (t % n_bins(T_SURFACE) > 0) then
do j = 1, t % n_bins(T_SURFACE)
uid = t % surface_bins(j) % scalar
if (dict_has_key(surface_dict, uid)) then
t % surface_bins(j) % scalar = dict_get_key(surface_dict, uid)
id = t % surface_bins(j) % scalar
if (dict_has_key(surface_dict, id)) then
t % surface_bins(j) % scalar = dict_get_key(surface_dict, id)
else
msg = "Could not find surface " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find surface " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end do
end if
@ -398,13 +395,13 @@ contains
if (t % n_bins(T_UNIVERSE) > 0) then
do j = 1, t % n_bins(T_UNIVERSE)
uid = t % universe_bins(j) % scalar
if (dict_has_key(universe_dict, uid)) then
t % universe_bins(j) % scalar = dict_get_key(universe_dict, uid)
id = t % universe_bins(j) % scalar
if (dict_has_key(universe_dict, id)) then
t % universe_bins(j) % scalar = dict_get_key(universe_dict, id)
else
msg = "Could not find universe " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find universe " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end do
end if
@ -414,13 +411,13 @@ contains
if (t % n_bins(T_MATERIAL) > 0) then
do j = 1, t % n_bins(T_MATERIAL)
uid = t % material_bins(j) % scalar
if (dict_has_key(material_dict, uid)) then
t % material_bins(j) % scalar = dict_get_key(material_dict, uid)
id = t % material_bins(j) % scalar
if (dict_has_key(material_dict, id)) then
t % material_bins(j) % scalar = dict_get_key(material_dict, id)
else
msg = "Could not find material " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find material " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end do
end if
@ -430,13 +427,13 @@ contains
if (t % n_bins(T_CELLBORN) > 0) then
do j = 1, t % n_bins(T_CELLBORN)
uid = t % cellborn_bins(j) % scalar
if (dict_has_key(cell_dict, uid)) then
t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, uid)
id = t % cellborn_bins(j) % scalar
if (dict_has_key(cell_dict, id)) then
t % cellborn_bins(j) % scalar = dict_get_key(cell_dict, id)
else
msg = "Could not find material " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find material " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end do
end if
@ -445,13 +442,13 @@ contains
! ADJUST MESH INDICES FOR EACH TALLY
if (t % n_bins(T_MESH) > 0) then
uid = t % mesh
if (dict_has_key(mesh_dict, uid)) then
t % mesh = dict_get_key(mesh_dict, uid)
id = t % mesh
if (dict_has_key(mesh_dict, id)) then
t % mesh = dict_get_key(mesh_dict, id)
else
msg = "Could not find mesh " // trim(int_to_str(uid)) // &
& " specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find mesh " // trim(int_to_str(id)) // &
& " specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
end if
end do
@ -534,10 +531,10 @@ contains
logical :: percent_in_atom ! nuclides specified in atom percent?
logical :: density_in_atom ! density specified in atom/b-cm?
character(10) :: key ! name of nuclide, e.g. 92235.03c
character(MAX_LINE_LEN) :: msg ! output/error message
type(Material), pointer :: mat => null()
! first find the index in the xsdata array for each nuclide in each material
! first find the index in the xs_listings array for each nuclide in each
! material
do i = 1, n_materials
mat => materials(i)
@ -545,9 +542,9 @@ contains
! given
if (.not. (all(mat%atom_percent > ZERO) .or. &
& all(mat%atom_percent < ZERO))) then
msg = "Cannot mix atom and weight percents in material " // &
& int_to_str(mat%uid)
call fatal_error(msg)
message = "Cannot mix atom and weight percents in material " // &
& int_to_str(mat % id)
call fatal_error()
end if
percent_in_atom = (mat%atom_percent(1) > ZERO)
@ -561,22 +558,22 @@ contains
! Check to make sure cross-section is continuous energy neutron table
n = len_trim(key)
if (key(n:n) /= 'c') then
msg = "Cross-section table " // trim(key) // " is not a " // &
"continuous-energy neutron table."
call fatal_error(msg)
message = "Cross-section table " // trim(key) // &
" is not a continuous-energy neutron table."
call fatal_error()
end if
if (dict_has_key(xsdata_dict, key)) then
index = dict_get_key(xsdata_dict, key)
mat % xsdata(j) = index
if (dict_has_key(xs_listing_dict, key)) then
index = dict_get_key(xs_listing_dict, key)
mat % xs_listing(j) = index
else
msg = "Cannot find cross-section " // trim(key) // " in specified &
&xsdata file."
call fatal_error(msg)
message = "Cannot find cross-section " // trim(key) // &
" in specified cross_sections.xml file."
call fatal_error()
end if
! determine atomic weight ratio
awr = xsdatas(index) % awr
awr = xs_listings(index) % awr
! if given weight percent, convert all values so that they are divided
! by awr. thus, when a sum is done over the values, it's actually
@ -597,8 +594,8 @@ contains
if (.not. density_in_atom) then
sum_percent = ZERO
do j = 1, mat % n_nuclides
index = mat % xsdata(j)
awr = xsdatas(index) % awr
index = mat % xs_listing(j)
awr = xs_listings(index) % awr
x = mat % atom_percent(j)
sum_percent = sum_percent + x*awr
end do

376
src/input_xml.f90
View file

@ -7,9 +7,9 @@ module input_xml
use geometry_header, only: Cell, Surface, Lattice
use global
use mesh_header, only: StructuredMesh
use output, only: message
use output, only: write_message
use string, only: lower_case, int_to_str, str_to_int, str_to_real, &
split_string
split_string, starts_with, ends_with
use tally_header, only: TallyObject
implicit none
@ -46,27 +46,48 @@ contains
integer :: n
integer :: coeffs_reqd
logical :: file_exists
character(MAX_FILE_LEN) :: env_variable
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: msg
character(MAX_LINE_LEN) :: filename
! Display output message
msg = "Reading settings XML file..."
call message(msg, 5)
message = "Reading settings XML file..."
call write_message(5)
! Check if settings.xml exists
filename = trim(path_input) // "settings.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
msg = "Settings XML file '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "Settings XML file '" // trim(filename) // "' does not exist!"
call fatal_error()
end if
! Initialize XML scalar variables
cross_sections_ = ""
verbosity_ = 0
! Parse settings.xml file
call read_xml_file_settings_t(filename)
! Cross-section library path
path_xsdata = trim(xslibrary % path)
! Find cross_sections.xml file -- the first place to look is the
! settings.xml file. If no file is found there, then we check the
! CROSS_SECTIONS environment variable
if (len_trim(cross_sections_) == 0) then
! No cross_sections.xml file specified in settings.xml, check environment
! variable
call get_environment_variable("CROSS_SECTIONS", env_variable)
if (len_trim(env_variable) == 0) then
message = "No cross_sections.xml file was specified in " // &
"settings.xml or in the CROSS_SECTIONS environment " // &
"variable."
call fatal_error()
else
path_cross_sections = trim(env_variable)
end if
else
path_cross_sections = trim(cross_sections_)
end if
! Criticality information
if (criticality % cycles > 0) then
@ -77,7 +98,7 @@ contains
end if
! Verbosity
verbosity = verbosity_
if (verbosity_ > 0) verbosity = verbosity_
if (associated(source_ % coeffs)) then
! Determine external source type
@ -88,19 +109,18 @@ contains
external_source % type = SRC_BOX
coeffs_reqd = 6
case default
msg = "Invalid source type: " // trim(source_ % type)
call fatal_error(msg)
message = "Invalid source type: " // trim(source_ % type)
call fatal_error()
end select
! Coefficients for external surface
n = size(source_ % coeffs)
if (n < coeffs_reqd) then
msg = "Not enough coefficients specified for external source."
print *, n, coeffs_reqd
call fatal_error(msg)
message = "Not enough coefficients specified for external source."
call fatal_error()
elseif (n > coeffs_reqd) then
msg = "Too many coefficients specified for external source."
call fatal_error(msg)
message = "Too many coefficients specified for external source."
call fatal_error()
else
allocate(external_source % values(n))
external_source % values = source_ % coeffs
@ -135,15 +155,14 @@ contains
integer :: coeffs_reqd
logical :: file_exists
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: msg
character(MAX_WORD_LEN) :: word
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
type(Lattice), pointer :: l => null()
! Display output message
msg = "Reading geometry XML file..."
call message(msg, 5)
message = "Reading geometry XML file..."
call write_message(5)
! ==========================================================================
! READ CELLS FROM GEOMETRY.XML
@ -152,8 +171,8 @@ contains
filename = trim(path_input) // "geometry.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
msg = "Geometry XML file '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "Geometry XML file '" // trim(filename) // "' does not exist!"
call fatal_error()
end if
! Parse geometry.xml file
@ -168,29 +187,30 @@ contains
c => cells(i)
! Copy data into cells
c % uid = cell_(i) % uid
c % id = cell_(i) % id
c % universe = cell_(i) % universe
c % material = cell_(i) % material
c % fill = cell_(i) % fill
! Check to make sure that either material or fill was specified
if (c % material == 0 .and. c % fill == 0) then
msg = "Neither material nor fill was specified for cell " // &
trim(int_to_str(c % uid))
call fatal_error(msg)
message = "Neither material nor fill was specified for cell " // &
trim(int_to_str(c % id))
call fatal_error()
end if
! Check to make sure that both material and fill haven't been
! specified simultaneously
if (c % material /= 0 .and. c % fill /= 0) then
msg = "Cannot specify material and fill simultaneously"
call fatal_error(msg)
message = "Cannot specify material and fill simultaneously"
call fatal_error()
end if
! Check to make sure that surfaces were specified
if (.not. associated(cell_(i) % surfaces)) then
msg = "No surfaces specified for cell " // trim(int_to_str(c % uid))
call fatal_error(msg)
message = "No surfaces specified for cell " // &
trim(int_to_str(c % id))
call fatal_error()
end if
! Allocate array for surfaces and copy
@ -200,7 +220,7 @@ contains
c % surfaces = cell_(i) % surfaces
! Add cell to dictionary
call dict_add_key(cell_dict, c % uid, i)
call dict_add_key(cell_dict, c % id, i)
! For cells, we also need to check if there's a new universe --
! also for every cell add 1 to the count of cells for the
@ -228,7 +248,7 @@ contains
s => surfaces(i)
! Copy data into cells
s % uid = surface_(i) % uid
s % id = surface_(i) % id
! Copy and interpret surface type
word = surface_(i) % type
@ -274,8 +294,8 @@ contains
s % type = SURF_GQ
coeffs_reqd = 10
case default
msg = "Invalid surface type: " // trim(surface_(i) % type)
call fatal_error(msg)
message = "Invalid surface type: " // trim(surface_(i) % type)
call fatal_error()
end select
! Check to make sure that the proper number of coefficients
@ -284,14 +304,14 @@ contains
n = size(surface_(i) % coeffs)
if (n < coeffs_reqd) then
msg = "Not enough coefficients specified for surface: " // &
trim(int_to_str(s % uid))
message = "Not enough coefficients specified for surface: " // &
trim(int_to_str(s % id))
print *, n, coeffs_reqd
call fatal_error(msg)
call fatal_error()
elseif (n > coeffs_reqd) then
msg = "Too many coefficients specified for surface: " // &
trim(int_to_str(s % uid))
call fatal_error(msg)
message = "Too many coefficients specified for surface: " // &
trim(int_to_str(s % id))
call fatal_error()
else
allocate(s % coeffs(n))
s % coeffs = surface_(i) % coeffs
@ -310,13 +330,13 @@ contains
case ('periodic')
s % bc = BC_PERIODIC
case default
msg = "Unknown boundary condition '" // trim(word) // &
"' specified on surface " // trim(int_to_str(s % uid))
call fatal_error(msg)
message = "Unknown boundary condition '" // trim(word) // &
"' specified on surface " // trim(int_to_str(s % id))
call fatal_error()
end select
! Add surface to dictionary
call dict_add_key(surface_dict, s % uid, i)
call dict_add_key(surface_dict, s % id, i)
end do
@ -330,8 +350,8 @@ contains
do i = 1, n_lattices
l => lattices(i)
! UID of lattice
l % uid = lattice_(i) % uid
! ID of lattice
l % id = lattice_(i) % id
! Read lattice type
word = lattice_(i) % type
@ -342,15 +362,15 @@ contains
case ('hex', 'hexagon', 'hexagonal')
l % type = LATTICE_HEX
case default
msg = "Invalid lattice type: " // trim(lattice_(i) % type)
call fatal_error(msg)
message = "Invalid lattice type: " // trim(lattice_(i) % type)
call fatal_error()
end select
! Read number of lattice cells in each dimension
n = size(lattice_(i) % dimension)
if (n /= 2 .and. n /= 3) then
msg = "Lattice must be two or three dimensions."
call fatal_error(msg)
message = "Lattice must be two or three dimensions."
call fatal_error()
end if
n_x = lattice_(i) % dimension(1)
n_y = lattice_(i) % dimension(2)
@ -359,18 +379,18 @@ contains
! Read lattice origin location
if (size(lattice_(i) % dimension) /= size(lattice_(i) % origin)) then
msg = "Number of entries on <origin> must be the same as the " // &
"number of entries on <dimension>."
call fatal_error(msg)
message = "Number of entries on <origin> must be the same as " // &
"the number of entries on <dimension>."
call fatal_error()
end if
l % x0 = lattice_(i) % origin(1)
l % y0 = lattice_(i) % origin(2)
! Read lattice widths
if (size(lattice_(i) % width) /= size(lattice_(i) % origin)) then
msg = "Number of entries on <width> must be the same as the " // &
"number of entries on <origin>."
call fatal_error(msg)
message = "Number of entries on <width> must be the same as " // &
"the number of entries on <origin>."
call fatal_error()
end if
l % width_x = lattice_(i) % width(1)
l % width_y = lattice_(i) % width(2)
@ -384,7 +404,7 @@ contains
end do
! Add lattice to dictionary
call dict_add_key(lattice_dict, l % uid, i)
call dict_add_key(lattice_dict, l % id, i)
end do
@ -406,21 +426,20 @@ contains
character(MAX_WORD_LEN) :: units
character(MAX_WORD_LEN) :: name
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: msg
type(Material), pointer :: m => null()
type(nuclide_xml), pointer :: nuc => null()
type(sab_xml), pointer :: sab => null()
! Display output message
msg = "Reading materials XML file..."
call message(msg, 5)
message = "Reading materials XML file..."
call write_message(5)
! Check is materials.xml exists
filename = trim(path_input) // "materials.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
msg = "Material XML file '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "Material XML file '" // trim(filename) // "' does not exist!"
call fatal_error()
end if
! Parse materials.xml file
@ -433,8 +452,8 @@ contains
do i = 1, n_materials
m => materials(i)
! Copy material uid
m % uid = material_(i) % uid
! Copy material id
m % id = material_(i) % id
! Copy density -- the default value for the units is given in the
! material_t.xml file and doesn't need to be specified here, hence case
@ -452,15 +471,16 @@ contains
case ('atom/cm3', 'atom/cc')
m % density = 1.0e-24 * val
case default
msg = "Unkwown units '" // trim(material_(i) % density % units) // &
"' specified on material " // trim(int_to_str(m % uid))
call fatal_error(msg)
message = "Unkwown units '" // trim(material_(i) % density % units) &
// "' specified on material " // trim(int_to_str(m % id))
call fatal_error()
end select
! Check to ensure material has at least one nuclide
if (.not. associated(material_(i) % nuclides)) then
msg = "No nuclides specified on material " // trim(int_to_str(m % uid))
call fatal_error(msg)
message = "No nuclides specified on material " // &
trim(int_to_str(m % id))
call fatal_error()
end if
! allocate arrays in Material object
@ -468,7 +488,7 @@ contains
m % n_nuclides = n
allocate(m % names(n))
allocate(m % nuclide(n))
allocate(m % xsdata(n))
allocate(m % xs_listing(n))
allocate(m % atom_density(n))
allocate(m % atom_percent(n))
@ -481,13 +501,13 @@ contains
! Check if no atom/weight percents were specified or if both atom and
! weight percents were specified
if (nuc % ao == ZERO .and. nuc % wo == ZERO) then
msg = "No atom or weight percent specified for nuclide " // &
message = "No atom or weight percent specified for nuclide " // &
trim(name)
call fatal_error(msg)
call fatal_error()
elseif (nuc % ao /= ZERO .and. nuc % wo /= ZERO) then
msg = "Cannot specify both atom and weight percents for a nuclide: " &
// trim(name)
call fatal_error(msg)
message = "Cannot specify both atom and weight percents for a " &
// "nuclide: " // trim(name)
call fatal_error()
end if
! Copy atom/weight percents
@ -504,13 +524,13 @@ contains
m % sab_name = name
m % has_sab_table = .true.
elseif (size(material_(i) % sab) > 1) then
msg = "Cannot have multiple S(a,b) tables on a single material."
call fatal_error(msg)
message = "Cannot have multiple S(a,b) tables on a single material."
call fatal_error()
end if
end do
! Add material to dictionary
call dict_add_key(material_dict, m % uid, i)
call dict_add_key(material_dict, m % id, i)
end do
@ -527,13 +547,12 @@ contains
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
integer :: uid ! user-specified identifier
integer :: id ! user-specified identifier
integer :: index ! index in meshes array
integer :: n ! size of arrays in mesh specification
integer :: n_words ! number of words read
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: msg
character(MAX_WORD_LEN) :: word
character(MAX_WORD_LEN) :: words(MAX_WORDS)
type(TallyObject), pointer :: t => null()
@ -549,8 +568,8 @@ contains
end if
! Display output message
msg = "Reading tallies XML file..."
call message(msg, 5)
message = "Reading tallies XML file..."
call write_message(5)
! Parse tallies.xml file
call read_xml_file_tallies_t(filename)
@ -569,8 +588,8 @@ contains
! Allocate tallies array
if (.not. associated(tally_)) then
n_tallies = 0
msg = "No tallies present in tallies.xml file!"
call warning(msg)
message = "No tallies present in tallies.xml file!"
call warning()
else
n_tallies = size(tally_)
allocate(tallies(n_tallies))
@ -582,8 +601,8 @@ contains
do i = 1, n_meshes
m => meshes(i)
! copy mesh uid
m % uid = mesh_(i) % id
! copy mesh id
m % id = mesh_(i) % id
! Read mesh type
word = mesh_(i) % type
@ -594,15 +613,15 @@ contains
case ('hex', 'hexagon', 'hexagonal')
m % type = LATTICE_HEX
case default
msg = "Invalid mesh type: " // trim(mesh_(i) % type)
call fatal_error(msg)
message = "Invalid mesh type: " // trim(mesh_(i) % type)
call fatal_error()
end select
! Determine number of dimensions for mesh
n = size(mesh_(i) % dimension)
if (n /= 2 .and. n /= 3) then
msg = "Mesh must be two or three dimensions."
call fatal_error(msg)
message = "Mesh must be two or three dimensions."
call fatal_error()
end if
m % n_dimension = n
@ -616,22 +635,22 @@ contains
! Read mesh origin location
if (m % n_dimension /= size(mesh_(i) % origin)) then
msg = "Number of entries on <origin> must be the same as the " // &
"number of entries on <dimension>."
call fatal_error(msg)
message = "Number of entries on <origin> must be the same as " // &
"the number of entries on <dimension>."
call fatal_error()
end if
m % origin = mesh_(i) % origin
! Read mesh widths
if (size(mesh_(i) % width) /= size(mesh_(i) % origin)) then
msg = "Number of entries on <width> must be the same as the " // &
"number of entries on <origin>."
call fatal_error(msg)
message = "Number of entries on <width> must be the same as " // &
"the number of entries on <origin>."
call fatal_error()
end if
m % width = mesh_(i) % width
! Add mesh to dictionary
call dict_add_key(mesh_dict, m % uid, i)
call dict_add_key(mesh_dict, m % id, i)
end do
! ==========================================================================
@ -648,15 +667,15 @@ contains
t % n_bins = 0
t % stride = 0
! Copy material uid
t % uid = tally_(i) % id
! Copy material id
t % id = tally_(i) % id
! Check to make sure that both cells and surfaces were not specified
if (len_trim(tally_(i) % filters % cell) > 0 .and. &
len_trim(tally_(i) % filters % surface) > 0) then
msg = "Cannot specify both cell and surface filters for tally " &
// trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Cannot specify both cell and surface filters for tally " &
// trim(int_to_str(t % id))
call fatal_error()
end if
! TODO: Parse logical expressions instead of just each word
@ -705,14 +724,14 @@ contains
t % mesh = tally_(i) % filters % mesh
if (t % mesh > 0) then
! Determine index in mesh array for this bin
uid = t % mesh
if (dict_has_key(mesh_dict, uid)) then
index = dict_get_key(mesh_dict, uid)
id = t % mesh
if (dict_has_key(mesh_dict, id)) then
index = dict_get_key(mesh_dict, id)
m => meshes(index)
else
msg = "Could not find mesh " // trim(int_to_str(uid)) // &
" specified on tally " // trim(int_to_str(t % uid))
call fatal_error(msg)
message = "Could not find mesh " // trim(int_to_str(id)) // &
" specified on tally " // trim(int_to_str(t % id))
call fatal_error()
end if
t % n_bins(T_MESH) = t % n_bins(T_MESH) + product(m % dimension)
@ -759,15 +778,15 @@ contains
case ('flux')
t % macro_bins(j) % scalar = MACRO_FLUX
if (t % n_bins(T_ENERGYOUT) > 0) then
msg = "Cannot tally flux with an outgoing energy filter."
call fatal_error(msg)
message = "Cannot tally flux with an outgoing energy filter."
call fatal_error()
end if
case ('total')
t % macro_bins(j) % scalar = MACRO_TOTAL
if (t % n_bins(T_ENERGYOUT) > 0) then
msg = "Cannot tally total reaction rate with an outgoing " &
// "energy filter."
call fatal_error(msg)
message = "Cannot tally total reaction rate with an " &
// "outgoing energy filter."
call fatal_error()
end if
case ('scatter')
t % macro_bins(j) % scalar = MACRO_SCATTER
@ -790,16 +809,16 @@ contains
case ('absorption')
t % macro_bins(j) % scalar = MACRO_ABSORPTION
if (t % n_bins(T_ENERGYOUT) > 0) then
msg = "Cannot tally absorption rate with an outgoing " &
message = "Cannot tally absorption rate with an outgoing " &
// "energy filter."
call fatal_error(msg)
call fatal_error()
end if
case ('fission')
t % macro_bins(j) % scalar = MACRO_FISSION
if (t % n_bins(T_ENERGYOUT) > 0) then
msg = "Cannot tally fission rate with an outgoing " &
message = "Cannot tally fission rate with an outgoing " &
// "energy filter."
call fatal_error(msg)
call fatal_error()
end if
case ('nu-fission')
t % macro_bins(j) % scalar = MACRO_NU_FISSION
@ -810,17 +829,17 @@ contains
! Check to make sure that current is the only desired response
! for this tally
if (n_words > 1) then
msg = "Cannot tally other macro reactions in the same " &
message = "Cannot tally other macro reactions in the same " &
// "tally as surface currents. Separate other macro " &
// "reactions into a distinct tally."
call fatal_error(msg)
call fatal_error()
end if
! Check to make sure that only the mesh filter was specified
if (t % mesh == 0 .or. t % n_bins(T_MESH) /= &
product(t % n_bins, t % n_bins > 0)) then
msg = "Surface currents must be used with a mesh filter only."
call fatal_error(msg)
message = "Surface currents must be used with a mesh filter only."
call fatal_error()
end if
! Since the number of bins for the mesh filter was already set
@ -829,8 +848,8 @@ contains
t % n_bins(T_MESH) = t % n_bins(T_MESH) - product(m % dimension)
! Get pointer to mesh
uid = t % mesh
index = dict_get_key(mesh_dict, uid)
id = t % mesh
index = dict_get_key(mesh_dict, id)
m => meshes(index)
! We need to increase the dimension by one since we also need
@ -844,8 +863,8 @@ contains
end if
case default
msg = "Unknown macro reaction: " // trim(words(j))
call fatal_error(msg)
message = "Unknown macro reaction: " // trim(words(j))
call fatal_error()
end select
end do
t % n_macro_bins = n_words
@ -866,19 +885,18 @@ contains
logical :: file_exists ! does tallies.xml file exist?
character(MAX_LINE_LEN) :: filename
character(MAX_LINE_LEN) :: msg
! Check if plot.xml exists
filename = trim(path_input) // "plot.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
msg = "Plot XML file '" // trim(filename) // "' does not exist!"
call fatal_error(msg)
message = "Plot XML file '" // trim(filename) // "' does not exist!"
call fatal_error()
end if
! Display output message
msg = "Reading plot XML file..."
call message(msg, 5)
message = "Reading plot XML file..."
call write_message(5)
! Parse plot.xml file
call read_xml_file_plot_t(filename)
@ -908,4 +926,102 @@ contains
end subroutine read_plot_xml
!===============================================================================
! READ_CROSS_SECTIONS_XML reads information from a cross_sections.xml file. This
! file contains a listing of the ACE cross sections that may be used.
!===============================================================================
subroutine read_cross_sections_xml()
use xml_data_cross_sections_t
integer :: i ! loop index
integer :: n_listings ! number of listings in cross_sections.xml
logical :: file_exists ! does cross_sections.xml exist?
character(MAX_WORD_LEN) :: directory
type(XsListing), pointer :: listing => null()
! 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
message = "Cross sections XML file '" // trim(path_cross_sections) // &
"' does not exist!"
call fatal_error()
end if
message = "Reading cross sections XML file..."
call write_message(5)
! Initialize directory_ variable
directory_ = ""
! Parse cross_sections.xml file
call read_xml_file_cross_sections_t(path_cross_sections)
! Copy directory information if present
directory = trim(directory_)
! Allocate xs_listings array
if (.not. associated(ace_tables_)) then
message = "No ACE table listings present in cross_sections.xml file!"
call fatal_error()
else
n_listings = size(ace_tables_)
allocate(xs_listings(n_listings))
end if
do i = 1, n_listings
listing => xs_listings(i)
! copy a number of attributes
listing % name = trim(ace_tables_(i) % name)
listing % alias = trim(ace_tables_(i) % alias)
listing % zaid = ace_tables_(i) % zaid
listing % awr = ace_tables_(i) % awr
listing % temp = ace_tables_(i) % temperature
! determine type of cross section
select case(ace_tables_(i) % type)
case ('neutron')
listing % type = ACE_NEUTRON
case ('thermal')
listing % type = ACE_THERMAL
case ('dosimetry')
listing % type = ACE_DOSIMETRY
end select
! determine metastable state
if (ace_tables_(i) % metastable == 0) then
listing % metastable = .false.
else
listing % metastable = .true.
end if
! determine whether binary/ascii
if (ace_tables_(i) % binary == 0) then
listing % binary = .false.
else
listing % binary = .true.
end if
! determine path of cross section table
if (starts_with(ace_tables_(i) % path, '/')) then
listing % path = ace_tables_(i) % path
else
if (ends_with(directory,'/')) then
listing % path = trim(directory) // trim(ace_tables_(i) % path)
else
listing % path = trim(directory) // '/' // trim(ace_tables_(i) % path)
end if
end if
! create dictionary entry for both name and alias
call dict_add_key(xs_listing_dict, listing % name, i)
call dict_add_key(xs_listing_dict, listing % alias, i)
end do
end subroutine read_cross_sections_xml
end module input_xml

6
src/interpolation.f90
View file

@ -2,6 +2,7 @@ module interpolation
use constants
use error, only: fatal_error
use global, only: message
use search, only: binary_search
implicit none
@ -34,7 +35,6 @@ contains
real(8) :: r ! interpolation factor
real(8) :: x0, x1 ! bounding x values
real(8) :: y0, y1 ! bounding y values
character(MAX_LINE_LEN) :: msg
! determine starting location
if (present(loc_start)) then
@ -76,8 +76,8 @@ contains
elseif (n_regions == 1) then
interp = data(loc_interp + 1)
elseif (n_regions > 1) then
msg = "Multiple interpolation regions not yet supported."
call fatal_error(msg)
message = "Multiple interpolation regions not yet supported."
call fatal_error()
end if
! handle special case of histogram interpolation

6
src/logging.f90
View file

@ -15,7 +15,7 @@ contains
character(MAX_WORD_LEN) :: path_log ! path of log file
logical :: file_exists ! does log file already exist?
integer :: ioError ! error status for file access
! integer :: ioError ! error status for file access
! Create filename for log file
path_log = trim(path_input) // ".log"
@ -27,8 +27,8 @@ contains
end if
! Open log file for writing
open(FILE=path_log, UNIT=UNIT_LOG, STATUS='replace', &
& ACTION='write', IOSTAT=ioError)
! open(FILE=path_log, UNIT=UNIT_LOG, STATUS='replace', &
! & ACTION='write', IOSTAT=ioError)
end subroutine create_log

6
src/main.f90
View file

@ -6,7 +6,7 @@ program main
use initialize, only: initialize_run
use mcnp_random, only: RN_init_particle
use mpi_routines, only: synchronize_bank
use output, only: message, header, print_runtime
use output, only: write_message, header, print_runtime
use particle_header, only: Particle
use plot, only: run_plot
use physics, only: transport
@ -71,7 +71,7 @@ contains
call timer_start(time_compute)
msg = "Simulating cycle " // trim(int_to_str(i_cycle)) // "..."
call message(msg, 8)
call write_message(8)
! Set all tallies to zero
n_bank = 0
@ -88,7 +88,7 @@ contains
end if
! set random number seed
i_particle = (i_cycle-1)*n_particles + p % uid
i_particle = (i_cycle-1)*n_particles + p % id
call RN_init_particle(i_particle)
! transport particle

5
src/material_header.f90
View file

@ -7,15 +7,14 @@ module material_header
!===============================================================================
type Material
integer :: uid ! unique identifier
integer :: id ! unique identifier
integer :: n_nuclides ! number of nuclides
character(10), allocatable :: names(:) ! isotope names
integer, allocatable :: xsdata(:) ! index in xsdata list
integer, allocatable :: xs_listing(:) ! index in xs_listings list
integer, allocatable :: nuclide(:) ! index in nuclides array
real(8) :: density ! total atom density in atom/b-cm
real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm
real(8), allocatable :: atom_percent(:) ! atom/weight percent (negative for weight)
real(8), allocatable :: total_xs(:) ! macroscopic cross-section
! S(a,b) data references
logical :: has_sab_table = .false.

2
src/mesh_header.f90
View file

@ -8,7 +8,7 @@ module mesh_header
!===============================================================================
type StructuredMesh
integer :: uid
integer :: id
integer :: type
integer :: n_dimension
integer, allocatable :: dimension(:)

42
src/mpi_routines.f90
View file

@ -4,7 +4,7 @@ module mpi_routines
use error, only: fatal_error
use global
use mcnp_random, only: rang, RN_init_particle, RN_skip
use output, only: message
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use tally_header, only: TallyObject
@ -35,7 +35,6 @@ contains
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer(MPI_ADDRESS_KIND) :: base
character(MAX_LINE_LEN) :: msg ! Error message
type(Bank) :: b
mpi_enabled = .true.
@ -43,22 +42,22 @@ contains
! Initialize MPI
call MPI_INIT(ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Failed to initialize MPI."
call fatal_error(msg)
message = "Failed to initialize MPI."
call fatal_error()
end if
! Determine number of processors
call MPI_COMM_SIZE(MPI_COMM_WORLD, n_procs, ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Could not determine number of processors."
call fatal_error(msg)
message = "Could not determine number of processors."
call fatal_error()
end if
! Determine rank of each processor
call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)
if (ierr /= MPI_SUCCESS) then
msg = "Could not determine MPI rank."
call fatal_error(msg)
message = "Could not determine MPI rank."
call fatal_error()
end if
! Determine master
@ -69,7 +68,7 @@ contains
end if
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % uid, bank_disp(1), ierr)
call MPI_GET_ADDRESS(b % id, bank_disp(1), ierr)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), ierr)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), ierr)
call MPI_GET_ADDRESS(b % E, bank_disp(4), ierr)
@ -113,7 +112,7 @@ contains
integer(8) :: finish ! ending index in local fission bank
integer(8) :: total ! total sites in global fission bank
integer(8) :: count ! index for source bank
integer(8) :: index ! index for uid -- accounts for all nodes
integer(8) :: index ! index for id -- accounts for all nodes
integer :: send_to_left ! # of bank sites to send/recv to or from left
integer :: send_to_right ! # of bank sites to send/recv to or from right
integer(8) :: sites_needed ! # of sites to be sampled
@ -122,7 +121,6 @@ contains
& temp_sites(:), & ! local array of extra sites on each node
& left_bank(:), & ! bank sites to send/recv to or from left node
& right_bank(:) ! bank sites to send/recv to or fram right node
character(MAX_LINE_LEN) :: msg
#ifdef MPI
integer :: ierr
@ -133,8 +131,8 @@ contains
real(8) :: t0, t1, t2, t3, t4
#endif
msg = "Collecting number of fission sites..."
call message(msg, 8)
message = "Collecting number of fission sites..."
call write_message(8)
#ifdef MPI
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
@ -159,8 +157,8 @@ contains
! Check if there are no fission sites
if (total == 0) then
msg = "No fission sites banked!"
call fatal_error(msg)
message = "No fission sites banked!"
call fatal_error()
end if
! Make sure all processors start at the same point for random sampling
@ -171,7 +169,7 @@ contains
allocate(temp_sites(2*work))
count = 0_8 ! Index for local source_bank
index = 0_8 ! Index for global source uid -- must account for all nodes
index = 0_8 ! Index for global source id -- must account for all nodes
if (total < n_particles) then
sites_needed = mod(n_particles,total)
@ -180,8 +178,8 @@ contains
end if
p_sample = real(sites_needed,8)/real(total,8)
msg = "Sampling fission sites..."
call message(msg, 8)
message = "Sampling fission sites..."
call write_message(8)
! ==========================================================================
! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
@ -252,8 +250,8 @@ contains
t2 = MPI_WTIME()
t_sync(2) = t_sync(2) + (t2 - t1)
msg = "Sending fission sites..."
call message(msg, 8)
message = "Sending fission sites..."
call write_message(8)
! ==========================================================================
! SEND BANK SITES TO NEIGHBORS
@ -281,8 +279,8 @@ contains
t_sync(3) = t_sync(3) + (t3 - t2)
#endif
msg = "Constructing source bank..."
call message(msg, 8)
message = "Constructing source bank..."
call write_message(8)
! ==========================================================================
! RECONSTRUCT SOURCE BANK

89
src/output.f90
View file

@ -152,14 +152,13 @@ contains
end subroutine header
!===============================================================================
! MESSAGE displays an informational message to the log file and the standard
! output stream.
! WRITE_MESSAGE displays an informational message to the log file and the
! standard output stream.
!===============================================================================
subroutine message(msg, level)
subroutine write_message(level)
character(*), intent(in) :: msg
integer, optional :: level
integer, optional :: level
integer :: n_lines
integer :: i
@ -168,13 +167,13 @@ contains
if (.not. master .and. present(level)) return
if (.not. present(level) .or. level <= verbosity) then
n_lines = (len_trim(msg)-1)/79 + 1
n_lines = (len_trim(message)-1)/79 + 1
do i = 1, n_lines
write(ou, fmt='(1X,A)') trim(msg(79*(i-1)+1:79*i))
write(ou, fmt='(1X,A)') trim(message(79*(i-1)+1:79*i))
end do
end if
end subroutine message
end subroutine write_message
!===============================================================================
! GET_TODAY determines the date and time at which the program began execution
@ -232,13 +231,13 @@ contains
select case (p % type)
case (NEUTRON)
write(ou,*) 'Neutron ' // int_to_str(p % uid)
write(ou,*) 'Neutron ' // int_to_str(p % id)
case (PHOTON)
write(ou,*) 'Photon ' // int_to_str(p % uid)
write(ou,*) 'Photon ' // int_to_str(p % id)
case (ELECTRON)
write(ou,*) 'Electron ' // int_to_str(p % uid)
write(ou,*) 'Electron ' // int_to_str(p % id)
case default
write(ou,*) 'Unknown Particle ' // int_to_str(p % uid)
write(ou,*) 'Unknown Particle ' // int_to_str(p % id)
end select
write(ou,*) ' x = ' // real_to_str(p % xyz(1))
write(ou,*) ' y = ' // real_to_str(p % xyz(2))
@ -258,20 +257,20 @@ contains
if (p % cell > 0) then
c => cells(p % cell)
write(ou,*) ' Cell = ' // int_to_str(c % uid)
write(ou,*) ' Cell = ' // int_to_str(c % id)
else
write(ou,*) ' Cell not determined'
end if
if (p % surface > 0) then
s => surfaces(p % surface)
write(ou,*) ' Surface = ' // int_to_str(s % uid)
write(ou,*) ' Surface = ' // int_to_str(s % id)
else
write(ou,*) ' Surface = None'
end if
u => universes(p % universe)
write(ou,*) ' Universe = ' // int_to_str(u % uid)
write(ou,*) ' Universe = ' // int_to_str(u % id)
write(ou,*)
nullify(c)
@ -315,26 +314,26 @@ contains
type(Lattice), pointer :: l => null()
type(Material), pointer :: m => null()
write(ou,*) 'Cell ' // int_to_str(c % uid)
temp = dict_get_key(cell_dict, c % uid)
write(ou,*) 'Cell ' // int_to_str(c % id)
temp = dict_get_key(cell_dict, c % id)
write(ou,*) ' Array Index = ' // int_to_str(temp)
u => universes(c % universe)
write(ou,*) ' Universe = ' // int_to_str(u % uid)
write(ou,*) ' Universe = ' // int_to_str(u % id)
select case (c % type)
case (CELL_NORMAL)
write(ou,*) ' Fill = NONE'
case (CELL_FILL)
u => universes(c % fill)
write(ou,*) ' Fill = Universe ' // int_to_str(u % uid)
write(ou,*) ' Fill = Universe ' // int_to_str(u % id)
case (CELL_LATTICE)
l => lattices(c % fill)
write(ou,*) ' Fill = Lattice ' // int_to_str(l % uid)
write(ou,*) ' Fill = Lattice ' // int_to_str(l % id)
end select
if (c % material == 0) then
write(ou,*) ' Material = NONE'
else
m => materials(c % material)
write(ou,*) ' Material = ' // int_to_str(m % uid)
write(ou,*) ' Material = ' // int_to_str(m % id)
end if
write(ou,*) ' Parent Cell = ' // int_to_str(c % parent)
string = ""
@ -373,12 +372,12 @@ contains
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
write(ou,*) 'Universe ' // int_to_str(univ % uid)
write(ou,*) 'Universe ' // int_to_str(univ % id)
write(ou,*) ' Level = ' // int_to_str(univ % level)
string = ""
do i = 1, univ % n_cells
c => cells(univ % cells(i))
string = trim(string) // ' ' // int_to_str(c % uid)
string = trim(string) // ' ' // int_to_str(c % id)
end do
write(ou,*) ' Cells =' // trim(string)
write(ou,*)
@ -395,7 +394,7 @@ contains
type(Lattice), pointer :: lat
write(ou,*) 'Lattice ' // int_to_str(lat % uid)
write(ou,*) 'Lattice ' // int_to_str(lat % id)
write(ou,*) ' n_x = ' // int_to_str(lat % n_x)
write(ou,*) ' n_y = ' // int_to_str(lat % n_y)
write(ou,*) ' x0 = ' // real_to_str(lat % x0)
@ -417,7 +416,7 @@ contains
integer :: i
character(MAX_LINE_LEN) :: string
write(ou,*) 'Surface ' // int_to_str(surf % uid)
write(ou,*) 'Surface ' // int_to_str(surf % id)
select case (surf % type)
case (SURF_PX)
string = "X Plane"
@ -493,7 +492,7 @@ contains
type(Nuclide), pointer :: nuc => null()
! Write identifier for material
write(ou,*) 'Material ' // int_to_str(mat % uid)
write(ou,*) 'Material ' // int_to_str(mat % id)
! Write total atom density in atom/b-cm
write(ou,*) ' Atom Density = ' // trim(real_to_str(mat % density)) &
@ -528,7 +527,7 @@ contains
type(TallyObject), pointer :: t
integer :: i
integer :: uid
integer :: id
character(MAX_LINE_LEN) :: string
type(Cell), pointer :: c => null()
type(Surface), pointer :: s => null()
@ -536,14 +535,14 @@ contains
type(Material), pointer :: m => null()
type(StructuredMesh), pointer :: sm => null()
write(ou,*) 'Tally ' // int_to_str(t % uid)
write(ou,*) 'Tally ' // int_to_str(t % id)
if (t % n_bins(T_CELL) > 0) then
string = ""
do i = 1, t % n_bins(T_CELL)
uid = t % cell_bins(i) % scalar
c => cells(uid)
string = trim(string) // ' ' // trim(int_to_str(c % uid))
id = t % cell_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(int_to_str(c % id))
end do
write(ou, *) ' Cell Bins:' // trim(string)
end if
@ -551,9 +550,9 @@ contains
if (t % n_bins(T_SURFACE) > 0) then
string = ""
do i = 1, t % n_bins(T_SURFACE)
uid = t % surface_bins(i) % scalar
s => surfaces(uid)
string = trim(string) // ' ' // trim(int_to_str(s % uid))
id = t % surface_bins(i) % scalar
s => surfaces(id)
string = trim(string) // ' ' // trim(int_to_str(s % id))
end do
write(ou, *) ' Surface Bins:' // trim(string)
end if
@ -561,9 +560,9 @@ contains
if (t % n_bins(T_UNIVERSE) > 0) then
string = ""
do i = 1, t % n_bins(T_UNIVERSE)
uid = t % universe_bins(i) % scalar
u => universes(uid)
string = trim(string) // ' ' // trim(int_to_str(u % uid))
id = t % universe_bins(i) % scalar
u => universes(id)
string = trim(string) // ' ' // trim(int_to_str(u % id))
end do
write(ou, *) ' Material Bins:' // trim(string)
end if
@ -571,17 +570,17 @@ contains
if (t % n_bins(T_MATERIAL) > 0) then
string = ""
do i = 1, t % n_bins(T_MATERIAL)
uid = t % material_bins(i) % scalar
m => materials(uid)
string = trim(string) // ' ' // trim(int_to_str(m % uid))
id = t % material_bins(i) % scalar
m => materials(id)
string = trim(string) // ' ' // trim(int_to_str(m % id))
end do
write(ou, *) ' Material Bins:' // trim(string)
end if
if (t % n_bins(T_MESH) > 0) then
string = ""
uid = t % mesh
sm => meshes(uid)
id = t % mesh
sm => meshes(id)
string = trim(string) // ' ' // trim(int_to_str(sm % dimension(1)))
do i = 2, sm % n_dimension
string = trim(string) // ' x ' // trim(int_to_str(sm % dimension(i)))
@ -592,9 +591,9 @@ contains
if (t % n_bins(T_CELLBORN) > 0) then
string = ""
do i = 1, t % n_bins(T_CELLBORN)
uid = t % cellborn_bins(i) % scalar
c => cells(uid)
string = trim(string) // ' ' // trim(int_to_str(c % uid))
id = t % cellborn_bins(i) % scalar
c => cells(id)
string = trim(string) // ' ' // trim(int_to_str(c % id))
end do
write(ou, *) ' Birth Region Bins:' // trim(string)
end if

3
src/particle_header.f90
View file

@ -11,7 +11,7 @@ module particle_header
type Particle
! Basic data
integer(8) :: uid ! Unique ID
integer(8) :: id ! Unique ID
integer :: type ! Particle type (n, p, e, etc)
! Physical data
@ -68,6 +68,7 @@ contains
p % type = NEUTRON
p % wgt = ONE
p % last_wgt = ONE
p % alive = .true.
p % n_bank = 0
p % cell = 0

146
src/physics.f90
View file

@ -11,9 +11,9 @@ module physics
use global
use interpolation, only: interpolate_tab1
use mcnp_random, only: rang
use output, only: message, print_particle
use output, only: write_message, print_particle
use particle_header, only: Particle
use tally, only: score_tally
use tally, only: score_tally, score_surface_current
use search, only: binary_search
use string, only: int_to_str
@ -36,7 +36,6 @@ contains
real(8) :: distance ! distance particle travels
logical :: found_cell ! found cell which particle is in?
logical :: in_lattice ! is surface crossing in lattice?
character(MAX_LINE_LEN) :: msg ! output/error message
type(Universe), pointer :: univ
if (p % cell == 0) then
@ -45,9 +44,9 @@ contains
! if particle couldn't be located, print error
if (.not. found_cell) then
write(msg, '(A,3ES11.3)') &
write(message, '(A,3ES11.3)') &
"Could not locate cell for particle at: ", p % xyz
call fatal_error(msg)
call fatal_error()
end if
! set birth cell attribute
@ -55,13 +54,13 @@ contains
end if
if (verbosity >= 9) then
msg = "Simulating Particle " // trim(int_to_str(p % uid))
call message(msg)
message = "Simulating Particle " // trim(int_to_str(p % id))
call write_message()
end if
if (verbosity >= 10) then
msg = " Born in cell " // trim(int_to_str(cells(p%cell)%uid))
call message(msg)
message = " Born in cell " // trim(int_to_str(cells(p%cell)%id))
call write_message()
end if
! find energy index, interpolation factor
@ -79,9 +78,6 @@ contains
! Select smaller of the two distances
distance = min(d_to_boundary, d_to_collision)
! Save original coordinates of particle
p % last_xyz = p % xyz
! Advance particle
p % xyz = p % xyz + distance * p % uvw
p % xyz_local = p % xyz_local + distance * p % uvw
@ -100,6 +96,9 @@ contains
! collision
p % surface = 0
call collision(p)
! Save coordinates at collision for tallying purposes
p % last_xyz = p % xyz
end if
end do
@ -393,8 +392,8 @@ contains
! check for very low energy
if (p % E < 1.0e-100_8) then
p % alive = .false.
! msg = "Killing neutron with extremely low energy"
! call warning(msg)
message = "Killing neutron with extremely low energy"
call warning()
end if
! Score collision estimator tallies for any macro tallies -- this is done
@ -413,6 +412,7 @@ contains
if (tallies_on) then
call score_tally(p, scattered, fissioned)
call score_surface_current(p)
end if
! Reset number of particles banked during collision
@ -444,7 +444,6 @@ contains
real(8) :: prob ! cumulative probability
real(8) :: cutoff ! random number
real(8) :: atom_density ! atom density of nuclide in atom/b-cm
character(MAX_LINE_LEN) :: msg ! output/error message
type(Material), pointer :: mat => null()
type(Nuclide), pointer :: nuc => null()
type(Reaction), pointer :: rxn => null()
@ -464,8 +463,8 @@ contains
! Check to make sure that a nuclide was sampled
if (i > mat % n_nuclides) then
msg = "Did not sample any nuclide during collision."
call fatal_error(msg)
message = "Did not sample any nuclide during collision."
call fatal_error()
end if
! Find atom density and microscopic total cross section
@ -528,9 +527,9 @@ contains
! Check to make sure partial fission reaction sampled
if (i > nuc % n_fission) then
msg = "Did not sample any partial fission reaction for " // &
"survival biasing in " // trim(nuc % name)
call fatal_error(msg)
message = "Did not sample any partial fission reaction " // &
"for survival biasing in " // trim(nuc % name)
call fatal_error()
end if
rxn => nuc % reactions(nuc % index_fission(i))
@ -636,10 +635,10 @@ contains
! Check to make sure inelastic scattering reaction sampled
if (i > nuc % n_reaction) then
msg = "Did not sample any reaction for nuclide " // &
message = "Did not sample any reaction for nuclide " // &
trim(nuc % name) // " on material " // &
trim(int_to_str(mat % uid))
call fatal_error(msg)
trim(int_to_str(mat % id))
call fatal_error()
end if
rxn => nuc % reactions(i)
@ -767,7 +766,6 @@ contains
real(8) :: mu_i1jk ! outgoing cosine k for E_in(i+1) and E_out(j)
real(8) :: prob ! probability for sampling Bragg edge
real(8) :: u, v, w ! directional cosines
character(MAX_LINE_LEN) :: msg
type(SAB_Table), pointer :: sab => null()
! Get pointer to S(a,b) table
@ -865,7 +863,7 @@ contains
j = 1
end if
else
msg = "Invalid secondary energy mode on S(a,b) table " // &
message = "Invalid secondary energy mode on S(a,b) table " // &
trim(sab % name)
end if
@ -1027,7 +1025,6 @@ contains
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
logical :: actual_event ! did fission actually occur? (no survival biasing)
character(MAX_LINE_LEN) :: msg ! error message
type(Nuclide), pointer :: nuc
! Get pointer to nuclide
@ -1079,7 +1076,7 @@ contains
if (nu == 0 .or. n_bank == 3*n_particles) return
do i = n_bank + 1, min(n_bank + nu, 3*n_particles)
! Bank source neutrons by copying particle data
fission_bank(i) % uid = p % uid
fission_bank(i) % id = p % id
fission_bank(i) % xyz = p % xyz
! sample cosine of angle
@ -1099,9 +1096,9 @@ contains
NR = nuc % nu_d_precursor_data(loc + 1)
NE = nuc % nu_d_precursor_data(loc + 2 + 2*NR)
if (NR > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&sampling delayed neutron precursor yield."
call fatal_error(msg)
call fatal_error()
end if
! interpolate on energy grid
@ -1254,27 +1251,6 @@ contains
end subroutine inelastic_scatter
!===============================================================================
! N_ABSORPTION handles all absorbing reactions, i.e. (n,gamma), (n,p), (n,a),
! etc.
!===============================================================================
subroutine n_absorption(p)
type(Particle), pointer :: p
integer :: cell_num ! user-specified cell number
character(MAX_LINE_LEN) :: msg ! output/error message
p % alive = .false.
if (verbosity >= 10) then
cell_num = cells(p % cell)%uid
msg = " Absorbed in cell " // trim(int_to_str(cell_num))
call message(msg)
end if
end subroutine n_absorption
!===============================================================================
! SAMPLE_ANGLE samples the cosine of the angle between incident and exiting
! particle directions either from 32 equiprobable bins or from a tabular
@ -1302,7 +1278,6 @@ contains
real(8) :: c_k ! cumulative frequency at k
real(8) :: c_k1 ! cumulative frequency at k+1
real(8) :: p0,p1 ! probability distribution
character(MAX_LINE_LEN) :: msg ! error message
! check if reaction has angular distribution -- if not, sample outgoing
! angle isotropically
@ -1379,20 +1354,20 @@ contains
mu = mu0 + (sqrt(p0*p0 + 2*frac*(xi - c_k))-p0)/frac
end if
else
msg = "Unknown interpolation type: " // trim(int_to_str(interp))
call fatal_error(msg)
message = "Unknown interpolation type: " // trim(int_to_str(interp))
call fatal_error()
end if
if (abs(mu) > ONE) then
msg = "Sampled cosine of angle outside [-1, 1)."
call warning(msg)
message = "Sampled cosine of angle outside [-1, 1)."
call warning()
mu = sign(ONE,mu)
end if
else
msg = "Unknown angular distribution type: " // trim(int_to_str(type))
call fatal_error(msg)
message = "Unknown angular distribution type: " // trim(int_to_str(type))
call fatal_error()
end if
end function sample_angle
@ -1501,15 +1476,14 @@ contains
real(8) :: E_max ! parameter for n-body dist
real(8) :: x, y, v ! intermediate variables for n-body dist
real(8) :: r1, r2, r3, r4, r5, r6
character(MAX_LINE_LEN) :: msg ! error message
! TODO: If there are multiple scattering laws, sample scattering law
! Check for multiple interpolation regions
if (edist % n_interp > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&attempting to sample secondary energy distribution."
call fatal_error(msg)
call fatal_error()
end if
! Determine which secondary energy distribution law to use
@ -1524,9 +1498,9 @@ contains
NE = edist % data(2 + 2*NR)
NET = edist % data(3 + 2*NR + NE)
if (NR > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&attempting to sample equiprobable energy bins."
call fatal_error(msg)
call fatal_error()
end if
! determine index on incoming energy grid and interpolation factor
@ -1591,9 +1565,9 @@ contains
NR = edist % data(1)
NE = edist % data(2 + 2*NR)
if (NR > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&attempting to sample continuous tabular distribution."
call fatal_error(msg)
call fatal_error()
end if
! find energy bin and calculate interpolation factor -- if the energy is
@ -1650,9 +1624,9 @@ contains
if (ND > 0) then
! discrete lines present
msg = "Discrete lines in continuous tabular distributed not &
message = "Discrete lines in continuous tabular distributed not &
&yet supported"
call fatal_error(msg)
call fatal_error()
end if
! determine outgoing energy bin
@ -1685,8 +1659,8 @@ contains
& p_l_k)/frac
end if
else
msg = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error(msg)
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error()
end if
! Now interpolate between incident energy bins i and i + 1
@ -1762,17 +1736,17 @@ contains
! KALBACH-MANN CORRELATED SCATTERING
if (.not. present(mu_out)) then
msg = "Law 44 called without giving mu_out as argument."
call fatal_error(msg)
message = "Law 44 called without giving mu_out as argument."
call fatal_error()
end if
! read number of interpolation regions and incoming energies
NR = edist % data(1)
NE = edist % data(2 + 2*NR)
if (NR > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&attempting to sample Kalbach-Mann distribution."
call fatal_error(msg)
call fatal_error()
end if
! find energy bin and calculate interpolation factor -- if the energy is
@ -1830,9 +1804,9 @@ contains
if (ND > 0) then
! discrete lines present
msg = "Discrete lines in continuous tabular distributed not &
message = "Discrete lines in continuous tabular distributed not &
&yet supported"
call fatal_error(msg)
call fatal_error()
end if
! determine outgoing energy bin
@ -1879,8 +1853,8 @@ contains
KM_R = R_k + (R_k1 - R_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
KM_A = A_k + (A_k1 - A_k)*(E_out - E_l_k)/(E_l_k1 - E_l_k)
else
msg = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error(msg)
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error()
end if
! Now interpolate between incident energy bins i and i + 1
@ -1905,17 +1879,17 @@ contains
! CORRELATED ENERGY AND ANGLE DISTRIBUTION
if (.not. present(mu_out)) then
msg = "Law 44 called without giving mu_out as argument."
call fatal_error(msg)
message = "Law 44 called without giving mu_out as argument."
call fatal_error()
end if
! read number of interpolation regions and incoming energies
NR = edist % data(1)
NE = edist % data(2 + 2*NR)
if (NR > 0) then
msg = "Multiple interpolation regions not supported while &
message = "Multiple interpolation regions not supported while &
&attempting to sample correlated energy-angle distribution."
call fatal_error(msg)
call fatal_error()
end if
! find energy bin and calculate interpolation factor -- if the energy is
@ -1973,9 +1947,9 @@ contains
if (ND > 0) then
! discrete lines present
msg = "Discrete lines in continuous tabular distributed not &
message = "Discrete lines in continuous tabular distributed not &
&yet supported"
call fatal_error(msg)
call fatal_error()
end if
! determine outgoing energy bin
@ -2009,8 +1983,8 @@ contains
& p_l_k)/frac
end if
else
msg = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error(msg)
message = "Unknown interpolation type: " // trim(int_to_str(INTT))
call fatal_error()
end if
! Now interpolate between incident energy bins i and i + 1
@ -2062,8 +2036,8 @@ contains
mu_out = mu_k + (sqrt(p_k*p_k + 2*frac*(r3 - c_k))-p_k)/frac
end if
else
msg = "Unknown interpolation type: " // trim(int_to_str(JJ))
call fatal_error(msg)
message = "Unknown interpolation type: " // trim(int_to_str(JJ))
call fatal_error()
end if
case (66)

6
src/search.f90
View file

@ -2,6 +2,7 @@ module search
use constants, only: ONE, MAX_LINE_LEN
use error, only: fatal_error
use global, only: message
contains
@ -20,14 +21,13 @@ contains
integer :: L
integer :: R
real(8) :: testval
character(MAX_LINE_LEN) :: msg
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
msg = "Value outside of array during binary search"
call fatal_error(msg)
message = "Value outside of array during binary search"
call fatal_error()
end if
do while (R - L > 1)

12
src/source.f90
View file

@ -5,7 +5,7 @@ module source
use cross_section_header, only: Nuclide
use global
use mcnp_random, only: rang, RN_init_particle
use output, only: message
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use physics, only: watt_spectrum
@ -33,7 +33,7 @@ contains
character(MAX_LINE_LEN) :: msg ! error message
msg = 'Initializing source particles...'
call message(msg, 6)
call write_message(6)
! Determine maximum amount of particles to simulate on each processor
maxwork = ceiling(real(n_particles)/n_procs,8)
@ -51,7 +51,7 @@ contains
! Initialize first cycle source bank
do i = 0, n_procs - 1
if (rank == i) then
! UID's of first and last source particles
! ID's of first and last source particles
bank_first = i*maxwork + 1
bank_last = min((i+1)*maxwork, n_particles)
@ -66,7 +66,7 @@ contains
! sample position
r = (/ (rang(), k = 1,3) /)
p % uid = j
p % id = j
p % xyz = p_min + r*(p_max - p_min)
p % xyz_local = p % xyz
p % last_xyz = p % xyz
@ -120,8 +120,8 @@ contains
! point to next source particle
p => source_bank(source_index)
! set uid
p % uid = bank_first + source_index - 1
! set id
p % id = bank_first + source_index - 1
end function get_source_particle

8
src/string.f90
View file

@ -2,6 +2,7 @@ module string
use constants, only: MAX_WORDS, MAX_LINE_LEN, ERROR_INT, ERROR_REAL
use error, only: warning
use global, only: message
implicit none
@ -31,7 +32,6 @@ contains
integer :: i ! current index
integer :: i_start ! starting index of word
integer :: i_end ! ending index of word
character(MAX_LINE_LEN) :: msg
i_start = 0
i_end = 0
@ -49,9 +49,9 @@ contains
if (i_end > 0) then
n = n + 1
if (i_end - i_start + 1 > len(words(n))) then
msg = "The word '" // string(i_start:i_end) // "' is longer than " &
& // "the space allocated for it."
call warning(msg)
message = "The word '" // string(i_start:i_end) // &
"' is longer than the space allocated for it."
call warning()
end if
words(n) = string(i_start:i_end)
! reset indices

410
src/tally.f90
View file

@ -1,12 +1,11 @@
module tally
use constants
use cross_section, only: get_macro_xs
use error, only: fatal_error
use global
use mesh, only: get_mesh_bin, bin_to_mesh_indices, get_mesh_indices
use mesh_header, only: StructuredMesh
use output, only: message, header
use output, only: write_message, header
use search, only: binary_search
use string, only: int_to_str, real_to_str
use tally_header, only: TallyScore, TallyMapItem, TallyMapElement
@ -36,13 +35,12 @@ contains
real(8), save :: k1 = 0. ! accumulated keff
real(8), save :: k2 = 0. ! accumulated keff**2
real(8) :: std ! stdev of keff over active cycles
character(MAX_LINE_LEN) :: msg ! output/error message
#ifdef MPI
integer :: ierr
#endif
msg = "Calculate cycle keff..."
call message(msg, 8)
message = "Calculate cycle keff..."
call write_message(8)
! set k1 and k2 at beginning of run
if (i_cycle == 1) then
@ -103,7 +101,6 @@ contains
integer :: n ! number of bins
integer :: filter_bins ! running total of number of filter bins
integer :: score_bins ! number of scoring bins
character(MAX_LINE_LEN) :: msg ! output/error message
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
@ -243,8 +240,8 @@ contains
if (n > 0) then
score_bins = n
else
msg = "Must have macro tally bins!"
call fatal_error(msg)
message = "Must have macro tally bins!"
call fatal_error()
end if
! Allocate scores for tally
@ -319,7 +316,6 @@ contains
logical :: in_mesh
logical :: has_energyout_bin
logical :: analog
character(MAX_LINE_LEN) :: msg
type(TallyObject), pointer :: t
type(StructuredMesh), pointer :: m
@ -334,11 +330,8 @@ contains
do i = 1, n_tallies
t => tallies(i)
! Handle surface current tallies separately
if (t % surface_current) then
call score_surface_current(p, t)
cycle
end if
! Surface current tallies are treated separately
if (t % surface_current) cycle
! =======================================================================
! DETERMINE SCORING BIN COMBINATION
@ -391,12 +384,7 @@ contains
call get_mesh_bin(m, p % xyz, mesh_bin, in_mesh)
if (.not. in_mesh) cycle
if (t % surface_current) then
msg = "Surface current mesh tally not yet implemented."
call fatal_error(msg)
else
bins(T_MESH) = mesh_bin
end if
bins(T_MESH) = mesh_bin
else
bins(T_MESH) = 1
end if
@ -500,8 +488,8 @@ contains
case (MACRO_NU_FISSION)
cycle
case default
msg = "Invalid macro reaction on analog tally."
call fatal_error(msg)
message = "Invalid macro reaction on analog tally."
call fatal_error()
end select
end if
@ -582,16 +570,17 @@ contains
end subroutine score_tally
!===============================================================================
! SCORE_SURFACE_CURRENT
! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually
! determining which mesh surfaces were crossed
!===============================================================================
subroutine score_surface_current(p, t)
subroutine score_surface_current(p)
type(Particle), pointer :: p
type(TallyObject), pointer :: t
integer :: i ! loop indices
integer :: j ! loop indices
integer :: k ! loop indices
integer :: ijk0(3) ! indices of starting coordinates
integer :: ijk1(3) ! indices of ending coordinates
integer :: n_cross ! number of surface crossings
@ -607,206 +596,211 @@ contains
logical :: x_same ! same starting/ending x index (i)
logical :: y_same ! same starting/ending y index (j)
logical :: z_same ! same starting/ending z index (k)
type(TallyObject), pointer :: t => null()
type(StructuredMesh), pointer :: m => null()
! Get pointer to mesh
m => meshes(t % mesh)
do i = 1, n_tallies
! Copy starting and ending location of particle
xyz0 = p % last_xyz
xyz1 = p % xyz
! Copy starting and ending location of particle
xyz0 = p % last_xyz
xyz1 = p % xyz
! Get pointer to tally and mesh
t => tallies(i)
m => meshes(t % mesh)
! Determine indices for starting and ending location
call get_mesh_indices(m, xyz0, ijk0, start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1, end_in_mesh)
! Skip non-surface-current tallies
if (.not. t % surface_current) cycle
! Check to make sure start or end is in mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) return
! Determine indices for starting and ending location
call get_mesh_indices(m, xyz0, ijk0, start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1, end_in_mesh)
! Calculate number of surface crossings
n_cross = sum(abs(ijk1 - ijk0))
if (n_cross == 0) return
! Check to make sure start or end is in mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) cycle
! Copy particle's direction
uvw = p % uvw
! Calculate number of surface crossings
n_cross = sum(abs(ijk1 - ijk0))
if (n_cross == 0) cycle
! ==========================================================================
! SPECIAL CASES WHERE TWO INDICES ARE THE SAME
! Copy particle's direction
uvw = p % uvw
x_same = (ijk0(1) == ijk1(1))
y_same = (ijk0(2) == ijk1(2))
z_same = (ijk0(3) == ijk1(3))
! =======================================================================
! SPECIAL CASES WHERE TWO INDICES ARE THE SAME
if (x_same .and. y_same) then
! Only z crossings
if (uvw(3) > 0) then
do i = ijk0(3), ijk1(3) - 1
ijk0(3) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
do i = ijk0(3) - 1, ijk1(3), -1
ijk0(3) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
return
elseif (x_same .and. z_same) then
! Only y crossings
if (uvw(2) > 0) then
do i = ijk0(2), ijk1(2) - 1
ijk0(2) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
do i = ijk0(2) - 1, ijk1(2), -1
ijk0(2) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
return
elseif (y_same .and. z_same) then
! Only x crossings
if (uvw(1) > 0) then
do i = ijk0(1), ijk1(1) - 1
ijk0(1) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
do i = ijk0(1) - 1, ijk1(1), -1
ijk0(1) = i
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
return
end if
x_same = (ijk0(1) == ijk1(1))
y_same = (ijk0(2) == ijk1(2))
z_same = (ijk0(3) == ijk1(3))
! ==========================================================================
! GENERIC CASE
! Bounding coordinates
do i = 1, 3
if (uvw(i) > 0) then
xyz_cross(i) = m % origin(i) + ijk0(i) * m % width(i)
else
xyz_cross(i) = m % origin(i) + (ijk0(i) - 1) * m % width(i)
end if
end do
do i = 1, n_cross
! Reset scoring bin index
score_index = 0
! Calculate distance to each bounding surface. We need to treat special
! case where the cosine of the angle is zero since this would result in a
! divide-by-zero.
do j = 1, 3
if (uvw(j) == 0) then
d(j) = INFINITY
if (x_same .and. y_same) then
! Only z crossings
if (uvw(3) > 0) then
do j = ijk0(3), ijk1(3) - 1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
d(j) = (xyz_cross(j) - xyz0(j))/uvw(j)
do j = ijk0(3) - 1, ijk1(3), -1
ijk0(3) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
cycle
elseif (x_same .and. z_same) then
! Only y crossings
if (uvw(2) > 0) then
do j = ijk0(2), ijk1(2) - 1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
do j = ijk0(2) - 1, ijk1(2), -1
ijk0(2) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
cycle
elseif (y_same .and. z_same) then
! Only x crossings
if (uvw(1) > 0) then
do j = ijk0(1), ijk1(1) - 1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
else
do j = ijk0(1) - 1, ijk1(1), -1
ijk0(1) = j
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end do
end if
cycle
end if
! =======================================================================
! GENERIC CASE
! Bounding coordinates
do j = 1, 3
if (uvw(j) > 0) then
xyz_cross(j) = m % origin(j) + ijk0(j) * m % width(j)
else
xyz_cross(j) = m % origin(j) + (ijk0(j) - 1) * m % width(j)
end if
end do
! Determine the closest bounding surface of the mesh cell by calculating
! the minimum distance
do k = 1, n_cross
! Reset scoring bin index
score_index = 0
distance = minval(d)
! Now use the minimum distance and diretion of the particle to determine
! which surface was crossed
if (distance == d(1)) then
if (uvw(1) > 0) then
! Crossing into right mesh cell -- this is treated as outgoing
! current from (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
! Calculate distance to each bounding surface. We need to treat
! special case where the cosine of the angle is zero since this would
! result in a divide-by-zero.
do j = 1, 3
if (uvw(j) == 0) then
d(j) = INFINITY
else
d(j) = (xyz_cross(j) - xyz0(j))/uvw(j)
end if
ijk0(1) = ijk0(1) + 1
xyz_cross(1) = xyz_cross(1) + m % width(1)
else
! Crossing into left mesh cell -- this is treated as incoming
! current in (i-1,j,k)
ijk0(1) = ijk0(1) - 1
xyz_cross(1) = xyz_cross(1) - m % width(1)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
end do
! Determine the closest bounding surface of the mesh cell by
! calculating the minimum distance
distance = minval(d)
! Now use the minimum distance and diretion of the particle to
! determine which surface was crossed
if (distance == d(1)) then
if (uvw(1) > 0) then
! Crossing into right mesh cell -- this is treated as outgoing
! current from (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_RIGHT
end if
ijk0(1) = ijk0(1) + 1
xyz_cross(1) = xyz_cross(1) + m % width(1)
else
! Crossing into left mesh cell -- this is treated as incoming
! current in (i-1,j,k)
ijk0(1) = ijk0(1) - 1
xyz_cross(1) = xyz_cross(1) - m % width(1)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_RIGHT
end if
end if
elseif (distance == d(2)) then
if (uvw(2) > 0) then
! Crossing into front mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
end if
ijk0(2) = ijk0(2) + 1
xyz_cross(2) = xyz_cross(2) + m % width(2)
else
! Crossing into back mesh cell -- this is treated as incoming
! current in (i,j-1,k)
ijk0(2) = ijk0(2) - 1
xyz_cross(2) = xyz_cross(2) - m % width(2)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
end if
end if
else if (distance == d(3)) then
if (uvw(3) > 0) then
! Crossing into top mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
end if
ijk0(3) = ijk0(3) + 1
xyz_cross(3) = xyz_cross(3) + m % width(3)
else
! Crossing into bottom mesh cell -- this is treated as incoming
! current in (i,j,k-1)
ijk0(3) = ijk0(3) - 1
xyz_cross(3) = xyz_cross(3) - m % width(3)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
end if
end if
end if
elseif (distance == d(2)) then
if (uvw(2) > 0) then
! Crossing into front mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_FRONT
end if
ijk0(2) = ijk0(2) + 1
xyz_cross(2) = xyz_cross(2) + m % width(2)
else
! Crossing into back mesh cell -- this is treated as incoming
! current in (i,j-1,k)
ijk0(2) = ijk0(2) - 1
xyz_cross(2) = xyz_cross(2) - m % width(2)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_FRONT
end if
! Check for errors
if (score_index < 0 .or. score_index > t % n_total_bins) then
message = "Score index outside range."
call fatal_error()
end if
else if (distance == d(3)) then
if (uvw(3) > 0) then
! Crossing into top mesh cell -- this is treated as outgoing
! current in (i,j,k)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + OUT_TOP
end if
ijk0(3) = ijk0(3) + 1
xyz_cross(3) = xyz_cross(3) + m % width(3)
else
! Crossing into bottom mesh cell -- this is treated as incoming
! current in (i,j,k-1)
ijk0(3) = ijk0(3) - 1
xyz_cross(3) = xyz_cross(3) - m % width(3)
if (all(ijk0 >= 0) .and. all(ijk0 <= m % dimension)) then
score_index = sum(t % stride(1:3) * ijk0) + IN_TOP
end if
! Add to surface current tally
if (score_index > 0) then
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
end if
! Check for errors
if (score_index < 0 .or. score_index > t % n_total_bins) then
print *, ijk0
print *, ijk1
print *, score_index
print *, t % stride(1:3), t % n_total_bins
call fatal_error("Score_index outside range.")
end if
! Calculate new coordinates
xyz0 = xyz0 + distance * uvw
end do
! Add to surface current tally
if (score_index > 0) then
call add_to_score(t % scores(score_index, 1), p % last_wgt)
end if
! Calculate new coordinates
xyz0 = xyz0 + distance * uvw
end do
end subroutine score_surface_current
@ -995,7 +989,7 @@ contains
t => tallies(i)
! Write header block
call header("TALLY " // trim(int_to_str(t % uid)), 3, UNIT_TALLY)
call header("TALLY " // trim(int_to_str(t % id)), 3, UNIT_TALLY)
! Handle surface current tallies separately
if (t % surface_current) then
@ -1227,19 +1221,19 @@ contains
select case(filter_type)
case (T_UNIVERSE)
index = t % universe_bins(bin) % scalar
label = int_to_str(universes(index) % uid)
label = int_to_str(universes(index) % id)
case (T_MATERIAL)
index = t % material_bins(bin) % scalar
label = int_to_str(materials(index) % uid)
label = int_to_str(materials(index) % id)
case (T_CELL)
index = t % cell_bins(bin) % scalar
label = int_to_str(cells(index) % uid)
label = int_to_str(cells(index) % id)
case (T_CELLBORN)
index = t % cellborn_bins(bin) % scalar
label = int_to_str(cells(index) % uid)
label = int_to_str(cells(index) % id)
case (T_SURFACE)
index = t % surface_bins(bin) % scalar
label = int_to_str(surfaces(index) % uid)
label = int_to_str(surfaces(index) % id)
case (T_MESH)
m => meshes(t % mesh)
allocate(ijk(m % n_dimension))

2
src/tally_header.f90
View file

@ -58,7 +58,7 @@ module tally_header
type TallyObject
! Basic data
integer :: uid
integer :: id
integer :: type
real(8) :: volume
logical :: surface_current = .false.

147
src/utils/convert_xsdata.py Executable file
View file

@ -0,0 +1,147 @@
#!/usr/bin/env python
import os
import sys
from xml.dom.minidom import getDOMImplementation
types = {1: "neutron", 2: "dosimetry", 3: "thermal"}
class Xsdata(object):
def __init__(self, filename):
self._table_dict = {}
self.tables = []
for line in open(filename, 'r'):
words = line.split()
# If this listing is just an alias listing, only assign the alias
# attribute
name = words[1]
alias = words[0]
table = self.find_table(name)
if table:
if name not in table.alias:
table.alias.append(alias)
continue
table = XsdataTable()
table.name = name
table.type = types[int(words[2])]
table.zaid = int(words[3])
table.metastable = int(words[4])
table.awr = float(words[5])
table.temperature = 8.6173423e-11 * float(words[6])
table.binary = int(words[7])
table.path = words[8]
self.tables.append(table)
self._table_dict[name] = table
# Check for common directory
self.directory = os.path.dirname(self.tables[0].path)
for table in self.tables:
if not table.path.startswith(self.directory):
self.directory = None
break
def to_xml(self):
# Create XML document
impl = getDOMImplementation()
doc = impl.createDocument(None, "cross_sections", None)
# Get root element
root = doc.documentElement
# Add a directory node
if self.directory:
directoryNode = doc.createElement("directory")
text = doc.createTextNode(self.directory)
directoryNode.appendChild(text)
root.appendChild(directoryNode)
for table in self.tables:
table.path = os.path.basename(table.path)
# Add a node for each table
for table in self.tables:
node = table.to_xml_node(doc)
root.appendChild(node)
return doc
def find_table(self, name):
if name in self._table_dict:
return self._table_dict[name]
else:
return None
class XsdataTable(object):
def __init__(self):
self.alias = []
def to_xml_node(self, doc):
node = doc.createElement("ace_table")
node.setAttribute("name", self.name)
for attribute in ["alias", "zaid", "type", "metastable",
"awr", "temperature", "binary", "path"]:
if hasattr(self, attribute):
# Join string for alias attribute
if attribute == "alias":
if not self.alias:
continue
string = " ".join(self.alias)
else:
string = "{0}".format(getattr(self,attribute))
# Skip metastable and binary if 0
if attribute == "metastable" and self.metastable == 0:
continue
if attribute == "binary" and self.binary == 0:
continue
# Create attribute node
# nodeAttr = doc.createElement(attribute)
# text = doc.createTextNode(string)
# nodeAttr.appendChild(text)
# node.appendChild(nodeAttr)
node.setAttribute(attribute, string)
return node
if __name__ == '__main__':
# Read command line arguments
if len(sys.argv) < 3:
sys.exit("Usage: convert_xsdata.py xsdataFile xmlFile")
xsdataFile = sys.argv[1]
xmlFile = sys.argv[2]
# Read xsdata and create XML document object
xsdataObject = Xsdata(xsdataFile)
doc = xsdataObject.to_xml()
# Reduce number of lines
lines = doc.toprettyxml(indent=' ')
lines = lines.replace('<alias>\n ','<alias>')
lines = lines.replace('\n </alias>','</alias>')
lines = lines.replace('<zaid>\n ','<zaid>')
lines = lines.replace('\n </zaid>','</zaid>')
lines = lines.replace('<type>\n ','<type>')
lines = lines.replace('\n </type>','</type>')
lines = lines.replace('<awr>\n ','<awr>')
lines = lines.replace('\n </awr>','</awr>')
lines = lines.replace('<temperature>\n ','<temperature>')
lines = lines.replace('\n </temperature>','</temperature>')
lines = lines.replace('<path>\n ','<path>')
lines = lines.replace('\n </path>','</path>')
lines = lines.replace('<metastable>\n ','<metastable>')
lines = lines.replace('\n </metastable>','</metastable>')
lines = lines.replace('<binary>\n ','<binary>')
lines = lines.replace('\n </binary>','</binary>')
# Write document in pretty XML to specified file
f = open(xmlFile, 'w')
f.write(lines)
f.close()

212
src/utils/convert_xsdir.py Executable file
View file

@ -0,0 +1,212 @@
#!/usr/bin/env python
import os
import sys
from xml.dom.minidom import getDOMImplementation
class Xsdir(object):
def __init__(self, filename):
self.f = open(filename, 'r')
self.filename = os.path.abspath(filename)
self.directory = os.path.dirname(filename)
self.awr = {}
self.tables = []
# Read first section (DATAPATH)
line = self.f.readline()
words = line.split()
if words:
if words[0].lower().startswith('datapath'):
index = line.index('=')
self.datapath = line[index+1:].strip()
# Read second section
line = self.f.readline()
words = line.split()
assert len(words) == 3
assert words[0].lower() == 'atomic'
assert words[1].lower() == 'weight'
assert words[2].lower() == 'ratios'
while True:
line = self.f.readline()
words = line.split()
# Check for end of second section
if len(words) % 2 != 0 or words[0] == 'directory':
break
for zaid, awr in zip(words[::2], words[1::2]):
self.awr[zaid] = awr
# Read third section
while words[0] != 'directory':
words = self.f.readline().split()
while True:
words = self.f.readline().split()
if not words:
break
# Handle continuation lines
while words[-1] == '+':
extraWords = self.f.readline().split()
words = words + extraWords
assert len(words) >= 7
# Create XsdirTable object and add to line
table = XsdirTable(self.directory)
self.tables.append(table)
# All tables have at least 7 attributes
table.name = words[0]
table.awr = float(words[1])
table.filename = words[2]
table.access = words[3]
table.filetype = int(words[4])
table.address = int(words[5])
table.tablelength = int(words[6])
if len(words) > 7:
table.recordlength = int(words[7])
if len(words) > 8:
table.entries = int(words[8])
if len(words) > 9:
table.temperature = float(words[9])
if len(words) > 10:
table.ptable = (words[10] == 'ptable')
def to_xml(self):
# Create XML document
impl = getDOMImplementation()
doc = impl.createDocument(None, "cross_sections", None)
# Get root element
root = doc.documentElement
# Add a directory node
if self.directory:
directoryNode = doc.createElement("directory")
text = doc.createTextNode(self.directory)
directoryNode.appendChild(text)
root.appendChild(directoryNode)
for table in self.tables:
table.path = os.path.basename(table.path)
# Add a node for each table
for table in self.tables:
if table.name[-1] in ['e', 'p', 'u', 'h', 'g' ,'m', 'd']:
continue
node = table.to_xml_node(doc)
root.appendChild(node)
return doc
class XsdirTable(object):
def __init__(self, directory=None):
self.directory = None
self.name = None
self.awr = None
self.filename = None
self.access = None
self.filetype = None
self.address = None
self.tablelength = None
self.recordlength = None
self.entries = None
self.temperature = None
self.ptable = False
@property
def path(self):
if self.directory:
return os.path.join(self.directory, self.filename)
else:
return self.filename
@path.setter
def path(self, value):
self.diretory = ''
self.filename = value
@property
def metastable(self):
# Only valid for neutron cross-sections
if not self.name.endswith('c'):
return
# Handle special case of Am-242 and Am-242m
if self.zaid == '95242':
return 1
elif self.zaid == '95642':
return 0
# All other cases
A = int(self.zaid) % 1000
if A > 600:
return 1
else:
return 0
@property
def zaid(self):
if self.name.endswith('c'):
return self.name[:self.name.find('.')]
else:
return 0
def to_xml_node(self, doc):
node = doc.createElement("ace_table")
node.setAttribute("name", self.name)
for attribute in ["alias", "zaid", "type", "metastable",
"awr", "temperature", "binary", "path"]:
if hasattr(self, attribute):
# Join string for alias attribute
if attribute == "alias":
if not self.alias:
continue
string = " ".join(self.alias)
else:
string = "{0}".format(getattr(self,attribute))
# Skip metastable and binary if 0
if attribute == "metastable" and self.metastable == 0:
continue
if attribute == "binary" and self.binary == 0:
continue
# Skip any attribute that is none
if getattr(self, attribute) is None:
continue
# Create attribute node
# nodeAttr = doc.createElement(attribute)
# text = doc.createTextNode(string)
# nodeAttr.appendChild(text)
# node.appendChild(nodeAttr)
node.setAttribute(attribute, string)
return node
if __name__ == '__main__':
# Read command line arguments
if len(sys.argv) < 3:
sys.exit("Usage: convert_xsdir.py xsdirFile xmlFile")
xsdirFile = sys.argv[1]
xmlFile = sys.argv[2]
# Read xsdata and create XML document object
xsdirObject = Xsdir(xsdirFile)
doc = xsdirObject.to_xml()
# Reduce number of lines
lines = doc.toprettyxml(indent=' ')
# Write document in pretty XML to specified file
f = open(xmlFile, 'w')
f.write(lines)
f.close()

23
src/xml-fortran/templates/cross_sections_t.xml Normal file
View file

@ -0,0 +1,23 @@
<?xml version="1.0"?>
<template>
<!-- This is the template for reading cross section listings in OpenMC -->
<options rootname="cross_sections" />
<typedef name="ace_table_xml">
<component name="name" type="word" length="15" />
<component name="alias" type="word" length="15" />
<component name="type" type="word" length="10" default="'neutron'" />
<component name="zaid" type="integer" default="0" />
<component name="metastable" type="integer" default="0" />
<component name="awr" type="double" default="0.0" />
<component name="temperature" type="double" default="0.0" />
<component name="binary" type="integer" default="0" />
<component name="path" type="word" length="255" />
</typedef>
<variable name="ace_tables_" tag="ace_table" type="ace_table_xml" dimension="1" />
<variable name="directory_" tag="directory" type="word" length="255" />
</template>

6
src/xml-fortran/templates/geometry_t.xml
View file

@ -6,7 +6,7 @@
<options rootname="geometry" />
<typedef name="cell_xml">
<component name="uid" type="integer" />
<component name="id" type="integer" />
<component name="universe" type="integer" default="0" />
<component name="material" type="integer" default="0" />
<component name="fill" type="integer" default="0" />
@ -14,14 +14,14 @@
</typedef>
<typedef name="surface_xml">
<component name="uid" type="integer" />
<component name="id" type="integer" />
<component name="type" type="word" length="15" />
<component name="coeffs" type="double-array" />
<component name="boundary" type="word" length="12" default="'transmit'" />
</typedef>
<typedef name="lattice_xml">
<component name="uid" type="integer" />
<component name="id" type="integer" />
<component name="type" type="word" length="12" />
<component name="dimension" type="integer-array" />
<component name="origin" type="double-array" />

2
src/xml-fortran/templates/materials_t.xml
View file

@ -31,7 +31,7 @@
<!-- Type for specifying a material -->
<typedef name="material_xml">
<component name="uid" type="integer" />
<component name="id" type="integer" />
<component name="density" type="density_xml" />
<component name="nuclides" tag="nuclide" type="nuclide_xml" dimension="1" />
<component name="sab" type="sab_xml" dimension="1" />

6
src/xml-fortran/templates/settings_t.xml
View file

@ -3,10 +3,6 @@
<options rootname="settings" />
<typedef name="xslibrary_xml">
<component name="path" type="word" length="250" />
</typedef>
<typedef name="criticality_xml">
<component name="cycles" type="integer" />
<component name="inactive" type="integer" />
@ -23,11 +19,11 @@
<component name="weight_avg" type="double" default="1.0" />
</typedef>
<variable name="xslibrary" type="xslibrary_xml" />
<variable name="criticality" type="criticality_xml" />
<variable name="verbosity_" tag="verbosity" type="integer" />
<variable name="source_" tag="source" type="source_xml" />
<variable name="survival_" tag="survival_biasing" type="word" length="3" />
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
</template>