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

Merge branch 'mcpl_output_dev' into mcpl_input

Browse source
This commit is contained in:
Erik B Knudsen 2022-10-14 01:43:40 +02:00
commit 309893230c
12 changed files with 337 additions and 84 deletions
Download patch file Download diff file Expand all files Collapse all files

20
CMakeLists.txt
View file

@ -159,6 +159,13 @@ if(${HDF5_VERSION} VERSION_GREATER_EQUAL 1.12.0)
list(APPEND cxxflags -DH5Oget_info_by_idx_vers=1 -DH5O_info_t_vers=1)
endif()
#===============================================================================
# MCPL
#===============================================================================
if (OPENMC_USE_MCPL)
find_package(MCPL REQUIRED)
endif()
#===============================================================================
# Set compile/link flags based on which compiler is being used
#===============================================================================
@ -186,10 +193,6 @@ if(OPENMC_ENABLE_COVERAGE)
list(APPEND ldflags --coverage)
endif()
if(OPENMC_USE_MCPL)
list(APPEND ldflags -lmcpl)
endif()
# Show flags being used
message(STATUS "OpenMC C++ flags: ${cxxflags}")
message(STATUS "OpenMC Linker flags: ${ldflags}")
@ -444,10 +447,6 @@ endif()
if (OPENMC_USE_MPI)
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MPI)
endif()
if(OPENMC_USE_MCPL)
target_compile_definitions(libopenmc PUBLIC -DOPENMC_MCPL)
endif()
# Set git SHA1 hash as a compile definition
if(GIT_FOUND)
@ -491,6 +490,11 @@ if (OPENMC_USE_MPI)
target_link_libraries(libopenmc MPI::MPI_CXX)
endif()
if (OPENMC_USE_MCPL)
target_compile_definitions(libopenmc PUBLIC OPENMC_MCPL)
target_link_libraries(libopenmc MCPL::mcpl)
endif()
#===============================================================================
# openmc executable
#===============================================================================

2
include/openmc/settings.h
View file

@ -47,7 +47,9 @@ extern "C" bool run_CE; //!< run with continuous-energy data?
extern bool source_latest; //!< write latest source at each batch?
extern bool source_separate; //!< write source to separate file?
extern bool source_write; //!< write source in HDF5 files?
extern bool source_mcpl_write; //!< write source in mcpl files?
extern bool surf_source_write; //!< write surface source file?
extern bool surf_mcpl_write; //!< write surface mcpl file?
extern bool surf_source_read; //!< read surface source file?
extern bool survival_biasing; //!< use survival biasing?
extern bool temperature_multipole; //!< use multipole data?

9
include/openmc/state_point.h
View file

@ -9,6 +9,10 @@
#include "openmc/particle.h"
#include "openmc/vector.h"
#ifdef OPENMC_MCPL
#include <mcpl.h>
#endif
namespace openmc {
void load_state_point();
@ -21,5 +25,10 @@ void write_tally_results_nr(hid_t file_id);
void restart_set_keff();
void write_unstructured_mesh_results();
#ifdef OPENMC_MCPL
void write_mcpl_source_point(const char *filename_, bool surf_source_bank = false);
void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank);
#endif
} // namespace openmc
#endif // OPENMC_STATE_POINT_H

9
openmc/settings.py
View file

@ -171,6 +171,7 @@ class Settings:
banked (int)
:max_particles: Maximum number of particles to be banked on
surfaces per process (int)
:mcpl: Output in the form of an MCPL-file (bool)
survival_biasing : bool
Indicate whether survival biasing is to be used
tabular_legendre : dict
@ -652,7 +653,7 @@ class Settings:
cv.check_type('surface source writing options', surf_source_write, Mapping)
for key, value in surf_source_write.items():
cv.check_value('surface source writing key', key,
('surface_ids', 'max_particles'))
('surface_ids', 'max_particles', 'mcpl'))
if key == 'surface_ids':
cv.check_type('surface ids for source banking', value,
Iterable, Integral)
@ -664,6 +665,9 @@ class Settings:
value, Integral)
cv.check_greater_than('maximum particle banks on surfaces per process',
value, 0)
elif key == 'mcpl':
cv.check_type('write to an MCPL-format file', value, bool)
self._surf_source_write = surf_source_write
@confidence_intervals.setter
@ -1032,6 +1036,9 @@ class Settings:
if 'max_particles' in self._surf_source_write:
subelement = ET.SubElement(element, "max_particles")
subelement.text = str(self._surf_source_write['max_particles'])
if 'mcpl' in self._surf_source_write:
subelement = ET.SubElement(element, "mcpl")
subelement.text = str(self._surf_source_write['mcpl']).lower()
def _create_confidence_intervals(self, root):
if self._confidence_intervals is not None:

5
openmc/source.py
View file

@ -223,7 +223,10 @@ class Source:
if self.particle != 'neutron':
element.set("particle", self.particle)
if self.file is not None:
element.set("file", self.file)
if (self.file.endswith('.mcpl')):
element.set("mcpl", self.file)
else:
element.set("file", self.file)
if self.library is not None:
element.set("library", self.library)
if self.parameters is not None:

12
src/settings.cpp
View file

@ -60,7 +60,9 @@ bool run_CE {true};
bool source_latest {false};
bool source_separate {false};
bool source_write {true};
bool source_mcpl_write {true};
bool surf_source_write {false};
bool surf_mcpl_write {false};
bool surf_source_read {false};
bool survival_biasing {false};
bool temperature_multipole {false};
@ -652,6 +654,9 @@ void read_settings_xml()
if (check_for_node(node_sp, "write")) {
source_write = get_node_value_bool(node_sp, "write");
}
if (check_for_node(node_sp, "mcpl")) {
source_write = get_node_value_bool(node_sp, "mcpl");
}
if (check_for_node(node_sp, "overwrite_latest")) {
source_latest = get_node_value_bool(node_sp, "overwrite_latest");
source_separate = source_latest;
@ -682,9 +687,14 @@ void read_settings_xml()
max_surface_particles =
std::stoll(get_node_value(node_ssw, "max_particles"));
}
#ifdef OPENMC_MCPL
if(check_for_node(node_ssw, "mcpl")){
surf_mcpl_write=true;
}
#endif
}
// If source is not seperate and is to be written out in the statepoint file,
// If source is not separate and is to be written out in the statepoint file,
// make sure that the sourcepoint batch numbers are contained in the
// statepoint list
if (!source_separate) {

40
src/simulation.cpp
View file

@ -389,16 +389,33 @@ void finalize_batch()
if (settings::run_mode == RunMode::EIGENVALUE) {
// Write out a separate source point if it's been specified for this batch
if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
settings::source_write && settings::source_separate) {
write_source_point(nullptr);
}
#ifdef OPENMC_MCPL
if(! settings::source_mcpl_write) {
#endif
if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
settings::source_write && settings::source_separate) {
write_source_point(nullptr);
}
// Write a continously-overwritten source point if requested.
if (settings::source_latest) {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
// Write a continously-overwritten source point if requested.
if (settings::source_latest) {
auto filename = settings::path_output + "source.h5";
write_source_point(filename.c_str());
}
#ifdef OPENMC_MCPL
} else {
if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
settings::source_mcpl_write && settings::source_separate) {
write_mcpl_source_point(nullptr);
}
// Write a continously-overwritten source point if requested.
if (settings::source_latest && settings::source_mcpl_write) {
auto filename = settings::path_output + "source.mcpl";
write_mcpl_source_point(filename.c_str());
}
}
#endif
}
// Write out surface source if requested.
@ -407,6 +424,13 @@ void finalize_batch()
auto filename = settings::path_output + "surface_source.h5";
write_source_point(filename.c_str(), true);
}
#ifdef OPENMC_MCPL
if(settings::surf_mcpl_write && simulation::current_batch == settings::n_batches){
auto filename = settings::path_output + "surface_source.mcpl";
write_mcpl_source_point(filename.c_str(), true);
}
#endif
}
void initialize_generation()

27
src/source.cpp
View file

@ -343,17 +343,28 @@ FileSource::FileSource(mcpl_file_t mcpl_file)
const mcpl_particle_t *mcpl_particle;
//extract particle from mcpl-file
mcpl_particle=mcpl_read(mcpl_file);
// check if it is a neutron or a photon. otherwise skip
while ( mcpl_particle->pdgcode!=2112 && mcpl_particle->pdgcode!=22 ) {
// check if it is a neutron, photon, electron, or positron. Otherwise skip.
int pdg=mcpl_particle->pdgcode;
while ( pdg!=2112 && pdg!=22 && pdg!=11 && pdg!=-11) {
mcpl_particle=mcpl_read(mcpl_file);
//should check for file exhaustion This could happen if particles are other than
//neutrons or photons
pdg=mcpl_particle->pdgcode;
//should check for file exhaustion. This could happen if particles are other than
//neutrons, photons, electrons, or positrons.
}
if(mcpl_particle->pdgcode==2112) {
site_.particle=ParticleType::neutron;
} else if (mcpl_particle->pdgcode==22) {
site_.particle=ParticleType::photon;
switch(pdg){
case 2112:
site_.particle=ParticleType::neutron;
break;
case 22:
site_.particle=ParticleType::photon;
break;
case 11:
site_.particle=ParticleType::electron;
break;
case -11:
site_.particle=ParticleType::positron;
break;
}
//particle is good, convert to openmc-formalism

143
src/state_point.cpp
View file

@ -28,6 +28,10 @@
#include "openmc/timer.h"
#include "openmc/vector.h"
#ifdef OPENMC_MCPL
#include <mcpl.h>
#endif
namespace openmc {
extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
@ -599,6 +603,46 @@ void write_source_point(const char* filename, bool surf_source_bank)
file_close(file_id);
}
#ifdef OPENMC_MCPL
void write_mcpl_source_point(const char *filename, bool surf_source_bank)
{
std::string filename_;
if (filename) {
filename_ = filename;
} else {
// Determine width for zero padding
int w = std::to_string(settings::n_max_batches).size();
filename_ = fmt::format("{0}source.{1:0{2}}.mcpl", settings::path_output,
simulation::current_batch, w);
}
mcpl_outfile_t file_id;
std::string line;
if (mpi::master) {
// this must be rewritten: file_open is h5-specific
file_id = mcpl_create_outfile(filename_.c_str());
//write_attribute(file_id, "filetype", "source");
//write header stuff (oopy in xml-files as binary blobs for instance))
if (VERSION_DEV){
line=fmt::format("OpenMC {0}.{1}.{2}-development",VERSION_MAJOR,VERSION_MINOR,VERSION_RELEASE);
} else {
line=fmt::format("OpenMC {0}.{1}.{2}",VERSION_MAJOR,VERSION_MINOR,VERSION_RELEASE);
}
mcpl_hdr_set_srcname(file_id,line.c_str());
}
write_mcpl_source_bank(file_id, surf_source_bank);
if (mpi::master) {
//change this - this is h5 specific
mcpl_close_outfile(file_id);
}
}
#endif
void write_source_bank(hid_t group_id, bool surf_source_bank)
{
hid_t banktype = h5banktype();
@ -715,6 +759,105 @@ void write_source_bank(hid_t group_id, bool surf_source_bank)
H5Tclose(banktype);
}
#ifdef OPENMC_MCPL
void write_mcpl_source_bank(mcpl_outfile_t file_id, bool surf_source_bank)
{
int64_t dims_size = settings::n_particles;
int64_t count_size = simulation::work_per_rank;
// Set vectors for source bank and starting bank index of each process
vector<int64_t>* bank_index = &simulation::work_index;
vector<SourceSite>* source_bank = &simulation::source_bank;
vector<int64_t> surf_source_index_vector;
vector<SourceSite> surf_source_bank_vector;
if(surf_source_bank) {
surf_source_index_vector = calculate_surf_source_size();
dims_size = surf_source_index_vector[mpi::n_procs];
count_size = simulation::surf_source_bank.size();
bank_index = &surf_source_index_vector;
// Copy data in a SharedArray into a vector.
surf_source_bank_vector.resize(count_size);
surf_source_bank_vector.assign(simulation::surf_source_bank.data(),
simulation::surf_source_bank.data() + count_size);
source_bank = &surf_source_bank_vector;
}
if (mpi::master) {
// Particles are writeen to disk from the master node only
// Save source bank sites since the array is overwritten below
#ifdef OPENMC_MPI
vector<SourceSite> temp_source {source_bank->begin(), source_bank->end()};
#endif
//loop over the other nodes and receive data - then write those.
for (int i = 0; i < mpi::n_procs; ++i) {
// number of particles for node node i
size_t count[] {
static_cast<size_t>((*bank_index)[i + 1] - (*bank_index)[i])};
#ifdef OPENMC_MPI
if (i>0)
MPI_Recv(source_bank->data(), count[0], mpi::source_site, i, i,
mpi::intracomm, MPI_STATUS_IGNORE);
#endif
// now write the source_bank data again.
for (vector<SourceSite>::iterator _site=source_bank->begin(); _site!=source_bank->end();_site++){
// particle is now at the iterator
// write it to the mcpl-file
mcpl_particle_t p;
p.position[0]=_site->r.x;
p.position[1]=_site->r.y;
p.position[2]=_site->r.z;
// mcpl requires that the direction vector is unit length
// which is also the case in openmc
p.direction[0]=_site->u.x;
p.direction[1]=_site->u.y;
p.direction[2]=_site->u.z;
// mcpl stores kinetic energy in MeV
p.ekin=_site->E*1e-6;
p.time=_site->time*1e3;
p.weight=_site->wgt;
switch(_site->particle){
case ParticleType::neutron:
p.pdgcode=2112;
break;
case ParticleType::photon:
p.pdgcode=22;
break;
case ParticleType::electron:
p.pdgcode=11;
break;
case ParticleType::positron:
p.pdgcode=-11;
break;
}
mcpl_add_particle(file_id,&p);
}
}
#ifdef OPENMC_MPI
// Restore state of source bank
std::copy(temp_source.begin(), temp_source.end(), source_bank->begin());
#endif
} else {
#ifdef OPENMC_MPI
MPI_Send(source_bank->data(), count_size, mpi::source_site, 0, mpi::rank,
mpi::intracomm);
#endif
}
}
#endif
// Determine member names of a compound HDF5 datatype
std::string dtype_member_names(hid_t dtype_id)
{

2
tests/regression_tests/source_mcpl_file/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
3.017557E-01 3.398770E-03
3.039964E-01 3.654869E-04

6
tests/regression_tests/source_mcpl_file/settings.xml
View file

@ -1,11 +1,15 @@
<?xml version="1.0"?>
<settings>
<state_point batches="10" />
<source_point mcpl="true" separate="true" />
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<mcpl>source.10.mcpl</mcpl>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>

146
tests/regression_tests/source_mcpl_file/test.py
View file

@ -1,65 +1,101 @@
import os
import subprocess
#!/usr/bin/env python
import glob
import os
from tests.testing_harness import TestHarness
from tests.testing_harness import *
class SourceMCPLFileTestHarness(TestHarness):
def execute_test(self):
"""Run OpenMC with the appropriate arguments and check the outputs."""
try:
self._create_input()
self._test_input_created()
self._run_openmc()
self._test_output_created()
results = self._get_results()
self._write_results(results)
self._compare_results()
finally:
self._cleanup()
settings1="""<?xml version="1.0"?>
<settings>
<state_point batches="10" />
<source_point mcpl="true" separate="true" />
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<space type="box">
<parameters>-4 -4 -4 4 4 4</parameters>
</space>
</source>
</settings>
"""
def _create_input(self):
compiled=subprocess.run(['gcc','-o','gen_dummy_mcpl.out','gen_dummy_mcpl.c','-lm','-lmcpl'])
assert compiled==0, 'Could not compile mcpl-file generator code'
subprocess.run(['./gen_dummy_mcpl.out'])
settings2 = """<?xml version="1.0"?>
<settings>
<eigenvalue>
<batches>10</batches>
<inactive>5</inactive>
<particles>1000</particles>
</eigenvalue>
<source>
<mcpl> source.10.{0} </mcpl>
</source>
</settings>
"""
def update_results(self):
"""Update the results_true using the current version of OpenMC."""
try:
self._create_input()
self._test_input_created()
self._run_openmc()
self._test_output_created()
results = self._get_results()
self._write_results(results)
self._overwrite_results()
finally:
self._cleanup()
def _test_input_created(self):
"""Check that the input mcpl.file was generated as it should"""
mcplfile=glob.glob(os.path.join(os.getcwd(),'source.10.mcpl'))
assert len(mcplfile) == 1, 'Either multiple or no mcpl files ' \
'exist.'
assert mcplfile[0].endswith('mcpl'), \
'output file does not end with mcpl.'
class SourceFileTestHarness(TestHarness):
def execute_test(self):
"""Run OpenMC with the appropriate arguments and check the outputs."""
try:
self._run_openmc()
self._test_output_created()
self._run_openmc_restart()
results = self._get_results()
self._write_results(results)
self._compare_results()
finally:
self._cleanup()
def _test_output_created(self):
"""Check that the output files were created"""
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
'exist.'
assert statepoint[0].endswith('h5'), \
'statepoint file does not end with h5.'
def update_results(self):
"""Update the results_true using the current version of OpenMC."""
try:
self._run_openmc()
self._test_output_created()
self._run_openmc_restart()
results = self._get_results()
self._write_results(results)
self._overwrite_results()
finally:
self._cleanup()
def _cleanup(self):
super()._cleanup()
source_mcpl=glob.glob(os.path.join(os.getcwd(),'source*.mcpl'))
for f in source_mcpl:
if (os.path.exists(f)):
os.remove(f)
def _test_output_created(self):
"""Make sure statepoint and source files have been created."""
statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))
assert len(statepoint) == 1, 'Either multiple or no statepoint files ' \
'exist.'
assert statepoint[0].endswith('h5'), \
'Statepoint file is not a HDF5 file.'
def test_mcpl_source_file():
harness = SourceMCPLFileTestHarness('source.10.mcpl')
harness.main()
source = glob.glob(os.path.join(os.getcwd(), 'source.10.mcpl*'))
assert len(source) == 1, 'Either multiple or no source files exist.'
assert source[0].endswith('mcpl') or source[0].endswith('mcpl.gz'), \
'Source file is not a MCPL file.'
def _run_openmc_restart(self):
# Get the name of the source file.
source = glob.glob(os.path.join(os.getcwd(), 'source.10.*'))
# Write the new settings.xml file.
with open('settings.xml','w') as fh:
fh.write(settings2.format(source[0].split('.')[-1]))
# Run OpenMC.
self._run_openmc()
def _cleanup(self):
TestHarness._cleanup(self)
output = glob.glob(os.path.join(os.getcwd(), 'source.*'))
#for f in output:
# if os.path.exists(f):
# os.remove(f)
with open('settings.xml','w') as fh:
fh.write(settings1)
def test_source_file():
harness = SourceFileTestHarness('statepoint.10.h5')
harness.main()