diff --git a/include/openmc/settings.h b/include/openmc/settings.h index f800e39f11..be8661f4dd 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -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 //============================================================================== -extern "C" void read_settings(pugi::xml_node* root); +extern "C" void read_settings_xml(); } // namespace openmc diff --git a/src/api.F90 b/src/api.F90 index bfa1314c6d..33f464fa37 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -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. diff --git a/src/initialize.F90 b/src/initialize.F90 index fac8f918cc..0dbaa83bff 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -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 diff --git a/src/initialize.cpp b/src/initialize.cpp index d2cdc6c5a5..2c68401513 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -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; } diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1d7979041e..6c1452e122 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 & - ¬ 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 , , etc. directly - node_mode = root - else - call warning(" 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(" or 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 cutoff is deprecated and should & - &be replaced by .") - 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 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 diff --git a/src/output.F90 b/src/output.F90 index 49dc380a9c..56d1ecfcc2 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/settings.cpp b/src/settings.cpp index 2805a20988..860e576268 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1,11 +1,19 @@ #include "openmc/settings.h" +#include // for numeric_limits +#include +#include + +#include + #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::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(" 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 , , etc. directly + node_mode = root; } else { - fatal_error("Unknown temperature method: " + temp_str); + warning(" 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(" or 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(*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 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 cutoff is deprecated and should " + "be replaced by ."); + 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(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(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 + + + // 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 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(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 " + " 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 diff --git a/src/tallies/trigger_header.F90 b/src/tallies/trigger_header.F90 index 87507224ad..8b3bdbff2d 100644 --- a/src/tallies/trigger_header.F90 +++ b/src/tallies/trigger_header.F90 @@ -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