From f8925b25bee777a5b68919e57f9932da7f6ba539 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 7 Feb 2019 10:05:43 -0600 Subject: [PATCH 01/10] Convert transport() and cross_surface() to C++ --- CMakeLists.txt | 1 - include/openmc/cell.h | 8 + include/openmc/error.h | 2 +- include/openmc/particle.h | 6 + include/openmc/physics.h | 2 +- include/openmc/physics_mg.h | 2 +- include/openmc/tallies/derivative.h | 7 + include/openmc/tallies/tally_scoring.h | 10 + src/cell.cpp | 6 +- src/geometry.F90 | 24 -- src/material.cpp | 5 - src/material_header.F90 | 8 - src/particle.cpp | 509 ++++++++++++++++++++++- src/particle_restart.F90 | 3 +- src/simulation.cpp | 7 +- src/surface.cpp | 6 +- src/tallies/derivative.cpp | 9 +- src/tallies/tally.F90 | 71 ---- src/tallies/tally_scoring.cpp | 31 +- src/track_output.F90 | 10 +- src/tracking.F90 | 539 ------------------------- 21 files changed, 566 insertions(+), 700 deletions(-) delete mode 100644 src/tracking.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 5a204fb66..6bc51e812 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -337,7 +337,6 @@ add_library(libopenmc SHARED src/stl_vector.F90 src/string.F90 src/surface_header.F90 - src/tracking.F90 src/track_output.F90 src/vector_header.F90 src/xml_interface.F90 diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c90cb7f69..ba202c343 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -195,5 +195,13 @@ public: }; #endif +//============================================================================== +// Non-member functions +//============================================================================== + +#ifdef DAGMC +int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed); +#endif + } // namespace openmc #endif // OPENMC_CELL_H diff --git a/include/openmc/error.h b/include/openmc/error.h index ee7168a4e..206c9b5b2 100644 --- a/include/openmc/error.h +++ b/include/openmc/error.h @@ -49,7 +49,7 @@ void warning(const std::stringstream& message) warning(message.str()); } -void write_message(const std::string& message, int level); +void write_message(const std::string& message, int level=0); inline void write_message(const std::stringstream& message, int level) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 84e628fe2..a81b3e77a 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -148,6 +148,12 @@ extern "C" { //! \param src Source site data void from_source(const Bank* src); + //! Transport the particle + void transport(); + + //! Cross a surface + void cross_surface(); + //! mark a particle as lost and create a particle restart file //! \param message A warning message to display void mark_as_lost(const char* message); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 5d480e0bf..0c2b399d1 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -14,7 +14,7 @@ namespace openmc { //============================================================================== //! Sample a nuclide and reaction and then calls the appropriate routine -extern "C" void collision(Particle* p); +void collision(Particle* p); //! Samples an incident neutron reaction void sample_neutron_reaction(Particle* p); diff --git a/include/openmc/physics_mg.h b/include/openmc/physics_mg.h index 8febad19a..511a0658b 100644 --- a/include/openmc/physics_mg.h +++ b/include/openmc/physics_mg.h @@ -12,7 +12,7 @@ namespace openmc { //! \brief samples particle behavior after a collision event. //! \param p Particle to operate on -extern "C" void +void collision_mg(Particle* p); //! \brief samples a reaction type. diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 3237fb3ba..21df00624 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -35,6 +35,13 @@ void apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, double atom_density, int score_bin, double& score); +void score_collision_derivative(const Particle* p); + +void score_track_derivative(const Particle* p, double distance); + +//! Set the flux derivatives on differential tallies to zero. +void zero_flux_derivs(); + } // namespace openmc //============================================================================== diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 8aa659325..12016956c 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -50,6 +50,16 @@ private: // Non-member functions //============================================================================== +void score_collision_tally(const Particle* p); + +void score_analog_tally_ce(const Particle* p); + +void score_analog_tally_mg(const Particle* p); + +void score_tracklength_tally(const Particle* p, double distance); + +void score_surface_tally(const Particle* p, const std::vector& tallies); + } // namespace openmc #endif // OPENMC_TALLIES_TALLY_SCORING_H diff --git a/src/cell.cpp b/src/cell.cpp index 642059506..56635309f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -805,12 +805,12 @@ extern "C" { } #ifdef DAGMC - int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed ) + int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed) { moab::EntityHandle surf = - surf_xed->dagmc_ptr_->entity_by_id(2,surf_xed->id_); + surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_); moab::EntityHandle vol = - cur_cell->dagmc_ptr_->entity_by_id(3,cur_cell->id_); + cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_); moab::EntityHandle new_vol; cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol); diff --git a/src/geometry.F90 b/src/geometry.F90 index 7456644e6..02ababae6 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -50,34 +50,10 @@ module geometry integer(C_INT), intent(out) :: lattice_translation(3) integer(C_INT), intent(out) :: next_level end subroutine distance_to_boundary - -#ifdef DAGMC - - function next_cell_c(current_cell, surface_crossed) & - bind(C, name="next_cell") result(new_cell) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: current_cell - type(C_PTR), intent(in), value :: surface_crossed - integer(C_INT32_T) :: new_cell - end function next_cell_c - -#endif - end interface contains -#ifdef DAGMC - - function next_cell(c, s) result(new_cell) - type(Cell), intent(in) :: c - type(Surface), intent(in) :: s - integer :: new_cell - new_cell = next_cell_c(c%ptr, s%ptr) - end function next_cell - -#endif - !=============================================================================== ! FIND_CELL determines what cell a source particle is in within a particular ! universe. If the base universe is passed, the particle should be found as long diff --git a/src/material.cpp b/src/material.cpp index 510cea966..41c9d0a01 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -1135,11 +1135,6 @@ extern "C" { return model::materials[i_mat - 1]->density_gpcc_; } - void material_calculate_xs(const Particle* p) - { - model::materials[p->material - 1]->calculate_xs(*p); - } - void free_memory_material() { for (Material *mat : model::materials) {delete mat;} diff --git a/src/material_header.F90 b/src/material_header.F90 index cda79ec60..dc1e68944 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -2,12 +2,9 @@ module material_header use, intrinsic :: ISO_C_BINDING - use particle_header, only: Particle - implicit none private - public :: material_calculate_xs public :: material_id public :: material_nuclide public :: material_nuclide_size @@ -21,11 +18,6 @@ module material_header integer(C_INT32_T) :: id end function - subroutine material_calculate_xs(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - function material_nuclide(i_mat, idx) bind(C) result(nuc) import C_INT32_T, C_INT integer(C_INT32_T), value :: i_mat diff --git a/src/particle.cpp b/src/particle.cpp index e5a38a983..26c4e8a7f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,16 +1,29 @@ #include "openmc/particle.h" -#include +#include // copy, min +#include // log, abs, copysign #include #include "openmc/bank.h" #include "openmc/capi.h" +#include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/error.h" +#include "openmc/geometry.h" #include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" +#include "openmc/nuclide.h" +#include "openmc/physics.h" +#include "openmc/physics_mg.h" +#include "openmc/random_lcg.h" #include "openmc/settings.h" +#include "openmc/surface.h" #include "openmc/simulation.h" +#include "openmc/tallies/derivative.h" +#include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" namespace openmc { @@ -119,6 +132,500 @@ Particle::from_source(const Bank* src) last_E = E; } +extern "C" void initialize_particle_track(); +extern "C" void write_particle_track(const Particle*); +extern "C" void add_particle_track(); +extern "C" void finalize_particle_track(const Particle*); + +void +Particle::transport() +{ + // Display message if high verbosity or trace is on + if (settings::verbosity >= 9 || simulation::trace) { + write_message("Simulating Particle " + std::to_string(id)); + } + + // Initialize number of events to zero + int n_event = 0; + + // Add paricle's starting weight to count for normalizing tallies later +#pragma omp atomic + total_weight += wgt; + + // Force calculation of cross-sections by setting last energy to zero + if (settings::run_CE) { + for (int i = 0; i < data::nuclides.size(); ++i) { + simulation::micro_xs[i].last_E = 0.0; + } + } + + // Prepare to write out particle track. + if (write_track) initialize_particle_track(); + + // Every particle starts with no accumulated flux derivative. + if (!model::active_tallies.empty()) zero_flux_derivs(); + + while (true) { + // Set the random number stream + if (type == static_cast(ParticleType::neutron)) { + prn_set_stream(STREAM_TRACKING); + } else { + prn_set_stream(STREAM_PHOTON); + } + + // Store pre-collision particle properties + last_wgt = wgt; + last_E = E; + std::copy(coord[0].uvw, coord[0].uvw + 3, last_uvw); + std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz); + + // If the cell hasn't been determined based on the particle's location, + // initiate a search for the current cell. This generally happens at the + // beginning of the history and again for any secondary particles + if (coord[n_coord - 1].cell == C_NONE) { + if (!find_cell(this, false)) { + this->mark_as_lost("Could not find the cell containing particle " + + std::to_string(id)); + return; + } + + // set birth cell attribute + if (cell_born == C_NONE) cell_born = coord[n_coord - 1].cell; + } + + // Write particle track. + if (write_track) write_particle_track(this); + + if (settings::check_overlaps) check_cell_overlap(this); + + // Calculate microscopic and macroscopic cross sections + if (material != MATERIAL_VOID) { + if (settings::run_CE) { + if (material != last_material || sqrtkT != last_sqrtkT) { + // If the material is the same as the last material and the + // temperature hasn't changed, we don't need to lookup cross + // sections again. + model::materials[material - 1]->calculate_xs(*this); + } + } else { + // Get the MG data + calculate_xs_c(material, g, sqrtkT, coord[n_coord-1].uvw, + simulation::material_xs.total, simulation::material_xs.absorption, + simulation::material_xs.nu_fission); + + // Finally, update the particle group while we have already checked + // for if multi-group + last_g = g; + } + } else { + simulation::material_xs.total = 0.0; + simulation::material_xs.absorption = 0.0; + simulation::material_xs.fission = 0.0; + simulation::material_xs.nu_fission = 0.0; + } + + // Find the distance to the nearest boundary + double d_boundary; + int surface_crossed; + int lattice_translation[3]; + int next_level; + distance_to_boundary(this, &d_boundary, &surface_crossed, + lattice_translation, &next_level); + + // Sample a distance to collision + double d_collision; + if (type == static_cast(ParticleType::electron) || + type == static_cast(ParticleType::positron)) { + d_collision = 0.0; + } else if (simulation::material_xs.total == 0.0) { + d_collision = INFINITY; + } else { + d_collision = -std::log(prn()) / simulation::material_xs.total; + } + + // Select smaller of the two distances + double distance = std::min(d_boundary, d_collision); + + // Advance particle + for (int j = 0; j < n_coord; ++j) { + // TODO: use Position + coord[j].xyz[0] += distance * coord[j].uvw[0]; + coord[j].xyz[1] += distance * coord[j].uvw[1]; + coord[j].xyz[2] += distance * coord[j].uvw[2]; + } + + // Score track-length tallies + if (!model::active_tracklength_tallies.empty()) { + score_tracklength_tally(this, distance); + } + + // Score track-length estimate of k-eff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type == static_cast(ParticleType::neutron)) { + global_tally_tracklength += wgt * distance * simulation::material_xs.nu_fission; + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) { + score_track_derivative(this, distance); + } + + if (d_collision > d_boundary) { + // ==================================================================== + // PARTICLE CROSSES SURFACE + + if (next_level > 0) n_coord = next_level; + + // Saving previous cell data + for (int j = 0; j < n_coord; ++j) { + last_cell[j] = coord[j].cell; + } + last_n_coord = n_coord; + + if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || + lattice_translation[2] != 0) { + // Particle crosses lattice boundary + surface = ERROR_INT; + cross_lattice(this, lattice_translation); + event = EVENT_LATTICE; + } else { + // Particle crosses surface + surface = surface_crossed; + this->cross_surface(); + event = EVENT_SURFACE; + } + // Score cell to cell partial currents + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } + } else { + // ==================================================================== + // PARTICLE HAS COLLISION + + // Score collision estimate of keff + if (settings::run_mode == RUN_MODE_EIGENVALUE && + type == static_cast(ParticleType::neutron)) { + global_tally_collision += wgt * simulation::material_xs.nu_fission + / simulation::material_xs.total; + } + + // Score surface current tallies -- this has to be done before the collision + // since the direction of the particle will change and we need to use the + // pre-collision direction to figure out what mesh surfaces were crossed + + if (!model::active_meshsurf_tallies.empty()) + score_surface_tally(this, model::active_meshsurf_tallies); + + // Clear surface component + surface = ERROR_INT; + + if (settings::run_CE) { + collision(this); + } else { + collision_mg(this); + } + + // Score collision estimator tallies -- this is done after a collision + // has occurred rather than before because we need information on the + // outgoing energy for any tallies with an outgoing energy filter + if (!model::active_collision_tallies.empty()) score_collision_tally(this); + if (!model::active_analog_tallies.empty()) { + if (settings::run_CE) { + score_analog_tally_ce(this); + } else { + score_analog_tally_mg(this); + } + } + + // Reset banked weight during collision + n_bank = 0; + wgt_bank = 0.0; + for (int& v : n_delayed_bank) v = 0; + + // Reset fission logical + fission = false; + + // Save coordinates for tallying purposes + std::copy(coord[0].xyz, coord[0].xyz + 3, last_xyz_current); + + // Set last material to none since cross sections will need to be + // re-evaluated + last_material = F90_NONE; + + // Set all uvws to base level -- right now, after a collision, only the + // base level uvws are changed + for (int j = 0; j < n_coord - 1; ++j) { + if (coord[j + 1].rotated) { + // If next level is rotated, apply rotation matrix + // FIXME + // coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % & + // rotation_matrix, coord[j].uvw) + } else { + // Otherwise, copy this level's direction + std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw); + } + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) score_collision_derivative(this); + } + + // If particle has too many events, display warning and kill it + ++n_event; + if (n_event == MAX_EVENTS) { + if (mpi::master) warning("Particle " + std::to_string(id) + + " underwent maximum number of events."); + alive = false; + } + + // Check for secondary particles if this particle is dead + if (!alive) { + // If no secondary particles, break out of event loop + if (n_secondary == 0) break; + + this->from_source(&secondary_bank[n_secondary - 1]); + --n_secondary; + n_event = 0; + + // Enter new particle in particle track file + if (write_track) add_particle_track(); + } + } + + // Finish particle track output. + if (write_track) { + write_particle_track(this); + finalize_particle_track(this); + } +} + +void +Particle::cross_surface() +{ + int i_surface = std::abs(surface); + // TODO: off-by-one + const auto& surf {model::surfaces[i_surface - 1]}; + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Crossing surface " + std::to_string(surf->id_)); + } + + if (surf->bc_ == BC_VACUUM && (settings::run_mode != RUN_MODE_PLOTTING)) { + // ======================================================================= + // PARTICLE LEAKS OUT OF PROBLEM + + // Kill particle + alive = false; + + // Score any surface current tallies -- note that the particle is moved + // forward slightly so that if the mesh boundary is on the surface, it is + // still processed + + if (!model::active_meshsurf_tallies.empty()) { + // TODO: Find a better solution to score surface currents than + // physically moving the particle forward slightly + + // TODO: Use Position + coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + } + + // Score to global leakage tally + global_tally_leakage += wgt; + + // Display message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Leaked out of surface " + std::to_string(surf->id_)); + } + return; + + } else if (surf->bc_ == BC_REFLECT && (settings::run_mode != RUN_MODE_PLOTTING)) { + // ======================================================================= + // PARTICLE REFLECTS FROM SURFACE + + // Do not handle reflective boundary conditions on lower universes + if (n_coord != 1) { + this->mark_as_lost("Cannot reflect particle " + std::to_string(id) + + " off surface in a lower universe."); + return; + } + + // Score surface currents since reflection causes the direction of the + // particle to change. For surface filters, we need to score the tallies + // twice, once before the particle's surface attribute has changed and + // once after. For mesh surface filters, we need to artificially move + // the particle slightly back in case the surface crossing is coincident + // with a mesh boundary + + if (!model::active_surface_tallies.empty()) { + score_surface_tally(this, model::active_surface_tallies); + } + + + if (!model::active_meshsurf_tallies.empty()) { + Position r {coord[0].xyz}; + coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + } + + // Reflect particle off surface + Direction u = surf->reflect(coord[0].xyz, coord[0].uvw); + + // Make sure new particle direction is normalized + double norm = u.norm(); + coord[0].uvw[0] = u.x/norm; + coord[0].uvw[1] = u.y/norm; + coord[0].uvw[2] = u.z/norm; + + // Reassign particle's cell and surface + coord[0].cell = last_cell[last_n_coord - 1]; + surface = -surface; + + // If a reflective surface is coincident with a lattice or universe + // boundary, it is necessary to redetermine the particle's coordinates in + // the lower universes. + + n_coord = 1; + if (!find_cell(this, true)) { + this->mark_as_lost("Couldn't find particle after reflecting from surface " + + std::to_string(surf->id_) + "."); + return; + } + + // Set previous coordinate going slightly past surface crossing + last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT*coord[0].uvw[0]; + last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT*coord[0].uvw[1]; + last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT*coord[0].uvw[2]; + + // Diagnostic message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Reflected from surface " + std::to_string(surf->id_)); + } + return; + + } else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) { + // ======================================================================= + // PERIODIC BOUNDARY + + // Do not handle periodic boundary conditions on lower universes + if (n_coord != 1) { + this->mark_as_lost("Cannot transfer particle " + std::to_string(id) + + " across surface in a lower universe. Boundary conditions must be " + "applied to root universe."); + return; + } + + // Score surface currents since reflection causes the direction of the + // particle to change -- artificially move the particle slightly back in + // case the surface crossing is coincident with a mesh boundary + if (!model::active_meshsurf_tallies.empty()) { + Position r {coord[0].xyz}; + coord[0].xyz[0] -= TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] -= TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] -= TINY_BIT * coord[0].uvw[2]; + score_surface_tally(this, model::active_meshsurf_tallies); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + } + + // Get a pointer to the partner periodic surface. Offset the index to + // correct for C vs. Fortran indexing. + auto surf_p = dynamic_cast(surf); + if (surf_p) { + auto other = dynamic_cast( + model::surfaces[surf_p->i_periodic_]); + + // Adjust the particle's location and direction. + Position r {coord[0].xyz}; + Direction u {coord[0].uvw}; + bool rotational = other->periodic_translate(surf_p, r, u); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + std::copy(&u.x, &u.x + 3, coord[0].uvw); + + // Reassign particle's surface + // TODO: off-by-one + surface = rotational ? + surf_p->i_periodic_ + 1 : + std::copysign(surf_p->i_periodic_ + 1, surface); + } + + // Figure out what cell particle is in now + n_coord = 1; + + if (!find_cell(this, true)) { + this->mark_as_lost("Couldn't find particle after hitting periodic " + "boundary on surface " + std::to_string(surf->id_) + "."); + return; + } + + // Set previous coordinate going slightly past surface crossing + last_xyz_current[0] = coord[0].xyz[0] + TINY_BIT * coord[0].uvw[0]; + last_xyz_current[1] = coord[0].xyz[1] + TINY_BIT * coord[0].uvw[1]; + last_xyz_current[2] = coord[0].xyz[2] + TINY_BIT * coord[0].uvw[2]; + + // Diagnostic message + if (settings::verbosity >= 10 || simulation::trace) { + write_message(" Hit periodic boundary on surface " + + std::to_string(surf->id_)); + } + return; + } + + // ========================================================================== + // SEARCH NEIGHBOR LISTS FOR NEXT CELL + +#ifdef DAGMC + if (settings::dagmc) { + int32_t i_cell = next_cell(model::cells[last_cell[0]], + model::surfaces[std::abs(surface)]); + // save material and temp + last_material = material; + last_sqrtkT = sqrtKT; + // set new cell value + coord[0].cell = i_cell + cell_instance = 0; + material = model::cells[i_cell]->material[0]; + sqrtKT = model::cells[i_cell]->sqrtKT[0]; + return + } +#endif + + if (find_cell(this, true)) return; + + // ========================================================================== + // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS + + // Remove lower coordinate levels and assignment of surface + surface = ERROR_INT; + n_coord = 1; + bool found = find_cell(this, false); + + if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) { + // If a cell is still not found, there are two possible causes: 1) there is + // a void in the model, and 2) the particle hit a surface at a tangent. If + // the particle is really traveling tangent to a surface, if we move it + // forward a tiny bit it should fix the problem. + + n_coord = 1; + coord[0].xyz[0] += TINY_BIT * coord[0].uvw[0]; + coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; + coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; + + // Couldn't find next cell anywhere// This probably means there is an actual + // undefined region in the geometry. + + if (!find_cell(this, false)) { + this->mark_as_lost("After particle " + std::to_string(id) + + " crossed surface " + std::to_string(surf->id_) + + " it could not be located in any cell and it did not leak."); + return; + } + } +} + void Particle::mark_as_lost(const char* message) { diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index dd1c89669..31ee5ff66 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -74,7 +74,8 @@ contains call set_particle_seed(particle_seed) ! Transport neutron - call transport(p) + ! FIXME + !call transport(p) ! Write output if particle made it call print_particle(p) diff --git a/src/simulation.cpp b/src/simulation.cpp index 7c7bbc4e6..acc488353 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -29,11 +29,9 @@ namespace openmc { // data/functions from Fortran side extern "C" void accumulate_tallies(); extern "C" void allocate_tally_results(); -extern "C" void init_tally_routines(); extern "C" void setup_active_tallies(); extern "C" void simulation_init_f(); extern "C" void simulation_finalize_f(); -extern "C" void transport(Particle* p); extern "C" void write_tallies(); } // namespace openmc @@ -76,9 +74,6 @@ int openmc_simulation_init() simulation::filter_matches.resize(model::tally_filters.size()); } - // Set up tally procedure pointers - init_tally_routines(); - // Allocate source bank, and for eigenvalue simulations also allocate the // fission bank allocate_banks(); @@ -213,7 +208,7 @@ int openmc_next_batch(int* status) initialize_history(&p, simulation::current_work); // transport particle - transport(&p); + p.transport(); } // Accumulate time for transport diff --git a/src/surface.cpp b/src/surface.cpp index 3ffbb3277..8b9b20ea1 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -340,7 +340,7 @@ bool SurfaceXPlane::periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const { Direction other_n = other->normal(r); - if (other_n.x == 1 and other_n.y == 0 and other_n.z == 0) { + if (other_n.x == 1 && other_n.y == 0 && other_n.z == 0) { r.x = x0_; return false; } else { @@ -397,11 +397,11 @@ void SurfaceYPlane::to_hdf5_inner(hid_t group_id) const write_dataset(group_id, "coefficients", coeffs); } -bool SurfaceYPlane::periodic_translate(const PeriodicSurface *other, +bool SurfaceYPlane::periodic_translate(const PeriodicSurface* other, Position& r, Direction& u) const { Direction other_n = other->normal(r); - if (other_n.x == 0 and other_n.y == 1 and other_n.z == 0) { + if (other_n.x == 0 && other_n.y == 1 && other_n.z == 0) { // The periodic partner is also aligned along y. Just change the y coord. r.y = y0_; return false; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index 600cd2876..c3bf8f1bf 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -571,7 +571,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, //! Adjust diff tally flux derivatives for a particle tracking event. -extern "C" void +void score_track_derivative(const Particle* p, double distance) { // A void material cannot be perturbed so it will not affect flux derivatives. @@ -627,8 +627,7 @@ score_track_derivative(const Particle* p, double distance) //! absorption will never be tallied. The paricle will be killed before any //! further tallies are scored. -extern "C" void -score_collision_derivative(const Particle* p) +void score_collision_derivative(const Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. if (p->material == MATERIAL_VOID) return; @@ -698,9 +697,7 @@ score_collision_derivative(const Particle* p) } } -//! Set the flux derivatives on differential tallies to zero. -extern "C" void -zero_flux_derivs() +void zero_flux_derivs() { for (auto& deriv : model::tally_derivs) deriv.flux_deriv = 0.; } diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 96b2fff06..003eb438d 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -11,87 +11,16 @@ module tally use message_passing use mgxs_interface use nuclide_header - use particle_header, only: LocalCoord, Particle use settings use simulation_header - use string, only: to_str use tally_derivative_header use tally_filter use tally_header implicit none - procedure(score_analog_tally_), pointer :: score_analog_tally => null() - - abstract interface - subroutine score_analog_tally_(p) - import Particle - type(Particle), intent(in) :: p - end subroutine score_analog_tally_ - end interface - - interface - subroutine score_analog_tally_ce(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - - subroutine score_analog_tally_mg(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - - subroutine score_tracklength_tally(p, distance) bind(C) - import Particle, C_DOUBLE - type(Particle) :: p - real(C_DOUBLE), value :: distance - end subroutine - - subroutine score_collision_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_meshsurface_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_surface_tally(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine score_track_derivative(p, distance) bind(C) - import Particle, C_DOUBLE - type(Particle) :: p - real(C_DOUBLE), value :: distance - end subroutine - - subroutine score_collision_derivative(p) bind(C) - import Particle - type(Particle) :: p - end subroutine - - subroutine zero_flux_derivs() bind(C) - end subroutine - end interface - contains -!=============================================================================== -! INIT_TALLY_ROUTINES Sets the procedure pointers needed for minimizing code -! with the CE and MG modes. -!=============================================================================== - - subroutine init_tally_routines() bind(C) - if (run_CE) then - score_analog_tally => score_analog_tally_ce - else - score_analog_tally => score_analog_tally_mg - end if - end subroutine init_tally_routines - !=============================================================================== ! ACCUMULATE_TALLIES accumulates the sum of the contributions from each history ! within the batch to a new random variable diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 8c509a0aa..fd96d9d85 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1976,8 +1976,7 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, // //! Analog tallies ar etriggered at every collision, not every event. -extern "C" void -score_analog_tally_ce(const Particle* p) +void score_analog_tally_ce(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2040,8 +2039,7 @@ score_analog_tally_ce(const Particle* p) // //! Analog tallies ar etriggered at every collision, not every event. -extern "C" void -score_analog_tally_mg(const Particle* p) +void score_analog_tally_mg(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2096,7 +2094,7 @@ score_analog_tally_mg(const Particle* p) //! collision) and thus cannot be done for tallies that require post-collision //! information. -extern "C" void +void score_tracklength_tally(const Particle* p, double distance) { // Determine the tracklength estimate of the flux @@ -2172,8 +2170,7 @@ score_tracklength_tally(const Particle* p, double distance) //! reactions that we didn't sample. It is assumed the material is not void //! since collisions do not occur in voids. -extern "C" void -score_collision_tally(const Particle* p) +void score_collision_tally(const Particle* p) { // Determine the collision estimate of the flux double flux; @@ -2243,8 +2240,8 @@ score_collision_tally(const Particle* p) //! Score surface or mesh-surface tallies for particle currents. -static void -score_surface_tally_inner(const Particle* p, const std::vector& tallies) +void +score_surface_tally(const Particle* p, const std::vector& tallies) { // No collision, so no weight change when survival biasing double flux = p->wgt; @@ -2289,20 +2286,4 @@ score_surface_tally_inner(const Particle* p, const std::vector& tallies) match.bins_present_ = false; } -//! Score mesh-surface tallies for particle currents. - -extern "C" void -score_meshsurface_tally(const Particle* p) -{ - score_surface_tally_inner(p, model::active_meshsurf_tallies); -} - -//! Score surface tallies for particle currents. - -extern "C" void -score_surface_tally(const Particle* p) -{ - score_surface_tally_inner(p, model::active_surface_tallies); -} - } // namespace openmc diff --git a/src/track_output.F90 b/src/track_output.F90 index 95af9e6b0..cf1625ca3 100644 --- a/src/track_output.F90 +++ b/src/track_output.F90 @@ -5,6 +5,8 @@ module track_output + use, intrinsic :: ISO_C_BINDING + use constants use hdf5_interface use particle_header, only: Particle @@ -34,7 +36,7 @@ contains ! information !=============================================================================== - subroutine initialize_particle_track() + subroutine initialize_particle_track() bind(C) allocate(tracks(1)) end subroutine initialize_particle_track @@ -42,7 +44,7 @@ contains ! WRITE_PARTICLE_TRACK copies particle position to an array. !=============================================================================== - subroutine write_particle_track(p) + subroutine write_particle_track(p) bind(C) type(Particle), intent(in) :: p real(8), allocatable :: new_coords(:, :) @@ -71,7 +73,7 @@ contains ! secondary particle !=============================================================================== - subroutine add_particle_track() + subroutine add_particle_track() bind(C) type(TrackCoordinates), allocatable :: new_tracks(:) integer :: i @@ -92,7 +94,7 @@ contains ! FINALIZE_PARTICLE_TRACK writes the particle track array to disk. !=============================================================================== - subroutine finalize_particle_track(p) + subroutine finalize_particle_track(p) bind(C) type(Particle), intent(in) :: p integer :: i diff --git a/src/tracking.F90 b/src/tracking.F90 deleted file mode 100644 index d3c85bd79..000000000 --- a/src/tracking.F90 +++ /dev/null @@ -1,539 +0,0 @@ -module tracking - - use, intrinsic :: ISO_C_BINDING - - use constants - use error, only: warning, write_message - use geometry_header, only: cells - use geometry, only: find_cell, distance_to_boundary, cross_lattice,& - check_cell_overlap -#ifdef DAGMC - use geometry, only: next_cell -#endif - - use material_header, only: material_calculate_xs - use message_passing - use mgxs_interface - use nuclide_header - use particle_header - use random_lcg, only: prn, prn_set_stream - use settings - use simulation_header - use string, only: to_str - use surface_header - use tally_header - use tally, only: score_analog_tally, score_tracklength_tally, & - score_collision_tally, score_surface_tally, & - score_meshsurface_tally, & - score_track_derivative, zero_flux_derivs, & - score_collision_derivative - use track_output, only: initialize_particle_track, write_particle_track, & - add_particle_track, finalize_particle_track - - implicit none - - interface - subroutine collision(p) bind(C) - import Particle - type(Particle), intent(inout) :: p - end subroutine - - subroutine collision_mg(p) bind(C) - import Particle, C_DOUBLE - type(Particle), intent(inout) :: p - end subroutine collision_mg - end interface - -contains - -!=============================================================================== -! TRANSPORT encompasses the main logic for moving a particle through geometry. -!=============================================================================== - - subroutine transport(p) bind(C) - - type(Particle), intent(inout) :: p - - integer :: j ! coordinate level - integer :: next_level ! next coordinate level to check - integer :: surface_crossed ! surface which particle is on - integer :: lattice_translation(3) ! in-lattice translation vector - integer :: n_event ! number of collisions/crossings - real(8) :: d_boundary ! distance to nearest boundary - real(8) :: d_collision ! sampled distance to collision - real(8) :: distance ! distance particle travels - logical :: found_cell ! found cell which particle is in? - - ! Display message if high verbosity or trace is on - if (verbosity >= 9 .or. trace) then - call write_message("Simulating Particle " // trim(to_str(p % id))) - end if - - ! Initialize number of events to zero - n_event = 0 - - ! Add paricle's starting weight to count for normalizing tallies later -!$omp atomic - total_weight = total_weight + p % wgt - - ! Force calculation of cross-sections by setting last energy to zero - if (run_CE) then - micro_xs % last_E = ZERO - end if - - ! Prepare to write out particle track. - if (p % write_track) then - call initialize_particle_track() - endif - - ! Every particle starts with no accumulated flux derivative. - if (active_tallies_size() > 0) call zero_flux_derivs() - - EVENT_LOOP: do - ! Set the random number stream - if (p % type == NEUTRON) then - call prn_set_stream(STREAM_TRACKING) - else - call prn_set_stream(STREAM_PHOTON) - end if - - ! Store pre-collision particle properties - p % last_wgt = p % wgt - p % last_E = p % E - p % last_uvw = p % coord(1) % uvw - p % last_xyz = p % coord(1) % xyz - - ! If the cell hasn't been determined based on the particle's location, - ! initiate a search for the current cell. This generally happens at the - ! beginning of the history and again for any secondary particles - if (p % coord(p % n_coord) % cell == C_NONE) then - call find_cell(p, found_cell) - if (.not. found_cell) then - call particle_mark_as_lost(p, "Could not find the cell containing" & - // " particle " // trim(to_str(p %id))) - return - end if - - ! set birth cell attribute - if (p % cell_born == C_NONE) p % cell_born = p % coord(p % n_coord) % cell - end if - - ! Write particle track. - if (p % write_track) call write_particle_track(p) - - if (check_overlaps) call check_cell_overlap(p) - - ! Calculate microscopic and macroscopic cross sections - if (p % material /= MATERIAL_VOID) then - if (run_CE) then - if (p % material /= p % last_material .or. & - p % sqrtkT /= p % last_sqrtkT) then - ! If the material is the same as the last material and the - ! temperature hasn't changed, we don't need to lookup cross - ! sections again. - call material_calculate_xs(p) - end if - else - ! Get the MG data - call calculate_xs_c(p % material, p % g, p % sqrtkT, & - p % coord(p % n_coord) % uvw, material_xs % total, & - material_xs % absorption, material_xs % nu_fission) - - - ! Finally, update the particle group while we have already checked - ! for if multi-group - p % last_g = p % g - end if - else - material_xs % total = ZERO - material_xs % absorption = ZERO - material_xs % fission = ZERO - material_xs % nu_fission = ZERO - end if - - ! Find the distance to the nearest boundary - call distance_to_boundary(p, d_boundary, surface_crossed, & - lattice_translation, next_level) - - ! Sample a distance to collision - if (p % type == ELECTRON .or. p % type == POSITRON) then - d_collision = ZERO - else if (material_xs % total == ZERO) then - d_collision = INFINITY - else - d_collision = -log(prn()) / material_xs % total - end if - - ! Select smaller of the two distances - distance = min(d_boundary, d_collision) - - ! Advance particle - do j = 1, p % n_coord - p % coord(j) % xyz = p % coord(j) % xyz + distance * p % coord(j) % uvw - end do - - ! Score track-length tallies - if (active_tracklength_tallies_size() > 0) then - call score_tracklength_tally(p, distance) - end if - - ! Score track-length estimate of k-eff - if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then - global_tally_tracklength = global_tally_tracklength + p % wgt * & - distance * material_xs % nu_fission - end if - - ! Score flux derivative accumulators for differential tallies. - if (active_tallies_size() > 0) call score_track_derivative(p, distance) - - if (d_collision > d_boundary) then - ! ==================================================================== - ! PARTICLE CROSSES SURFACE - - if (next_level > 0) p % n_coord = next_level - - ! Saving previous cell data - do j = 1, p % n_coord - p % last_cell(j) = p % coord(j) % cell - end do - p % last_n_coord = p % n_coord - - if (any(lattice_translation /= 0)) then - ! Particle crosses lattice boundary - p % surface = ERROR_INT - call cross_lattice(p, lattice_translation) - p % event = EVENT_LATTICE - else - ! Particle crosses surface - p % surface = surface_crossed - - call cross_surface(p) - p % event = EVENT_SURFACE - end if - ! Score cell to cell partial currents - if (active_surface_tallies_size() > 0) call score_surface_tally(p) - else - ! ==================================================================== - ! PARTICLE HAS COLLISION - - ! Score collision estimate of keff - if (run_mode == MODE_EIGENVALUE .and. p % type == NEUTRON) then - global_tally_collision = global_tally_collision + p % wgt * & - material_xs % nu_fission / material_xs % total - end if - - ! Score surface current tallies -- this has to be done before the collision - ! since the direction of the particle will change and we need to use the - ! pre-collision direction to figure out what mesh surfaces were crossed - - if (active_meshsurf_tallies_size() > 0) call score_meshsurface_tally(p) - - ! Clear surface component - p % surface = ERROR_INT - - if (run_CE) then - call collision(p) - else - call collision_mg(p) - end if - - ! Score collision estimator tallies -- this is done after a collision - ! has occurred rather than before because we need information on the - ! outgoing energy for any tallies with an outgoing energy filter - if (active_collision_tallies_size() > 0) call score_collision_tally(p) - if (active_analog_tallies_size() > 0) call score_analog_tally(p) - - ! Reset banked weight during collision - p % n_bank = 0 - p % wgt_bank = ZERO - p % n_delayed_bank(:) = 0 - - ! Reset fission logical - p % fission = .false. - - ! Save coordinates for tallying purposes - p % last_xyz_current = p % coord(1) % xyz - - ! Set last material to none since cross sections will need to be - ! re-evaluated - p % last_material = NONE - - ! Set all uvws to base level -- right now, after a collision, only the - ! base level uvws are changed - do j = 1, p % n_coord - 1 - if (p % coord(j + 1) % rotated) then - ! If next level is rotated, apply rotation matrix - p % coord(j + 1) % uvw = matmul(cells(p % coord(j) % cell + 1) % & - rotation_matrix, p % coord(j) % uvw) - else - ! Otherwise, copy this level's direction - p % coord(j + 1) % uvw = p % coord(j) % uvw - end if - end do - - ! Score flux derivative accumulators for differential tallies. - if (active_tallies_size() > 0) call score_collision_derivative(p) - end if - - ! If particle has too many events, display warning and kill it - n_event = n_event + 1 - if (n_event == MAX_EVENTS) then - if (master) call warning("Particle " // trim(to_str(p%id)) & - &// " underwent maximum number of events.") - p % alive = .false. - end if - - ! Check for secondary particles if this particle is dead - if (.not. p % alive) then - if (p % n_secondary > 0) then - call particle_from_source(p, p % secondary_bank(p % n_secondary)) - p % n_secondary = p % n_secondary - 1 - n_event = 0 - - ! Enter new particle in particle track file - if (p % write_track) call add_particle_track() - else - exit EVENT_LOOP - end if - end if - end do EVENT_LOOP - - ! Finish particle track output. - if (p % write_track) then - call write_particle_track(p) - call finalize_particle_track(p) - endif - - end subroutine transport - -!=============================================================================== -! CROSS_SURFACE handles all surface crossings, whether the particle leaks out of -! the geometry, is reflected, or crosses into a new lattice or cell -!=============================================================================== - - subroutine cross_surface(p) - type(Particle), intent(inout) :: p - - real(8) :: u ! x-component of direction - real(8) :: v ! y-component of direction - real(8) :: w ! z-component of direction - real(8) :: norm ! "norm" of surface normal - real(8) :: xyz(3) ! Saved global coordinate - integer :: i_surface ! index in surfaces -#ifdef DAGMC - integer :: i_cell ! index of new cell -#endif - logical :: rotational ! if rotational periodic BC applied - logical :: found ! particle found in universe? - class(Surface), pointer :: surf - class(Surface), pointer :: surf2 ! periodic partner surface - - i_surface = abs(p % surface) - surf => surfaces(i_surface) - if (verbosity >= 10 .or. trace) then - call write_message(" Crossing surface " // trim(to_str(surf % id()))) - end if - - if (surf % bc() == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE LEAKS OUT OF PROBLEM - - ! Kill particle - p % alive = .false. - - ! Score any surface current tallies -- note that the particle is moved - ! forward slightly so that if the mesh boundary is on the surface, it is - ! still processed - - if (active_meshsurf_tallies_size() > 0) then - ! TODO: Find a better solution to score surface currents than - ! physically moving the particle forward slightly - - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - end if - - ! Score to global leakage tally - global_tally_leakage = global_tally_leakage + p % wgt - - ! Display message - if (verbosity >= 10 .or. trace) then - call write_message(" Leaked out of surface " & - &// trim(to_str(surf % id()))) - end if - return - - elseif (surf % bc() == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then - ! ======================================================================= - ! PARTICLE REFLECTS FROM SURFACE - - ! Do not handle reflective boundary conditions on lower universes - if (p % n_coord /= 1) then - call particle_mark_as_lost(p, "Cannot reflect particle " & - // trim(to_str(p % id)) // " off surface in a lower universe.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change. For surface filters, we need to score the tallies - ! twice, once before the particle's surface attribute has changed and - ! once after. For mesh surface filters, we need to artificially move - ! the particle slightly back in case the surface crossing is coincident - ! with a mesh boundary - - if (active_surface_tallies_size() > 0) call score_surface_tally(p) - - - if (active_meshsurf_tallies_size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - p % coord(1) % xyz = xyz - end if - - ! Reflect particle off surface - call surf % reflect(p%coord(1)%xyz, p%coord(1)%uvw) - - ! Make sure new particle direction is normalized - u = p%coord(1)%uvw(1) - v = p%coord(1)%uvw(2) - w = p%coord(1)%uvw(3) - norm = sqrt(u*u + v*v + w*w) - p%coord(1)%uvw(:) = [u, v, w] / norm - - ! Reassign particle's cell and surface - p % coord(1) % cell = p % last_cell(p % last_n_coord) - p % surface = -p % surface - - ! If a reflective surface is coincident with a lattice or universe - ! boundary, it is necessary to redetermine the particle's coordinates in - ! the lower universes. - - p % n_coord = 1 - call find_cell(p, found, .true.) - if (.not. found) then - call particle_mark_as_lost(p, "Couldn't find particle after reflecting& - & from surface " // trim(to_str(surf % id())) // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Reflected from surface " & - &// trim(to_str(surf%id()))) - end if - return - elseif (surf % bc() == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then - ! ======================================================================= - ! PERIODIC BOUNDARY - - ! Do not handle periodic boundary conditions on lower universes - if (p % n_coord /= 1) then - call particle_mark_as_lost(p, "Cannot transfer particle " & - // trim(to_str(p % id)) // " across surface in a lower universe.& - & Boundary conditions must be applied to universe 0.") - return - end if - - ! Score surface currents since reflection causes the direction of the - ! particle to change -- artificially move the particle slightly back in - ! case the surface crossing is coincident with a mesh boundary - if (active_meshsurf_tallies_size() > 0) then - xyz = p % coord(1) % xyz - p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - call score_meshsurface_tally(p) - p % coord(1) % xyz = xyz - end if - - ! Get a pointer to the partner periodic surface. Offset the index to - ! correct for C vs. Fortran indexing. - surf2 => surfaces(surf % i_periodic() + 1) - - ! Adjust the particle's location and direction. - rotational = surf2 % periodic_translate(surf, p % coord(1) % xyz, & - p % coord(1) % uvw) - - ! Reassign particle's surface - if (rotational) then - p % surface = surf % i_periodic() + 1 - else - p % surface = sign(surf % i_periodic() + 1, p % surface) - end if - - ! Figure out what cell particle is in now - p % n_coord = 1 - call find_cell(p, found, .true.) - if (.not. found) then - call particle_mark_as_lost(p, "Couldn't find particle after hitting & - &periodic boundary on surface " // trim(to_str(surf % id())) & - // ".") - return - end if - - ! Set previous coordinate going slightly past surface crossing - p % last_xyz_current = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - - ! Diagnostic message - if (verbosity >= 10 .or. trace) then - call write_message(" Hit periodic boundary on surface " & - // trim(to_str(surf%id()))) - end if - return - end if - - ! ========================================================================== - ! SEARCH NEIGHBOR LISTS FOR NEXT CELL - -#ifdef DAGMC - if (dagmc) then - i_cell = next_cell(cells(p % last_cell(1) + 1), surfaces(abs(p % surface))) - ! save material and temp - p % last_material = p % material - p % last_sqrtkT = p % sqrtKT - ! set new cell value - p % coord(1) % cell = i_cell-1 ! decrement for C++ indexing - p % cell_instance = 1 - p % material = cells(i_cell) % material(1) - p % sqrtKT = cells(i_cell) % sqrtKT(0) - return - end if -#endif - - call find_cell(p, found, .true.) - if (found) return - - ! ========================================================================== - ! COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS - - ! Remove lower coordinate levels and assignment of surface - p % surface = ERROR_INT - p % n_coord = 1 - call find_cell(p, found) - - if (run_mode /= MODE_PLOTTING .and. (.not. found)) then - ! If a cell is still not found, there are two possible causes: 1) there is - ! a void in the model, and 2) the particle hit a surface at a tangent. If - ! the particle is really traveling tangent to a surface, if we move it - ! forward a tiny bit it should fix the problem. - - p % n_coord = 1 - p % coord(1) % xyz = p % coord(1) % xyz + TINY_BIT * p % coord(1) % uvw - call find_cell(p, found) - - ! Couldn't find next cell anywhere! This probably means there is an actual - ! undefined region in the geometry. - - if (.not. found) then - call particle_mark_as_lost(p, "After particle " // trim(to_str(p % id)) & - // " crossed surface " // trim(to_str(surf % id())) & - // " it could not be located in any cell and it did not leak.") - return - end if - end if - - end subroutine cross_surface - -end module tracking From a1de6b82a27b08ed9c157a7978ddd7b6029ba9c0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 23:20:19 -0600 Subject: [PATCH 02/10] Remove surfaces on Fortran side --- CMakeLists.txt | 1 - src/api.F90 | 4 +- src/geometry.F90 | 2 - src/input_xml.F90 | 35 ----------- src/surface.cpp | 68 ++++++--------------- src/surface_header.F90 | 135 ----------------------------------------- 6 files changed, 23 insertions(+), 222 deletions(-) delete mode 100644 src/surface_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 6bc51e812..1c0cab1dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -336,7 +336,6 @@ add_library(libopenmc SHARED src/state_point.F90 src/stl_vector.F90 src/string.F90 - src/surface_header.F90 src/track_output.F90 src/vector_header.F90 src/xml_interface.F90 diff --git a/src/api.F90 b/src/api.F90 index c1d934e36..c4c47ee00 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -64,7 +64,6 @@ contains use photon_header use settings use simulation_header - use surface_header use tally_filter_header use tally_header @@ -90,6 +89,9 @@ contains subroutine free_memory_volume() bind(C) end subroutine + subroutine free_memory_surfaces() bind(C) + end subroutine + subroutine sab_clear() bind(C) end subroutine end interface diff --git a/src/geometry.F90 b/src/geometry.F90 index 02ababae6..0abe96ee4 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -5,8 +5,6 @@ module geometry use particle_header use simulation_header use settings - use surface_header - use string, only: to_str use, intrinsic :: ISO_C_BINDING diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 238bb7302..ccdf142f2 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -19,7 +19,6 @@ module input_xml use nuclide_header use photon_header use random_lcg, only: prn - use surface_header use settings use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, & @@ -116,7 +115,6 @@ contains call write_message("Reading DAGMC geometry...", 5) call load_dagmc_geometry() - call allocate_surfaces() call allocate_cells() ! setup universe data structs @@ -153,7 +151,6 @@ contains integer :: n_cells_in_univ real(8) :: phi, theta, psi logical :: file_exists - logical :: boundary_exists character(MAX_LINE_LEN) :: filename type(Cell), pointer :: c type(XMLDocument) :: doc @@ -188,28 +185,8 @@ contains ! ========================================================================== ! READ SURFACES FROM GEOMETRY.XML - ! This variable is used to check whether at least one boundary condition was - ! applied to a surface - boundary_exists = .false. - call read_surfaces(root % ptr) - ! Allocate surfaces array - allocate(surfaces(surfaces_size())) - do i = 1, size(surfaces) - surfaces(i) % ptr = surface_pointer(i - 1); - - if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true. - end do - - ! Check to make sure a boundary condition was applied to at least one - ! surface - if (run_mode /= MODE_PLOTTING) then - if (.not. boundary_exists) then - call fatal_error("No boundary conditions were applied to any surfaces!") - end if - end if - ! ========================================================================== ! READ CELLS FROM GEOMETRY.XML @@ -316,18 +293,6 @@ contains end subroutine read_geometry_xml - subroutine allocate_surfaces() - integer :: i - - ! Allocate surfaces array - allocate(surfaces(surfaces_size())) - - do i = 1, size(surfaces) - surfaces(i) % ptr = surface_pointer(i - 1); - end do - - end subroutine allocate_surfaces - subroutine allocate_cells() integer :: i type(Cell), pointer :: c diff --git a/src/surface.cpp b/src/surface.cpp index 8b9b20ea1..0cceee9f4 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -7,8 +7,9 @@ #include "openmc/error.h" #include "openmc/hdf5_interface.h" -#include "openmc/xml_interface.h" +#include "openmc/settings.h" #include "openmc/string_utils.h" +#include "openmc/xml_interface.h" namespace openmc { @@ -1243,59 +1244,30 @@ read_surfaces(pugi::xml_node* node) } } } + + // Check to make sure a boundary condition was applied to at least one + // surface + bool boundary_exists = false; + for (const auto& surf : model::surfaces) { + if (surf->bc_ != BC_TRANSMIT) { + boundary_exists = true; + break; + } + } + if (settings::run_mode != RUN_MODE_PLOTTING && !boundary_exists) { + fatal_error("No boundary conditions were applied to any surfaces!"); + } } //============================================================================== // Fortran compatibility functions //============================================================================== -extern "C" { - Surface* surface_pointer(int surf_ind) {return model::surfaces[surf_ind];} - - int surface_id(Surface* surf) {return surf->id_;} - - int surface_bc(Surface* surf) {return surf->bc_;} - - void surface_reflect(Surface* surf, double xyz[3], double uvw[3]) - { - Position r {xyz}; - Direction u {uvw}; - u = surf->reflect(r, u); - - uvw[0] = u.x; - uvw[1] = u.y; - uvw[2] = u.z; - } - - int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic_;} - - bool - surface_periodic(PeriodicSurface* surf, PeriodicSurface* other, double xyz[3], - double uvw[3]) - { - Position r {xyz}; - Direction u {uvw}; - bool rotational = surf->periodic_translate(other, r, u); - - // Copy back to arrays - xyz[0] = r.x; - xyz[1] = r.y; - xyz[2] = r.z; - uvw[0] = u.x; - uvw[1] = u.y; - uvw[2] = u.z; - - return rotational; - } - - void free_memory_surfaces_c() - { - for (Surface* surf : model::surfaces) {delete surf;} - model::surfaces.clear(); - model::surface_map.clear(); - } - - int surfaces_size() { return model::surfaces.size(); } +extern "C" void free_memory_surfaces() +{ + for (Surface* surf : model::surfaces) {delete surf;} + model::surfaces.clear(); + model::surface_map.clear(); } } // namespace openmc diff --git a/src/surface_header.F90 b/src/surface_header.F90 deleted file mode 100644 index 17f9515f3..000000000 --- a/src/surface_header.F90 +++ /dev/null @@ -1,135 +0,0 @@ -module surface_header - - use, intrinsic :: ISO_C_BINDING - - use hdf5_interface - - implicit none - - interface - pure function surface_pointer(surf_ind) bind(C) result(ptr) - use ISO_C_BINDING - implicit none - integer(C_INT), intent(in), value :: surf_ind - type(C_PTR) :: ptr - end function surface_pointer - - pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - integer(C_INT) :: id - end function surface_id_c - - pure function surface_bc_c(surf_ptr) bind(C, name='surface_bc') result(bc) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - integer(C_INT) :: bc - end function surface_bc_c - - pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) & - bind(C, name='surface_reflect') - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(inout) :: uvw(3); - end subroutine surface_reflect_c - - pure function surface_i_periodic_c(surf_ptr) & - bind(C, name="surface_i_periodic") result(i_periodic) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - integer(C_INT) :: i_periodic - end function surface_i_periodic_c - - function surface_periodic_c(surf_ptr, other_ptr, xyz, uvw) & - bind(C, name="surface_periodic") result(rotational) - use ISO_C_BINDING - implicit none - type(C_PTR), intent(in), value :: surf_ptr - type(C_PTR), intent(in), value :: other_ptr - real(C_DOUBLE), intent(inout) :: xyz(3); - real(C_DOUBLE), intent(inout) :: uvw(3); - logical(C_BOOL) :: rotational - end function surface_periodic_c - - subroutine free_memory_surfaces_c() bind(C) - end subroutine free_memory_surfaces_c - end interface - -!=============================================================================== -! SURFACE type defines a first- or second-order surface that can be used to -! construct closed volumes (cells) -!=============================================================================== - - type :: Surface - type(C_PTR) :: ptr - - contains - - procedure :: id => surface_id - procedure :: bc => surface_bc - procedure :: reflect => surface_reflect - procedure :: i_periodic => surface_i_periodic - procedure :: periodic_translate => surface_periodic - - end type Surface - - type(Surface), allocatable, target :: surfaces(:) - - interface - function surfaces_size() result(sz) bind(C) - import C_INT - integer(C_INT) :: sz - end function - end interface - -contains - - pure function surface_id(this) result(id) - class(Surface), intent(in) :: this - integer(C_INT) :: id - id = surface_id_c(this % ptr) - end function surface_id - - pure function surface_bc(this) result(bc) - class(Surface), intent(in) :: this - integer(C_INT) :: bc - bc = surface_bc_c(this % ptr) - end function surface_bc - - pure subroutine surface_reflect(this, xyz, uvw) - class(Surface), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3); - real(C_DOUBLE), intent(inout) :: uvw(3); - call surface_reflect_c(this % ptr, xyz, uvw) - end subroutine surface_reflect - - pure function surface_i_periodic(this) result(i_periodic) - class(Surface), intent(in) :: this - integer(C_INT) :: i_periodic - i_periodic = surface_i_periodic_c(this % ptr) - end function surface_i_periodic - - function surface_periodic(this, other, xyz, uvw) result(rotational) - class(Surface), intent(in) :: this - class(Surface), intent(in) :: other - real(C_DOUBLE), intent(inout) :: xyz(3); - real(C_DOUBLE), intent(inout) :: uvw(3); - logical(C_BOOL) :: rotational - rotational = surface_periodic_c(this % ptr, other % ptr, xyz, uvw) - end function surface_periodic - -!=============================================================================== -! FREE_MEMORY_SURFACES deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_surfaces() - if (allocated(surfaces)) deallocate(surfaces) - call free_memory_surfaces_c() - end subroutine free_memory_surfaces - -end module surface_header From 2591a7054a15f0d462f726ce48406290826884a2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 13 Feb 2019 23:57:13 -0600 Subject: [PATCH 03/10] Remove nuclides array on Fortran side --- CMakeLists.txt | 3 - include/openmc/capi.h | 3 +- openmc/capi/material.py | 2 +- openmc/capi/nuclide.py | 12 +- openmc/capi/tally.py | 2 +- src/algorithm.F90 | 141 -------- src/api.F90 | 1 + src/cross_sections.cpp | 7 - src/endf.F90 | 268 --------------- src/geometry_header.F90 | 3 +- src/input_xml.F90 | 9 - src/nuclide.cpp | 82 ++--- src/nuclide_header.F90 | 613 +---------------------------------- src/particle_restart.F90 | 5 +- src/reaction_header.F90 | 270 --------------- src/simulation.F90 | 4 +- src/tallies/tally.F90 | 3 - src/tallies/tally_header.F90 | 1 - 18 files changed, 42 insertions(+), 1387 deletions(-) delete mode 100644 src/algorithm.F90 delete mode 100644 src/endf.F90 delete mode 100644 src/reaction_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1c0cab1dc..96547ebbf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -303,14 +303,12 @@ target_compile_options(faddeeva PRIVATE ${cxxflags}) #=============================================================================== add_library(libopenmc SHARED - src/algorithm.F90 src/bank_header.F90 src/api.F90 src/dagmc_header.F90 src/constants.F90 src/dict_header.F90 src/eigenvalue.F90 - src/endf.F90 src/error.F90 src/geometry.F90 src/geometry_header.F90 @@ -328,7 +326,6 @@ add_library(libopenmc SHARED src/photon_header.F90 src/pugixml/pugixml_f.F90 src/random_lcg.F90 - src/reaction_header.F90 src/relaxng src/settings.F90 src/simulation_header.F90 diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 66d718ff5..80d431063 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -83,7 +83,7 @@ extern "C" { int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_next_batch(int* status); - int openmc_nuclide_name(int index, char** name); + int openmc_nuclide_name(int index, const char** name); int openmc_particle_restart(); int openmc_plot_geometry(); int openmc_reset(); @@ -165,7 +165,6 @@ extern "C" { // Global variables extern char openmc_err_msg[256]; - extern int n_nuclides; extern int32_t n_realizations; extern int32_t n_sab_tables; extern int32_t n_tallies; diff --git a/openmc/capi/material.py b/openmc/capi/material.py index 5057c8ff4..4bf19a1e0 100644 --- a/openmc/capi/material.py +++ b/openmc/capi/material.py @@ -163,7 +163,7 @@ class Material(_FortranObjectWithID): _dll.openmc_material_get_densities(self._index, nuclides, densities, n) # Convert to appropriate types and return - nuclide_list = [Nuclide(nuclides[i] + 1).name for i in range(n.value)] + nuclide_list = [Nuclide(nuclides[i]).name for i in range(n.value)] density_array = as_array(densities, (n.value,)) return nuclide_list, density_array diff --git a/openmc/capi/nuclide.py b/openmc/capi/nuclide.py index e57653c64..824f27e32 100644 --- a/openmc/capi/nuclide.py +++ b/openmc/capi/nuclide.py @@ -23,6 +23,7 @@ _dll.openmc_load_nuclide.errcheck = _error_handler _dll.openmc_nuclide_name.argtypes = [c_int, POINTER(c_char_p)] _dll.openmc_nuclide_name.restype = c_int _dll.openmc_nuclide_name.errcheck = _error_handler +_dll.nuclides_size.restype = c_int def load_nuclide(name): @@ -70,12 +71,7 @@ class Nuclide(_FortranObject): def name(self): name = c_char_p() _dll.openmc_nuclide_name(self._index, name) - - # Find blank in name - i = 0 - while name.value[i:i+1] != b' ': - i += 1 - return name.value[:i].decode() + return name.value.decode() class _NuclideMapping(Mapping): @@ -91,10 +87,10 @@ class _NuclideMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield Nuclide(i + 1).name + yield Nuclide(i).name def __len__(self): - return c_int.in_dll(_dll, 'n_nuclides').value + return _dll.nuclides_size() def __repr__(self): return repr(dict(self)) diff --git a/openmc/capi/tally.py b/openmc/capi/tally.py index df8ef9d08..029baf4bb 100644 --- a/openmc/capi/tally.py +++ b/openmc/capi/tally.py @@ -278,7 +278,7 @@ class Tally(_FortranObjectWithID): nucs = POINTER(c_int)() n = c_int() _dll.openmc_tally_get_nuclides(self._index, nucs, n) - return [Nuclide(nucs[i]+1).name if nucs[i] >= 0 else 'total' + return [Nuclide(nucs[i]).name if nucs[i] >= 0 else 'total' for i in range(n.value)] @nuclides.setter diff --git a/src/algorithm.F90 b/src/algorithm.F90 deleted file mode 100644 index d9f9e2e22..000000000 --- a/src/algorithm.F90 +++ /dev/null @@ -1,141 +0,0 @@ -module algorithm - - use constants - use stl_vector, only: VectorInt, VectorReal - - implicit none - - interface sort - module procedure sort_int, sort_real, sort_vector_int, sort_vector_real - end interface sort - - interface find - module procedure find_int, find_real, find_vector_int, find_vector_real - end interface find - -contains - -!=============================================================================== -! SORT sorts an array in place using an insertion sort. -!=============================================================================== - - pure subroutine sort_int(array) - integer, intent(inout) :: array(:) - - integer :: k, m - integer :: temp - - if (size(array) > 1) then - SORT: do k = 2, size(array) - ! Save value to move - m = k - temp = array(k) - - MOVE_OVER: do while (m > 1) - ! Check if insertion value is greater than (m-1)th value - if (temp >= array(m - 1)) exit - - ! Move values over until hitting one that's not larger - array(m) = array(m - 1) - m = m - 1 - end do MOVE_OVER - - ! Put the original value into its new position - array(m) = temp - end do SORT - end if - end subroutine sort_int - - pure subroutine sort_real(array) - real(8), intent(inout) :: array(:) - - integer :: k, m - real(8) :: temp - - if (size(array) > 1) then - SORT: do k = 2, size(array) - ! Save value to move - m = k - temp = array(k) - - MOVE_OVER: do while (m > 1) - ! Check if insertion value is greater than (m-1)th value - if (temp >= array(m - 1)) exit - - ! Move values over until hitting one that's not larger - array(m) = array(m - 1) - m = m - 1 - end do MOVE_OVER - - ! Put the original value into its new position - array(m) = temp - end do SORT - end if - end subroutine sort_real - - pure subroutine sort_vector_int(vec) - type(VectorInt), intent(inout) :: vec - - call sort_int(vec % data(1:vec%size())) - end subroutine sort_vector_int - - pure subroutine sort_vector_real(vec) - type(VectorReal), intent(inout) :: vec - - call sort_real(vec % data(1:vec%size())) - end subroutine sort_vector_real - -!=============================================================================== -! FIND determines the index of the first occurrence of a value in an array. If -! the value does not appear in the array, -1 is returned. -!=============================================================================== - - pure function find_int(array, val) result(index) - integer, intent(in) :: array(:) - integer, intent(in) :: val - integer :: index - - integer :: i - - index = -1 - do i = 1, size(array) - if (array(i) == val) then - index = i - exit - end if - end do - end function find_int - - pure function find_real(array, val) result(index) - real(8), intent(in) :: array(:) - real(8), intent(in) :: val - integer :: index - - integer :: i - - index = -1 - do i = 1, size(array) - if (array(i) == val) then - index = i - exit - end if - end do - end function find_real - - pure function find_vector_int(vec, val) result(index) - type(VectorInt), intent(in) :: vec - integer, intent(in) :: val - integer :: index - - index = find_int(vec % data(1:vec % size()), val) - end function find_vector_int - - pure function find_vector_real(vec, val) result(index) - type(VectorReal), intent(in) :: vec - real(8), intent(in) :: val - integer :: index - - index = find_real(vec % data(1:vec % size()), val) - end function find_vector_real - -end module algorithm diff --git a/src/api.F90 b/src/api.F90 index c4c47ee00..5200927cf 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -61,6 +61,7 @@ contains use bank_header use geometry_header use material_header + use nuclide_header use photon_header use settings use simulation_header diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 23af2f624..150eb51cb 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -179,9 +179,6 @@ void read_cross_sections_xml() } } -extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr, - const double* temps, int n, int n_nuclide); - void read_ce_cross_sections(const std::vector>& nuc_temps, const std::vector>& thermal_temps) @@ -226,10 +223,6 @@ read_ce_cross_sections(const std::vector>& nuc_temps, data::nuclides.push_back(std::make_unique( group, nuc_temps[i_nuc], i_nuclide)); - // Read from Fortran too - nuclide_from_hdf5(group, data::nuclides.back().get(), - &nuc_temps[i_nuc].front(), nuc_temps[i_nuc].size(), i_nuclide + 1); - close_group(group); file_close(file_id); diff --git a/src/endf.F90 b/src/endf.F90 deleted file mode 100644 index 58b554997..000000000 --- a/src/endf.F90 +++ /dev/null @@ -1,268 +0,0 @@ -module endf - - use constants - use string, only: to_str - - implicit none - -contains - -!=============================================================================== -! REACTION_NAME gives the name of the reaction for a given MT value -!=============================================================================== - - pure function reaction_name(MT) result(string) - - integer, intent(in) :: MT - character(MAX_WORD_LEN) :: string - - select case (MT) - ! Special reactions for tallies - case (SCORE_FLUX) - string = "flux" - case (SCORE_TOTAL) - string = "total" - case (SCORE_SCATTER) - string = "scatter" - case (SCORE_NU_SCATTER) - string = "nu-scatter" - case (SCORE_ABSORPTION) - string = "absorption" - case (SCORE_FISSION) - string = "fission" - case (SCORE_NU_FISSION) - string = "nu-fission" - case (SCORE_DECAY_RATE) - string = "decay-rate" - case (SCORE_DELAYED_NU_FISSION) - string = "delayed-nu-fission" - case (SCORE_PROMPT_NU_FISSION) - string = "prompt-nu-fission" - case (SCORE_KAPPA_FISSION) - string = "kappa-fission" - case (SCORE_CURRENT) - string = "current" - case (SCORE_EVENTS) - string = "events" - case (SCORE_INVERSE_VELOCITY) - string = "inverse-velocity" - case (SCORE_FISS_Q_PROMPT) - string = "fission-q-prompt" - case (SCORE_FISS_Q_RECOV) - string = "fission-q-recoverable" - - ! Normal ENDF-based reactions - case (TOTAL_XS) - string = '(n,total)' - case (ELASTIC) - string = '(n,elastic)' - case (N_LEVEL) - string = '(n,level)' - case (N_2ND) - string = '(n,2nd)' - case (N_2N) - string = '(n,2n)' - case (N_3N) - string = '(n,3n)' - case (N_FISSION) - string = '(n,fission)' - case (N_F) - string = '(n,f)' - case (N_NF) - string = '(n,nf)' - case (N_2NF) - string = '(n,2nf)' - case (N_NA) - string = '(n,na)' - case (N_N3A) - string = '(n,n3a)' - case (N_2NA) - string = '(n,2na)' - case (N_3NA) - string = '(n,3na)' - case (N_NP) - string = '(n,np)' - case (N_N2A) - string = '(n,n2a)' - case (N_2N2A) - string = '(n,2n2a)' - case (N_ND) - string = '(n,nd)' - case (N_NT) - string = '(n,nt)' - case (N_N3HE) - string = '(n,nHe-3)' - case (N_ND2A) - string = '(n,nd2a)' - case (N_NT2A) - string = '(n,nt2a)' - case (N_4N) - string = '(n,4n)' - case (N_3NF) - string = '(n,3nf)' - case (N_2NP) - string = '(n,2np)' - case (N_3NP) - string = '(n,3np)' - case (N_N2P) - string = '(n,n2p)' - case (N_NPA) - string = '(n,npa)' - case (N_N1 : N_N40) - string = '(n,n' // trim(to_str(MT-50)) // ')' - case (N_NC) - string = '(n,nc)' - case (N_DISAPPEAR) - string = '(n,disappear)' - case (N_GAMMA) - string = '(n,gamma)' - case (N_P) - string = '(n,p)' - case (N_D) - string = '(n,d)' - case (N_T) - string = '(n,t)' - case (N_3HE) - string = '(n,3He)' - case (N_A) - string = '(n,a)' - case (N_2A) - string = '(n,2a)' - case (N_3A) - string = '(n,3a)' - case (N_2P) - string = '(n,2p)' - case (N_PA) - string = '(n,pa)' - case (N_T2A) - string = '(n,t2a)' - case (N_D2A) - string = '(n,d2a)' - case (N_PD) - string = '(n,pd)' - case (N_PT) - string = '(n,pt)' - case (N_DA) - string = '(n,da)' - case (201) - string = '(n,Xn)' - case (202) - string = '(n,Xgamma)' - case (203) - string = '(n,Xp)' - case (204) - string = '(n,Xd)' - case (205) - string = '(n,Xt)' - case (206) - string = '(n,X3He)' - case (207) - string = '(n,Xa)' - case (444) - string = '(damage)' - case (COHERENT) - string = 'coherent scatter' - case (INCOHERENT) - string = 'incoherent scatter' - case (PAIR_PROD_ELEC) - string = 'pair production, electron' - case (PAIR_PROD) - string = 'pair production' - case (PAIR_PROD_NUC) - string = 'pair production, nuclear' - case (PHOTOELECTRIC) - string = 'photoelectric' - case (534 : 572) - string = 'photoelectric, ' // trim(SUBSHELLS(MT - 533)) // ' subshell' - case (600 : 648) - string = '(n,p' // trim(to_str(MT-600)) // ')' - case (649) - string = '(n,pc)' - case (650 : 698) - string = '(n,d' // trim(to_str(MT-650)) // ')' - case (699) - string = '(n,dc)' - case (700 : 748) - string = '(n,t' // trim(to_str(MT-700)) // ')' - case (749) - string = '(n,tc)' - case (750 : 798) - string = '(n,3He' // trim(to_str(MT-750)) // ')' - case (799) - string = '(n,3Hec)' - case (800 : 848) - string = '(n,a' // trim(to_str(MT-800)) // ')' - case (849) - string = '(n,ac)' - case default - string = 'MT=' // trim(to_str(MT)) - end select - - end function reaction_name - -!=============================================================================== -! IS_FISSION determines if a given MT number is that of a fission event. This -! accounts for aggregate fission (MT=18) as well as partial fission reactions. -!=============================================================================== - - function is_fission(MT) result(fission_event) - - integer, intent(in) :: MT - logical :: fission_event - - if (MT == N_FISSION .or. MT == N_F .or. MT == N_NF .or. MT == N_2NF & - .or. MT == N_3NF) then - fission_event = .true. - else - fission_event = .false. - end if - - end function is_fission - -!=============================================================================== -! IS_DISAPPEARANCE determines if a given MT number is that of a disappearance -! reaction, i.e. a reaction with no neutron in the exit channel -!=============================================================================== - - function is_disappearance(MT) result(dis) - - integer, intent(in) :: MT - logical :: dis - - if (MT >= N_DISAPPEAR .and. MT <= N_DA) then - dis = .true. - elseif (MT >= N_P0 .and. MT <= N_AC) then - dis = .true. - elseif (any(MT == [N_TA, N_DT, N_P3HE, N_D3HE, N_3HEA, N_3P])) then - dis = .true. - else - dis = .false. - end if - - end function is_disappearance - -!=============================================================================== -! IS_INELASTIC_SCATTER determines if a given MT number is that of an inelastic -! scattering event -!=============================================================================== - - function is_inelastic_scatter(MT) result(retval) - - integer, intent(in) :: MT - logical :: retval - - if (MT < 100) then - if (is_fission(MT)) then - retval = .false. - else - retval = (MT >= MISC .and. MT /= 27) - end if - elseif (MT <= 200) then - retval = (.not. is_disappearance(MT)) - else - retval = .false. - end if - - end function - -end module endf diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index b6c454b9d..563a60687 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -3,7 +3,8 @@ module geometry_header use, intrinsic :: ISO_C_BINDING use dict_header, only: DictIntInt - use nuclide_header + use error + use string, only: to_str implicit none diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ccdf142f2..94a1d6190 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2,10 +2,8 @@ module input_xml use, intrinsic :: ISO_C_BINDING - use algorithm, only: find use constants use dict_header, only: DictIntInt, DictCharInt, DictEntryCI - use endf, only: reaction_name use error, only: fatal_error, warning, write_message, openmc_err_msg use geometry_header #ifdef DAGMC @@ -320,11 +318,6 @@ contains type(XMLNode) :: root interface - function nuclides_size() bind(C) result(n) - import C_INT - integer(C_INT) :: n - end function - function elements_size() bind(C) result(n) import C_INT integer(C_INT) :: n @@ -360,9 +353,7 @@ contains call read_materials(root % ptr) ! Set total number of nuclides and elements - n_nuclides = nuclides_size() n_elements = elements_size() - allocate(nuclides(n_nuclides)) ! Close materials XML file call doc % clear() diff --git a/src/nuclide.cpp b/src/nuclide.cpp index ac02a9c71..9b9a9ef31 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -892,18 +892,11 @@ void check_data_version(hid_t file_id) // C API //============================================================================== -extern "C" void extend_nuclides(); -extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr, - const double* temps, int n, int n_nuclide); - extern "C" int openmc_load_nuclide(const char* name) { if (data::nuclide_map.find(name) == data::nuclide_map.end()) { const auto& it = data::library_map.find({Library::Type::neutron, name}); if (it != data::library_map.end()) { - // Extend nuclides array on Fortran side - extend_nuclides(); - // Get filename for library containing nuclide int idx = it->second; std::string& filename = data::libraries[idx].path_; @@ -920,10 +913,6 @@ extern "C" int openmc_load_nuclide(const char* name) data::nuclides.push_back(std::make_unique( group, temperature, i_nuclide)); - // Read from Fortran too - nuclide_from_hdf5(group, data::nuclides.back().get(), - &temperature.front(), temperature.size(), i_nuclide + 1); - close_group(group); file_close(file_id); @@ -943,6 +932,30 @@ extern "C" int openmc_load_nuclide(const char* name) return 0; } +extern "C" int +openmc_get_nuclide_index(const char* name, int* index) +{ + auto it = data::nuclide_map.find(name); + if (it == data::nuclide_map.end()) { + set_errmsg("No nuclide named '" + std::string{name} + "' has been loaded."); + return OPENMC_E_DATA; + } + *index = it->second; + return 0; +} + +extern "C" int +openmc_nuclide_name(int index, const char** name) +{ + if (index >= 0 && index < data::nuclides.size()) { + *name = data::nuclides[index]->name_.data(); + return 0; + } else { + set_errmsg("Index in nuclides vector is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + //============================================================================== // Fortran compatibility functions //============================================================================== @@ -956,47 +969,6 @@ set_particle_energy_bounds(int particle, double E_min, double E_max) extern "C" int nuclides_size() { return data::nuclide_map.size(); } -extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); } - -extern "C" double nuclide_awr(int i_nuc) { return data::nuclides[i_nuc - 1]->awr_; } - -extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx) -{ - return nuc->reactions_[i_rx-1].get(); -} - -extern "C" void nuclide_calculate_xs_c(Nuclide* nuc, int i_sab, double E, - int i_log_union, double sqrtkT, double sab_frac) -{ - nuc->calculate_xs(i_sab, E, i_log_union, sqrtkT, sab_frac); -} - -extern "C" void nuclide_calculate_elastic_xs_c(Nuclide* nuc) -{ - nuc->calculate_elastic_xs(); -} - -extern "C" double nuclide_nu_c(Nuclide* nuc, double E, int emission_mode, int group) -{ - return nuc->nu(E, static_cast(emission_mode - 1), group); -} - -extern "C" double nuclide_fission_q_prompt(Nuclide* nuc, double E) -{ - return nuc->fission_q_prompt_ ? (*nuc->fission_q_prompt_)(E) : 0.0; -} - -extern "C" double nuclide_fission_q_recov(Nuclide* nuc, double E) -{ - return nuc->fission_q_recov_ ? (*nuc->fission_q_recov_)(E) : 0.0; -} - -extern "C" void multipole_deriv_eval(Nuclide* nuc, double E, double sqrtkT, - double* sig_s, double* sig_a, double* sig_f) -{ - std::tie(*sig_s, *sig_a, *sig_f) = nuc->multipole_->evaluate_deriv(E, sqrtkT); -} - extern "C" bool multipole_in_range(const Nuclide* nuc, double E) { return nuc->multipole_ && E >= nuc->multipole_->E_min_&& @@ -1009,12 +981,6 @@ extern "C" void nuclides_clear() data::nuclide_map.clear(); } -extern "C" int nuclide_map_get(const char* name) -{ - auto it = data::nuclide_map.find(name); - return it == data::nuclide_map.end() ? -1 : it->second + 1; -} - extern "C" NuclideMicroXS* micro_xs_ptr(); extern "C" ElementMicroXS* micro_photon_xs_ptr(); diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index f84716964..b0a3f6444 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -3,82 +3,10 @@ module nuclide_header use, intrinsic :: ISO_FORTRAN_ENV use, intrinsic :: ISO_C_BINDING - use algorithm, only: sort, find use constants - use endf, only: is_fission, is_disappearance - use error - use hdf5_interface - use message_passing - use reaction_header, only: Reaction - use settings - use stl_vector, only: VectorInt, VectorReal - use string implicit none -!=============================================================================== -! Nuclide contains the base nuclidic data for a nuclide described as needed -! for continuous-energy neutron transport. -!=============================================================================== - - type EnergyGrid - real(8), allocatable :: energy(:) ! energy values corresponding to xs - end type EnergyGrid - - ! Positions for first dimension of Nuclide % xs - integer, parameter :: & - XS_TOTAL = 1, & - XS_ABSORPTION = 2, & - XS_FISSION = 3, & - XS_NU_FISSION = 4, & - XS_PHOTON_PROD = 5 - - ! The array within SumXS is of shape (5, n_energy) where the first dimension - ! corresponds to the following values: 1) total, 2) absorption (MT > 100), 3) - ! fission, 4) neutron production, 5) photon production - type SumXS - real(8), allocatable :: value(:,:) - end type SumXS - - type :: Nuclide - ! Nuclide meta-data - character(20) :: name ! name of nuclide, e.g. U235 - integer :: Z ! atomic number - integer :: A ! mass number - integer :: metastable ! metastable state - real(8) :: awr ! Atomic Weight Ratio - integer :: i_nuclide ! The nuclides index in the nuclides array - real(8), allocatable :: kTs(:) ! temperature in eV (k*T) - - ! Fission information - logical :: fissionable = .false. ! nuclide is fissionable? - - ! Energy grid for each temperature - type(EnergyGrid), allocatable :: grid(:) - - ! Microscopic cross sections - type(SumXS), allocatable :: xs(:) - - ! Fission information - integer, allocatable :: index_fission(:) ! indices in reactions - - ! Reactions - type(Reaction), allocatable :: reactions(:) - - ! Array that maps MT values to index in reactions; used at tally-time. Note - ! that ENDF-102 does not assign any MT values above 891. - integer :: reaction_index(891) - - type(C_PTR) :: ptr - - contains - procedure :: init_grid => nuclide_init_grid - procedure :: nu => nuclide_nu - procedure, private :: create_derived => nuclide_create_derived - procedure :: calculate_xs => nuclide_calculate_xs - procedure :: calculate_elastic_xs => nuclide_calculate_elastic_xs - end type Nuclide - !=============================================================================== ! NUCLIDEMICROXS contains cached microscopic cross sections for a particular ! nuclide at the current energy @@ -140,487 +68,25 @@ module nuclide_header real(C_DOUBLE) :: pair_production ! macroscopic pair production xs end type MaterialMacroXS - ! Nuclear data for each nuclide - type(Nuclide), allocatable, target :: nuclides(:) - integer(C_INT), bind(C) :: n_nuclides - ! Cross section caches type(NuclideMicroXS), allocatable, target :: micro_xs(:) ! Cache for each nuclide type(MaterialMacroXS), bind(C) :: material_xs ! Cache for current material !$omp threadprivate(micro_xs, material_xs) interface - function library_present_c(type, name) result(b) bind(C, name='library_present') - import C_INT, C_CHAR, C_BOOL - integer(C_INT), value :: type - character(kind=C_CHAR), intent(in) :: name(*) - logical(C_BOOL) :: b - end function - - function library_path_c(type, name) result(path) bind(C, name='library_path') - import C_INT, C_CHAR, C_PTR - integer(C_INT), value :: type - character(kind=C_CHAR), intent(in) :: name(*) - type(C_PTR) :: path - end function - - subroutine multipole_deriv_eval(ptr, E, sqrtkT, sig_s, sig_a, sig_f) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - real(C_DOUBLE), value :: sqrtkT - real(C_DOUBLE), intent(out) :: sig_s - real(C_DOUBLE), intent(out) :: sig_a - real(C_DOUBLE), intent(out) :: sig_f - end subroutine - - function multipole_in_range(ptr, E) result(b) bind(C) - import C_PTR, C_DOUBLE, C_BOOL - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - logical(C_BOOL) :: b - end function - - function nuclide_awr(i_nuc) result(awr) bind(C) - import C_INT, C_DOUBLE - integer(C_INT), value :: i_nuc - real(C_DOUBLE) :: awr - end function - - function nuclide_fission_q_prompt(ptr, E) result(q) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - real(C_DOUBLE) :: q - end function - - function nuclide_fission_q_recov(ptr, E) result(q) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - real(C_DOUBLE) :: q - end function - - function nuclide_map_get(name) result(idx) bind(C) - import C_CHAR, C_INT - character(kind=C_CHAR), intent(in) :: name(*) - integer(C_INT) :: idx + function nuclides_size() bind(C) result(n) + import C_INT + integer(C_INT) :: n end function end interface contains - function library_path(type, name) result(path) - integer, intent(in) :: type - character(len=*), intent(in) :: name - character(MAX_FILE_LEN) :: path - - type(C_PTR) :: ptr - character(kind=C_CHAR), pointer :: string(:) - - ptr = library_path_c(type, to_c_string(name)) - call c_f_pointer(ptr, string, [255]) - path = to_f_string(string) - end function - - function library_present(type, name) result(b) - integer, intent(in) :: type - character(len=*), intent(in) :: name - logical :: b - - b = library_present_c(type, to_c_string(name)) - end function - function micro_xs_ptr() result(ptr) bind(C) type(C_PTR) :: ptr ptr = C_LOC(micro_xs(1)) end function - subroutine nuclide_from_hdf5(group_id, ptr, temps, n, i_nuclide) bind(C) - integer(HID_T), value :: group_id - type(C_PTR), value :: ptr - type(C_PTR), value :: temps - integer(C_INT), value :: n - integer(C_INT), value :: i_nuclide - - real(C_DOUBLE), pointer :: temperature(:) ! list of desired temperatures - - integer :: i - integer :: i_closest - integer :: n_temperature - integer(HID_T) :: energy_group, energy_dset - integer(HID_T) :: kT_group - integer(HID_T) :: rxs_group - integer(HID_T) :: rx_group - integer(HSIZE_T) :: j - integer(HSIZE_T) :: dims(1) - character(MAX_WORD_LEN) :: temp_str - character(MAX_WORD_LEN), allocatable :: dset_names(:) - character(MAX_WORD_LEN), allocatable :: grp_names(:) - real(8), allocatable :: temps_available(:) ! temperatures available - real(8) :: temp_desired - real(8) :: temp_actual - type(VectorInt) :: MTs - type(VectorInt) :: temps_to_read - - ! Get array passed - call c_f_pointer(temps, temperature, [n]) - - associate (this => nuclides(i_nuclide)) - this % ptr = ptr - - ! Get name of nuclide from group - this % name = get_name(group_id) - - ! Get rid of leading '/' - this % name = trim(this % name(2:)) - - call read_attribute(this % Z, group_id, 'Z') - call read_attribute(this % A, group_id, 'A') - call read_attribute(this % metastable, group_id, 'metastable') - call read_attribute(this % awr, group_id, 'atomic_weight_ratio') - kT_group = open_group(group_id, 'kTs') - - ! Determine temperatures available - call get_datasets(kT_group, dset_names) - allocate(temps_available(size(dset_names))) - do i = 1, size(dset_names) - ! Read temperature value - call read_dataset(temps_available(i), kT_group, trim(dset_names(i))) - temps_available(i) = temps_available(i) / K_BOLTZMANN - end do - call sort(temps_available) - - ! If only one temperature is available, revert to nearest temperature - if (size(temps_available) == 1 .and. temperature_method == TEMPERATURE_INTERPOLATION) then - if (master) then - call warning("Cross sections for " // trim(this % name) // " are only & - &available at one temperature. Reverting to nearest temperature & - &method.") - end if - temperature_method = TEMPERATURE_NEAREST - end if - - ! Determine actual temperatures to read -- start by checking whether a - ! temperature range was given, in which case all temperatures in the range - ! are loaded irrespective of what temperatures actually appear in the model - if (temperature_range(2) > ZERO) then - do i = 1, size(temps_available) - temp_actual = temps_available(i) - if (temperature_range(1) <= temp_actual .and. temp_actual <= temperature_range(2)) then - call temps_to_read % push_back(nint(temp_actual)) - end if - end do - end if - - select case (temperature_method) - case (TEMPERATURE_NEAREST) - ! Find nearest temperatures - do i = 1, n - temp_desired = temperature(i) - i_closest = minloc(abs(temps_available - temp_desired), dim=1) - temp_actual = temps_available(i_closest) - if (abs(temp_actual - temp_desired) < temperature_tolerance) then - if (find(temps_to_read, nint(temp_actual)) == -1) then - call temps_to_read % push_back(nint(temp_actual)) - - ! Write warning for resonance scattering data if 0K is not available - if (abs(temp_actual - temp_desired) > 0 .and. temp_desired == 0 & - .and. master) then - call warning(trim(this % name) // " does not contain 0K data & - &needed for resonance scattering options selected. Using & - &data at " // trim(to_str(temp_actual)) & - // " K instead.") - end if - end if - else - call fatal_error("Nuclear data library does not contain cross & - §ions for " // trim(this % name) // " at or near " // & - trim(to_str(nint(temp_desired))) // " K.") - end if - end do - - case (TEMPERATURE_INTERPOLATION) - ! If temperature interpolation or multipole is selected, get a list of - ! bounding temperatures for each actual temperature present in the model - TEMP_LOOP: do i = 1, n - temp_desired = temperature(i) - - do j = 1, size(temps_available) - 1 - if (temps_available(j) <= temp_desired .and. & - temp_desired < temps_available(j + 1)) then - if (find(temps_to_read, nint(temps_available(j))) == -1) then - call temps_to_read % push_back(nint(temps_available(j))) - end if - if (find(temps_to_read, nint(temps_available(j + 1))) == -1) then - call temps_to_read % push_back(nint(temps_available(j + 1))) - end if - cycle TEMP_LOOP - end if - end do - - call fatal_error("Nuclear data library does not contain cross sections & - &for " // trim(this % name) // " at temperatures that bound " // & - trim(to_str(nint(temp_desired))) // " K.") - end do TEMP_LOOP - - end select - - ! Sort temperatures to read - call sort(temps_to_read) - - n_temperature = temps_to_read % size() - allocate(this % kTs(n_temperature)) - allocate(this % grid(n_temperature)) - - ! Get kT values - do i = 1, n_temperature - ! Get temperature as a string - temp_str = trim(to_str(temps_to_read % data(i))) // "K" - - ! Read exact temperature value - call read_dataset(this % kTs(i), kT_group, trim(temp_str)) - end do - call close_group(kT_group) - - ! Read energy grid - energy_group = open_group(group_id, 'energy') - do i = 1, n_temperature - temp_str = trim(to_str(temps_to_read % data(i))) // "K" - energy_dset = open_dataset(energy_group, temp_str) - call get_shape(energy_dset, dims) - allocate(this % grid(i) % energy(int(dims(1), 4))) - call read_dataset(this % grid(i) % energy, energy_dset) - call close_dataset(energy_dset) - end do - call close_group(energy_group) - - ! Get MT values based on group names - rxs_group = open_group(group_id, 'reactions') - call get_groups(rxs_group, grp_names) - do j = 1, size(grp_names) - if (starts_with(grp_names(j), "reaction_")) then - call MTs % push_back(int(str_to_int(grp_names(j)(10:12)))) - end if - end do - - ! Read reactions - allocate(this % reactions(MTs % size())) - do i = 1, size(this % reactions) - rx_group = open_group(rxs_group, 'reaction_' // trim(& - zero_padded(MTs % data(i), 3))) - - ! Set pointer for each reaction - call this % reactions(i) % init(this % ptr, i) - - call close_group(rx_group) - end do - call close_group(rxs_group) - - ! Create derived cross section data - call this % create_derived() - - ! Finalize with the nuclide index - this % i_nuclide = i_nuclide - end associate - - end subroutine nuclide_from_hdf5 - - subroutine nuclide_create_derived(this) - class(Nuclide), intent(inout) :: this - - integer :: i, j, k, l - integer :: t - integer :: n - integer :: n_grid - integer :: i_fission - integer :: n_temperature - type(VectorInt) :: MTs - - n_temperature = size(this % kTs) - allocate(this % xs(n_temperature)) - this % reaction_index(:) = 0 - do i = 1, n_temperature - ! Allocate and initialize derived cross sections - n_grid = size(this % grid(i) % energy) - allocate(this % xs(i) % value(5,n_grid)) - this % xs(i) % value(:,:) = ZERO - end do - - i_fission = 0 - - do i = 1, size(this % reactions) - call MTs % push_back(this % reactions(i) % MT) - this % reaction_index(this % reactions(i) % MT) = i - - associate (rx => this % reactions(i)) - do t = 1, n_temperature - j = rx % xs_threshold(t) - n = rx % xs_size(t) - - ! Calculate photon production cross section - do k = 1, rx % products_size() - if (rx % product_particle(k) == PHOTON) then - do l = 1, n - this % xs(t) % value(XS_PHOTON_PROD,l+j-1) = & - this % xs(t) % value(XS_PHOTON_PROD,l+j-1) + & - rx % xs(t, l) * rx % product_yield(k, & - this % grid(t) % energy(l+j-1)) - end do - end if - end do - - ! Skip gas production cross sections (MT=200+), etc. - if (rx % MT > N_5N2P .and. rx % MT < N_P0) cycle - - ! Skip any reaction that has been marked as redundant - if (rx % redundant) cycle - - ! Add contribution to total cross section - do k = j, j + n - 1 - this % xs(t) % value(XS_TOTAL,k) = this % xs(t) % & - value(XS_TOTAL,k) + rx % xs(t, k - j + 1) - end do - - ! Add contribution to absorption cross section - if (is_disappearance(rx % MT)) then - do k = j, j + n - 1 - this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % & - value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1) - end do - end if - - ! Information about fission reactions - if (t == 1) then - if (rx % MT == N_FISSION) then - allocate(this % index_fission(1)) - elseif (rx % MT == N_F) then - allocate(this % index_fission(PARTIAL_FISSION_MAX)) - end if - end if - - ! Add contribution to fission cross section - if (is_fission(rx % MT)) then - this % fissionable = .true. - do k = j, j + n - 1 - this % xs(t) % value(XS_FISSION,k) = this % xs(t) % & - value(XS_FISSION,k) + rx % xs(t, k - j + 1) - - ! Also need to add fission cross sections to absorption - this % xs(t) % value(XS_ABSORPTION,k) = this % xs(t) % & - value(XS_ABSORPTION,k) + rx % xs(t, k - j + 1) - end do - - ! Keep track of this reaction for easy searching later - if (t == 1) then - i_fission = i_fission + 1 - this % index_fission(i_fission) = i - end if - end if ! fission - end do ! temperature - end associate ! rx - end do ! reactions - - ! Calculate nu-fission cross section - do t = 1, n_temperature - if (this % fissionable) then - n_grid = size(this % grid(t) % energy) - do i = 1, n_grid - this % xs(t) % value(XS_NU_FISSION,i) = & - this % nu(this % grid(t) % energy(i), EMISSION_TOTAL) * & - this % xs(t) % value(XS_FISSION,i) - end do - end if - end do - end subroutine nuclide_create_derived - -!=============================================================================== -! NUCLIDE_NU is an interface to the number of fission neutrons produced -!=============================================================================== - - pure function nuclide_nu(this, E, emission_mode, group) result(nu) - class(Nuclide), intent(in) :: this - real(8), intent(in) :: E - integer, intent(in) :: emission_mode - integer, optional, intent(in) :: group - real(8) :: nu - - interface - pure function nuclide_nu_c(ptr, E, emission_mode, group) result(nu) bind(C) - import C_PTR, C_DOUBLE, C_INT - type(C_PTR), value, intent(in) :: ptr - real(C_DOUBLE), value, intent(in) :: E - integer(C_INT), value, intent(in) :: emission_mode - integer(C_INT), value, intent(in) :: group - real(C_DOUBLE) :: nu - end function - end interface - - if (present(group)) then - nu = nuclide_nu_c(this % ptr, E, emission_mode, group) - else - nu = nuclide_nu_c(this % ptr, E, emission_mode, 0) - end if - end function nuclide_nu - - subroutine nuclide_init_grid(this) - class(nuclide), intent(inout) :: this - interface - subroutine nuclide_init_grid_c(ptr) bind(C) - import C_PTR - type(C_PTR), value :: ptr - end subroutine - end interface - call nuclide_init_grid_c(this % ptr) - end subroutine - -!=============================================================================== -! NUCLIDE_CALCULATE_XS determines microscopic cross sections for the nuclide -! at the energy of the given particle -!=============================================================================== - - subroutine nuclide_calculate_xs(this, i_sab, E, i_log_union, sqrtkT, & - sab_frac) - class(Nuclide), intent(in) :: this ! Nuclide object - integer, intent(in) :: i_sab ! index into sab_tables array - real(8), intent(in) :: E ! energy - integer, intent(in) :: i_log_union ! index into logarithmic mapping array or - ! material union energy grid - real(8), intent(in) :: sqrtkT ! square root of kT, material dependent - real(8), intent(in) :: sab_frac ! fraction of atoms affected by S(a,b) - - interface - subroutine nuclide_calculate_xs_c(ptr, i_sab, E, i_log_union, sqrtkT, sab_frac) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i_sab - real(C_DOUBLE), value :: E - integer(C_INT), value :: i_log_union - real(C_DOUBLE), value :: sqrtkT - real(C_DOUBLE), value :: sab_frac - end subroutine - end interface - - call nuclide_calculate_xs_c(this % ptr, i_sab - 1, E, i_log_union, sqrtkT, sab_frac) - end subroutine nuclide_calculate_xs - -!=============================================================================== -! NUCLIDE_CALCULATE_ELASTIC_XS precalculates the free atom elastic scattering -! cross section. Normally it is not needed until a collision actually occurs in -! a material. However, in the thermal and unresolved resonance regions, we have -! to calculate it early to adjust the total cross section correctly. -!=============================================================================== - - subroutine nuclide_calculate_elastic_xs(this) - class(Nuclide), intent(in) :: this - interface - subroutine nuclide_calculate_elastic_xs_c(ptr) bind(C) - import C_PTR - type(C_PTR), value :: ptr - end subroutine - end interface - call nuclide_calculate_elastic_xs_c(this % ptr) - end subroutine nuclide_calculate_elastic_xs !=============================================================================== ! FREE_MEMORY_NUCLIDE deallocates global arrays defined in this module @@ -635,80 +101,9 @@ contains end subroutine end interface - ! Deallocate cross section data, listings, and cache - if (allocated(nuclides)) then - ! First call the clear routines - deallocate(nuclides) - call nuclides_clear() - end if - n_nuclides = 0 - + call nuclides_clear() call library_clear() end subroutine free_memory_nuclide -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_get_nuclide_index(name, index) result(err) bind(C) - ! Return the index in the nuclides array of a nuclide with a given name - character(kind=C_CHAR), intent(in) :: name(*) - integer(C_INT), intent(out) :: index - integer(C_INT) :: err - - character(:), allocatable :: name_ - - ! Copy array of C_CHARs to normal Fortran string - name_ = to_f_string(name) - - err = 0 - if (allocated(nuclides)) then - index = nuclide_map_get(name) - if (index == -1) then - err = E_DATA - call set_errmsg("No nuclide named '" // trim(name_) // & - "' has been loaded.") - end if - else - err = E_ALLOCATE - call set_errmsg("Memory for nuclides has not been allocated.") - end if - end function openmc_get_nuclide_index - - - function openmc_nuclide_name(index, name) result(err) bind(C) - ! Return the name of a nuclide with a given index - integer(C_INT), value, intent(in) :: index - type(c_ptr), intent(out) :: name - integer(C_INT) :: err - - character(C_CHAR), pointer :: name_ - - err = E_UNASSIGNED - if (allocated(nuclides)) then - if (index >= 1 .and. index <= size(nuclides)) then - name_ => nuclides(index) % name(1:1) - name = C_LOC(name_) - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in nuclides array is out of bounds.") - end if - else - err = E_ALLOCATE - call set_errmsg("Memory for nuclides has not been allocated yet.") - end if - end function openmc_nuclide_name - - subroutine extend_nuclides() bind(C) - type(Nuclide), allocatable :: new_nuclides(:) - - ! allocate extra space in nuclides array - allocate(new_nuclides(n_nuclides + 1)) - new_nuclides(1:n_nuclides) = nuclides(:) - call move_alloc(FROM=new_nuclides, TO=nuclides) - n_nuclides = n_nuclides + 1 - end subroutine - end module nuclide_header diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 index 31ee5ff66..c05f9f3df 100644 --- a/src/particle_restart.F90 +++ b/src/particle_restart.F90 @@ -7,14 +7,13 @@ module particle_restart use error, only: write_message use hdf5_interface, only: file_open, file_close, read_dataset, HID_T use mgxs_interface, only: energy_bin_avg - use nuclide_header, only: micro_xs, n_nuclides + use nuclide_header, only: micro_xs, nuclides_size use particle_header use photon_header, only: micro_photon_xs, n_elements use random_lcg, only: set_particle_seed use settings use simulation_header use tally_header, only: n_tallies - use tracking, only: transport implicit none private @@ -49,7 +48,7 @@ contains verbosity = 10 !$omp parallel - allocate(micro_xs(n_nuclides)) + allocate(micro_xs(nuclides_size())) allocate(micro_photon_xs(n_elements)) !$omp end parallel call set_micro_xs() diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 deleted file mode 100644 index bc0fc8f32..000000000 --- a/src/reaction_header.F90 +++ /dev/null @@ -1,270 +0,0 @@ -module reaction_header - - use, intrinsic :: ISO_C_BINDING - - use hdf5_interface - - implicit none - private - -!=============================================================================== -! REACTION contains the cross-section and secondary energy and angle -! distributions for a single reaction in a continuous-energy ACE-format table -!=============================================================================== - - type, public :: Reaction - type(C_PTR) :: ptr - integer(C_INT) :: MT ! ENDF MT value - real(C_DOUBLE) :: Q_value ! Reaction Q value - logical(C_BOOL) :: scatter_in_cm ! scattering system in center-of-mass? - logical(C_BOOL) :: redundant ! redundant reaction? - contains - procedure :: init - procedure :: mt_ - procedure :: q_value_ - procedure :: scatter_in_cm_ - procedure :: redundant_ - procedure :: product_decay_rate - procedure :: product_emission_mode - procedure :: product_particle - procedure :: product_sample - procedure :: product_yield - procedure :: products_size - procedure :: sample_elastic_mu - procedure :: xs - procedure :: xs_size - procedure :: xs_threshold - end type Reaction - - interface - function nuclide_reaction(nuc_ptr, i_rx) result(ptr) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: nuc_ptr - integer(C_INT), value :: i_rx - type(C_PTR) :: ptr - end function - - function reaction_mt(ptr) result(mt) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT) :: mt - end function - - function reaction_q_value(ptr) result(q_value) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE) :: q_value - end function - - function reaction_scatter_in_cm(ptr) result(b) bind(C) - import C_PTR, C_BOOL - type(C_PTR), value :: ptr - logical(C_BOOL) :: b - end function - - function reaction_redundant(ptr) result(b) bind(C) - import C_PTR, C_BOOL - type(C_PTR), value :: ptr - logical(C_BOOL) :: b - end function - - pure function reaction_product_decay_rate(ptr, product) result(rate) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: product - real(C_DOUBLE) :: rate - end function - - pure function reaction_product_emission_mode(ptr, product) result(m) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: product - integer(C_INT) :: m - end function - - pure function reaction_product_particle(ptr, product) result(particle) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: product - integer(C_INT) :: particle - end function - - subroutine reaction_product_sample(ptr, product, E_in, E_out, mu) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: product - real(C_DOUBLE), value :: E_in - real(C_DOUBLE), intent(out) :: E_out - real(C_DOUBLE), intent(out) :: mu - end subroutine - - pure function reaction_product_yield(ptr, product, E) result(val) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: product - real(C_DOUBLE), intent(in), value :: E - real(C_DOUBLE) :: val - end function - - pure function reaction_products_size(ptr) result(sz) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT) :: sz - end function - - function reaction_sample_elastic_mu(ptr, E) result(mu) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE), value :: E - real(C_DOUBLE) :: mu - end function - - function reaction_xs(ptr, temperature, energy) result(xs) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: temperature - integer(C_INT), value :: energy - real(C_DOUBLE) :: xs - end function - - function reaction_xs_size(ptr, temperature) result(sz) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), value :: temperature - integer(C_INT) :: sz - end function - - function reaction_xs_threshold(ptr, temperature) result(threshold) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), value :: temperature - integer(C_INT) :: threshold - end function - end interface - -contains - - subroutine init(this, nuc_ptr, i_rx) - class(Reaction), intent(inout) :: this - type(C_PTR), intent(in) :: nuc_ptr - integer(C_INT), intent(in) :: i_rx - - this % ptr = nuclide_reaction(nuc_ptr, i_rx) - this % MT = reaction_mt(this % ptr) - this % Q_value = reaction_q_value(this % ptr) - this % scatter_in_cm = reaction_scatter_in_cm(this % ptr) - this % redundant = reaction_redundant(this % ptr) - end subroutine - - function mt_(this) result(mt) - class(Reaction), intent(in) :: this - integer(C_INT) :: MT - - mt = reaction_mt(this % ptr) - end function - - function q_value_(this) result(q_value) - class(Reaction), intent(in) :: this - real(C_DOUBLE) :: q_value - - q_value = reaction_q_value(this % ptr) - end function - - function scatter_in_cm_(this) result(cm) - class (Reaction), intent(in) :: this - logical(C_BOOL) :: cm - - cm = reaction_scatter_in_cm(this % ptr) - end function - - function redundant_(this) result(redundant) - class (Reaction), intent(in) :: this - logical(C_BOOL) :: redundant - - redundant = reaction_redundant(this % ptr) - end function - - pure function product_decay_rate(this, product) result(rate) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: product - real(C_DOUBLE) :: rate - - rate = reaction_product_decay_rate(this % ptr, product) - end function - - pure function product_emission_mode(this, product) result(m) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: product - integer(C_INT) :: m - - m = reaction_product_emission_mode(this % ptr, product) - end function - - pure function product_particle(this, product) result(p) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: product - integer(C_INT) :: p - - p = reaction_product_particle(this % ptr, product) - end function - - subroutine product_sample(this, product, E_in, E_out, mu) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: product - real(C_DOUBLE), intent(in) :: E_in - real(C_DOUBLE), intent(out) :: E_out - real(C_DOUBLE), intent(out) :: mu - - call reaction_product_sample(this % ptr, product, E_in, E_out, mu) - end subroutine - - pure function product_yield(this, product, E) result(val) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: product - real(C_DOUBLE), intent(in) :: E - real(C_DOUBLE) :: val - - val = reaction_product_yield(this % ptr, product, E) - end function - - pure function products_size(this) result(sz) - class(Reaction), intent(in) :: this - integer(C_INT) :: sz - - sz = reaction_products_size(this % ptr) - end function - - function sample_elastic_mu(this, E) result(mu) - class(Reaction), intent(in) :: this - real(C_DOUBLE), intent(in) :: E - real(C_DOUBLE) :: mu - - mu = reaction_sample_elastic_mu(this % ptr, E) - end function - - function xs(this, temperature, energy) result(val) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: temperature - integer(C_INT), intent(in) :: energy - real(C_DOUBLE) :: val - - val = reaction_xs(this % ptr, temperature, energy) - end function - - function xs_size(this, temperature) result(sz) - class(Reaction), intent(in) :: this - integer(C_INT) :: temperature - integer(C_INT) :: sz - - sz = reaction_xs_size(this % ptr, temperature) - end function - - function xs_threshold(this, temperature) result(val) - class(Reaction), intent(in) :: this - integer(C_INT), intent(in) :: temperature - integer(C_INT) :: val - - val = reaction_xs_threshold(this % ptr, temperature) - end function - -end module reaction_header diff --git a/src/simulation.F90 b/src/simulation.F90 index f088481dd..2a521807d 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -2,7 +2,7 @@ module simulation use, intrinsic :: ISO_C_BINDING - use nuclide_header, only: micro_xs, n_nuclides + use nuclide_header, only: micro_xs, nuclides_size use photon_header, only: micro_photon_xs, n_elements implicit none @@ -18,7 +18,7 @@ contains !$omp parallel ! Allocate array for microscopic cross section cache - allocate(micro_xs(n_nuclides)) + allocate(micro_xs(nuclides_size())) allocate(micro_photon_xs(n_elements)) !$omp end parallel diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 003eb438d..e7cdcbc1d 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -4,13 +4,10 @@ module tally use bank_header use constants - use dict_header, only: EMPTY - use error, only: fatal_error use geometry_header use material_header use message_passing use mgxs_interface - use nuclide_header use settings use simulation_header use tally_derivative_header diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 14ac94107..f5e2c566d 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -7,7 +7,6 @@ module tally_header use dict_header, only: DictIntInt use hdf5_interface, only: HID_T, HSIZE_T use message_passing, only: n_procs, master - use nuclide_header, only: nuclide_map_get use settings, only: reduce_tallies, run_mode use stl_vector, only: VectorInt use string, only: to_lower, to_f_string, str_to_int, to_str, to_c_string From 7e1abc439bdbb167b303685f345bfcbcc9c77f55 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Feb 2019 06:27:56 -0600 Subject: [PATCH 04/10] Remove three Fortran modules and a bunch of constants --- CMakeLists.txt | 3 - src/constants.F90 | 109 +-------------- src/eigenvalue.F90 | 15 --- src/input_xml.F90 | 7 +- src/math.F90 | 20 --- src/mesh_header.F90 | 296 ----------------------------------------- src/nuclide_header.F90 | 2 +- 7 files changed, 7 insertions(+), 445 deletions(-) delete mode 100644 src/eigenvalue.F90 delete mode 100644 src/math.F90 delete mode 100644 src/mesh_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 96547ebbf..2c625517a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -308,7 +308,6 @@ add_library(libopenmc SHARED src/dagmc_header.F90 src/constants.F90 src/dict_header.F90 - src/eigenvalue.F90 src/error.F90 src/geometry.F90 src/geometry_header.F90 @@ -316,8 +315,6 @@ add_library(libopenmc SHARED src/initialize.F90 src/input_xml.F90 src/material_header.F90 - src/math.F90 - src/mesh_header.F90 src/message_passing.F90 src/mgxs_interface.F90 src/nuclide_header.F90 diff --git a/src/constants.F90 b/src/constants.F90 index e68f5a6c3..9061bb800 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -7,16 +7,8 @@ module constants ! ============================================================================ ! VERSIONING NUMBERS - ! OpenMC major, minor, and release numbers - integer, parameter :: VERSION_MAJOR = 0 - integer, parameter :: VERSION_MINOR = 10 - integer, parameter :: VERSION_RELEASE = 0 - integer, parameter :: & - VERSION(3) = [VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE] - ! Version numbers for binary files - integer, parameter :: VERSION_TRACK(2) = [2, 0] - integer, parameter :: VERSION_VOLUME(2) = [1, 0] + integer, parameter :: VERSION_TRACK(2) = [2, 0] ! ============================================================================ ! ADJUSTABLE PARAMETERS @@ -24,15 +16,6 @@ module constants ! NOTE: This is the only section of the constants module that should ever be ! adjusted. Modifying constants in other sections may cause the code to fail. - ! Significance level for confidence intervals - real(8), parameter :: CONFIDENCE_LEVEL = 0.95_8 - - ! Used for surface current tallies - real(8), parameter :: TINY_BIT = 1e-8_8 - - ! Maximum number of collisions/crossings - integer, parameter :: MAX_EVENTS = 1000000 - ! Maximum number of secondary particles created integer, parameter :: MAX_SECONDARY = 1000 @@ -50,24 +33,12 @@ module constants real(8), parameter :: & PI = 3.1415926535898_8, & ! pi - MASS_NEUTRON = 1.00866491588_8, & ! mass of a neutron in amu - MASS_NEUTRON_EV = 939.5654133e6_8, & ! mass of a neutron in eV/c^2 - C_LIGHT = 2.99792458e8_8, & ! speed of light in m/s - K_BOLTZMANN = 8.6173303e-5_8, & ! Boltzmann constant in eV/K INFINITY = huge(0.0_8), & ! positive infinity ZERO = 0.0_8, & HALF = 0.5_8, & ONE = 1.0_8, & TWO = 2.0_8 - ! Electron subshell labels - character(3), parameter :: SUBSHELLS(39) = [ & - 'K ', 'L1 ', 'L2 ', 'L3 ', 'M1 ', 'M2 ', 'M3 ', 'M4 ', 'M5 ', & - 'N1 ', 'N2 ', 'N3 ', 'N4 ', 'N5 ', 'N6 ', 'N7 ', 'O1 ', 'O2 ', & - 'O3 ', 'O4 ', 'O5 ', 'O6 ', 'O7 ', 'O8 ', 'O9 ', 'P1 ', 'P2 ', & - 'P3 ', 'P4 ', 'P5 ', 'P6 ', 'P7 ', 'P8 ', 'P9 ', 'P10', 'P11', & - 'Q1 ', 'Q2 ', 'Q3 '] - ! ============================================================================ ! GEOMETRY-RELATED CONSTANTS @@ -99,55 +70,11 @@ module constants ELECTRON = 2, & POSITRON = 3 - ! Reaction types - integer, parameter :: & - TOTAL_XS = 1, ELASTIC = 2, N_NONELASTIC = 3, N_LEVEL = 4, MISC = 5, & - N_2ND = 11, N_2N = 16, N_3N = 17, N_FISSION = 18, N_F = 19, & - N_NF = 20, N_2NF = 21, N_NA = 22, N_N3A = 23, N_2NA = 24, & - N_3NA = 25, N_NP = 28, N_N2A = 29, N_2N2A = 30, N_ND = 32, & - N_NT = 33, N_N3HE = 34, N_ND2A = 35, N_NT2A = 36, N_4N = 37, & - N_3NF = 38, N_2NP = 41, N_3NP = 42, N_N2P = 44, N_NPA = 45, & - N_N1 = 51, N_N40 = 90, N_NC = 91, N_DISAPPEAR = 101, N_GAMMA = 102, & - N_P = 103, N_D = 104, N_T = 105, N_3HE = 106, N_A = 107, & - N_2A = 108, N_3A = 109, N_2P = 111, N_PA = 112, N_T2A = 113, & - N_D2A = 114, N_PD = 115, N_PT = 116, N_DA = 117, N_5N = 152, & - N_6N = 153, N_2NT = 154, N_TA = 155, N_4NP = 156, N_3ND = 157, & - N_NDA = 158, N_2NPA = 159, N_7N = 160, N_8N = 161, N_5NP = 162, & - N_6NP = 163, N_7NP = 164, N_4NA = 165, N_5NA = 166, N_6NA = 167, & - N_7NA = 168, N_4ND = 169, N_5ND = 170, N_6ND = 171, N_3NT = 172, & - N_4NT = 173, N_5NT = 174, N_6NT = 175, N_2N3HE = 176, N_3N3HE = 177, & - N_4N3HE = 178, N_3N2P = 179, N_3N3A = 180, N_3NPA = 181, N_DT = 182, & - N_NPD = 183, N_NPT = 184, N_NDT = 185, N_NP3HE = 186, N_ND3HE = 187, & - N_NT3HE = 188, N_NTA = 189, N_2N2P = 190, N_P3HE = 191, N_D3HE = 192, & - N_3HEA = 193, N_4N2P = 194, N_4N2A = 195, N_4NPA = 196, N_3P = 197, & - N_N3P = 198, N_3N2PA = 199, N_5N2P = 200, N_P0 = 600, N_PC = 649, & - N_D0 = 650, N_DC = 699, N_T0 = 700, N_TC = 749, N_3HE0 = 750, & - N_3HEC = 799, N_A0 = 800, N_AC = 849, N_2N0 = 875, N_2NC = 891, & - COHERENT = 502, INCOHERENT = 504, PHOTOELECTRIC = 522, & - PAIR_PROD_ELEC = 515, PAIR_PROD = 516, PAIR_PROD_NUC = 517 - - ! Depletion reactions - integer, parameter :: DEPLETION_RX(6) = [N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N] - ! MGXS Table Types integer, parameter :: & MGXS_ISOTROPIC = 1, & ! Isotropically Weighted Data MGXS_ANGLE = 2 ! Data by Angular Bins - ! Secondary particle emission type - integer, parameter :: & - EMISSION_PROMPT = 1, & ! Prompt emission of secondary particle - EMISSION_DELAYED = 2, & ! Delayed emission of secondary particle - EMISSION_TOTAL = 3 ! Yield represents total emission (prompt + delayed) - - ! Maximum number of partial fission reactions - integer, parameter :: PARTIAL_FISSION_MAX = 4 - - ! Temperature treatment method - integer, parameter :: & - TEMPERATURE_NEAREST = 1, & - TEMPERATURE_INTERPOLATION = 2 - ! ============================================================================ ! TALLY-RELATED CONSTANTS @@ -206,21 +133,6 @@ module constants FILTER_MATERIAL = 2, & FILTER_CELL = 3 - ! Tally surface current directions - integer, parameter :: & - OUT_LEFT = 1, & ! x min - IN_LEFT = 2, & ! x min - OUT_RIGHT = 3, & ! x max - IN_RIGHT = 4, & ! x max - OUT_BACK = 5, & ! y min - IN_BACK = 6, & ! y min - OUT_FRONT = 7, & ! y max - IN_FRONT = 8, & ! y max - OUT_BOTTOM = 9, & ! z min - IN_BOTTOM = 10, & ! z min - OUT_TOP = 11, & ! z max - IN_TOP = 12 ! z max - ! Global tally parameters integer, parameter :: N_GLOBAL_TALLIES = 4 integer, parameter :: & @@ -235,25 +147,6 @@ module constants DIFF_NUCLIDE_DENSITY = 2, & DIFF_TEMPERATURE = 3 - - ! Mgxs::get_xs enumerated types - integer(C_INT), parameter :: & - MG_GET_XS_TOTAL = 0, & - MG_GET_XS_ABSORPTION = 1, & - MG_GET_XS_INVERSE_VELOCITY = 2, & - MG_GET_XS_DECAY_RATE = 3, & - MG_GET_XS_SCATTER = 4, & - MG_GET_XS_SCATTER_MULT = 5, & - MG_GET_XS_SCATTER_FMU_MULT = 6, & - MG_GET_XS_SCATTER_FMU = 7, & - MG_GET_XS_FISSION = 8, & - MG_GET_XS_KAPPA_FISSION = 9, & - MG_GET_XS_PROMPT_NU_FISSION = 10, & - MG_GET_XS_DELAYED_NU_FISSION = 11, & - MG_GET_XS_NU_FISSION = 12, & - MG_GET_XS_CHI_PROMPT = 13, & - MG_GET_XS_CHI_DELAYED = 14 - ! ============================================================================ ! RANDOM NUMBER STREAM CONSTANTS diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 deleted file mode 100644 index 4de3ae8a6..000000000 --- a/src/eigenvalue.F90 +++ /dev/null @@ -1,15 +0,0 @@ -module eigenvalue - - use, intrinsic :: ISO_C_BINDING - - implicit none - - interface - function openmc_get_keff(k_combined) result(err) bind(C) - import C_INT, C_DOUBLE - real(C_DOUBLE), intent(out) :: k_combined(2) - integer(C_INT) :: err - end function - end interface - -end module eigenvalue diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 94a1d6190..574d2949b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,12 +11,10 @@ module input_xml #endif use hdf5_interface use material_header - use mesh_header use message_passing use mgxs_interface use nuclide_header use photon_header - use random_lcg, only: prn use settings use stl_vector, only: VectorInt, VectorReal, VectorChar use string, only: to_lower, to_str, str_to_int, & @@ -424,6 +422,11 @@ contains import C_PTR type(C_PTR) :: node_ptr end subroutine + + subroutine read_meshes(node_ptr) bind(C) + import C_PTR + type(C_PTR) :: node_ptr + end subroutine end interface ! Check if tallies.xml exists diff --git a/src/math.F90 b/src/math.F90 deleted file mode 100644 index 81856c93c..000000000 --- a/src/math.F90 +++ /dev/null @@ -1,20 +0,0 @@ -module math - - use, intrinsic :: ISO_C_BINDING - - implicit none - private - public :: t_percentile - - interface - - pure function t_percentile(p, df) bind(C) result(t) - use ISO_C_BINDING - implicit none - real(C_DOUBLE), value, intent(in) :: p - integer(C_INT), value, intent(in) :: df - real(C_DOUBLE) :: t - end function t_percentile - end interface - -end module math diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 deleted file mode 100644 index a0c14077a..000000000 --- a/src/mesh_header.F90 +++ /dev/null @@ -1,296 +0,0 @@ -module mesh_header - - use, intrinsic :: ISO_C_BINDING - - implicit none - -!=============================================================================== -! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by -! congruent squares or cubes -!=============================================================================== - - type :: RegularMesh - type(C_PTR) :: ptr - contains - procedure :: id => regular_id - procedure :: volume_frac => regular_volume_frac - procedure :: n_dimension => regular_n_dimension - procedure :: dimension => regular_dimension - procedure :: lower_left => regular_lower_left - procedure :: upper_right => regular_upper_right - procedure :: width => regular_width - - procedure :: get_bin => regular_get_bin - procedure :: get_indices => regular_get_indices - procedure :: get_bin_from_indices => regular_get_bin_from_indices - procedure :: get_indices_from_bin => regular_get_indices_from_bin - end type RegularMesh - - interface - function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), optional, intent(out) :: index_start - integer(C_INT32_T), optional, intent(out) :: index_end - integer(C_INT) :: err - end function openmc_extend_meshes - - function openmc_get_mesh_index(id, index) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value :: id - integer(C_INT32_T), intent(out) :: index - integer(C_INT) :: err - end function openmc_get_mesh_index - - function openmc_mesh_get_id(index, id) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value :: index - integer(C_INT32_T), intent(out) :: id - integer(C_INT) :: err - end function openmc_mesh_get_id - - function openmc_mesh_set_id(index, id) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: id - integer(C_INT) :: err - end function openmc_mesh_set_id - - function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C) - import C_INT32_T, C_PTR, C_INT - integer(C_INT32_T), value, intent(in) :: index - type(C_PTR), intent(out) :: dims - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - end function openmc_mesh_get_dimension - - function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C) - import C_INT32_T, C_INT - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), value, intent(in) :: n - integer(C_INT), intent(in) :: dims(n) - integer(C_INT) :: err - end function openmc_mesh_set_dimension - - function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C) - import C_INT32_T, C_PTR, C_INT - integer(C_INT32_T), value, intent(in) :: index - type(C_PTR), intent(out) :: ll - type(C_PTR), intent(out) :: ur - type(C_PTR), intent(out) :: width - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - end function openmc_mesh_get_params - - function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C) - import C_INT32_T, C_INT, C_DOUBLE - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT), value, intent(in) :: n - real(C_DOUBLE), intent(in), optional :: ll(n) - real(C_DOUBLE), intent(in), optional :: ur(n) - real(C_DOUBLE), intent(in), optional :: width(n) - integer(C_INT) :: err - end function openmc_mesh_set_params - - function mesh_id(ptr) result(id) bind(C) - import C_PTR, C_INT32_T - type(C_PTR), value :: ptr - integer(C_INT32_T) :: id - end function - - function mesh_volume_frac(ptr) result(volume_frac) bind(C) - import C_PTR, C_DOUBLE - type(C_PTR), value :: ptr - real(C_DOUBLE) :: volume_frac - end function - - function mesh_n_dimension(ptr) result(n) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT) :: n - end function - - function mesh_dimension(ptr, i) result(d) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: ptr - integer(C_INT), value :: i - integer(C_INT) :: d - end function - - function mesh_lower_left(ptr, i) result(ll) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: ll - end function - - function mesh_upper_right(ptr, i) result(ur) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: ur - end function - - function mesh_width(ptr, i) result(w) bind(C) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), value :: ptr - integer(C_INT), value :: i - real(C_DOUBLE) :: w - end function - - pure function mesh_get_bin(ptr, xyz) result(bin) bind(C) - import C_PTR, C_DOUBLE, C_INT - type(C_PTR), intent(in), value :: ptr - real(C_DOUBLE), intent(in) :: xyz(*) - integer(C_INT) :: bin - end function - - pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in) :: ijk(*) - integer(C_INT) :: bin - end function - - pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) - import C_PTR, C_DOUBLE, C_INT, C_BOOL - type(C_PTR), intent(in), value :: ptr - real(C_DOUBLE), intent(in) :: xyz(*) - integer(C_INT), intent(out) :: ijk(*) - logical(C_BOOL), intent(out) :: in_mesh - end subroutine - - pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: ptr - integer(C_INT), intent(in), value :: bin - integer(C_INT), intent(out) :: ijk(*) - end subroutine - - function mesh_ptr(i) result(ptr) bind(C) - import C_INT, C_PTR - integer(C_INT), value :: i - type(C_PTR) :: ptr - end function - - subroutine read_meshes(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine - - function n_meshes() result(n) bind(C) - import C_INT - integer(C_INT) :: n - end function - end interface - -contains - - function meshes(i) result(m) - integer, intent(in) :: i - type(RegularMesh) :: m - - m % ptr = mesh_ptr(i) - end function - - function regular_id(this) result(id) - class(RegularMesh), intent(in) :: this - integer(C_INT32_T) :: id - id = mesh_id(this % ptr) - end function - - function regular_volume_frac(this) result(volume_frac) - class(RegularMesh), intent(in) :: this - real(C_DOUBLE) :: volume_frac - volume_frac = mesh_volume_frac(this % ptr) - end function - - function regular_n_dimension(this) result(n) - class(RegularMesh), intent(in) :: this - integer(C_INT) :: n - n = mesh_n_dimension(this % ptr) - end function - - function regular_dimension(this, i) result(d) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - integer(C_INT) :: d - d = mesh_dimension(this % ptr, i) - end function - - function regular_lower_left(this, i) result(ll) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: ll - ll = mesh_lower_left(this % ptr, i) - end function - - function regular_upper_right(this, i) result(ur) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: ur - ur = mesh_upper_right(this % ptr, i) - end function - - function regular_width(this, i) result(w) - class(RegularMesh), intent(in) :: this - integer(C_INT), intent(in) :: i - real(C_DOUBLE) :: w - w = mesh_width(this % ptr, i) - end function - -!=============================================================================== -! GET_MESH_BIN determines the tally bin for a particle in a structured mesh -!=============================================================================== - - pure subroutine regular_get_bin(this, xyz, bin) - class(RegularMesh), intent(in) :: this - real(8), intent(in) :: xyz(:) ! coordinates - integer, intent(out) :: bin ! tally bin - - bin = mesh_get_bin(this % ptr, xyz) - end subroutine - -!=============================================================================== -! GET_MESH_INDICES determines the indices of a particle in a structured mesh -!=============================================================================== - - pure subroutine regular_get_indices(this, xyz, ijk, in_mesh) - class(RegularMesh), intent(in) :: this - real(8), intent(in) :: xyz(:) ! coordinates to check - integer, intent(out) :: ijk(:) ! indices in mesh - logical, intent(out) :: in_mesh ! were given coords in mesh? - - logical(C_BOOL) :: in_mesh_ - - call mesh_get_indices(this % ptr, xyz, ijk, in_mesh_) - in_mesh = in_mesh_ - end subroutine regular_get_indices - -!=============================================================================== -! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number -! for use in a TallyObject results array -!=============================================================================== - - pure function regular_get_bin_from_indices(this, ijk) result(bin) - class(RegularMesh), intent(in) :: this - integer, intent(in) :: ijk(:) - integer :: bin - - bin = mesh_get_bin_from_indices(this % ptr, ijk) - end function regular_get_bin_from_indices - -!=============================================================================== -! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to -! (i), (i,j), or (i,j,k) indices -!=============================================================================== - - pure subroutine regular_get_indices_from_bin(this, bin, ijk) - class(RegularMesh), intent(in) :: this - integer, intent(in) :: bin - integer, intent(out) :: ijk(:) - - call mesh_get_indices_from_bin(this % ptr, bin, ijk) - end subroutine regular_get_indices_from_bin - -end module mesh_header diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index b0a3f6444..57a283c65 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -31,7 +31,7 @@ module nuclide_header ! Cross sections for depletion reactions (note that these are not stored in ! macroscopic cache) - real(C_DOUBLE) :: reaction(size(DEPLETION_RX)) + real(C_DOUBLE) :: reaction(6) ! Indicies and factors needed to compute cross sections from the data tables integer(C_INT) :: index_grid ! Index on nuclide energy grid From 24c84e394ca3b075bbb11c188d113470a833a2fc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Feb 2019 06:39:48 -0600 Subject: [PATCH 05/10] Make sure output() doesn't always write to stderr --- src/error.cpp | 4 ++-- src/output.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/error.cpp b/src/error.cpp index 38001b72a..cc63c9bc3 100644 --- a/src/error.cpp +++ b/src/error.cpp @@ -33,7 +33,7 @@ void output(const std::string& message, std::ostream& out, int indent) while (i_start < length) { if (length - i_start < line_len) { // Remainder of message will fit on line - std::cerr << message.substr(i_start) << '\n'; + out << message.substr(i_start) << '\n'; break; } else { @@ -42,7 +42,7 @@ void output(const std::string& message, std::ostream& out, int indent) auto pos = s.find_last_of(' '); // Write up to last space, or whole line if no space is present - std::cerr << s.substr(0, pos) << '\n' << std::setw(indent) << " "; + out << s.substr(0, pos) << '\n' << std::setw(indent) << " "; // Advance starting position i_start += (pos == std::string::npos) ? line_len : pos + 1; diff --git a/src/output.cpp b/src/output.cpp index dbfb2ab7b..5f4a53cdc 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -65,7 +65,7 @@ void title() " ############### %%%%%%%%%%%%%%%%\n" << " ############ %%%%%%%%%%%%%%%\n" << " ######## %%%%%%%%%%%%%%\n" << - " %%%%%%%%%%%\n"; + " %%%%%%%%%%%\n\n"; // Write version information std::cout << @@ -88,7 +88,7 @@ void title() #ifdef _OPENMP // Write number of OpenMP threads - std::cout << " OpenMC Threads | " << omp_get_max_threads() << '\n'; + std::cout << " OpenMP Threads | " << omp_get_max_threads() << '\n'; #endif std::cout << '\n'; } From 968e9a95d4393490c0aa7a5f98c4638b7fee8819 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Feb 2019 12:43:05 -0600 Subject: [PATCH 06/10] Convert particle restart to C++ --- CMakeLists.txt | 3 +- include/openmc/capi.h | 1 - include/openmc/hdf5_interface.h | 13 ++- include/openmc/particle_restart.h | 10 +++ include/openmc/simulation.h | 3 + src/constants.F90 | 12 --- src/main.cpp | 4 +- src/particle_restart.F90 | 142 ------------------------------ src/particle_restart.cpp | 109 +++++++++++++++++++++++ src/simulation.cpp | 2 - 10 files changed, 137 insertions(+), 162 deletions(-) create mode 100644 include/openmc/particle_restart.h delete mode 100644 src/particle_restart.F90 create mode 100644 src/particle_restart.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 2c625517a..43e250c34 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -319,10 +319,8 @@ add_library(libopenmc SHARED src/mgxs_interface.F90 src/nuclide_header.F90 src/particle_header.F90 - src/particle_restart.F90 src/photon_header.F90 src/pugixml/pugixml_f.F90 - src/random_lcg.F90 src/relaxng src/settings.F90 src/simulation_header.F90 @@ -370,6 +368,7 @@ add_library(libopenmc SHARED src/nuclide.cpp src/output.cpp src/particle.cpp + src/particle_restart.cpp src/photon.cpp src/physics.cpp src/physics_common.cpp diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 80d431063..5ec852d2d 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -84,7 +84,6 @@ extern "C" { int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh); int openmc_next_batch(int* status); int openmc_nuclide_name(int index, const char** name); - int openmc_particle_restart(); int openmc_plot_geometry(); int openmc_reset(); int openmc_run(); diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 1fc378bde..f55b14171 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -223,8 +223,9 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) // Templates/overloads for read_dataset and related methods //============================================================================== -template -void read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false) +template inline +std::enable_if_t>::value> +read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false) { read_dataset(obj_id, name, H5TypeMap::type_id, &buffer, indep); } @@ -242,6 +243,14 @@ read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false) str = std::string{buffer, n}; } +// array version +template inline void +read_dataset(hid_t dset, const char* name, std::array& buffer, bool indep=false) +{ + read_dataset(dset, name, H5TypeMap::type_id, buffer.data(), indep); +} + +// vector version template void read_dataset(hid_t dset, std::vector& vec, bool indep=false) { diff --git a/include/openmc/particle_restart.h b/include/openmc/particle_restart.h new file mode 100644 index 000000000..24ea237a4 --- /dev/null +++ b/include/openmc/particle_restart.h @@ -0,0 +1,10 @@ +#ifndef OPENMC_PARTICLE_RESTART_H +#define OPENMC_PARTICLE_RESTART_H + +namespace openmc { + +void run_particle_restart(); + +} // namespace openmc + +#endif // OPENMC_PARTICLE_RESTART_H diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index f9e1f7c96..af333f714 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -83,6 +83,9 @@ void finalize_generation(); //! Determine overall generation number extern "C" int overall_generation(); +extern "C" void simulation_init_f(); +extern "C" void simulation_finalize_f(); + #ifdef OPENMC_MPI void broadcast_results(); #endif diff --git a/src/constants.F90 b/src/constants.F90 index 9061bb800..c1566d776 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -147,18 +147,6 @@ module constants DIFF_NUCLIDE_DENSITY = 2, & DIFF_TEMPERATURE = 3 - ! ============================================================================ - ! RANDOM NUMBER STREAM CONSTANTS - - integer(C_INT), bind(C, name='N_STREAMS') :: N_STREAMS - integer(C_INT), bind(C, name='STREAM_TRACKING') :: STREAM_TRACKING - integer(C_INT), bind(C, name='STREAM_TALLIES') :: STREAM_TALLIES - integer(C_INT), bind(C, name='STREAM_SOURCE') :: STREAM_SOURCE - integer(C_INT), bind(C, name='STREAM_URR_PTABLE') :: STREAM_URR_PTABLE - integer(C_INT), bind(C, name='STREAM_VOLUME') :: STREAM_VOLUME - integer(C_INT), bind(C, name='STREAM_PHOTON') :: STREAM_PHOTON - integer(C_INT64_T), parameter :: DEFAULT_SEED = 1_8 - ! ============================================================================ ! MISCELLANEOUS CONSTANTS diff --git a/src/main.cpp b/src/main.cpp index b757c67ff..936c89b78 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -5,6 +5,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/message_passing.h" +#include "openmc/particle_restart.h" #include "openmc/settings.h" @@ -36,7 +37,8 @@ int main(int argc, char* argv[]) { err = openmc_plot_geometry(); break; case RUN_MODE_PARTICLE: - if (mpi::master) err = openmc_particle_restart(); + if (mpi::master) run_particle_restart(); + err = 0; break; case RUN_MODE_VOLUME: err = openmc_calculate_volumes(); diff --git a/src/particle_restart.F90 b/src/particle_restart.F90 deleted file mode 100644 index c05f9f3df..000000000 --- a/src/particle_restart.F90 +++ /dev/null @@ -1,142 +0,0 @@ -module particle_restart - - use, intrinsic :: ISO_FORTRAN_ENV - - use bank_header, only: Bank - use constants - use error, only: write_message - use hdf5_interface, only: file_open, file_close, read_dataset, HID_T - use mgxs_interface, only: energy_bin_avg - use nuclide_header, only: micro_xs, nuclides_size - use particle_header - use photon_header, only: micro_photon_xs, n_elements - use random_lcg, only: set_particle_seed - use settings - use simulation_header - use tally_header, only: n_tallies - - implicit none - private - public :: openmc_particle_restart - -contains - -!=============================================================================== -! OPENMC_PARTICLE_RESTART is the main routine that runs the particle restart -!=============================================================================== - - function openmc_particle_restart() result(err) bind(C) - integer(C_INT) :: err - - integer(8) :: particle_seed - integer :: previous_run_mode - type(Particle) :: p - - interface - subroutine set_micro_xs() bind(C) - end subroutine - - subroutine print_particle(p) bind(C) - import Particle - type(Particle), intent(in) :: p - end subroutine - end interface - - err = 0 - - ! Set verbosity high - verbosity = 10 - - !$omp parallel - allocate(micro_xs(nuclides_size())) - allocate(micro_photon_xs(n_elements)) - !$omp end parallel - call set_micro_xs() - - ! Initialize the particle to be tracked - call particle_initialize(p) - - ! Read in the restart information - call read_particle_restart(p, previous_run_mode) - - ! Set all tallies to 0 for now (just tracking errors) - n_tallies = 0 - - ! Compute random number seed - select case (previous_run_mode) - case (MODE_EIGENVALUE) - particle_seed = (total_gen + overall_generation() - 1)*n_particles + p % id - case (MODE_FIXEDSOURCE) - particle_seed = p % id - end select - - call set_particle_seed(particle_seed) - - ! Transport neutron - ! FIXME - !call transport(p) - - ! Write output if particle made it - call print_particle(p) - - deallocate(micro_xs) - - end function openmc_particle_restart - -!=============================================================================== -! READ_PARTICLE_RESTART reads the particle restart file -!=============================================================================== - - subroutine read_particle_restart(p, previous_run_mode) - type(Particle), intent(inout) :: p - integer, intent(inout) :: previous_run_mode - - integer(HID_T) :: file_id - character(MAX_WORD_LEN) :: tempstr - - ! Write meessage - call write_message("Loading particle restart file " & - // trim(path_particle_restart) // "...", 5) - - ! Open file - file_id = file_open(path_particle_restart, 'r') - - ! Read data from file - call read_dataset(current_batch, file_id, 'current_batch') - call read_dataset(gen_per_batch, file_id, 'generations_per_batch') - call read_dataset(current_gen, file_id, 'current_generation') - call read_dataset(n_particles, file_id, 'n_particles') - call read_dataset(tempstr, file_id, 'run_mode') - select case (tempstr) - case ('eigenvalue') - previous_run_mode = MODE_EIGENVALUE - case ('fixed source') - previous_run_mode = MODE_FIXEDSOURCE - end select - call read_dataset(p % id, file_id, 'id') - call read_dataset(p % type, file_id, 'type') - call read_dataset(p % wgt, file_id, 'weight') - call read_dataset(p % E, file_id, 'energy') - call read_dataset(p % coord(1) % xyz, file_id, 'xyz') - call read_dataset(p % coord(1) % uvw, file_id, 'uvw') - - ! Set energy group and average energy in multi-group mode - if (.not. run_CE) then - p % g = int(p % E) - p % E = energy_bin_avg(p % g) - end if - - ! Set particle last attributes - p % last_wgt = p % wgt - p % last_xyz_current = p % coord(1)%xyz - p % last_xyz = p % coord(1)%xyz - p % last_uvw = p % coord(1)%uvw - p % last_E = p % E - p % last_g = p % g - - ! Close hdf5 file - call file_close(file_id) - - end subroutine read_particle_restart - -end module particle_restart diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp new file mode 100644 index 000000000..89781bb6c --- /dev/null +++ b/src/particle_restart.cpp @@ -0,0 +1,109 @@ +#include "openmc/particle_restart.h" + +#include "openmc/constants.h" +#include "openmc/hdf5_interface.h" +#include "openmc/mgxs_interface.h" +#include "openmc/nuclide.h" +#include "openmc/output.h" +#include "openmc/particle.h" +#include "openmc/random_lcg.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/tallies/tally.h" + +#include // for copy +#include +#include + +namespace openmc { + +void read_particle_restart(Particle& p, int& previous_run_mode) +{ + // Write meessage + write_message("Loading particle restart file " + + settings::path_particle_restart + "...", 5); + + // Open file + hid_t file_id = file_open(settings::path_particle_restart, 'r'); + + // Read data from file + read_dataset(file_id, "current_batch", simulation::current_batch); + read_dataset(file_id, "generations_per_batch", settings::gen_per_batch); + read_dataset(file_id, "current_generation", simulation::current_gen); + read_dataset(file_id, "n_particles", settings::n_particles); + std::string mode; + read_dataset(file_id, "run_mode", mode); + if (mode == "eigenvalue") { + previous_run_mode = RUN_MODE_EIGENVALUE; + } else if (mode == "fixed source") { + previous_run_mode = RUN_MODE_FIXEDSOURCE; + } + read_dataset(file_id, "id", p.id); + read_dataset(file_id, "type", p.type); + read_dataset(file_id, "weight", p.wgt); + read_dataset(file_id, "energy", p.E); + std::array x; + read_dataset(file_id, "xyz", x); + std::copy(x.data(), x.data() + 3, p.coord[0].xyz); + read_dataset(file_id, "uvw", x); + std::copy(x.data(), x.data() + 3, p.coord[0].uvw); + + // Set energy group and average energy in multi-group mode + if (!settings::run_CE) { + p.g = p.E; + p.E = data::energy_bin_avg[p.g - 1]; + } + + // Set particle last attributes + p.last_wgt = p.wgt; + std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz_current); + std::copy(p.coord[0].xyz, p.coord[0].xyz + 3, p.last_xyz); + std::copy(p.coord[0].uvw, p.coord[0].uvw + 3, p.last_uvw); + p.last_E = p.E; + p.last_g = p.g; + + // Close hdf5 file + file_close(file_id); +} + +void run_particle_restart() +{ + // Set verbosity high + settings::verbosity = 10; + + simulation_init_f(); + set_micro_xs(); + + // Initialize the particle to be tracked + Particle p; + p.initialize(); + + // Read in the restart information + int previous_run_mode; + read_particle_restart(p, previous_run_mode); + + // Set all tallies to 0 for now (just tracking errors) + model::tallies.clear(); + + // Compute random number seed + int64_t particle_seed; + switch (previous_run_mode) { + case RUN_MODE_EIGENVALUE: + particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id; + break; + case RUN_MODE_FIXEDSOURCE: + particle_seed = p.id; + break; + } + set_particle_seed(particle_seed); + + // Transport neutron + p.transport(); + + // Write output if particle made it + print_particle(&p); + + simulation_finalize_f(); +} + +} // namespace openmc diff --git a/src/simulation.cpp b/src/simulation.cpp index acc488353..1fdda3061 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -30,8 +30,6 @@ namespace openmc { extern "C" void accumulate_tallies(); extern "C" void allocate_tally_results(); extern "C" void setup_active_tallies(); -extern "C" void simulation_init_f(); -extern "C" void simulation_finalize_f(); extern "C" void write_tallies(); } // namespace openmc From 47c19353e4f704038c34cf3562d90947c19e6fdb Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 14 Feb 2019 15:30:00 -0600 Subject: [PATCH 07/10] Fix rotation matrix in transport() --- src/geometry_header.F90 | 5 ----- src/input_xml.F90 | 35 ----------------------------------- src/particle.cpp | 8 +++++--- 3 files changed, 5 insertions(+), 43 deletions(-) diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 563a60687..afb9e2c0d 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -82,11 +82,6 @@ module geometry_header type Cell type(C_PTR) :: ptr - - ! Rotation matrix and translation vector - real(8), allocatable :: rotation(:) - real(8), allocatable :: rotation_matrix(:,:) - contains procedure :: id => cell_id diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 574d2949b..60c6485ac 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -217,41 +217,6 @@ contains // to_str(c % id())) end if - ! Rotation matrix - if (check_for_node(node_cell, "rotation")) then - ! Rotations can only be applied to cells that are being filled with - ! another universe - if (c % fill() == C_NONE) then - call fatal_error("Cannot apply a rotation to cell " // trim(to_str(& - &c % id())) // " because it is not filled with another universe") - end if - - ! Read number of rotation parameters - n = node_word_count(node_cell, "rotation") - if (n /= 3) then - call fatal_error("Incorrect number of rotation parameters on cell " & - // to_str(c % id())) - end if - - ! Copy rotation angles in x,y,z directions - allocate(c % rotation(3)) - call get_node_array(node_cell, "rotation", c % rotation) - phi = -c % rotation(1) * PI/180.0_8 - theta = -c % rotation(2) * PI/180.0_8 - psi = -c % rotation(3) * PI/180.0_8 - - ! Calculate rotation matrix based on angles given - allocate(c % rotation_matrix(3,3)) - c % rotation_matrix = reshape((/ & - cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), & - -cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), & - cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), & - sin(phi)*cos(theta), & - sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), & - -sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), & - cos(phi)*cos(theta) /), (/ 3,3 /)) - end if - ! Add cell to dictionary call cell_dict % set(c % id(), i) diff --git a/src/particle.cpp b/src/particle.cpp index 26c4e8a7f..73ec98ce6 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -357,9 +357,11 @@ Particle::transport() for (int j = 0; j < n_coord - 1; ++j) { if (coord[j + 1].rotated) { // If next level is rotated, apply rotation matrix - // FIXME - // coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % & - // rotation_matrix, coord[j].uvw) + const auto& m {model::cells[coord[j].cell]->rotation_}; + Direction u {coord[j].uvw}; + coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z; + coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z; + coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z; } else { // Otherwise, copy this level's direction std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw); From 964fdfde3001b7040eedfaddf1a3667dba6ca76e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Feb 2019 06:35:59 -0600 Subject: [PATCH 08/10] Remove cells array on Fortran side --- CMakeLists.txt | 1 - include/openmc/cell.h | 2 + include/openmc/dagmc.h | 1 + include/openmc/geometry_aux.h | 4 +- include/openmc/lattice.h | 6 + include/openmc/surface.h | 6 + openmc/capi/cell.py | 5 +- src/api.F90 | 6 +- src/cell.cpp | 127 ++++++++-------- src/dagmc.cpp | 19 ++- src/geometry.F90 | 1 - src/geometry_aux.cpp | 47 +++++- src/geometry_header.F90 | 279 ---------------------------------- src/initialize.cpp | 1 - src/input_xml.F90 | 226 --------------------------- src/lattice.cpp | 7 +- src/source.cpp | 2 +- src/surface.cpp | 7 +- src/tallies/tally.F90 | 1 - tests/unit_tests/test_capi.py | 2 +- 20 files changed, 162 insertions(+), 588 deletions(-) delete mode 100644 src/geometry_header.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index 43e250c34..65292367d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -310,7 +310,6 @@ add_library(libopenmc SHARED src/dict_header.F90 src/error.F90 src/geometry.F90 - src/geometry_header.F90 src/hdf5_interface.F90 src/initialize.F90 src/input_xml.F90 diff --git a/include/openmc/cell.h b/include/openmc/cell.h index ba202c343..cb77d7e41 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -199,6 +199,8 @@ public: // Non-member functions //============================================================================== +void read_cells(pugi::xml_node node); + #ifdef DAGMC int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed); #endif diff --git a/include/openmc/dagmc.h b/include/openmc/dagmc.h index 4a129ef09..5a944a42d 100644 --- a/include/openmc/dagmc.h +++ b/include/openmc/dagmc.h @@ -26,6 +26,7 @@ extern moab::DagMC* DAG; extern "C" void load_dagmc_geometry(); extern "C" void free_memory_dagmc(); +void read_geometry_dagmc(); extern "C" pugi::xml_document* read_uwuw_materials(); bool get_uwuw_materials_xml(std::string& s); diff --git a/include/openmc/geometry_aux.h b/include/openmc/geometry_aux.h index 3e5cb93e1..a6bfe06e2 100644 --- a/include/openmc/geometry_aux.h +++ b/include/openmc/geometry_aux.h @@ -10,6 +10,8 @@ namespace openmc { +void read_geometry_xml(); + //============================================================================== //! Replace Universe, Lattice, and Material IDs with indices. //============================================================================== @@ -111,7 +113,7 @@ extern "C" int maximum_levels(int32_t univ); //! Deallocates global vectors and maps for cells, universes, and lattices. //============================================================================== -extern "C" void free_memory_geometry_c(); +extern "C" void free_memory_geometry(); } // namespace openmc #endif // OPENMC_GEOMETRY_AUX_H diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h index 5a5eff0b4..68728de5f 100644 --- a/include/openmc/lattice.h +++ b/include/openmc/lattice.h @@ -273,5 +273,11 @@ private: std::array pitch_; //!< Lattice tile width and height }; +//============================================================================== +// Non-member functions +//============================================================================== + +void read_lattices(pugi::xml_node node); + } // namespace openmc #endif // OPENMC_LATTICE_H diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 3df791ba3..ce9f78399 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -395,6 +395,12 @@ public: void to_hdf5_inner(hid_t group_id) const; }; +//============================================================================== +// Non-member functions +//============================================================================== + +void read_surfaces(pugi::xml_node node); + //============================================================================== // Fortran compatibility functions //============================================================================== diff --git a/openmc/capi/cell.py b/openmc/capi/cell.py index 1a4c8996f..894e4c79e 100644 --- a/openmc/capi/cell.py +++ b/openmc/capi/cell.py @@ -38,6 +38,7 @@ _dll.openmc_cell_set_temperature.errcheck = _error_handler _dll.openmc_get_cell_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_cell_index.restype = c_int _dll.openmc_get_cell_index.errcheck = _error_handler +_dll.cells_size.restype = c_int class Cell(_FortranObjectWithID): @@ -156,10 +157,10 @@ class _CellMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield Cell(index=i + 1).id + yield Cell(index=i).id def __len__(self): - return c_int32.in_dll(_dll, 'n_cells').value + return _dll.cells_size() def __repr__(self): return repr(dict(self)) diff --git a/src/api.F90 b/src/api.F90 index 5200927cf..cb3856861 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -39,7 +39,7 @@ contains err = E_UNASSIGNED if (found) then - index = p % coord(p % n_coord) % cell + 1 + index = p % coord(p % n_coord) % cell instance = p % cell_instance err = 0 else @@ -59,7 +59,6 @@ contains subroutine free_memory() bind(C) use bank_header - use geometry_header use material_header use nuclide_header use photon_header @@ -93,6 +92,9 @@ contains subroutine free_memory_surfaces() bind(C) end subroutine + subroutine free_memory_geometry() bind(C) + end subroutine + subroutine sab_clear() bind(C) end subroutine end interface diff --git a/src/cell.cpp b/src/cell.cpp index 56635309f..2a2bd82d6 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -23,8 +23,6 @@ namespace openmc { namespace model { -int32_t n_cells {0}; - std::vector cells; std::unordered_map cell_map; @@ -643,18 +641,18 @@ void DAGCell::to_hdf5(hid_t group_id) const { return; } // Non-method functions //============================================================================== -extern "C" void -read_cells(pugi::xml_node* node) +void read_cells(pugi::xml_node node) { // Count the number of cells. - for (pugi::xml_node cell_node: node->children("cell")) {model::n_cells++;} - if (model::n_cells == 0) { + int n_cells = 0; + for (pugi::xml_node cell_node: node.children("cell")) {n_cells++;} + if (n_cells == 0) { fatal_error("No cells found in geometry.xml!"); } // Loop over XML cell elements and populate the array. - model::cells.reserve(model::n_cells); - for (pugi::xml_node cell_node : node->children("cell")) { + model::cells.reserve(n_cells); + for (pugi::xml_node cell_node : node.children("cell")) { model::cells.push_back(new CSGCell(cell_node)); } @@ -699,9 +697,8 @@ read_cells(pugi::xml_node* node) extern "C" int openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n) { - if (index >= 1 && index <= model::cells.size()) { - //TODO: off-by-one - Cell& c {*model::cells[index - 1]}; + if (index >= 0 && index < model::cells.size()) { + Cell& c {*model::cells[index]}; *type = c.type_; if (c.type_ == FILL_MATERIAL) { *indices = c.material_.data(); @@ -721,9 +718,8 @@ extern "C" int openmc_cell_set_fill(int32_t index, int type, int32_t n, const int32_t* indices) { - if (index >= 1 && index <= model::cells.size()) { - //TODO: off-by-one - Cell& c {*model::cells[index - 1]}; + if (index >= 0 && index < model::cells.size()) { + Cell& c {*model::cells[index]}; if (type == FILL_MATERIAL) { c.type_ = FILL_MATERIAL; c.material_.clear(); @@ -756,9 +752,9 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n, extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { - if (index >= 1 && index <= model::cells.size()) { + if (index >= 0 && index < model::cells.size()) { //TODO: off-by-one - Cell& c {*model::cells[index - 1]}; + Cell& c {*model::cells[index]}; if (instance) { if (*instance >= 0 && *instance < c.sqrtkT_.size()) { @@ -781,28 +777,66 @@ openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) return 0; } +//! Return the index in the cells array of a cell with a given ID +extern "C" int +openmc_get_cell_index(int32_t id, int32_t* index) +{ + auto it = model::cell_map.find(id); + if (it != model::cell_map.end()) { + *index = it->second; + return 0; + } else { + set_errmsg("No cell exists with ID=" + std::to_string(id) + "."); + return OPENMC_E_INVALID_ID; + } +} + +//! Return the ID of a cell +extern "C" int +openmc_cell_get_id(int32_t index, int32_t* id) +{ + if (index >= 0 && index < model::cells.size()) { + *id = model::cells[index]->id_; + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Set the ID of a cell +extern "C" int +openmc_cell_set_id(int32_t index, int32_t id) +{ + if (index >= 0 && index < model::cells.size()) { + model::cells[index]->id_ = id; + model::cell_map[id] = index; + return 0; + } else { + set_errmsg("Index in cells array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } +} + +//! Extend the cells array by n elements +extern "C" int +openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end) +{ + if (index_start) *index_start = model::cells.size(); + if (index_end) *index_end = model::cells.size() + n - 1; + for (int32_t i = 0; i < n; i++) { + model::cells.push_back(new CSGCell()); + } + return 0; +} + + //============================================================================== // Fortran compatibility functions //============================================================================== extern "C" { - Cell* cell_pointer(int32_t cell_ind) {return model::cells[cell_ind];} - - int32_t cell_id(Cell* c) {return c->id_;} - - void - cell_set_id(Cell* c, int32_t id) - { - c->id_ = id; - - // Find the index of this cell and update the cell map. - for (int i = 0; i < model::cells.size(); i++) { - if (model::cells[i] == c) { - model::cell_map[id] = i; - break; - } - } - } + int cells_size() { return model::cells.size(); } #ifdef DAGMC int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed) @@ -818,33 +852,6 @@ extern "C" { return cur_cell->dagmc_ptr_->index_by_handle(new_vol); } #endif - - int32_t cell_universe(Cell* c) {return c->universe_;} - - int32_t cell_fill(Cell* c) {return c->fill_;} - - int cell_material_size(Cell* c) {return c->material_.size();} - - //TODO: off-by-one - int32_t cell_material(Cell* c, int i) - { - int32_t mat = c->material_[i-1]; - if (mat == MATERIAL_VOID) return MATERIAL_VOID; - return mat + 1; - } - - int cell_sqrtkT_size(Cell* c) {return c->sqrtkT_.size();} - - double cell_sqrtkT(Cell* c, int i) {return c->sqrtkT_[i];} - - void extend_cells_c(int32_t n) - { - model::cells.reserve(model::cells.size() + n); - for (int32_t i = 0; i < n; i++) { - model::cells.push_back(new CSGCell()); - } - model::n_cells = model::cells.size(); - } } diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 22f9d392a..3af23bc46 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -4,6 +4,7 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/geometry_aux.h" #include "openmc/string_utils.h" #include "openmc/settings.h" #include "openmc/geometry.h" @@ -133,9 +134,9 @@ void load_dagmc_geometry() /// Cells (Volumes) \\\ // initialize cell objects - model::n_cells = model::DAG->num_entities(3); + int n_cells = model::DAG->num_entities(3); moab::EntityHandle graveyard = 0; - for (int i = 0; i < model::n_cells; i++) { + for (int i = 0; i < n_cells; i++) { moab::EntityHandle vol_handle = model::DAG->entity_by_index(3, i+1); // set cell ids using global IDs @@ -312,6 +313,20 @@ void load_dagmc_geometry() return; } +void read_geometry_dagmc() +{ + // Check if dagmc.h5m exists + std::string filename = settings::path_input + "dagmc.h5m" + if (!file_exists(filename)) { + fatal_error("Geometry DAGMC file '" + filename + "' does not exist!"); + } + + write_message("Reading DAGMC geometry...", 5); + load_dagmc_geometry(); + + model::root_universe = find_root_universe() +} + void free_memory_dagmc() { delete model::DAG; diff --git a/src/geometry.F90 b/src/geometry.F90 index 0abe96ee4..1b7ead90a 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1,7 +1,6 @@ module geometry use constants - use geometry_header use particle_header use simulation_header use settings diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 9073c511d..f69a63d9b 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -4,10 +4,14 @@ #include #include +#include "pugixml.hpp" + #include "openmc/cell.h" #include "openmc/constants.h" #include "openmc/container_util.h" +#include "openmc/dagmc.h" #include "openmc/error.h" +#include "openmc/file_utils.h" #include "openmc/geometry.h" #include "openmc/lattice.h" #include "openmc/material.h" @@ -19,6 +23,46 @@ namespace openmc { +void read_geometry_xml() +{ +#ifdef DAGMC + if (settings::dagmc) { + read_geometry_dagmc(); + return; + } +#endif + + // Display output message + write_message("Reading geometry XML file...", 5); + + // Check if geometry.xml exists + std::string filename = settings::path_input + "geometry.xml"; + if (!file_exists(filename)) { + fatal_error("Geometry XML file '" + filename + "' does not exist!"); + } + + // Parse settings.xml file + pugi::xml_document doc; + auto result = doc.load_file(filename.c_str()); + if (!result) { + fatal_error("Error processing geometry.xml file."); + } + + // Get root element + pugi::xml_node root = doc.document_element(); + + // Read surfaces, cells, lattice + read_surfaces(root); + read_cells(root); + read_lattices(root); + + // ========================================================================== + // SETUP UNIVERSES + + // Allocate universes, universe cell arrays, and assign base universe + model::root_universe = find_root_universe(); +} + //============================================================================== void @@ -473,12 +517,11 @@ maximum_levels(int32_t univ) //============================================================================== void -free_memory_geometry_c() +free_memory_geometry() { for (Cell* c : model::cells) {delete c;} model::cells.clear(); model::cell_map.clear(); - model::n_cells = 0; for (Universe* u : model::universes) {delete u;} model::universes.clear(); diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 deleted file mode 100644 index afb9e2c0d..000000000 --- a/src/geometry_header.F90 +++ /dev/null @@ -1,279 +0,0 @@ -module geometry_header - - use, intrinsic :: ISO_C_BINDING - - use dict_header, only: DictIntInt - use error - use string, only: to_str - - implicit none - - interface - function cell_pointer(cell_ind) bind(C) result(ptr) - import C_PTR, C_INT32_T - integer(C_INT32_T), intent(in), value :: cell_ind - type(C_PTR) :: ptr - end function cell_pointer - - function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T) :: id - end function cell_id_c - - subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id') - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T), intent(in), value :: id - end subroutine cell_set_id_c - - function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & - result(universe) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T) :: universe - end function cell_universe_c - - function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT32_T) :: fill - end function cell_fill_c - - function cell_material_size_c(cell_ptr) bind(C, name='cell_material_size') & - result(n) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT) :: n - end function cell_material_size_c - - function cell_material_c(cell_ptr, i) bind(C, name='cell_material') & - result(mat) - import C_PTR, C_INT, C_INT32_T - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT), intent(in), value :: i - integer(C_INT32_T) :: mat - end function cell_material_c - - function cell_sqrtkT_size_c(cell_ptr) bind(C, name='cell_sqrtkT_size') & - result(n) - import C_PTR, C_INT - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT) :: n - end function cell_sqrtkT_size_c - - function cell_sqrtkT_c(cell_ptr, i) bind(C, name='cell_sqrtkT') & - result(sqrtkT) - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: cell_ptr - integer(C_INT), intent(in), value :: i - real(C_DOUBLE) :: sqrtkT - end function cell_sqrtkT_c - - subroutine extend_cells_c(n) bind(C) - import C_INT32_t - integer(C_INT32_T), intent(in), value :: n - end subroutine extend_cells_c - end interface - -!=============================================================================== -! CELL defines a closed volume by its bounding surfaces -!=============================================================================== - - type Cell - type(C_PTR) :: ptr - contains - - procedure :: id => cell_id - procedure :: set_id => cell_set_id - procedure :: universe => cell_universe - procedure :: fill => cell_fill - procedure :: material_size => cell_material_size - procedure :: material => cell_material - procedure :: sqrtkT_size => cell_sqrtkT_size - procedure :: sqrtkT => cell_sqrtkT - - end type Cell - - ! array index of the root universe - integer(C_INT), bind(C) :: root_universe - - integer(C_INT32_T), bind(C) :: n_cells ! # of cells - integer(C_INT32_T), bind(C) :: n_universes ! # of universes - - type(Cell), allocatable, target :: cells(:) - - ! Dictionaries which map user IDs to indices in the global arrays - type(DictIntInt) :: cell_dict - -contains - -!=============================================================================== - - function cell_id(this) result(id) - class(Cell), intent(in) :: this - integer(C_INT32_T) :: id - id = cell_id_c(this % ptr) - end function cell_id - - subroutine cell_set_id(this, id) - class(Cell), intent(in) :: this - integer(C_INT32_T), intent(in) :: id - call cell_set_id_c(this % ptr, id) - end subroutine cell_set_id - - function cell_universe(this) result(universe) - class(Cell), intent(in) :: this - integer(C_INT32_T) :: universe - universe = cell_universe_c(this % ptr) - end function cell_universe - - function cell_fill(this) result(fill) - class(Cell), intent(in) :: this - integer(C_INT32_T) :: fill - fill = cell_fill_c(this % ptr) - end function cell_fill - - function cell_material_size(this) result(n) - class(Cell), intent(in) :: this - integer(C_INT) :: n - n = cell_material_size_c(this % ptr) - end function cell_material_size - - function cell_material(this, i) result(mat) - class(Cell), intent(in) :: this - integer, intent(in) :: i - integer(C_INT32_T) :: mat - mat = cell_material_c(this % ptr, i) - end function cell_material - - function cell_sqrtkT_size(this) result(n) - class(Cell), intent(in) :: this - integer :: n - n = cell_sqrtkT_size_c(this % ptr) - end function cell_sqrtkT_size - - function cell_sqrtkT(this, i) result(sqrtkT) - class(Cell), intent(in) :: this - integer, intent(in) :: i - real(C_DOUBLE) :: sqrtkT - sqrtkT = cell_sqrtkT_c(this % ptr, i) - end function cell_sqrtkT - -!=============================================================================== -! FREE_MEMORY_GEOMETRY deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_geometry() - interface - subroutine free_memory_geometry_c() bind(C) - end subroutine free_memory_geometry_c - end interface - - call free_memory_geometry_c() - - n_cells = 0 - n_universes = 0 - - if (allocated(cells)) deallocate(cells) - - call cell_dict % clear() - - end subroutine free_memory_geometry - -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_extend_cells(n, index_start, index_end) result(err) bind(C) - ! Extend the cells array by n elements - integer(C_INT32_T), value, intent(in) :: n - integer(C_INT32_T), optional, intent(out) :: index_start - integer(C_INT32_T), optional, intent(out) :: index_end - integer(C_INT32_T) :: i - integer(C_INT) :: err - - type(Cell), allocatable :: temp(:) ! temporary cells array - - if (n_cells == 0) then - ! Allocate cells array - allocate(cells(n)) - else - ! Allocate cells array with increased size - allocate(temp(n_cells + n)) - - ! Copy original cells to temporary array - temp(1:n_cells) = cells - - ! Move allocation from temporary array - call move_alloc(FROM=temp, TO=cells) - end if - - ! Return indices in cells array - if (present(index_start)) index_start = n_cells + 1 - if (present(index_end)) index_end = n_cells + n - - ! Extend the C++ cells array and get pointers to the C++ objects - call extend_cells_c(n) - do i = n_cells - n, n_cells - cells(i) % ptr = cell_pointer(i - 1) - end do - - err = 0 - end function openmc_extend_cells - - - function openmc_get_cell_index(id, index) result(err) bind(C) - ! Return the index in the cells array of a cell with a given ID - integer(C_INT32_T), value :: id - integer(C_INT32_T), intent(out) :: index - integer(C_INT) :: err - - if (allocated(cells)) then - if (cell_dict % has(id)) then - index = cell_dict % get(id) - err = 0 - else - err = E_INVALID_ID - call set_errmsg("No cell exists with ID=" // trim(to_str(id)) // ".") - end if - else - err = E_ALLOCATE - call set_errmsg("Memory has not been allocated for cells.") - end if - end function openmc_get_cell_index - - - function openmc_cell_get_id(index, id) result(err) bind(C) - ! Return the ID of a cell - integer(C_INT32_T), value :: index - integer(C_INT32_T), intent(out) :: id - integer(C_INT) :: err - - if (index >= 1 .and. index <= size(cells)) then - id = cells(index) % id() - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in cells array is out of bounds.") - end if - end function openmc_cell_get_id - - - function openmc_cell_set_id(index, id) result(err) bind(C) - ! Set the ID of a cell - integer(C_INT32_T), value, intent(in) :: index - integer(C_INT32_T), value, intent(in) :: id - integer(C_INT) :: err - - if (index >= 1 .and. index <= n_cells) then - call cells(index) % set_id(id) - call cell_dict % set(id, index) - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in cells array is out of bounds.") - end if - end function openmc_cell_set_id - -end module geometry_header diff --git a/src/initialize.cpp b/src/initialize.cpp index e454917d0..59662d577 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -31,7 +31,6 @@ // data/functions from Fortran side extern "C" void read_command_line(); -extern "C" void read_geometry_xml(); extern "C" void read_materials_xml(); extern "C" void read_plots_xml(); extern "C" void read_tallies_xml(); diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 60c6485ac..8da108e35 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3,9 +3,7 @@ module input_xml use, intrinsic :: ISO_C_BINDING use constants - use dict_header, only: DictIntInt, DictCharInt, DictEntryCI use error, only: fatal_error, warning, write_message, openmc_err_msg - use geometry_header #ifdef DAGMC use dagmc_header #endif @@ -30,42 +28,11 @@ module input_xml save interface - subroutine count_cell_instances(univ_indx) bind(C) - import C_INT32_T - integer(C_INT32_T), intent(in), value :: univ_indx - end subroutine count_cell_instances - - subroutine read_surfaces(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine read_surfaces - - subroutine read_cells(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine read_cells - - subroutine read_lattices(node_ptr) bind(C) - import C_PTR - type(C_PTR) :: node_ptr - end subroutine read_lattices - subroutine read_materials(node_ptr) bind(C) import C_PTR type(C_PTR) :: node_ptr end subroutine read_materials - function find_root_universe() bind(C) result(root) - import C_INT32_T - integer(C_INT32_T) :: root - end function find_root_universe - - function maximum_levels(univ) bind(C) result(n) - import C_INT32_T, C_INT - integer(C_INT32_T), intent(in), value :: univ - integer(C_INT) :: n - end function maximum_levels - subroutine read_plots(node_ptr) bind(C) import C_PTR type(C_PTR) :: node_ptr @@ -81,199 +48,6 @@ module input_xml contains -#ifdef DAGMC - -!=============================================================================== -! READ_GEOMETRY_DAGMC reads data from a DAGMC .h5m file, checking -! for material properties and surface boundary conditions -! some universe information is spoofed for now -!=============================================================================== - - subroutine read_geometry_dagmc() - - integer :: i - integer :: univ_id - integer :: n_cells_in_univ - logical :: file_exists - character(MAX_LINE_LEN) :: filename - type(Cell), pointer :: c - type(VectorInt) :: univ_ids ! List of all universe IDs - type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each - ! universe contains - - ! Check if dagmc.h5m exists - filename = trim(path_input) // "dagmc.h5m" - inquire(FILE=filename, EXIST=file_exists) - if (.not. file_exists) then - call fatal_error("Geometry DAGMC file '" // trim(filename) // "' does not & - &exist!") - end if - - call write_message("Reading DAGMC geometry...", 5) - call load_dagmc_geometry() - call allocate_cells() - - ! setup universe data structs - do i = 1, n_cells - c => cells(i) - ! additional metadata spoofing - univ_id = c % universe() - - if (.not. cells_in_univ_dict % has(univ_id)) then - n_universes = n_universes + 1 - n_cells_in_univ = 1 - call univ_ids % push_back(univ_id) - else - n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id) - end if - call cells_in_univ_dict % set(univ_id, n_cells_in_univ) - end do - - root_universe = find_root_universe() - - end subroutine read_geometry_dagmc - -#endif - -!=============================================================================== -! READ_GEOMETRY_XML reads data from a geometry.xml file and parses it, checking -! for errors and placing properly-formatted data in the right data structures -!=============================================================================== - - subroutine read_geometry_xml() bind(C) - - integer :: i, n - integer :: univ_id - integer :: n_cells_in_univ - real(8) :: phi, theta, psi - logical :: file_exists - character(MAX_LINE_LEN) :: filename - type(Cell), pointer :: c - type(XMLDocument) :: doc - type(XMLNode) :: root - type(XMLNode) :: node_cell - type(XMLNode), allocatable :: node_cell_list(:) - type(VectorInt) :: univ_ids ! List of all universe IDs - type(DictIntInt) :: cells_in_univ_dict ! Used to count how many cells each - ! universe contains -#ifdef DAGMC - if (dagmc) then - call read_geometry_dagmc() - return - end if -#endif - - ! Display output message - call write_message("Reading geometry XML file...", 5) - - ! Check if geometry.xml exists - filename = trim(path_input) // "geometry.xml" - inquire(FILE=filename, EXIST=file_exists) - if (.not. file_exists) then - call fatal_error("Geometry XML file '" // trim(filename) // "' does not & - &exist!") - end if - - ! Parse geometry.xml file - call doc % load_file(filename) - root = doc % document_element() - - ! ========================================================================== - ! READ SURFACES FROM GEOMETRY.XML - - call read_surfaces(root % ptr) - - ! ========================================================================== - ! READ CELLS FROM GEOMETRY.XML - - call read_cells(root % ptr) - - ! Get pointer to list of XML - call get_node_list(root, "cell", node_cell_list) - - ! Get number of tags - n_cells = size(node_cell_list) - - ! Check for no cells - if (n_cells == 0) then - call fatal_error("No cells found in geometry.xml!") - end if - - ! Allocate cells array - allocate(cells(n_cells)) - - n_universes = 0 - do i = 1, n_cells - c => cells(i) - - c % ptr = cell_pointer(i - 1) - - ! Get pointer to i-th cell node - node_cell = node_cell_list(i) - - ! Check to make sure 'id' hasn't been used - if (cell_dict % has(c % id())) then - call fatal_error("Two or more cells use the same unique ID: " & - // to_str(c % id())) - end if - - ! Add cell to dictionary - call cell_dict % set(c % id(), i) - - ! For cells, we also need to check if there's a new universe -- - ! also for every cell add 1 to the count of cells for the - ! specified universe - univ_id = c % universe() - if (.not. cells_in_univ_dict % has(univ_id)) then - n_universes = n_universes + 1 - n_cells_in_univ = 1 - call univ_ids % push_back(univ_id) - else - n_cells_in_univ = 1 + cells_in_univ_dict % get(univ_id) - end if - call cells_in_univ_dict % set(univ_id, n_cells_in_univ) - - end do - - ! ========================================================================== - ! READ LATTICES FROM GEOMETRY.XML - - call read_lattices(root % ptr) - - ! ========================================================================== - ! SETUP UNIVERSES - - ! Allocate universes, universe cell arrays, and assign base universe - root_universe = find_root_universe() - - ! Clear dictionary - call cells_in_univ_dict%clear() - - ! Close geometry XML file - call doc % clear() - - end subroutine read_geometry_xml - - subroutine allocate_cells() - integer :: i - type(Cell), pointer :: c - - ! Allocate cells array - allocate(cells(n_cells)) - - do i = 1, n_cells - c => cells(i) - c % ptr = cell_pointer(i - 1) - ! Check to make sure 'id' hasn't been used - if (cell_dict % has(c % id())) then - call fatal_error("Two or more cells use the same unique ID: " & - // to_str(c % id())) - end if - ! Add cell to dictionary - call cell_dict % set(c % id(), i) - end do - end subroutine allocate_cells - subroutine read_materials_xml() bind(C) logical :: file_exists ! does materials.xml exist? character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml diff --git a/src/lattice.cpp b/src/lattice.cpp index d42440439..d8185de62 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -864,13 +864,12 @@ HexLattice::to_hdf5_inner(hid_t lat_group) const // Non-method functions //============================================================================== -extern "C" void -read_lattices(pugi::xml_node *node) +void read_lattices(pugi::xml_node node) { - for (pugi::xml_node lat_node : node->children("lattice")) { + for (pugi::xml_node lat_node : node.children("lattice")) { model::lattices.push_back(new RectLattice(lat_node)); } - for (pugi::xml_node lat_node : node->children("hex_lattice")) { + for (pugi::xml_node lat_node : node.children("hex_lattice")) { model::lattices.push_back(new HexLattice(lat_node)); } diff --git a/src/source.cpp b/src/source.cpp index 0e23bb3dd..231a95e21 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -172,7 +172,7 @@ Bank SourceDistribution::sample() const if (space_box) { if (space_box->only_fissionable()) { // Determine material - auto c = model::cells[cell_index - 1]; + auto c = model::cells[cell_index]; auto mat_index = c->material_.size() == 1 ? c->material_[0] : c->material_[instance]; diff --git a/src/surface.cpp b/src/surface.cpp index 0cceee9f4..cd1f4fc2b 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1061,12 +1061,11 @@ void SurfaceQuadric::to_hdf5_inner(hid_t group_id) const //============================================================================== -extern "C" void -read_surfaces(pugi::xml_node* node) +void read_surfaces(pugi::xml_node node) { // Count the number of surfaces. int n_surfaces = 0; - for (pugi::xml_node surf_node : node->children("surface")) {n_surfaces++;} + for (pugi::xml_node surf_node : node.children("surface")) {n_surfaces++;} if (n_surfaces == 0) { fatal_error("No surfaces found in geometry.xml!"); } @@ -1076,7 +1075,7 @@ read_surfaces(pugi::xml_node* node) { pugi::xml_node surf_node; int i_surf; - for (surf_node = node->child("surface"), i_surf = 0; surf_node; + for (surf_node = node.child("surface"), i_surf = 0; surf_node; surf_node = surf_node.next_sibling("surface"), i_surf++) { std::string surf_type = get_node_value(surf_node, "type", true, true); diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index e7cdcbc1d..823ef1acc 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -4,7 +4,6 @@ module tally use bank_header use constants - use geometry_header use material_header use message_passing use mgxs_interface diff --git a/tests/unit_tests/test_capi.py b/tests/unit_tests/test_capi.py index 89cf5d708..d5943b7be 100644 --- a/tests/unit_tests/test_capi.py +++ b/tests/unit_tests/test_capi.py @@ -72,7 +72,7 @@ def test_cell(capi_init): cell = openmc.capi.cells[1] assert isinstance(cell.fill, openmc.capi.Material) cell.fill = openmc.capi.materials[1] - assert str(cell) == 'Cell[1]' + assert str(cell) == 'Cell[0]' def test_new_cell(capi_init): From 39961dc0a41dae9791735a2ac24c09550e4435a0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 15 Feb 2019 11:42:49 -0600 Subject: [PATCH 09/10] Fix most DAGMC-related errors --- include/openmc/cell.h | 8 ++++---- src/cell.cpp | 34 +++++++++++++++------------------- src/dagmc.cpp | 4 ++-- src/particle.cpp | 18 +++++++++++------- 4 files changed, 32 insertions(+), 32 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index cb77d7e41..e2c3f400a 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -9,14 +9,14 @@ #include "hdf5.h" #include "pugixml.hpp" +#ifdef DAGMC +#include "DagMC.hpp" +#endif #include "openmc/constants.h" #include "openmc/neighbor_list.h" #include "openmc/position.h" - -#ifdef DAGMC -#include "DagMC.hpp" -#endif +#include "openmc/surface.h" namespace openmc { diff --git a/src/cell.cpp b/src/cell.cpp index 2a2bd82d6..af77094c9 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -830,29 +830,25 @@ openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end) return 0; } +#ifdef DAGMC +int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed) +{ + moab::EntityHandle surf = + surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_); + moab::EntityHandle vol = + cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_); + + moab::EntityHandle new_vol; + cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol); + + return cur_cell->dagmc_ptr_->index_by_handle(new_vol); +} +#endif //============================================================================== // Fortran compatibility functions //============================================================================== -extern "C" { - int cells_size() { return model::cells.size(); } - - #ifdef DAGMC - int32_t next_cell(DAGCell* cur_cell, DAGSurface* surf_xed) - { - moab::EntityHandle surf = - surf_xed->dagmc_ptr_->entity_by_id(2, surf_xed->id_); - moab::EntityHandle vol = - cur_cell->dagmc_ptr_->entity_by_id(3, cur_cell->id_); - - moab::EntityHandle new_vol; - cur_cell->dagmc_ptr_->next_vol(surf, vol, new_vol); - - return cur_cell->dagmc_ptr_->index_by_handle(new_vol); - } - #endif -} - +extern "C" int cells_size() { return model::cells.size(); } } // namespace openmc diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 3af23bc46..a398c9f75 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -316,7 +316,7 @@ void load_dagmc_geometry() void read_geometry_dagmc() { // Check if dagmc.h5m exists - std::string filename = settings::path_input + "dagmc.h5m" + std::string filename = settings::path_input + "dagmc.h5m"; if (!file_exists(filename)) { fatal_error("Geometry DAGMC file '" + filename + "' does not exist!"); } @@ -324,7 +324,7 @@ void read_geometry_dagmc() write_message("Reading DAGMC geometry...", 5); load_dagmc_geometry(); - model::root_universe = find_root_universe() + model::root_universe = find_root_universe(); } void free_memory_dagmc() diff --git a/src/particle.cpp b/src/particle.cpp index 73ec98ce6..4c39087ce 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -581,17 +581,21 @@ Particle::cross_surface() #ifdef DAGMC if (settings::dagmc) { - int32_t i_cell = next_cell(model::cells[last_cell[0]], - model::surfaces[std::abs(surface)]); + auto cellp = dynamic_cast(model::cells[last_cell[0]]); + // TODO: off-by-one + auto surfp = dynamic_cast(model::surfaces[std::abs(surface) - 1]); + int32_t i_cell = next_cell(cellp, surfp) - 1; // save material and temp last_material = material; - last_sqrtkT = sqrtKT; + last_sqrtkT = sqrtkT; // set new cell value - coord[0].cell = i_cell + coord[0].cell = i_cell; cell_instance = 0; - material = model::cells[i_cell]->material[0]; - sqrtKT = model::cells[i_cell]->sqrtKT[0]; - return + // TODO: off-by-one + int mat = model::cells[i_cell]->material_[0]; + material = (mat == MATERIAL_VOID) ? mat : mat + 1; + sqrtkT = model::cells[i_cell]->sqrtkT_[0]; + return; } #endif From 7325897ee9f4a0b80d215ae47ef51fc4005ea623 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 19 Feb 2019 21:48:20 -0600 Subject: [PATCH 10/10] Address @smharper comments on #1169 --- include/openmc/hdf5_interface.h | 4 +++ include/openmc/particle.h | 4 +-- include/openmc/simulation.h | 1 + include/openmc/tallies/derivative.h | 12 ++++++++ include/openmc/tallies/tally.h | 2 -- include/openmc/tallies/tally_scoring.h | 31 +++++++++++++++++++ src/cell.cpp | 1 - src/geometry_aux.cpp | 3 -- src/particle.cpp | 41 ++++++++++++-------------- src/simulation.cpp | 2 +- src/tallies/derivative.cpp | 9 ------ src/tallies/tally.cpp | 4 +-- src/tallies/tally_scoring.cpp | 24 --------------- 13 files changed, 72 insertions(+), 66 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index f55b14171..dcf7ab8dd 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -223,6 +223,10 @@ read_attribute(hid_t obj_id, const char* name, std::vector& vec) // Templates/overloads for read_dataset and related methods //============================================================================== +// Template for scalars. We need to be careful that the compiler does not use +// this version of read_dataset for vectors, arrays, or other non-scalar types. +// enable_if_t allows us to conditionally remove the function from overload +// resolution when the type T doesn't meet a certain criterion. template inline std::enable_if_t>::value> read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index a81b3e77a..01c7eaefc 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -148,10 +148,10 @@ extern "C" { //! \param src Source site data void from_source(const Bank* src); - //! Transport the particle + //! Transport a particle from birth to death void transport(); - //! Cross a surface + //! Cross a surface and handle boundary conditions void cross_surface(); //! mark a particle as lost and create a particle restart file diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index af333f714..c7053c608 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -36,6 +36,7 @@ extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? extern "C" int total_gen; //!< total number of generations simulated +extern "C" double total_weight; //!< Total source weight in a batch extern "C" int64_t work; //!< number of particles per process extern std::vector k_generation; diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h index 21df00624..02532afc0 100644 --- a/include/openmc/tallies/derivative.h +++ b/include/openmc/tallies/derivative.h @@ -35,8 +35,20 @@ void apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, double atom_density, int score_bin, double& score); +//! Adjust diff tally flux derivatives for a particle scattering event. +// +//! Note that this subroutine will be called after absorption events in +//! addition to scattering events, but any flux derivatives scored after an +//! absorption will never be tallied. The paricle will be killed before any +//! further tallies are scored. +// +//! \param p The particle being tracked void score_collision_derivative(const Particle* p); +//! Adjust diff tally flux derivatives for a particle tracking event. +// +//! \param p The particle being tracked +//! \param distance The distance in [cm] traveled by the particle void score_track_derivative(const Particle* p, double distance); //! Set the flux derivatives on differential tallies to zero. diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 464ec9f3d..7767110c2 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -101,8 +101,6 @@ private: // Global variable declarations //============================================================================== -extern "C" double total_weight; - namespace model { extern std::vector> tallies; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 12016956c..af03b1aa2 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -50,14 +50,45 @@ private: // Non-member functions //============================================================================== +//! Score tallies using a 1 / Sigma_t estimate of the flux. +// +//! This is triggered after every collision. It is invalid for tallies that +//! require post-collison information because it can score reactions that didn't +//! actually occur, and we don't a priori know what the outcome will be for +//! reactions that we didn't sample. It is assumed the material is not void +//! since collisions do not occur in voids. +// +//! \param p The particle being tracked void score_collision_tally(const Particle* p); +//! Score tallies based on a simple count of events (for continuous energy). +// +//! Analog tallies are triggered at every collision, not every event. +// +//! \param p The particle being tracked void score_analog_tally_ce(const Particle* p); +//! Score tallies based on a simple count of events (for multigroup). +// +//! Analog tallies are triggered at every collision, not every event. +// +//! \param p The particle being tracked void score_analog_tally_mg(const Particle* p); +//! Score tallies using a tracklength estimate of the flux. +// +//! This is triggered at every event (surface crossing, lattice crossing, or +//! collision) and thus cannot be done for tallies that require post-collision +//! information. +// +//! \param p The particle being tracked +//! \param distance The distance in [cm] traveled by the particle void score_tracklength_tally(const Particle* p, double distance); +//! Score surface or mesh-surface tallies for particle currents. +// +//! \param p The particle being tracked +//! \param tallies A vector of tallies to score to void score_surface_tally(const Particle* p, const std::vector& tallies); } // namespace openmc diff --git a/src/cell.cpp b/src/cell.cpp index af77094c9..6a962eced 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -753,7 +753,6 @@ extern "C" int openmc_cell_set_temperature(int32_t index, double T, const int32_t* instance) { if (index >= 0 && index < model::cells.size()) { - //TODO: off-by-one Cell& c {*model::cells[index]}; if (instance) { diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index f69a63d9b..d38b5f8f5 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -56,9 +56,6 @@ void read_geometry_xml() read_cells(root); read_lattices(root); - // ========================================================================== - // SETUP UNIVERSES - // Allocate universes, universe cell arrays, and assign base universe model::root_universe = find_root_universe(); } diff --git a/src/particle.cpp b/src/particle.cpp index 4c39087ce..db2421d78 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -149,8 +149,8 @@ Particle::transport() int n_event = 0; // Add paricle's starting weight to count for normalizing tallies later -#pragma omp atomic - total_weight += wgt; + #pragma omp atomic + simulation::total_weight += wgt; // Force calculation of cross-sections by setting last energy to zero if (settings::run_CE) { @@ -375,8 +375,8 @@ Particle::transport() // If particle has too many events, display warning and kill it ++n_event; if (n_event == MAX_EVENTS) { - if (mpi::master) warning("Particle " + std::to_string(id) - + " underwent maximum number of events."); + warning("Particle " + std::to_string(id) + + " underwent maximum number of events."); alive = false; } @@ -533,26 +533,23 @@ Particle::cross_surface() std::copy(&r.x, &r.x + 3, coord[0].xyz); } - // Get a pointer to the partner periodic surface. Offset the index to - // correct for C vs. Fortran indexing. + // Get a pointer to the partner periodic surface auto surf_p = dynamic_cast(surf); - if (surf_p) { - auto other = dynamic_cast( - model::surfaces[surf_p->i_periodic_]); + auto other = dynamic_cast( + model::surfaces[surf_p->i_periodic_]); - // Adjust the particle's location and direction. - Position r {coord[0].xyz}; - Direction u {coord[0].uvw}; - bool rotational = other->periodic_translate(surf_p, r, u); - std::copy(&r.x, &r.x + 3, coord[0].xyz); - std::copy(&u.x, &u.x + 3, coord[0].uvw); + // Adjust the particle's location and direction. + Position r {coord[0].xyz}; + Direction u {coord[0].uvw}; + bool rotational = other->periodic_translate(surf_p, r, u); + std::copy(&r.x, &r.x + 3, coord[0].xyz); + std::copy(&u.x, &u.x + 3, coord[0].uvw); - // Reassign particle's surface - // TODO: off-by-one - surface = rotational ? - surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface); - } + // Reassign particle's surface + // TODO: off-by-one + surface = rotational ? + surf_p->i_periodic_ + 1 : + std::copysign(surf_p->i_periodic_ + 1, surface); // Figure out what cell particle is in now n_coord = 1; @@ -620,7 +617,7 @@ Particle::cross_surface() coord[0].xyz[1] += TINY_BIT * coord[0].uvw[1]; coord[0].xyz[2] += TINY_BIT * coord[0].uvw[2]; - // Couldn't find next cell anywhere// This probably means there is an actual + // Couldn't find next cell anywhere! This probably means there is an actual // undefined region in the geometry. if (!find_cell(this, false)) { diff --git a/src/simulation.cpp b/src/simulation.cpp index 1fdda3061..096ea87ff 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -312,7 +312,7 @@ void initialize_batch() } // Reset total starting particle weight used for normalizing tallies - total_weight = 0.0; + simulation::total_weight = 0.0; // Determine if this batch is the first inactive or active batch. bool first_inactive = false; diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp index c3bf8f1bf..b28b78d0a 100644 --- a/src/tallies/derivative.cpp +++ b/src/tallies/derivative.cpp @@ -569,8 +569,6 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide, } } -//! Adjust diff tally flux derivatives for a particle tracking event. - void score_track_derivative(const Particle* p, double distance) { @@ -620,13 +618,6 @@ score_track_derivative(const Particle* p, double distance) } } -//! Adjust diff tally flux derivatives for a particle scattering event. -// -//! Note that this subroutine will be called after absorption events in -//! addition to scattering events, but any flux derivatives scored after an -//! absorption will never be tallied. The paricle will be killed before any -//! further tallies are scored. - void score_collision_derivative(const Particle* p) { // A void material cannot be perturbed so it will not affect flux derivatives. diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 4bfb15eff..93ba72a19 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -585,9 +585,9 @@ void reduce_tally_results() // We also need to determine the total starting weight of particles from the // last realization double weight_reduced; - MPI_Reduce(&total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM, + MPI_Reduce(&simulation::total_weight, &weight_reduced, 1, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); - if (mpi::master) total_weight = weight_reduced; + if (mpi::master) simulation::total_weight = weight_reduced; } #endif diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index fd96d9d85..690f61f81 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -1972,10 +1972,6 @@ score_all_nuclides(const Particle* p, int i_tally, double flux, } } -//! Score tallies based on a simple count of events (for continuous energy). -// -//! Analog tallies ar etriggered at every collision, not every event. - void score_analog_tally_ce(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { @@ -2035,10 +2031,6 @@ void score_analog_tally_ce(const Particle* p) match.bins_present_ = false; } -//! Score tallies based on a simple count of events (for multigroup). -// -//! Analog tallies ar etriggered at every collision, not every event. - void score_analog_tally_mg(const Particle* p) { for (auto i_tally : model::active_analog_tallies) { @@ -2088,12 +2080,6 @@ void score_analog_tally_mg(const Particle* p) match.bins_present_ = false; } -//! Score tallies using a tracklength estimate of the flux. -// -//! This is triggered at every event (surface crossing, lattice crossing, or -//! collision) and thus cannot be done for tallies that require post-collision -//! information. - void score_tracklength_tally(const Particle* p, double distance) { @@ -2162,14 +2148,6 @@ score_tracklength_tally(const Particle* p, double distance) match.bins_present_ = false; } -//! Score tallies using a 1 / Sigma_t estimate of the flux. -// -//! This is triggered after every collision. It is invalid for tallies that -//! require post-collison information because it can score reactions that didn't -//! actually occur, and we don't a priori know what the outcome will be for -//! reactions that we didn't sample. It is assumed the material is not void -//! since collisions do not occur in voids. - void score_collision_tally(const Particle* p) { // Determine the collision estimate of the flux @@ -2238,8 +2216,6 @@ void score_collision_tally(const Particle* p) match.bins_present_ = false; } -//! Score surface or mesh-surface tallies for particle currents. - void score_surface_tally(const Particle* p, const std::vector& tallies) {