Convert most of read_settings_xml to C++

This commit is contained in:
Paul Romano 2018-08-24 15:47:47 -05:00
parent 93f36b574d
commit f20496c906
8 changed files with 610 additions and 600 deletions

View file

@ -47,15 +47,14 @@ extern "C" bool write_all_tracks; //!< write track files for every partic
extern "C" bool write_initial_source; //!< write out initial source file?
// Paths to various files
// TODO: Make strings instead of char* once Fortran is gone
extern "C" char* path_input;
extern "C" char* path_statepoint;
extern "C" char* path_sourcepoint;
extern "C" char* path_particle_restart;
extern std::string path_cross_sections;
extern std::string path_multipole;
extern std::string path_output;
extern std::string path_cross_sections; //!< path to cross_sections.xml
extern std::string path_input; //!< directory where main .xml files resides
extern std::string path_multipole; //!< directory containing multipole files
extern std::string path_output; //!< directory where output files are written
extern std::string path_particle_restart; //!< path to a particle restart file
extern std::string path_source;
extern std::string path_sourcepoint; //!< path to a source file
extern std::string path_statepoint; //!< path to a statepoint file
extern "C" int32_t index_entropy_mesh; //!< Index of entropy mesh in global mesh array
extern "C" int32_t index_ufs_mesh; //!< Index of UFS mesh in global mesh array
@ -95,7 +94,7 @@ extern "C" double weight_survive; //!< Survival weight after Russian roul
//! \param[in] root XML node for <settings>
//==============================================================================
extern "C" void read_settings(pugi::xml_node* root);
extern "C" void read_settings_xml();
} // namespace openmc

View file

@ -133,7 +133,7 @@ contains
legendre_to_tabular_points = C_NONE
n_batch_interval = 1
n_lost_particles = 0
n_particles = 0
n_particles = -1
n_source_points = 0
n_state_points = 0
n_tallies = 0
@ -150,7 +150,7 @@ contains
restart_run = .false.
root_universe = -1
run_CE = .true.
run_mode = NONE
run_mode = -1
satisfy_triggers = .false.
call openmc_set_seed(DEFAULT_SEED)
source_latest = .false.

View file

@ -17,10 +17,28 @@ module initialize
implicit none
type(C_PTR), bind(C, name='path_input') :: openmc_path_input
type(C_PTR), bind(C, name='path_statepoint') :: openmc_path_statepoint
type(C_PTR), bind(C, name='path_sourcepoint') :: openmc_path_sourcepoint
type(C_PTR), bind(C, name='path_particle_restart') :: openmc_path_particle_restart
interface
function openmc_path_input() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_output() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_particle_restart() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_statepoint() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
function openmc_path_sourcepoint() result(ptr) bind(C)
import C_PTR
type(C_PTR) :: ptr
end function
end interface
contains
@ -164,29 +182,24 @@ contains
end function is_null
end interface
if (.not. is_null(openmc_path_input)) then
call c_f_pointer(openmc_path_input, string, [255])
if (.not. is_null(openmc_path_input())) then
call c_f_pointer(openmc_path_input(), string, [255])
path_input = to_f_string(string)
else
path_input = ''
end if
if (.not. is_null(openmc_path_statepoint)) then
call c_f_pointer(openmc_path_statepoint, string, [255])
if (.not. is_null(openmc_path_statepoint())) then
call c_f_pointer(openmc_path_statepoint(), string, [255])
path_state_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_sourcepoint)) then
call c_f_pointer(openmc_path_sourcepoint, string, [255])
if (.not. is_null(openmc_path_sourcepoint())) then
call c_f_pointer(openmc_path_sourcepoint(), string, [255])
path_source_point = to_f_string(string)
end if
if (.not. is_null(openmc_path_particle_restart)) then
call c_f_pointer(openmc_path_particle_restart, string, [255])
if (.not. is_null(openmc_path_particle_restart())) then
call c_f_pointer(openmc_path_particle_restart(), string, [255])
path_particle_restart = to_f_string(string)
end if
! Add slash at end of directory if it isn't there
if (len_trim(path_input) > 0 .and. .not. ends_with(path_input, "/")) then
path_input = trim(path_input) // "/"
end if
end subroutine read_command_line
end module initialize

View file

@ -15,6 +15,7 @@
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
// data/functions from Fortran side
extern "C" void print_usage();
@ -93,13 +94,6 @@ void initialize_mpi(MPI_Comm intracomm)
#endif // OPENMC_MPI
inline bool ends_with(std::string const& value, std::string const& ending)
{
if (ending.size() > value.size()) return false;
return std::equal(ending.rbegin(), ending.rend(), value.rbegin());
}
int
parse_command_line(int argc, char* argv[])
{
@ -214,7 +208,14 @@ parse_command_line(int argc, char* argv[])
}
// Determine directory where XML input files are
if (argc > 1 && last_flag < argc) settings::path_input = argv[last_flag + 1];
if (argc > 1 && last_flag < argc - 1) {
settings::path_input = std::string(argv[last_flag + 1]);
// Add slash at end of directory if it isn't there
if (!ends_with(settings::path_input, "/")) {
settings::path_input += "/";
}
}
return 0;
}

View file

@ -22,7 +22,7 @@ module input_xml
use output, only: title, header, print_plot
use photon_header
use plot_header
use random_lcg, only: prn, openmc_set_seed
use random_lcg, only: prn
use surface_header
use set_header, only: SetChar
use settings
@ -79,10 +79,8 @@ module input_xml
type(C_PTR) :: node_ptr
end subroutine read_lattices
subroutine read_settings(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
end subroutine read_settings
subroutine read_settings_xml() bind(C)
end subroutine read_settings_xml
subroutine read_materials(node_ptr) bind(C)
import C_PTR
@ -203,7 +201,8 @@ contains
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_settings_xml()
subroutine read_settings_xml_f(root_ptr) bind(C)
type(C_PTR), value :: root_ptr
character(MAX_LINE_LEN) :: temp_str
integer :: i
@ -211,303 +210,25 @@ contains
integer :: temp_int
integer :: temp_int_array3(3)
integer(C_INT32_T) :: i_start, i_end
integer(C_INT64_T) :: seed
integer(C_INT) :: err
integer, allocatable :: temp_int_array(:)
integer :: n_tracks
logical :: file_exists
character(MAX_LINE_LEN) :: filename
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_mode
type(XMLNode) :: node_cutoff
type(XMLNode) :: node_entropy
type(XMLNode) :: node_ufs
type(XMLNode) :: node_sp
type(XMLNode) :: node_output
type(XMLNode) :: node_res_scat
type(XMLNode) :: node_trigger
type(XMLNode) :: node_vol
type(XMLNode) :: node_tab_leg
type(XMLNode), allocatable :: node_mesh_list(:)
type(XMLNode), allocatable :: node_vol_list(:)
! Check if settings.xml exists
filename = trim(path_input) // "settings.xml"
inquire(FILE=filename, EXIST=file_exists)
if (.not. file_exists) then
if (run_mode /= MODE_PLOTTING) then
call fatal_error("Settings XML file '" // trim(filename) // "' does &
&not exist! In order to run OpenMC, you first need a set of input &
&files; at a minimum, this includes settings.xml, geometry.xml, &
&and materials.xml. Please consult the user's guide at &
&http://openmc.readthedocs.io for further information.")
else
! The settings.xml file is optional if we just want to make a plot.
return
end if
end if
! Get proper XMLNode type given pointer
root % ptr = root_ptr
! Parse settings.xml file
call doc % load_file(filename)
root = doc % document_element()
! Read settings from C++ side
call read_settings(root % ptr)
! Verbosity
if (check_for_node(root, "verbosity")) then
call get_node_value(root, "verbosity", verbosity)
end if
! To this point, we haven't displayed any output since we didn't know what
! the verbosity is. Now that we checked for it, show the title if necessary
if (master) then
if (verbosity >= 2) call title()
end if
call write_message("Reading settings XML file...", 5)
! Find if a multi-group or continuous-energy simulation is desired
if (check_for_node(root, "energy_mode")) then
call get_node_value(root, "energy_mode", temp_str)
temp_str = trim(to_lower(temp_str))
if (temp_str == "mg" .or. temp_str == "multi-group") then
run_CE = .false.
else if (temp_str == "ce" .or. temp_str == "continuous-energy") then
run_CE = .true.
end if
end if
! Look for deprecated cross_sections.xml file in settings.xml
if (check_for_node(root, "cross_sections")) then
call warning("Setting cross_sections in settings.xml has been deprecated.&
& The cross_sections are now set in materials.xml and the &
&cross_sections input to materials.xml and the OPENMC_CROSS_SECTIONS&
& environment variable will take precendent over setting &
&cross_sections in settings.xml.")
call get_node_value(root, "cross_sections", path_cross_sections)
end if
! Look for deprecated windowed_multipole file in settings.xml
if (run_mode /= MODE_PLOTTING) then
if (check_for_node(root, "multipole_library")) then
call warning("Setting multipole_library in settings.xml has been &
&deprecated. The multipole_library is now set in materials.xml and&
& the multipole_library input to materials.xml and the &
&OPENMC_MULTIPOLE_LIBRARY environment variable will take &
&precendent over setting multipole_library in settings.xml.")
call get_node_value(root, "multipole_library", path_multipole)
end if
if (.not. ends_with(path_multipole, "/")) &
path_multipole = trim(path_multipole) // "/"
end if
if (.not. run_CE) then
! Scattering Treatments
if (check_for_node(root, "max_order")) then
call get_node_value(root, "max_order", max_order)
else
! Set to default of largest int - 1, which means to use whatever is
! contained in library.
! This is largest int - 1 because for legendre scattering, a value of
! 1 is added to the order; adding 1 to huge(0) gets you the largest
! negative integer, which is not what we want.
max_order = huge(0) - 1
end if
else
max_order = 0
end if
! Check for a trigger node and get trigger information
if (check_for_node(root, "trigger")) then
node_trigger = root % child("trigger")
! Check if trigger(s) are to be turned on
call get_node_value(node_trigger, "active", trigger_on)
if (trigger_on) then
if (check_for_node(node_trigger, "max_batches") )then
call get_node_value(node_trigger, "max_batches", n_max_batches)
else
call fatal_error("The max_batches must be specified with triggers")
end if
! Get the batch interval to check triggers
if (.not. check_for_node(node_trigger, "batch_interval"))then
pred_batches = .true.
else
call get_node_value(node_trigger, "batch_interval", temp_int)
n_batch_interval = temp_int
if (n_batch_interval <= 0) then
call fatal_error("The batch interval must be greater than zero")
end if
end if
end if
end if
! Check run mode if it hasn't been set from the command line
if (run_mode == NONE) then
if (check_for_node(root, "run_mode")) then
call get_node_value(root, "run_mode", temp_str)
select case (to_lower(temp_str))
case ("eigenvalue")
run_mode = MODE_EIGENVALUE
case ("fixed source")
run_mode = MODE_FIXEDSOURCE
case ("plot")
run_mode = MODE_PLOTTING
case ("particle restart")
run_mode = MODE_PARTICLE
case ("volume")
run_mode = MODE_VOLUME
case default
call fatal_error("Unrecognized run mode: " // &
trim(temp_str) // ".")
end select
! Assume XML specifics <particles>, <batches>, etc. directly
node_mode = root
else
call warning("<run_mode> should be specified.")
! Make sure that either eigenvalue or fixed source was specified
node_mode = root % child("eigenvalue")
if (node_mode % associated()) then
if (run_mode == NONE) run_mode = MODE_EIGENVALUE
else
node_mode = root % child("fixed_source")
if (node_mode % associated()) then
if (run_mode == NONE) run_mode = MODE_FIXEDSOURCE
else
call fatal_error("<eigenvalue> or <fixed_source> not specified.")
end if
end if
end if
end if
if (run_mode == MODE_EIGENVALUE .or. run_mode == MODE_FIXEDSOURCE) then
! Read run parameters
call get_run_parameters(node_mode)
! Check number of active batches, inactive batches, and particles
if (n_batches <= n_inactive) then
call fatal_error("Number of active batches must be greater than zero.")
elseif (n_inactive < 0) then
call fatal_error("Number of inactive batches must be non-negative.")
elseif (n_particles <= 0) then
call fatal_error("Number of particles must be greater than zero.")
end if
end if
! Copy random number seed if specified
if (check_for_node(root, "seed")) then
call get_node_value(root, "seed", seed)
call openmc_set_seed(seed)
end if
! Check for electron treatment
if (check_for_node(root, "electron_treatment")) then
call get_node_value(root, "electron_treatment", temp_str)
select case (to_lower(temp_str))
case ("led")
electron_treatment = ELECTRON_LED
case ("ttb")
electron_treatment = ELECTRON_TTB
case default
call fatal_error("Unrecognized electron treatment: " // &
trim(temp_str) // ".")
end select
end if
! Check for photon transport
if (check_for_node(root, "photon_transport")) then
call get_node_value(root, "photon_transport", photon_transport)
if (.not. run_CE .and. photon_transport) then
call fatal_error("Photon transport is not currently supported &
&in Multi-group mode")
end if
end if
! Number of bins for logarithmic grid
if (check_for_node(root, "log_grid_bins")) then
call get_node_value(root, "log_grid_bins", n_log_bins)
if (n_log_bins < 1) then
call fatal_error("Number of bins for logarithmic grid must be &
&greater than zero.")
end if
else
n_log_bins = 8000
end if
! Number of OpenMP threads
if (check_for_node(root, "threads")) then
#ifdef _OPENMP
if (n_threads == NONE) then
call get_node_value(root, "threads", n_threads)
if (n_threads < 1) then
call fatal_error("Invalid number of threads: " // to_str(n_threads))
end if
call omp_set_num_threads(n_threads)
end if
#else
if (master) call warning("Ignoring number of threads.")
#endif
end if
! ==========================================================================
! EXTERNAL SOURCE
! Handled on C++ side
! Check if we want to write out source
if (check_for_node(root, "write_initial_source")) then
call get_node_value(root, "write_initial_source", write_initial_source)
end if
! Survival biasing
if (check_for_node(root, "survival_biasing")) then
call get_node_value(root, "survival_biasing", survival_biasing)
end if
! Probability tables
if (check_for_node(root, "ptables")) then
call get_node_value(root, "ptables", urr_ptables_on)
end if
! Cutoffs
if (check_for_node(root, "cutoff")) then
node_cutoff = root % child("cutoff")
if (check_for_node(node_cutoff, "weight")) then
call get_node_value(node_cutoff, "weight", weight_cutoff)
end if
if (check_for_node(node_cutoff, "weight_avg")) then
call get_node_value(node_cutoff, "weight_avg", weight_survive)
end if
if (check_for_node(node_cutoff, "energy_neutron")) then
call get_node_value(node_cutoff, "energy_neutron", energy_cutoff(1))
elseif (check_for_node(node_cutoff, "energy")) then
call warning("The use of an <energy> cutoff is deprecated and should &
&be replaced by <energy_neutron>.")
call get_node_value(node_cutoff, "energy", energy_cutoff(1))
end if
if (check_for_node(node_cutoff, "energy_photon")) then
call get_node_value(node_cutoff, "energy_photon", energy_cutoff(2))
end if
if (check_for_node(node_cutoff, "energy_electron")) then
call get_node_value(node_cutoff, "energy_electron", energy_cutoff(3))
end if
if (check_for_node(node_cutoff, "energy_positron")) then
call get_node_value(node_cutoff, "energy_positron", energy_cutoff(4))
end if
end if
! Particle trace
if (check_for_node(root, "trace")) then
call get_node_array(root, "trace", temp_int_array3)
trace_batch = temp_int_array3(1)
trace_gen = temp_int_array3(2)
trace_particle = int(temp_int_array3(3), 8)
if (run_mode == MODE_EIGENVALUE) then
! Preallocate space for keff and entropy by generation
call k_generation % reserve(n_max_batches*gen_per_batch)
call entropy % reserve(n_max_batches*gen_per_batch)
end if
! Particle tracks
@ -694,22 +415,9 @@ contains
call sourcepoint_batch % add(statepoint_batch % get_item(i))
end do
end if
! Check if the user has specified to write binary source file
if (check_for_node(node_sp, "separate")) then
call get_node_value(node_sp, "separate", source_separate)
end if
if (check_for_node(node_sp, "write")) then
call get_node_value(node_sp, "write", source_write)
end if
if (check_for_node(node_sp, "overwrite_latest")) then
call get_node_value(node_sp, "overwrite_latest", source_latest)
source_separate = source_latest
end if
else
! If no <source_point> tag was present, by default we keep source bank in
! statepoint file and write it out at statepoints intervals
source_separate = .false.
n_source_points = n_state_points
do i = 1, n_state_points
call sourcepoint_batch % add(statepoint_batch % get_item(i))
@ -729,91 +437,10 @@ contains
end do
end if
! Check if the user has specified to not reduce tallies at the end of every
! batch
if (check_for_node(root, "no_reduce")) then
call get_node_value(root, "no_reduce", reduce_tallies)
end if
! Check if the user has specified to use confidence intervals for
! uncertainties rather than standard deviations
if (check_for_node(root, "confidence_intervals")) then
call get_node_value(root, "confidence_intervals", confidence_intervals)
end if
! Check for output options
if (check_for_node(root, "output")) then
! Get pointer to output node
node_output = root % child("output")
! Check for summary option
if (check_for_node(node_output, "summary")) then
call get_node_value(node_output, "summary", output_summary)
end if
! Check for ASCII tallies output option
if (check_for_node(node_output, "tallies")) then
call get_node_value(node_output, "tallies", output_tallies)
end if
! Set output directory if a path has been specified
if (check_for_node(node_output, "path")) then
call get_node_value(node_output, "path", path_output)
if (.not. ends_with(path_output, "/")) &
path_output = trim(path_output) // "/"
end if
end if
! Check for cmfd run
if (check_for_node(root, "run_cmfd")) then
call get_node_value(root, "run_cmfd", cmfd_run)
end if
! Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) then
node_res_scat = root % child("resonance_scattering")
! See if resonance scattering is enabled
if (check_for_node(node_res_scat, "enable")) then
call get_node_value(node_res_scat, "enable", res_scat_on)
else
res_scat_on = .true.
end if
! Determine what method is used
if (check_for_node(node_res_scat, "method")) then
call get_node_value(node_res_scat, "method", temp_str)
select case(to_lower(temp_str))
case ('ares')
res_scat_method = RES_SCAT_ARES
case ('dbrc')
res_scat_method = RES_SCAT_DBRC
case ('wcm')
res_scat_method = RES_SCAT_WCM
case default
call fatal_error("Unrecognized resonance elastic scattering method: " &
// trim(temp_str) // ".")
end select
end if
! Minimum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_min")) then
call get_node_value(node_res_scat, "energy_min", res_scat_energy_min)
end if
if (res_scat_energy_min < ZERO) then
call fatal_error("Lower resonance scattering energy bound is negative")
end if
! Maximum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_max")) then
call get_node_value(node_res_scat, "energy_max", res_scat_energy_max)
end if
if (res_scat_energy_max < res_scat_energy_min) then
call fatal_error("Upper resonance scattering energy bound is below the &
&lower resonance scattering energy bound.")
end if
! Get nuclides that resonance scattering should be applied to
if (check_for_node(node_res_scat, "nuclides")) then
n = node_word_count(node_res_scat, "nuclides")
@ -832,138 +459,7 @@ contains
call volume_calcs(i) % from_xml(node_vol)
end do
! Get temperature settings
if (check_for_node(root, "temperature_default")) then
call get_node_value(root, "temperature_default", temperature_default)
end if
if (check_for_node(root, "temperature_method")) then
call get_node_value(root, "temperature_method", temp_str)
select case (to_lower(temp_str))
case ('nearest')
temperature_method = TEMPERATURE_NEAREST
case ('interpolation')
temperature_method = TEMPERATURE_INTERPOLATION
case default
call fatal_error("Unknown temperature method: " // trim(temp_str))
end select
end if
if (check_for_node(root, "temperature_tolerance")) then
call get_node_value(root, "temperature_tolerance", temperature_tolerance)
end if
if (check_for_node(root, "temperature_multipole")) then
call get_node_value(root, "temperature_multipole", temperature_multipole)
end if
if (check_for_node(root, "temperature_range")) then
call get_node_array(root, "temperature_range", temperature_range)
end if
! Check for tabular_legendre options
if (check_for_node(root, "tabular_legendre")) then
! Get pointer to tabular_legendre node
node_tab_leg = root % child("tabular_legendre")
! Check for enable option
if (check_for_node(node_tab_leg, "enable")) then
call get_node_value(node_tab_leg, "enable", legendre_to_tabular)
end if
! Check for the number of points
if (check_for_node(node_tab_leg, "num_points")) then
call get_node_value(node_tab_leg, "num_points", &
legendre_to_tabular_points)
if (legendre_to_tabular_points <= 1 .and. (.not. run_CE)) then
call fatal_error("The 'num_points' subelement/attribute of the &
&'tabular_legendre' element must contain a value greater than 1")
end if
end if
end if
! Check whether create fission sites
if (run_mode == MODE_FIXEDSOURCE) then
if (check_for_node(root, "create_fission_neutrons")) then
call get_node_value(root, "create_fission_neutrons", &
create_fission_neutrons)
end if
end if
! Close settings XML file
call doc % clear()
end subroutine read_settings_xml
!===============================================================================
! GET_RUN_PARAMETERS
!===============================================================================
subroutine get_run_parameters(node_base)
type(XMLNode), intent(in) :: node_base
character(MAX_LINE_LEN) :: temp_str
type(XMLNode) :: node_keff_trigger
! Check number of particles
if (.not. check_for_node(node_base, "particles")) then
call fatal_error("Need to specify number of particles.")
end if
! Get number of particles if it wasn't specified as a command-line argument
if (n_particles == 0) then
call get_node_value(node_base, "particles", n_particles)
end if
! Get number of basic batches
call get_node_value(node_base, "batches", n_batches)
if (.not. trigger_on) then
n_max_batches = n_batches
end if
n_inactive = 0
gen_per_batch = 1
! Get number of inactive batches
if (run_mode == MODE_EIGENVALUE) then
call get_node_value(node_base, "inactive", n_inactive)
if (check_for_node(node_base, "generations_per_batch")) then
call get_node_value(node_base, "generations_per_batch", gen_per_batch)
end if
! Preallocate space for keff and entropy by generation
call k_generation % reserve(n_max_batches*gen_per_batch)
call entropy % reserve(n_max_batches*gen_per_batch)
! Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) then
node_keff_trigger = node_base % child("keff_trigger")
if (check_for_node(node_keff_trigger, "type")) then
call get_node_value(node_keff_trigger, "type", temp_str)
temp_str = trim(to_lower(temp_str))
select case (temp_str)
case ('std_dev')
keff_trigger % trigger_type = STANDARD_DEVIATION
case ('variance')
keff_trigger % trigger_type = VARIANCE
case ('rel_err')
keff_trigger % trigger_type = RELATIVE_ERROR
case default
call fatal_error("Unrecognized keff trigger type " // temp_str)
end select
else
call fatal_error("Specify keff trigger type in settings XML")
end if
if (check_for_node(node_keff_trigger, "threshold")) then
call get_node_value(node_keff_trigger, "threshold", &
keff_trigger % threshold)
else
call fatal_error("Specify keff trigger threshold in settings XML")
end if
end if
end if
end subroutine get_run_parameters
end subroutine read_settings_xml_f
!===============================================================================
! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking

View file

@ -41,7 +41,7 @@ contains
! developers, version, and date/time which the problem was run.
!===============================================================================
subroutine title()
subroutine title() bind(C)
#ifdef _OPENMP
use omp_lib

View file

@ -1,11 +1,19 @@
#include "openmc/settings.h"
#include <limits> // for numeric_limits
#include <sstream>
#include <string>
#include <omp.h>
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/distribution.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/random_lcg.h"
#include "openmc/source.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
@ -46,33 +54,33 @@ bool urr_ptables_on {true};
bool write_all_tracks {false};
bool write_initial_source {false};
char* path_input;
char* path_statepoint;
char* path_sourcepoint;
char* path_particle_restart;
std::string path_cross_sections;
std::string path_input;
std::string path_multipole;
std::string path_output;
std::string path_particle_restart;
std::string path_source;
std::string path_sourcepoint;
std::string path_statepoint;
int32_t index_entropy_mesh {-1};
int32_t index_ufs_mesh {-1};
int32_t n_batches;
int32_t n_inactive;
int32_t n_inactive {0};
int32_t gen_per_batch {1};
int64_t n_particles {0};
int64_t n_particles {-1};
int electron_treatment {ELECTRON_TTB};
double energy_cutoff[4] {0.0, 1000.0, 0.0, 0.0};
int legendre_to_tabular_points {C_NONE};
int max_order;
int n_log_bins;
int max_order {0};
int n_log_bins {8000};
int n_max_batches;
int res_scat_method {RES_SCAT_ARES};
double res_scat_energy_min {0.01};
double res_scat_energy_max {1000.0};
int run_mode;
int run_mode {-1};
int temperature_method {TEMPERATURE_NEAREST};
double temperature_tolerance {10.0};
double temperature_default {293.6};
@ -85,87 +93,316 @@ int verbosity {7};
double weight_cutoff {0.25};
double weight_survive {1.0};
// TODO: Move to separate file
struct KTrigger {
int type;
double threshold;
};
extern "C" KTrigger keff_trigger;
} // namespace settings
//==============================================================================
// Functions
//==============================================================================
void read_settings(pugi::xml_node* root)
void get_run_parameters(pugi::xml_node node_base)
{
using namespace settings;
using namespace pugi;
// Check number of particles
if (!check_for_node(node_base, "particles")) {
fatal_error("Need to specify number of particles.");
}
// Get number of particles if it wasn't specified as a command-line argument
if (n_particles == -1) {
n_particles = std::stoll(get_node_value(node_base, "particles"));
}
// Get number of basic batches
if (check_for_node(node_base, "batches")) {
n_batches = std::stoi(get_node_value(node_base, "batches"));
}
if (!trigger_on) n_max_batches = n_batches;
// Get number of inactive batches
if (run_mode == RUN_MODE_EIGENVALUE) {
if (check_for_node(node_base, "inactive")) {
n_inactive = std::stoi(get_node_value(node_base, "inactive"));
}
if (check_for_node(node_base, "generations_per_batch")) {
gen_per_batch = std::stoi(get_node_value(node_base, "generations_per_batch"));
}
// TODO: Preallocate space for keff and entropy by generation
// TODO: Read keff_trigger information
// Get the trigger information for keff
if (check_for_node(node_base, "keff_trigger")) {
xml_node node_keff_trigger = node_base.child("keff_trigger");
if (check_for_node(node_keff_trigger, "type")) {
auto temp = get_node_value(node_keff_trigger, "type", true, true);
if (temp == "std_dev") {
keff_trigger.type = STANDARD_DEVIATION;
} else if (temp == "variance") {
keff_trigger.type = VARIANCE;
} else if ( temp == "rel_err") {
keff_trigger.type = RELATIVE_ERROR;
} else {
fatal_error("Unrecognized keff trigger type " + temp);
}
} else {
fatal_error("Specify keff trigger type in settings XML");
}
if (check_for_node(node_keff_trigger, "threshold")) {
keff_trigger.threshold = std::stod(get_node_value(
node_keff_trigger, "threshold"));
} else {
fatal_error("Specify keff trigger threshold in settings XML");
}
}
}
}
extern "C" void title();
extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr);
extern "C" void
read_settings_xml()
{
using namespace settings;
using namespace pugi;
// Check if settings.xml exists
std::string filename = std::string(path_input) + "settings.xml";
if (!file_exists(filename)) {
if (run_mode != RUN_MODE_PLOTTING) {
std::stringstream msg;
msg << "Settings XML file '" << filename << "' does not exist! In order "
"to run OpenMC, you first need a set of input files; at a minimum, this "
"includes settings.xml, geometry.xml, and materials.xml. Please consult "
"the user's guide at http://openmc.readthedocs.io for further "
"information.";
fatal_error(msg);
} else {
// The settings.xml file is optional if we just want to make a plot.
return;
}
}
// Parse settings.xml file
xml_document doc;
auto result = doc.load_file("settings.xml");
if (!result) {
fatal_error("Error processing settings.xml file.");
}
// Get root element
xml_node root = doc.document_element();
// Verbosity
if (check_for_node(root, "verbosity")) {
verbosity = std::stoi(get_node_value(root, "verbosity"));
}
// To this point, we haven't displayed any output since we didn't know what
// the verbosity is. Now that we checked for it, show the title if necessary
if (openmc_master) {
if (verbosity >= 2) title();
}
write_message("Reading settings XML file...", 5);
// Find if a multi-group or continuous-energy simulation is desired
if (check_for_node(root, "energy_mode")) {
std::string temp_str = get_node_value(root, "energy_mode", true, true);
if (temp_str == "mg" || temp_str == "multi-group") {
run_CE = false;
} else if (temp_str == "ce" || temp_str == "continuous-energy") {
run_CE = true;
}
}
// Look for deprecated cross_sections.xml file in settings.xml
if (check_for_node(*root, "cross_sections")) {
if (check_for_node(root, "cross_sections")) {
warning("Setting cross_sections in settings.xml has been deprecated."
" The cross_sections are now set in materials.xml and the "
"cross_sections input to materials.xml and the OPENMC_CROSS_SECTIONS"
" environment variable will take precendent over setting "
"cross_sections in settings.xml.");
path_cross_sections = get_node_value(*root, "cross_sections");
path_cross_sections = get_node_value(root, "cross_sections");
}
// Look for deprecated windowed_multipole file in settings.xml
if (run_mode != RUN_MODE_PLOTTING) {
if (check_for_node(*root, "multipole_library")) {
if (check_for_node(root, "multipole_library")) {
warning("Setting multipole_library in settings.xml has been "
"deprecated. The multipole_library is now set in materials.xml and"
" the multipole_library input to materials.xml and the "
"OPENMC_MULTIPOLE_LIBRARY environment variable will take "
"precendent over setting multipole_library in settings.xml.");
path_multipole = get_node_value(*root, "multipole_library");
path_multipole = get_node_value(root, "multipole_library");
}
if (!ends_with(path_multipole, "/")) {
path_multipole += "/";
}
}
// Check for output options
if (check_for_node(*root, "output")) {
if (!run_CE) {
// Scattering Treatments
if (check_for_node(root, "max_order")) {
max_order = std::stoi(get_node_value(root, "max_order"));
} else {
// Set to default of largest int - 1, which means to use whatever is
// contained in library.
// This is largest int - 1 because for legendre scattering, a value of
// 1 is added to the order; adding 1 to huge(0) gets you the largest
// negative integer, which is not what we want.
max_order = std::numeric_limits<int>::max() - 1;
}
}
// Get pointer to output node
pugi::xml_node node_output = root->child("output");
// Check for a trigger node and get trigger information
if (check_for_node(root, "trigger")) {
xml_node node_trigger = root.child("trigger");
// Set output directory if a path has been specified
if (check_for_node(node_output, "path")) {
path_output = get_node_value(node_output, "path");
if (!ends_with(path_output, "/")) {
path_output += "/";
// Check if trigger(s) are to be turned on
trigger_on = get_node_value_bool(node_trigger, "active");
if (trigger_on) {
if (check_for_node(node_trigger, "max_batches") ){
n_max_batches = std::stoi(get_node_value(node_trigger, "max_batches"));
} else {
fatal_error("<max_batches> must be specified with triggers");
}
// Get the batch interval to check triggers
if (!check_for_node(node_trigger, "batch_interval")){
trigger_predict = true;
} else {
trigger_batch_interval = std::stoi(get_node_value(node_trigger, "batch_interval"));
if (trigger_batch_interval <= 0) {
fatal_error("Trigger batch interval must be greater than zero");
}
}
}
}
// Get temperature settings
if (check_for_node(*root, "temperature_default")) {
temperature_default = std::stod(get_node_value(*root, "temperature_default"));
}
if (check_for_node(*root, "temperature_method")) {
auto temp_str = get_node_value(*root, "temperature_method", true, true);
if (temp_str == "nearest") {
temperature_method = TEMPERATURE_NEAREST;
} else if (temp_str == "interpolation") {
temperature_method = TEMPERATURE_INTERPOLATION;
// Check run mode if it hasn't been set from the command line
xml_node node_mode;
if (run_mode == -1) {
if (check_for_node(root, "run_mode")) {
std::string temp_str = get_node_value(root, "run_mode", true, true);
if (temp_str == "eigenvalue") {
run_mode = RUN_MODE_EIGENVALUE;
} else if (temp_str == "fixed source") {
run_mode = RUN_MODE_FIXEDSOURCE;
} else if (temp_str == "plot") {
run_mode = RUN_MODE_PLOTTING;
} else if (temp_str == "particle restart") {
run_mode = RUN_MODE_PARTICLE;
} else if (temp_str == "volume") {
run_mode = RUN_MODE_VOLUME;
} else {
fatal_error("Unrecognized run mode: " + temp_str);
}
// Assume XML specifies <particles>, <batches>, etc. directly
node_mode = root;
} else {
fatal_error("Unknown temperature method: " + temp_str);
warning("<run_mode> should be specified.");
// Make sure that either eigenvalue or fixed source was specified
node_mode = root.child("eigenvalue");
if (node_mode) {
if (run_mode == -1) run_mode = RUN_MODE_EIGENVALUE;
} else {
node_mode = root.child("fixed_source");
if (node_mode) {
if (run_mode == -1) run_mode = RUN_MODE_FIXEDSOURCE;
} else {
fatal_error("<eigenvalue> or <fixed_source> not specified.");
}
}
}
}
if (check_for_node(*root, "temperature_tolerance")) {
temperature_tolerance = std::stod(get_node_value(*root, "temperature_tolerance"));
if (run_mode == RUN_MODE_EIGENVALUE || run_mode == RUN_MODE_FIXEDSOURCE) {
// Read run parameters
get_run_parameters(node_mode);
// Check number of active batches, inactive batches, and particles
if (n_batches <= n_inactive) {
fatal_error("Number of active batches must be greater than zero.");
} else if (n_inactive < 0) {
fatal_error("Number of inactive batches must be non-negative.");
} else if (n_particles <= 0) {
fatal_error("Number of particles must be greater than zero.");
}
}
if (check_for_node(*root, "temperature_multipole")) {
temperature_multipole = get_node_value_bool(*root, "temperature_multipole");
// Copy random number seed if specified
if (check_for_node(root, "seed")) {
auto seed = std::stoll(get_node_value(root, "seed"));
openmc_set_seed(seed);
}
if (check_for_node(*root, "temperature_range")) {
auto range = get_node_array<double>(*root, "temperature_range");
temperature_range[0] = range[0];
temperature_range[1] = range[1];
// Check for electron treatment
if (check_for_node(root, "electron_treatment")) {
auto temp_str = get_node_value(root, "electron_treatment", true, true);
if (temp_str == "led") {
electron_treatment = ELECTRON_LED;
} else if (temp_str == "ttb") {
electron_treatment = ELECTRON_TTB;
} else {
fatal_error("Unrecognized electron treatment: " + temp_str + ".");
}
}
// Check for photon transport
if (check_for_node(root, "photon_transport")) {
photon_transport = get_node_value_bool(root, "photon_transport");
if (!run_CE && photon_transport) {
fatal_error("Photon transport is not currently supported in "
"multigroup mode");
}
}
// Number of bins for logarithmic grid
if (check_for_node(root, "log_grid_bins")) {
n_log_bins = std::stoi(get_node_value(root, "log_grid_bins"));
if (n_log_bins < 1) {
fatal_error("Number of bins for logarithmic grid must be greater "
"than zero.");
}
}
// Number of OpenMP threads
if (check_for_node(root, "threads")) {
#ifdef _OPENMP
if (openmc_n_threads == 0) {
openmc_n_threads = std::stoi(get_node_value(root, "threads"));
if (openmc_n_threads < 1) {
std::stringstream msg;
msg << "Invalid number of threads: " << openmc_n_threads;
fatal_error(msg);
}
omp_set_num_threads(openmc_n_threads);
}
#else
if (openmc_master) warning("Ignoring number of threads.");
#endif
}
// ==========================================================================
// EXTERNAL SOURCE
// Get point to list of <source> elements and make sure there is at least one
for (pugi::xml_node node : root->children("source")) {
for (pugi::xml_node node : root.children("source")) {
external_sources.emplace_back(node);
}
@ -178,6 +415,268 @@ void read_settings(pugi::xml_node* root)
};
external_sources.push_back(std::move(source));
}
// Check if we want to write out source
if (check_for_node(root, "write_initial_source")) {
write_initial_source = get_node_value_bool(root, "write_initial_source");
}
// Survival biasing
if (check_for_node(root, "survival_biasing")) {
survival_biasing = get_node_value_bool(root, "survival_biasing");
}
// Probability tables
if (check_for_node(root, "ptables")) {
urr_ptables_on = get_node_value_bool(root, "ptables");
}
// Cutoffs
if (check_for_node(root, "cutoff")) {
xml_node node_cutoff = root.child("cutoff");
if (check_for_node(node_cutoff, "weight")) {
weight_cutoff = std::stod(get_node_value(node_cutoff, "weight"));
}
if (check_for_node(node_cutoff, "weight_avg")) {
weight_survive = std::stod(get_node_value(node_cutoff, "weight_avg"));
}
if (check_for_node(node_cutoff, "energy_neutron")) {
energy_cutoff[0] = std::stod(get_node_value(node_cutoff, "energy_neutron"));
} else if (check_for_node(node_cutoff, "energy")) {
warning("The use of an <energy> cutoff is deprecated and should "
"be replaced by <energy_neutron>.");
energy_cutoff[0] = std::stod(get_node_value(node_cutoff, "energy"));
}
if (check_for_node(node_cutoff, "energy_photon")) {
energy_cutoff[1] = std::stod(get_node_value(node_cutoff, "energy_photon"));
}
if (check_for_node(node_cutoff, "energy_electron")) {
energy_cutoff[2] = std::stof(get_node_value(node_cutoff, "energy_electron"));
}
if (check_for_node(node_cutoff, "energy_positron")) {
energy_cutoff[3] = std::stod(get_node_value(node_cutoff, "energy_positron"));
}
}
// Particle trace
if (check_for_node(root, "trace")) {
auto temp = get_node_array<int64_t>(root, "trace");
trace_batch = temp.at(0);
trace_gen = temp.at(1);
trace_particle = temp.at(2);
}
// Particle tracks
if (check_for_node(root, "track")) {
// Get values and make sure there are three per particle
auto temp = get_node_array<int64_t>(root, "track");
if (temp.size() % 3 != 0) {
fatal_error("Number of integers specified in 'track' is not "
"divisible by 3. Please provide 3 integers per particle to be "
"tracked.");
}
// Reshape into track_identifiers
//allocate(track_identifiers(3, n_tracks/3))
//track_identifiers = reshape(temp_int_array, [3, n_tracks/3])
}
// TODO: Read meshes
// TODO: Read <state_point>
// Check if the user has specified to write source points
if (check_for_node(root, "source_point")) {
// Get source_point node
xml_node node_sp = root.child("source_point");
// TODO: Read source point batches
// Check if the user has specified to write binary source file
if (check_for_node(node_sp, "separate")) {
source_separate = get_node_value_bool(node_sp, "separate");
}
if (check_for_node(node_sp, "write")) {
source_write = get_node_value_bool(node_sp, "write");
}
if (check_for_node(node_sp, "overwrite_latest")) {
source_latest = get_node_value_bool(node_sp, "overwrite_latest");
source_separate = source_latest;
}
} else {
// If no <source_point> tag was present, by default we keep source bank in
// statepoint file and write it out at statepoints intervals
source_separate = false;
// TODO: add defaults
}
// TODO: Check source points are subset
// Check if the user has specified to not reduce tallies at the end of every
// batch
if (check_for_node(root, "no_reduce")) {
reduce_tallies = get_node_value_bool(root, "no_reduce");
}
// Check if the user has specified to use confidence intervals for
// uncertainties rather than standard deviations
if (check_for_node(root, "confidence_intervals")) {
confidence_intervals = get_node_value_bool(root, "confidence_intervals");
}
// Check for output options
if (check_for_node(root, "output")) {
// Get pointer to output node
pugi::xml_node node_output = root.child("output");
// Check for summary option
if (check_for_node(node_output, "summary")) {
output_summary = get_node_value_bool(node_output, "summary");
}
// Check for ASCII tallies output option
if (check_for_node(node_output, "tallies")) {
output_tallies = get_node_value_bool(node_output, "tallies");
}
// Set output directory if a path has been specified
if (check_for_node(node_output, "path")) {
path_output = get_node_value(node_output, "path");
if (!ends_with(path_output, "/")) {
path_output += "/";
}
}
}
// Check for cmfd run
if (check_for_node(root, "run_cmfd")) {
cmfd_run = get_node_value_bool(root, "run_cmfd");
}
// Resonance scattering parameters
if (check_for_node(root, "resonance_scattering")) {
xml_node node_res_scat = root.child("resonance_scattering");
// See if resonance scattering is enabled
if (check_for_node(node_res_scat, "enable")) {
res_scat_on = get_node_value_bool(node_res_scat, "enable");
} else {
res_scat_on = true;
}
// Determine what method is used
if (check_for_node(node_res_scat, "method")) {
auto temp = get_node_value(node_res_scat, "method", true, true);
if (temp == "ares") {
res_scat_method = RES_SCAT_ARES;
} else if (temp == "dbrc") {
res_scat_method = RES_SCAT_DBRC;
} else if (temp == "wcm") {
res_scat_method = RES_SCAT_WCM;
} else {
fatal_error("Unrecognized resonance elastic scattering method: "
+ temp + ".");
}
}
// Minimum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_min")) {
res_scat_energy_min = std::stod(get_node_value(node_res_scat, "energy_min"));
}
if (res_scat_energy_min < 0.0) {
fatal_error("Lower resonance scattering energy bound is negative");
}
// Maximum energy for resonance scattering
if (check_for_node(node_res_scat, "energy_max")) {
res_scat_energy_max = std::stod(get_node_value(node_res_scat, "energy_max"));
}
if (res_scat_energy_max < res_scat_energy_min) {
fatal_error("Upper resonance scattering energy bound is below the"
"lower resonance scattering energy bound.");
}
// TODO: Get resonance scattering nuclides
}
// TODO: Get volume calculations
// Get temperature settings
if (check_for_node(root, "temperature_default")) {
temperature_default = std::stod(get_node_value(root, "temperature_default"));
}
if (check_for_node(root, "temperature_method")) {
auto temp = get_node_value(root, "temperature_method", true, true);
if (temp == "nearest") {
temperature_method = TEMPERATURE_NEAREST;
} else if (temp == "interpolation") {
temperature_method = TEMPERATURE_INTERPOLATION;
} else {
fatal_error("Unknown temperature method: " + temp);
}
}
if (check_for_node(root, "temperature_tolerance")) {
temperature_tolerance = std::stod(get_node_value(root, "temperature_tolerance"));
}
if (check_for_node(root, "temperature_multipole")) {
temperature_multipole = get_node_value_bool(root, "temperature_multipole");
}
if (check_for_node(root, "temperature_range")) {
auto range = get_node_array<double>(root, "temperature_range");
temperature_range[0] = range.at(0);
temperature_range[1] = range.at(1);
}
// Check for tabular_legendre options
if (check_for_node(root, "tabular_legendre")) {
// Get pointer to tabular_legendre node
xml_node node_tab_leg = root.child("tabular_legendre");
// Check for enable option
if (check_for_node(node_tab_leg, "enable")) {
legendre_to_tabular = get_node_value_bool(node_tab_leg, "enable");
}
// Check for the number of points
if (check_for_node(node_tab_leg, "num_points")) {
legendre_to_tabular_points = std::stoi(get_node_value(
node_tab_leg, "num_points"));
if (legendre_to_tabular_points <= 1 && !run_CE) {
fatal_error("The 'num_points' subelement/attribute of the "
"<tabular_legendre> element must contain a value greater than 1");
}
}
}
// Check whether create fission sites
if (run_mode == RUN_MODE_FIXEDSOURCE) {
if (check_for_node(root, "create_fission_neutrons")) {
create_fission_neutrons = get_node_value_bool(root, "create_fission_neutrons");
}
}
// Read remaining settings from Fortran side
read_settings_xml_f(root.internal_object());
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" {
const char* openmc_path_input() {
return settings::path_input.c_str();
}
const char* openmc_path_statepoint() {
return settings::path_statepoint.c_str();
}
const char* openmc_path_sourcepoint() {
return settings::path_sourcepoint.c_str();
}
const char* openmc_path_particle_restart() {
return settings::path_particle_restart.c_str();
}
}
} // namespace openmc

View file

@ -1,5 +1,7 @@
module trigger_header
use, intrinsic :: ISO_C_BINDING
use constants, only: NONE, N_FILTER_TYPES, ZERO
implicit none
@ -22,11 +24,11 @@ module trigger_header
!===============================================================================
! KTRIGGER describes a user-specified precision trigger for k-effective
!===============================================================================
type, public :: KTrigger
integer :: trigger_type = 0
real(8) :: threshold = ZERO
type, public, bind(C) :: KTrigger
integer(C_INT) :: trigger_type = 0
real(C_DOUBLE) :: threshold = ZERO
end type KTrigger
type(KTrigger), public :: keff_trigger ! trigger for k-effective
type(KTrigger), public, bind(C) :: keff_trigger ! trigger for k-effective
end module trigger_header