From 9e0e8685776333b4d5b289f48aa994a7420d4f1c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 28 Aug 2018 06:59:39 -0500 Subject: [PATCH 01/21] Add C++ version of RegularMesh class --- CMakeLists.txt | 1 + include/openmc/hdf5_interface.h | 9 + include/openmc/mesh.h | 58 +++++ include/openmc/xml_interface.h | 12 ++ src/mesh.cpp | 368 ++++++++++++++++++++++++++++++++ 5 files changed, 448 insertions(+) create mode 100644 include/openmc/mesh.h create mode 100644 src/mesh.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 57d2c6dacb..78d1caf4d8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -394,6 +394,7 @@ add_library(libopenmc SHARED src/lattice.cpp src/material.cpp src/math_functions.cpp + src/mesh.cpp src/message_passing.cpp src/mgxs.cpp src/mgxs_interface.cpp diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 48ce41d0e9..944fe75cb5 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -347,6 +347,15 @@ write_dataset(hid_t obj_id, const char* name, const std::vector& buffer) write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } +template inline void +write_dataset(hid_t obj_id, const char* name, const xt::xarray& arr) +{ + auto s = arr.shape(); + std::vector dims {s.cbegin(), s.cend()}; + write_dataset(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, + arr.data(), false); +} + inline void write_dataset(hid_t obj_id, const char* name, Position r) { diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h new file mode 100644 index 0000000000..151701e0b0 --- /dev/null +++ b/include/openmc/mesh.h @@ -0,0 +1,58 @@ +#ifndef OPENMC_MESH_H +#define OPENMC_MESH_H + +#include +#include +#include + +#include "hdf5.h" +#include "pugixml.hpp" +#include "xtensor/xarray.hpp" + +#include "openmc/position.h" + +namespace openmc { + +//============================================================================== +//! Tessellation of n-dimensional Euclidean space by congruent squares or cubes +//============================================================================== + +class RegularMesh { +public: + // Constructors + RegularMesh(pugi::xml_node node); + + // Methods + int get_bin(Position r); + int get_bin_from_indices(const int* ijk); + void get_indices(Position r, int* ijk, bool* in_mesh); + void get_indices_from_bin(int bin, int* ijk); + bool intersects(Position r0, Position r1); + void to_hdf5(hid_t group); + + int id_ {-1}; //!< User-specified ID + int n_dimension_; + double volume_frac_; + xt::xarray shape_; + xt::xarray lower_left_; + xt::xarray upper_right_; + xt::xarray width_; + +private: + bool intersects_1d(Position r0, Position r1); + bool intersects_2d(Position r0, Position r1); + bool intersects_3d(Position r0, Position r1); +}; + +//============================================================================== +// Global variables +//============================================================================== + +extern std::vector global_meshes; + +extern std::unordered_map mesh_map; + + +} // namespace openmc + +#endif // OPENMC_MESH_H diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h index fb87ffdae1..85db26faa3 100644 --- a/include/openmc/xml_interface.h +++ b/include/openmc/xml_interface.h @@ -1,11 +1,14 @@ #ifndef OPENMC_XML_INTERFACE_H #define OPENMC_XML_INTERFACE_H +#include // for size_t #include // for stringstream #include #include #include "pugixml.hpp" +#include "xtensor/xarray.hpp" +#include "xtensor/xadapt.hpp" namespace openmc { @@ -38,5 +41,14 @@ std::vector get_node_array(pugi::xml_node node, const char* name, return values; } +template +xt::xarray get_node_xarray(pugi::xml_node node, const char* name, + bool lowercase=false) +{ + std::vector v = get_node_array(node, name, lowercase); + std::vector shape = {v.size()}; + return xt::adapt(v, shape); +} + } // namespace openmc #endif // OPENMC_XML_INTERFACE_H diff --git a/src/mesh.cpp b/src/mesh.cpp new file mode 100644 index 0000000000..9e4aae0969 --- /dev/null +++ b/src/mesh.cpp @@ -0,0 +1,368 @@ +#include "openmc/mesh.h" + +#include // for ceil +#include + +#include "xtensor/xeval.hpp" +#include "xtensor/xmath.hpp" + +#include "openmc/error.h" +#include "openmc/hdf5_interface.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector global_meshes; + +std::unordered_map mesh_map; + +//============================================================================== +// RegularMesh implementation +//============================================================================== + +RegularMesh::RegularMesh(pugi::xml_node node) +{ + // Copy mesh id + if (check_for_node(node, "id")) { + id_ = std::stoi(get_node_value(node, "id")); + + // Check to make sure 'id' hasn't been used + if (mesh_map.find(id_) != mesh_map.end()) { + fatal_error("Two or more meshes use the same unique ID: " + + std::to_string(id_)); + } + } + + // Read mesh type + if (check_for_node(node, "type")) { + auto temp = get_node_value(node, "type", true, true); + if (temp == "regular") { + // TODO: move elsewhere + } else { + fatal_error("Invalid mesh type: " + temp); + } + } + + // Determine number of dimensions for mesh + if (check_for_node(node, "dimension")) { + shape_ = get_node_xarray(node, "dimension"); + int n = n_dimension_ = shape_.size(); + if (n != 1 && n != 2 && n != 3) { + fatal_error("Mesh must be one, two, or three dimensions."); + } + + // Check that dimensions are all greater than zero + if (xt::any(shape_ <= 0)) { + fatal_error("All entries on the element for a tally " + "mesh must be positive."); + } + } + + // Check for lower-left coordinates + if (check_for_node(node, "lower_left")) { + // Read mesh lower-left corner location + lower_left_ = get_node_xarray(node, "lower_left"); + } else { + fatal_error("Must specify on a mesh."); + } + + if (check_for_node(node, "width")) { + // Make sure both upper-right or width were specified + if (check_for_node(node, "upper_right")) { + fatal_error("Cannot specify both and on a mesh."); + } + + width_ = get_node_xarray(node, "width"); + + // Check to ensure width has same dimensions + auto n = width_.size(); + if (n != lower_left_.size()) { + fatal_error("Number of entries on must be the same as " + "the number of entries on ."); + } + + // Check for negative widths + if (xt::any(width_ < 0.0)) { + fatal_error("Cannot have a negative on a tally mesh."); + } + + // Set width and upper right coordinate + upper_right_ = xt::eval(lower_left_ + shape_ * width_); + + } else if (check_for_node(node, "upper_right")) { + upper_right_ = get_node_xarray(node, "upper_right"); + + // Check to ensure width has same dimensions + auto n = upper_right_.size(); + if (n != lower_left_.size()) { + fatal_error("Number of entries on must be the " + "same as the number of entries on ."); + } + + // Check that upper-right is above lower-left + if (xt::any(upper_right_ < lower_left_)) { + fatal_error("The coordinates must be greater than " + "the coordinates on a tally mesh."); + } + + // Set width and upper right coordinate + width_ = xt::eval((upper_right_ - lower_left_) / shape_); + } else { + fatal_error("Must specify either and on a mesh."); + } + + if (shape_.dimension() > 0) { + if (shape_.size() != lower_left_.size()) { + fatal_error("Number of entries on must be the same " + "as the number of entries on ."); + } + + // Set volume fraction + volume_frac_ = 1.0/xt::prod(shape_)(); + } +} + +int RegularMesh::get_bin(Position r) +{ + // Loop over the dimensions of the mesh + for (int i = 0; i < n_dimension_; ++i) { + // Check for cases where particle is outside of mesh + if (r[i] < lower_left_[i]) { + return -1; + } else if (r[i] > upper_right_[i]) { + return -1; + } + } + + // Determine indices + int ijk[n_dimension_]; + bool in_mesh; + get_indices(r, ijk, &in_mesh); + if (!in_mesh) return -1; + + // Convert indices to bin + return get_bin_from_indices(ijk); +} + +int RegularMesh::get_bin_from_indices(const int* ijk) +{ + if (n_dimension_ == 1) { + return ijk[0]; + } else if (n_dimension_ == 2) { + return (ijk[1] - 1)*shape_[0] + ijk[0]; + } else if (n_dimension_ == 3) { + return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0]; + } +} + +void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) +{ + // Find particle in mesh + *in_mesh = true; + for (int i = 0; i < n_dimension_; ++i) { + ijk[i] = std::ceil((r[i] - lower_left_[i]) / width_[i]); + + // Check if indices are within bounds + if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false; + } +} + +void RegularMesh::get_indices_from_bin(int bin, int* ijk) +{ + if (n_dimension_ == 1) { + ijk[0] = bin; + } else if (n_dimension_ == 2) { + ijk[0] = (bin - 1) % shape_[0] + 1; + ijk[1] = (bin - 1) / shape_[0] + 1; + } else if (n_dimension_ == 3) { + ijk[0] = (bin - 1) % shape_[0] + 1; + ijk[1] = (bin - 1) % (shape_[0] * shape_[1]) / shape_[0] + 1; + ijk[2] = (bin - 1) / (shape_[0] * shape_[1]) + 1; + } +} + +bool RegularMesh::intersects(Position r0, Position r1) +{ + switch(n_dimension_) { + case 1: + return intersects_1d(r0, r1); + case 2: + return intersects_2d(r0, r1); + case 3: + return intersects_3d(r0, r1); + } +} + +bool RegularMesh::intersects_1d(Position r0, Position r1) +{ + // Copy coordinates of mesh lower_left and upper_right + double left = lower_left_[0]; + double right = upper_right_[0]; + + // Check if line intersects either left or right surface + if (r0.x < left) { + return r1.x > left; + } else if (r0.x < right) { + return r1.x < left || r1.x > right; + } else { + return r1.x < right; + } +} + +bool RegularMesh::intersects_2d(Position r0, Position r1) +{ + // Copy coordinates of starting point + double x0 = r0.x; + double y0 = r0.y; + + // Copy coordinates of ending point + double x1 = r1.x; + double y1 = r1.y; + + // Copy coordinates of mesh lower_left + double xm0 = lower_left_[0]; + double ym0 = lower_left_[1]; + + // Copy coordinates of mesh upper_right + double xm1 = upper_right_[0]; + double ym1 = upper_right_[1]; + + // Check if line intersects left surface -- calculate the intersection point y + if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { + double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); + if (yi >= ym0 && yi < ym1) { + return true; + } + } + + // Check if line intersects back surface -- calculate the intersection point + // x + if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) { + double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0); + if (xi >= xm0 && xi < xm1) { + return true; + } + } + + // Check if line intersects right surface -- calculate the intersection + // point y + if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { + double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); + if (yi >= ym0 && yi < ym1) { + return true; + } + } + + // Check if line intersects front surface -- calculate the intersection point + // x + if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) { + double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0); + if (xi >= xm0 && xi < xm1) { + return true; + } + } + return false; +} + +bool RegularMesh::intersects_3d(Position r0, Position r1) +{ + // Copy coordinates of starting point + double x0 = r0.x; + double y0 = r0.y; + double z0 = r0.z; + + // Copy coordinates of ending point + double x1 = r1.x; + double y1 = r1.y; + double z1 = r1.z; + + // Copy coordinates of mesh lower_left + double xm0 = lower_left_[0]; + double ym0 = lower_left_[1]; + double zm0 = lower_left_[2]; + + // Copy coordinates of mesh upper_right + double xm1 = upper_right_[0]; + double ym1 = upper_right_[1]; + double zm1 = upper_right_[2]; + + // Check if line intersects left surface -- calculate the intersection point + // (y,z) + if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { + double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); + double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0); + if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { + return true; + } + } + + // Check if line intersects back surface -- calculate the intersection point + // (x,z) + if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) { + double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0); + double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0); + if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { + return true; + } + } + + // Check if line intersects bottom surface -- calculate the intersection + // point (x,y) + if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) { + double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0); + double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0); + if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { + return true; + } + } + + // Check if line intersects right surface -- calculate the intersection point + // (y,z) + if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { + double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); + double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0); + if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { + return true; + } + } + + // Check if line intersects front surface -- calculate the intersection point + // (x,z) + if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) { + double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0); + double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0); + if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { + return true; + } + } + + // Check if line intersects top surface -- calculate the intersection point + // (x,y) + if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) { + double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0); + double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0); + if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { + return true; + } + } + return false; +} + +void RegularMesh::to_hdf5(hid_t group) +{ + hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); + + write_dataset(mesh_group, "type", "regular"); + write_dataset(mesh_group, "dimension", shape_); + write_dataset(mesh_group, "lower_left", lower_left_); + write_dataset(mesh_group, "upper_right", upper_right_); + write_dataset(mesh_group, "width", width_); + + close_group(mesh_group); +} + +} // namespace openmc From 19de269ae65dd70bd9f900f018e55262947b01a2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Aug 2018 10:31:13 -0500 Subject: [PATCH 02/21] Populate meshes vector from settings.xml --- include/openmc/mesh.h | 5 ++- src/cmfd_execute.F90 | 1 - src/mesh.cpp | 2 +- src/output.F90 | 1 - src/settings.cpp | 98 ++++++++++++++++++++++++++++++++++++++++++- 5 files changed, 101 insertions(+), 6 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 151701e0b0..1d75d778d5 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -1,7 +1,8 @@ #ifndef OPENMC_MESH_H #define OPENMC_MESH_H -#include +#include // for size_t +#include // for unique_ptr #include #include @@ -48,7 +49,7 @@ private: // Global variables //============================================================================== -extern std::vector global_meshes; +extern std::vector> meshes; extern std::unordered_map mesh_map; diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 7fe153d273..74264a036e 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -214,7 +214,6 @@ contains use bank_header, only: source_bank use constants, only: ZERO, ONE use error, only: warning, fatal_error - use mesh_header, only: RegularMesh use mesh, only: count_bank_sites use message_passing use string, only: to_str diff --git a/src/mesh.cpp b/src/mesh.cpp index 9e4aae0969..fb9ce2d0dd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -16,7 +16,7 @@ namespace openmc { // Global variables //============================================================================== -std::vector global_meshes; +std::vector> meshes; std::unordered_map mesh_map; diff --git a/src/output.F90 b/src/output.F90 index b4ecf497ea..b7b418a39b 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -10,7 +10,6 @@ module output use error, only: fatal_error, warning use geometry_header use math, only: t_percentile - use mesh_header, only: RegularMesh, meshes use message_passing, only: master, n_procs use mgxs_interface use nuclide_header diff --git a/src/settings.cpp b/src/settings.cpp index b7310c37f8..7e727178a4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -1,7 +1,9 @@ #include "openmc/settings.h" +#include // for ceil, pow #include // for numeric_limits #include +#include // for out_of_range #include #include @@ -13,6 +15,7 @@ #include "openmc/distribution_spatial.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/mesh.h" #include "openmc/random_lcg.h" #include "openmc/source.h" #include "openmc/string_utils.h" @@ -484,7 +487,100 @@ read_settings_xml() //track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) } - // TODO: Read meshes + // Read meshes + for (auto node : root.children("mesh")) { + // Read mesh and add to vector + meshes.emplace_back(new RegularMesh{node}); + + // Map ID to position in vector + mesh_map[meshes.back()->id_] = meshes.size() - 1; + } + + // Shannon Entropy mesh + if (check_for_node(root, "entropy_mesh")) { + int temp = std::stoi(get_node_value(root, "entropy_mesh")); + try { + index_entropy_mesh = mesh_map.at(temp); + } catch (std::out_of_range) { + std::stringstream msg; + msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; + fatal_error(msg); + } + } else if (check_for_node(root, "entropy")) { + warning("Specifying a Shannon entropy mesh via the element " + "is deprecated. Please create a mesh using and then reference " + "it by specifying its ID in an element."); + + // Read entropy mesh from + auto node_entropy = root.child("entropy"); + meshes.emplace_back(new RegularMesh{node_entropy}); + + // Set entropy mesh index + index_entropy_mesh = meshes.size() - 1; + + // Assign ID and set mapping + meshes.back()->id_ = 10000; + mesh_map[10000] = index_entropy_mesh; + } + + if (index_entropy_mesh >= 0) { + auto& m = *meshes[index_entropy_mesh]; + if (m.shape_.dimension() == 0) { + // If the user did not specify how many mesh cells are to be used in + // each direction, we automatically determine an appropriate number of + // cells + int n = std::ceil(std::pow(settings::n_particles / 20, 1.0/3.0)); + m.shape_ = {n, n, n}; + m.n_dimension_ = 3; + + // Calculate width + m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_; + } + + // TODO: Allocate entropy_P + // Allocate space for storing number of fission sites in each mesh cell + //allocate(entropy_p(1, product(m % dimension))) + + // Turn on Shannon entropy calculation + settings::entropy_on = true; + } + + // Uniform fission source weighting mesh + if (check_for_node(root, "ufs_mesh")) { + auto temp = std::stoi(get_node_value(root, "ufs_mesh")); + try { + index_ufs_mesh = mesh_map.at(temp); + } catch (std::out_of_range) { + std::stringstream msg; + msg << "Mesh " << temp << " specified for uniform fission site method " + "does not exist."; + fatal_error(msg); + } + } else if (check_for_node(root, "uniform_fs")) { + warning("Specifying a UFS mesh via the element " + "is deprecated. Please create a mesh using and then reference " + "it by specifying its ID in a element."); + + // Read entropy mesh from + auto node_ufs = root.child("uniform_fs"); + meshes.emplace_back(new RegularMesh{node_ufs}); + + // Set entropy mesh index + index_ufs_mesh = meshes.size() - 1; + + // Assign ID and set mapping + meshes.back()->id_ = 10001; + mesh_map[10001] = index_entropy_mesh; + } + + if (index_ufs_mesh >= 0) { + // Allocate array to store source fraction for UFS + // TODO: Allocate source_frac + //allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension))) + + // Turn on uniform fission source weighting + settings::ufs_on = true; + } // TODO: Read From b7feb47971837ba1dd5a230d8bddcedc61975ae1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Aug 2018 11:30:21 -0500 Subject: [PATCH 03/21] Add norm method on Position --- include/openmc/position.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/include/openmc/position.h b/include/openmc/position.h index 97ab6f0e66..d7a13c01e7 100644 --- a/include/openmc/position.h +++ b/include/openmc/position.h @@ -1,6 +1,7 @@ #ifndef OPENMC_POSITION_H #define OPENMC_POSITION_H +#include #include namespace openmc { @@ -46,6 +47,9 @@ struct Position { inline double dot(Position other) { return x*other.x + y*other.y + z*other.z; } + inline double norm() { + return std::sqrt(x*x + y*y + z*z); + } // Data members double x = 0.; From fe76ad695f1d57d7caf03e6f2be915b281726da8 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Aug 2018 11:30:32 -0500 Subject: [PATCH 04/21] Add bins_crossed method for RegularMesh --- include/openmc/mesh.h | 3 + src/mesh.cpp | 154 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 157 insertions(+) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 1d75d778d5..7bd5a2c290 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -10,6 +10,7 @@ #include "pugixml.hpp" #include "xtensor/xarray.hpp" +#include "openmc/particle.h" #include "openmc/position.h" namespace openmc { @@ -24,6 +25,8 @@ public: RegularMesh(pugi::xml_node node); // Methods + void bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths); int get_bin(Position r); int get_bin_from_indices(const int* ijk); void get_indices(Position r, int* ijk, bool* in_mesh); diff --git a/src/mesh.cpp b/src/mesh.cpp index fb9ce2d0dd..58a4e776bd 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -5,7 +5,10 @@ #include "xtensor/xeval.hpp" #include "xtensor/xmath.hpp" +#include "xtensor/xsort.hpp" +#include "xtensor/xtensor.hpp" +#include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" #include "openmc/xml_interface.h" @@ -352,6 +355,157 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) return false; } +void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, + std::vector& lengths) +{ + constexpr int MAX_SEARCH_ITER = 100; + + // ======================================================================== + // Determine if the track intersects the tally mesh. + + // Copy the starting and ending coordinates of the particle. Offset these + // just a bit for the purposes of determining if there was an intersection + // in case the mesh surfaces coincide with lattice/geometric surfaces which + // might produce finite-precision errors. + Position last_r {p->last_xyz}; + Position r {p->coord[0].xyz}; + Direction u {p->coord[0].uvw}; + + Position r0 = last_r + TINY_BIT*u; + Position r1 = r - TINY_BIT*u; + + // Determine indices for starting and ending location. + int n = n_dimension_; + xt::xtensor ijk0 = xt::empty({n}); + bool start_in_mesh; + get_indices(r0, ijk0.data(), &start_in_mesh); + xt::xtensor ijk1 = xt::empty({n}); + bool end_in_mesh; + get_indices(r1, ijk1.data(), &end_in_mesh); + + // If this is the first iteration of the filter loop, check if the track + // intersects any part of the mesh. + if ((!start_in_mesh) && (!end_in_mesh)) { + if (!intersects(r0, r1)) return; + } + + // ======================================================================== + // Figure out which mesh cell to tally. + + // Copy the un-modified coordinates the particle direction. + r0 = Position{p->last_xyz}; + r1 = Position{p->coord[0].xyz}; + + // Compute the length of the entire track. + double total_distance = (r1 - r0).norm(); + + // We are looking for the first valid mesh bin. Check to see if the + // particle starts inside the mesh. + if (!start_in_mesh) { + xt::xtensor d = xt::zeros({n}); + + // The particle does not start in the mesh. Note that we nudged the + // start and end coordinates by a TINY_BIT each so we will have + // difficulty resolving tracks that are less than 2*TINY_BIT in length. + // If the track is that short, it is also insignificant so we can + // safely ignore it in the tallies. + if (total_distance < 2*TINY_BIT) return; + + // The particle does not start in the mesh so keep iterating the ijk0 + // indices to cross the nearest mesh surface until we've found a valid + // bin. MAX_SEARCH_ITER prevents an infinite loop. + int search_iter = 0; + int j; + while (xt::any(ijk0 < 1) || xt::any(ijk0 > shape_)) { + if (search_iter == MAX_SEARCH_ITER) { + warning("Failed to find a mesh intersection on a tally mesh filter."); + return; + } + + for (j = 0; j < n; ++j) { + if (std::abs(u[j]) < FP_PRECISION) { + d(j) = INFTY; + } else if (u[j] > 0) { + double xyz_cross = lower_left_[j] + ijk0(j) * width_[j]; + d(j) = (xyz_cross - r0[j]) / u[j]; + } else { + double xyz_cross = lower_left_[j] + (ijk0(j) - 1) * width_[j]; + d(j) = (xyz_cross - r0[j]) / u[j]; + } + } + + int j = xt::argmin(d)(0); + if (u[j] > 0.0) { + ++ijk0(j); + } else { + --ijk0(j); + } + + ++search_iter; + } + + // Advance position + r0 += d(j) * u; + } + + while (true) { + // ======================================================================== + // Compute the length of the track segment in the appropiate mesh cell and + // return. + + double distance; + int j; + if (ijk0 == ijk1) { + // The track ends in this cell. Use the particle end location rather + // than the mesh surface. + distance = (r1 - r0).norm(); + } else { + // The track exits this cell. Determine the distance to the closest mesh + // surface. + xt::xtensor d = xt::zeros({n}); + for (int j = 0; j < n; ++j) { + if (abs(u[j]) < FP_PRECISION) { + d(j) = INFTY; + } else if (u[j] > 0) { + double xyz_cross = lower_left_[j] + ijk0(j) * width_[j]; + d(j) = (xyz_cross - r0[j]) / u[j]; + } else { + double xyz_cross = lower_left_[j] + (ijk0(j) - 1) * width_[j]; + d(j) = (xyz_cross - r0[j]) / u[j]; + } + } + j = xt::argmin(d)(0); + distance = d(j); + } + + // Assign the next tally bin and the score. + int bin = get_bin_from_indices(ijk0.data()); + bins.push_back(bin); + lengths.push_back(distance / total_distance); + + // Find the next mesh cell that the particle enters. + + // If the particle track ends in that bin, { we are done. + if (ijk0 == ijk1) break; + + // Translate the starting coordintes by the distance to that face. This + // should be the xyz that we computed the distance to in the last + // iteration of the filter loop. + r0 += distance * u; + + // Increment the indices into the next mesh cell. + if (u[j] > 0.0) { + ++ijk0(j); + } else { + --ijk0(j); + } + + // If the next indices are invalid, then the track has left the mesh and + // we are done. + if (xt::any(ijk0 < 1) || xt::any(ijk0 > shape_)) break; + } +} + void RegularMesh::to_hdf5(hid_t group) { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); From fb22413e8d2ce18130c2f51115e695a85877d90e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 29 Aug 2018 14:33:05 -0500 Subject: [PATCH 05/21] Move mesh C API functions to C++ --- include/openmc/capi.h | 2 +- include/openmc/mesh.h | 2 +- src/api.F90 | 2 +- src/mesh.cpp | 141 +++++++++++++++++++++ src/mesh_header.F90 | 278 ++++++++++-------------------------------- src/settings.cpp | 4 +- 6 files changed, 213 insertions(+), 216 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 7fd79fccfa..7ce42aa478 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -71,7 +71,7 @@ extern "C" { int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n); int openmc_mesh_set_id(int32_t index, int32_t id); int openmc_mesh_set_dimension(int32_t index, int n, const int* dims); - int openmc_mesh_set_params(int32_t index, const double* ll, const double* ur, const double* width, int n); + int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width); 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); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 7bd5a2c290..71649fa619 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -1,7 +1,6 @@ #ifndef OPENMC_MESH_H #define OPENMC_MESH_H -#include // for size_t #include // for unique_ptr #include #include @@ -22,6 +21,7 @@ namespace openmc { class RegularMesh { public: // Constructors + RegularMesh() = default; RegularMesh(pugi::xml_node node); // Methods diff --git a/src/api.F90 b/src/api.F90 index 75876080dc..2a6775a0c8 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -11,7 +11,7 @@ module openmc_api use hdf5_interface use material_header use math - use mesh_header + use mesh_header, only: free_memory_mesh use message_passing use nuclide_header use initialize, only: openmc_init_f diff --git a/src/mesh.cpp b/src/mesh.cpp index 58a4e776bd..8aecd1a22a 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,5 +1,6 @@ #include "openmc/mesh.h" +#include // for size_t #include // for ceil #include @@ -8,6 +9,7 @@ #include "xtensor/xsort.hpp" #include "xtensor/xtensor.hpp" +#include "openmc/capi.h" #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" @@ -519,4 +521,143 @@ void RegularMesh::to_hdf5(hid_t group) close_group(mesh_group); } +//============================================================================== +// C API functions +//============================================================================== + +//! Extend the meshes array by n elements +extern "C" int +openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end) +{ + if (!index_start) *index_start = meshes.size(); + for (int i = 0; i < n; ++i) { + meshes.emplace_back(new RegularMesh{}); + } + if (!index_end) *index_end = meshes.size() - 1; + + return 0; +} + +//! Return the index in the meshes array of a mesh with a given ID +extern "C" int +openmc_get_mesh_index(int32_t id, int32_t* index) +{ + auto pair = mesh_map.find(id); + if (pair == mesh_map.end()) { + set_errmsg("No mesh exists with ID=" + std::to_string(id) + "."); + return OPENMC_E_INVALID_ID; + } + *index = pair->second; + return 0; +} + +// Return the ID of a mesh +extern "C" int +openmc_mesh_get_id(int32_t index, int32_t* id) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + *id = meshes[index]->id_; + return 0; +} + +//! Set the ID of a mesh +extern "C" int +openmc_mesh_set_id(int32_t index, int32_t id) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + meshes[index]->id_ = id; + mesh_map[id] = index; + return 0; +} + +//! Get the dimension of a mesh +extern "C" int +openmc_mesh_get_dimension(int32_t index, int** dims, int* n) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + *dims = meshes[index]->shape_.data(); + *n = meshes[index]->n_dimension_; + return 0; +} + +//! Set the dimension of a mesh +extern "C" int +openmc_mesh_set_dimension(int32_t index, int n, const int* dims) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + // Copy dimension + std::vector shape = {static_cast(n)}; + auto& m = meshes[index]; + m->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape); + m->n_dimension_ = m->shape_.size(); + + return 0; +} + +//! Get the mesh parameters +extern "C" int +openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + auto& m = meshes[index]; + if (m->lower_left_.dimension() == 0) { + set_errmsg("Mesh parameters have not been set."); + return OPENMC_E_ALLOCATE; + } + + *ll = m->lower_left_.data(); + *ur = m->upper_right_.data(); + *n = m->n_dimension_; + return 0; +} + +//! Set the mesh parameters +extern "C" int +openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, + const double* width) +{ + if (index < 0 || index >= meshes.size()) { + set_errmsg("Index in meshes array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + + auto& m = meshes[index]; + std::vector shape = {static_cast(n)}; + if (ll && ur) { + m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); + m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); + m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_; + } else if (ll && width) { + m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); + m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); + m->upper_right_ = m->lower_left_ + m->shape_ * m->width_; + } else if (ur && width) { + m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); + m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); + m->lower_left_ = m->upper_right_ - m->shape_ * m->width_; + } else { + set_errmsg("At least two parameters must be specified."); + return OPENMC_E_INVALID_ARGUMENT; + } + + return 0; +} + } // namespace openmc diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index fb22d2ad31..52c220976c 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -52,6 +52,73 @@ module mesh_header ! Dictionary that maps user IDs to indices in 'meshes' type(DictIntInt), public :: mesh_dict + 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 + end interface + contains subroutine regular_from_xml(this, node) @@ -562,215 +629,4 @@ contains call mesh_dict % clear() end subroutine free_memory_mesh -!=============================================================================== -! C API FUNCTIONS -!=============================================================================== - - function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C) - ! Extend the meshes 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_INT) :: err - - type(RegularMesh), allocatable :: temp(:) ! temporary meshes array - - if (n_meshes == 0) then - ! Allocate meshes array - allocate(meshes(n)) - else - ! Allocate meshes array with increased size - allocate(temp(n_meshes + n)) - - ! Copy original meshes to temporary array - temp(1:n_meshes) = meshes - - ! Move allocation from temporary array - call move_alloc(FROM=temp, TO=meshes) - end if - - ! Return indices in meshes array - if (present(index_start)) index_start = n_meshes + 1 - if (present(index_end)) index_end = n_meshes + n - n_meshes = n_meshes + n - - err = 0 - end function openmc_extend_meshes - - - function openmc_get_mesh_index(id, index) result(err) bind(C) - ! Return the index in the meshes array of a mesh with a given ID - integer(C_INT32_T), value :: id - integer(C_INT32_T), intent(out) :: index - integer(C_INT) :: err - - if (allocated(meshes)) then - if (mesh_dict % has(id)) then - index = mesh_dict % get(id) - err = 0 - else - err = E_INVALID_ID - call set_errmsg("No mesh exists with ID=" // trim(to_str(id)) // ".") - end if - else - err = E_ALLOCATE - call set_errmsg("Memory has not been allocated for meshes.") - end if - end function openmc_get_mesh_index - - - function openmc_mesh_get_id(index, id) result(err) bind(C) - ! Return the ID of a mesh - integer(C_INT32_T), value :: index - integer(C_INT32_T), intent(out) :: id - integer(C_INT) :: err - - if (index >= 1 .and. index <= size(meshes)) then - id = meshes(index) % id - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_get_id - - - function openmc_mesh_set_id(index, id) result(err) bind(C) - ! Set the ID of a mesh - 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_meshes) then - meshes(index) % id = id - call mesh_dict % set(id, index) - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_set_id - - - function openmc_mesh_get_dimension(index, dims, n) result(err) bind(C) - ! Get the dimension of a mesh - integer(C_INT32_T), value, intent(in) :: index - type(C_PTR), intent(out) :: dims - integer(C_INT), intent(out) :: n - integer(C_INT) :: err - - if (index >= 1 .and. index <= n_meshes) then - dims = C_LOC(meshes(index) % dimension) - n = meshes(index) % n_dimension - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_get_dimension - - - function openmc_mesh_set_dimension(index, n, dims) result(err) bind(C) - ! Set the dimension of a mesh - 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 - - if (index >= 1 .and. index <= n_meshes) then - associate (m => meshes(index)) - if (allocated(m % dimension)) deallocate (m % dimension) - if (allocated(m % lower_left)) deallocate (m % lower_left) - if (allocated(m % upper_right)) deallocate (m % upper_right) - if (allocated(m % width)) deallocate (m % width) - - m % n_dimension = n - allocate(m % dimension(n)) - allocate(m % lower_left(n)) - allocate(m % upper_right(n)) - allocate(m % width(n)) - - ! Copy dimension - m % dimension(:) = dims - end associate - err = 0 - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_set_dimension - - - function openmc_mesh_get_params(index, ll, ur, width, n) result(err) bind(C) - ! Get the mesh parameters - 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 - - err = 0 - if (index >= 1 .and. index <= n_meshes) then - associate (m => meshes(index)) - if (allocated(m % lower_left)) then - ll = C_LOC(m % lower_left(1)) - ur = C_LOC(m % upper_right(1)) - width = C_LOC(m % width(1)) - n = m % n_dimension - else - err = E_ALLOCATE - call set_errmsg("Mesh parameters have not been set.") - end if - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_get_params - - - function openmc_mesh_set_params(index, n, ll, ur, width) result(err) bind(C) - ! Set the mesh parameters - 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 - - err = 0 - if (index >= 1 .and. index <= n_meshes) then - associate (m => meshes(index)) - if (allocated(m % lower_left)) deallocate (m % lower_left) - if (allocated(m % upper_right)) deallocate (m % upper_right) - if (allocated(m % width)) deallocate (m % width) - - allocate(m % lower_left(n)) - allocate(m % upper_right(n)) - allocate(m % width(n)) - - if (present(ll) .and. present(ur)) then - m % lower_left(:) = ll - m % upper_right(:) = ur - m % width(:) = (ur - ll) / m % dimension - elseif (present(ll) .and. present(width)) then - m % lower_left(:) = ll - m % width(:) = width - m % upper_right(:) = ll + width * m % dimension - elseif (present(ur) .and. present(width)) then - m % upper_right(:) = ur - m % width(:) = width - m % lower_left(:) = ur - width * m % dimension - else - err = E_INVALID_ARGUMENT - call set_errmsg("At least two parameters must be specified.") - end if - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in meshes array is out of bounds.") - end if - end function openmc_mesh_set_params - end module mesh_header diff --git a/src/settings.cpp b/src/settings.cpp index 7e727178a4..86f3c1fef4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -501,7 +501,7 @@ read_settings_xml() int temp = std::stoi(get_node_value(root, "entropy_mesh")); try { index_entropy_mesh = mesh_map.at(temp); - } catch (std::out_of_range) { + } catch (const std::out_of_range& e) { std::stringstream msg; msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; fatal_error(msg); @@ -550,7 +550,7 @@ read_settings_xml() auto temp = std::stoi(get_node_value(root, "ufs_mesh")); try { index_ufs_mesh = mesh_map.at(temp); - } catch (std::out_of_range) { + } catch (const std::out_of_range& e) { std::stringstream msg; msg << "Mesh " << temp << " specified for uniform fission site method " "does not exist."; From 6a1d653547b2e5ce65fd1b9f830d2509901b564b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 06:44:37 -0500 Subject: [PATCH 06/21] Convert count_bank_sites, Shannon entropy, UFS to C++ (incomplete) --- include/openmc/capi.h | 1 + include/openmc/eigenvalue.h | 37 ++++++++ include/openmc/mesh.h | 3 + include/openmc/message_passing.h | 1 + src/api.F90 | 4 +- src/bank_header.F90 | 18 +++- src/eigenvalue.F90 | 95 ------------------- src/eigenvalue.cpp | 151 +++++++++++++++++++++++++++++++ src/initialize.cpp | 16 ++-- src/input_xml.F90 | 104 --------------------- src/mesh.cpp | 123 +++++++++++++++++++++++++ src/mesh_header.F90 | 71 +++++++-------- src/message_passing.cpp | 1 + src/output.F90 | 12 ++- src/physics.F90 | 23 ++--- src/physics_mg.F90 | 24 ++--- src/settings.cpp | 8 -- src/simulation.F90 | 12 ++- src/simulation_header.F90 | 9 -- src/state_point.F90 | 48 ++++------ src/summary.F90 | 1 - src/tallies/tally.F90 | 1 - 22 files changed, 438 insertions(+), 325 deletions(-) create mode 100644 include/openmc/eigenvalue.h create mode 100644 src/eigenvalue.cpp diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 7ce42aa478..c155e7bd13 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -37,6 +37,7 @@ extern "C" { int openmc_filter_set_type(int32_t index, const char* type); int openmc_finalize(); int openmc_find_cell(double* xyz, int32_t* index, int32_t* instance); + int openmc_fission_bank(struct Bank** ptr, int64_t* n); int openmc_get_cell_index(int32_t id, int32_t* index); int openmc_get_filter_index(int32_t id, int32_t* index); void openmc_get_filter_next_id(int32_t* id); diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h new file mode 100644 index 0000000000..c2e76e26f6 --- /dev/null +++ b/include/openmc/eigenvalue.h @@ -0,0 +1,37 @@ +#ifndef OPENMC_EIGENVALUE_H +#define OPENMC_EIGENVALUE_H + +#include // for int64_t + +#include "openmc/particle.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +extern std::vector entropy; //!< Shannon entropy at each generation +extern xt::xtensor source_frac; //!< Source fraction for UFS + +extern "C" int64_t n_bank; + +//============================================================================== +// Non-member functions +//============================================================================== + +//! Calculates the Shannon entropy of the fission source distribution to assess +//! source convergence +extern "C" void shannon_entropy(); + +//! Determines the source fraction in each UFS mesh cell and reweights the +//! source bank so that the sum of the weights is equal to n_particles. The +//! 'source_frac' variable is used later to bias the production of fission sites +extern "C" void ufs_count_sites(); + +//! Get UFS weight corresponding to particle's location +extern "C" double ufs_get_weight(const Particle* p); + +} // namespace openmc + +#endif // OPENMC_EIGENVALUE_H diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 71649fa619..51a0a627f5 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -34,6 +34,9 @@ public: bool intersects(Position r0, Position r1); void to_hdf5(hid_t group); + xt::xarray count_sites(int64_t n, const Bank* bank, + int n_energy, const double* energies, bool* outside); + int id_ {-1}; //!< User-specified ID int n_dimension_; double volume_frac_; diff --git a/include/openmc/message_passing.h b/include/openmc/message_passing.h index 1dab43139b..c910210db5 100644 --- a/include/openmc/message_passing.h +++ b/include/openmc/message_passing.h @@ -10,6 +10,7 @@ namespace mpi { extern int rank; extern int n_procs; + extern bool master; #ifdef OPENMC_MPI extern MPI_Datatype bank; diff --git a/src/api.F90 b/src/api.F90 index 2a6775a0c8..b39a2f8fa6 100644 --- a/src/api.F90 +++ b/src/api.F90 @@ -11,7 +11,6 @@ module openmc_api use hdf5_interface use material_header use math - use mesh_header, only: free_memory_mesh use message_passing use nuclide_header use initialize, only: openmc_init_f @@ -305,6 +304,9 @@ contains interface subroutine free_memory_source() bind(C) end subroutine + + subroutine free_memory_mesh() bind(C) + end subroutine free_memory_mesh end interface call free_memory_geometry() diff --git a/src/bank_header.F90 b/src/bank_header.F90 index 2a998740eb..c39bf8b784 100644 --- a/src/bank_header.F90 +++ b/src/bank_header.F90 @@ -28,7 +28,7 @@ module bank_header type(Bank), allocatable, target :: master_fission_bank(:) #endif - integer(8) :: n_bank ! # of sites in fission bank + integer(C_INT64_T), bind(C) :: n_bank ! # of sites in fission bank !$omp threadprivate(fission_bank, n_bank) @@ -71,4 +71,20 @@ contains end if end function openmc_source_bank + function openmc_fission_bank(ptr, n) result(err) bind(C) + ! Return a pointer to the source bank + type(C_PTR), intent(out) :: ptr + integer(C_INT64_T), intent(out) :: n + integer(C_INT) :: err + + if (.not. allocated(fission_bank)) then + err = E_ALLOCATE + call set_errmsg("Fission bank has not been allocated.") + else + err = 0 + ptr = C_LOC(fission_bank) + n = size(fission_bank) + end if + end function openmc_fission_bank + end module bank_header diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index d47bc4dab8..38284c2798 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -294,46 +294,6 @@ contains end subroutine synchronize_bank -!=============================================================================== -! SHANNON_ENTROPY calculates the Shannon entropy of the fission source -! distribution to assess source convergence -!=============================================================================== - - subroutine shannon_entropy() - - integer :: i ! index for mesh elements - real(8) :: entropy_gen ! entropy at this generation - logical :: sites_outside ! were there sites outside entropy box? - - associate (m => meshes(index_entropy_mesh)) - ! count number of fission sites over mesh - call count_bank_sites(m, fission_bank, entropy_p, & - size_bank=n_bank, sites_outside=sites_outside) - - ! display warning message if there were sites outside entropy box - if (sites_outside) then - if (master) call warning("Fission source site(s) outside of entropy box.") - end if - - ! sum values to obtain shannon entropy - if (master) then - ! Normalize to total weight of bank sites - entropy_p = entropy_p / sum(entropy_p) - - entropy_gen = ZERO - do i = 1, size(entropy_p, 2) - if (entropy_p(1,i) > ZERO) then - entropy_gen = entropy_gen - & - entropy_p(1,i) * log(entropy_p(1,i))/log(TWO) - end if - end do - - ! Add value to vector - call entropy % push_back(entropy_gen) - end if - end associate - end subroutine shannon_entropy - !=============================================================================== ! CALCULATE_GENERATION_KEFF collects the single-processor tracklength k's onto ! the master processor and normalizes them. This should work whether or not the @@ -577,61 +537,6 @@ contains end function openmc_get_keff -!=============================================================================== -! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and -! reweights the source bank so that the sum of the weights is equal to -! n_particles. The 'source_frac' variable is used later to bias the production -! of fission sites -!=============================================================================== - - subroutine count_source_for_ufs() - - real(8) :: total ! total weight in source bank - logical :: sites_outside ! were there sites outside the ufs mesh? -#ifdef OPENMC_MPI - integer :: n ! total number of ufs mesh cells - integer :: mpi_err ! MPI error code -#endif - - associate (m => meshes(index_ufs_mesh)) - - if (current_batch == 1 .and. current_gen == 1) then - ! On the first generation, just assume that the source is already evenly - ! distributed so that effectively the production of fission sites is not - ! biased - - source_frac = m % volume_frac - - else - ! count number of source sites in each ufs mesh cell - call count_bank_sites(m, source_bank, source_frac, & - sites_outside=sites_outside, size_bank=work) - - ! Check for sites outside of the mesh - if (master .and. sites_outside) then - call fatal_error("Source sites outside of the UFS mesh!") - end if - -#ifdef OPENMC_MPI - ! Send source fraction to all processors - n = product(m % dimension) - call MPI_BCAST(source_frac, n, MPI_REAL8, 0, mpi_intracomm, mpi_err) -#endif - - ! Normalize to total weight to get fraction of source in each cell - total = sum(source_frac) - source_frac = source_frac / total - - ! Since the total starting weight is not equal to n_particles, we need to - ! renormalize the weight of the source sites - - source_bank % wgt = source_bank % wgt * n_particles / total - end if - - end associate - - end subroutine count_source_for_ufs - #ifdef _OPENMP !=============================================================================== ! JOIN_BANK_FROM_THREADS joins threadprivate fission banks into a single fission diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp new file mode 100644 index 0000000000..9ce5b9ef80 --- /dev/null +++ b/src/eigenvalue.cpp @@ -0,0 +1,151 @@ +#include "openmc/eigenvalue.h" + +#include "xtensor/xmath.hpp" +#include "xtensor/xtensor.hpp" +#include "xtensor/xview.hpp" + +#include "openmc/capi.h" +#include "openmc/error.h" +#include "openmc/hdf5_interface.h" +#include "openmc/mesh.h" +#include "openmc/message_passing.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +std::vector entropy; +xt::xtensor source_frac; + +//============================================================================== +// Non-member functions +//============================================================================== + +void shannon_entropy() +{ + // Get pointer to entropy mesh + auto& m = meshes[settings::index_entropy_mesh]; + + // Get pointer to fission bank + Bank* fission_bank; + int64_t n; + openmc_fission_bank(&fission_bank, &n); + + // Get source weight in each mesh bin + bool sites_outside; + xt::xtensor p = m->count_sites( + n_bank, fission_bank, 0, nullptr, &sites_outside); + + // display warning message if there were sites outside entropy box + if (sites_outside) { + if (mpi::master) warning("Fission source site(s) outside of entropy box."); + } + + // sum values to obtain shannon entropy + if (mpi::master) { + // Normalize to total weight of bank sites + p /= xt::sum(p); + + double H = 0.0; + for (auto p_i : p) { + if (p_i > 0.0) { + H -= p_i * std::log(p_i)/std::log(2.0); + } + } + + // Add value to vector + entropy.push_back(H); + } +} + +void ufs_count_sites() +{ + auto &m = meshes[settings::index_ufs_mesh]; + + if (openmc_current_batch == 1 && openmc_current_gen == 1) { + // On the first generation, just assume that the source is already evenly + // distributed so that effectively the production of fission sites is not + // biased + + auto s = xt::view(source_frac, xt::all()); + s = m->volume_frac_; + + } else { + // Get pointer to source bank + Bank* source_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + + + // count number of source sites in each ufs mesh cell + bool sites_outside; + source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr, + &sites_outside); + + // Check for sites outside of the mesh + if (mpi::master && sites_outside) { + fatal_error("Source sites outside of the UFS mesh!"); + } + +#ifdef OPENMC_MPI + // Send source fraction to all processors + int n = xt::prod(m->shape_)(); + MPI_Bcast(source_frac.data(), n, MPI_DOUBLE, 0, mpi::intracomm) +#endif + + // Normalize to total weight to get fraction of source in each cell + double total = xt::sum(source_frac)(); + source_frac /= total; + + // Since the total starting weight is not equal to n_particles, we need to + // renormalize the weight of the source sites + for (int i = 0; i < openmc_work; ++i) { + source_bank[i].wgt *= settings::n_particles / total; + } + } +} + +double ufs_get_weight(const Particle* p) +{ + auto& m = meshes[settings::index_entropy_mesh]; + + // Determine indices on ufs mesh for current location + // TODO: off by one + int mesh_bin = m->get_bin({p->coord[0].xyz}) - 1; + if (mesh_bin < 0) { + p->write_restart(); + fatal_error("Source site outside UFS mesh!"); + } + + if (source_frac(mesh_bin) != 0.0) { + return m->volume_frac_ / source_frac(mesh_bin); + } else { + return 1.0; + } +} + +extern "C" void entropy_to_hdf5(hid_t group) +{ + if (settings::entropy_on) { + write_dataset(group, "entropy", entropy); + } +} + +extern "C" void entropy_from_hdf5(hid_t group) +{ + if (settings::entropy_on) { + read_dataset(group, "entropy", entropy); + } +} + +extern "C" double entropy_c(int i) +{ + return entropy.at(i - 1); +} + + +} // namespace openmc diff --git a/src/initialize.cpp b/src/initialize.cpp index 2c68401513..19fc82c859 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -61,7 +61,7 @@ namespace openmc { #ifdef OPENMC_MPI void initialize_mpi(MPI_Comm intracomm) { - openmc::mpi::intracomm = intracomm; + mpi::intracomm = intracomm; // Initialize MPI int flag; @@ -69,13 +69,13 @@ void initialize_mpi(MPI_Comm intracomm) if (!flag) MPI_Init(nullptr, nullptr); // Determine number of processes and rank for each - MPI_Comm_size(intracomm, &openmc::mpi::n_procs); - MPI_Comm_rank(intracomm, &openmc::mpi::rank); + MPI_Comm_size(intracomm, &mpi::n_procs); + MPI_Comm_rank(intracomm, &mpi::rank); // Set variable for Fortran side - openmc_n_procs = openmc::mpi::n_procs; - openmc_rank = openmc::mpi::rank; - openmc_master = (openmc::mpi::rank == 0); + openmc_n_procs = mpi::n_procs; + openmc_rank = mpi::rank; + openmc_master = mpi::master = (mpi::rank == 0); // Create bank datatype Bank b; @@ -88,8 +88,8 @@ void initialize_mpi(MPI_Comm intracomm) }; int blocks[] {1, 3, 3, 1, 1}; MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_INT}; - MPI_Type_create_struct(5, blocks, disp, types, &openmc::mpi::bank); - MPI_Type_commit(&openmc::mpi::bank); + MPI_Type_create_struct(5, blocks, disp, types, &mpi::bank); + MPI_Type_commit(&mpi::bank); } #endif // OPENMC_MPI diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ffee754503..61b3208043 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -249,110 +249,6 @@ contains track_identifiers = reshape(temp_int_array, [3, n_tracks/3]) end if - ! Read meshes - call get_node_list(root, "mesh", node_mesh_list) - - ! Check for user meshes and allocate - n = size(node_mesh_list) - if (n > 0) then - err = openmc_extend_meshes(n, i_start, i_end) - end if - - do i = 1, n - associate (m => meshes(i_start + i - 1)) - ! Instantiate mesh from XML node - call m % from_xml(node_mesh_list(i)) - - ! Add mesh to dictionary - call mesh_dict % set(m % id, i_start + i - 1) - end associate - end do - - ! Shannon Entropy mesh - if (check_for_node(root, "entropy_mesh")) then - call get_node_value(root, "entropy_mesh", temp_int) - if (mesh_dict % has(temp_int)) then - index_entropy_mesh = mesh_dict % get(temp_int) - else - call fatal_error("Mesh " // to_str(temp_int) // " specified for & - &Shannon entropy does not exist.") - end if - elseif (check_for_node(root, "entropy")) then - call warning("Specifying a Shannon entropy mesh via the element & - &is deprecated. Please create a mesh using and then reference & - &it by specifying its ID in an element.") - - ! Get pointer to entropy node - node_entropy = root % child("entropy") - - err = openmc_extend_meshes(1, index_entropy_mesh) - - associate (m => meshes(index_entropy_mesh)) - ! Assign ID - m % id = 10000 - - call m % from_xml(node_entropy) - end associate - end if - - if (index_entropy_mesh > 0) then - associate(m => meshes(index_entropy_mesh)) - if (.not. allocated(m % dimension)) then - ! If the user did not specify how many mesh cells are to be used in - ! each direction, we automatically determine an appropriate number of - ! cells - m % n_dimension = 3 - allocate(m % dimension(3)) - m % dimension = ceiling((n_particles/20)**(ONE/THREE)) - - ! Calculate width - m % width = (m % upper_right - m % lower_left) / m % dimension - end if - - ! Allocate space for storing number of fission sites in each mesh cell - allocate(entropy_p(1, product(m % dimension))) - end associate - - ! Turn on Shannon entropy calculation - entropy_on = .true. - end if - - ! Uniform fission source weighting mesh - if (check_for_node(root, "ufs_mesh")) then - call get_node_value(root, "ufs_mesh", temp_int) - if (mesh_dict % has(temp_int)) then - index_ufs_mesh = mesh_dict % get(temp_int) - else - call fatal_error("Mesh " // to_str(temp_int) // " specified for & - &uniform fission site method does not exist.") - end if - elseif (check_for_node(root, "uniform_fs")) then - call warning("Specifying a UFS mesh via the element & - &is deprecated. Please create a mesh using and then reference & - &it by specifying its ID in a element.") - - ! Get pointer to ufs node - node_ufs = root % child("uniform_fs") - - err = openmc_extend_meshes(1, index_ufs_mesh) - - ! Allocate mesh object and coordinates on mesh - associate (m => meshes(index_ufs_mesh)) - ! Assign ID - m % id = 10001 - - call m % from_xml(node_ufs) - end associate - end if - - if (index_ufs_mesh > 0) then - ! Allocate array to store source fraction for UFS - allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension))) - - ! Turn on uniform fission source weighting - ufs = .true. - end if - ! Check if the user has specified to write state points if (check_for_node(root, "state_point")) then diff --git a/src/mesh.cpp b/src/mesh.cpp index 8aecd1a22a..47f03ad3de 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,9 +1,14 @@ #include "openmc/mesh.h" +#include // for copy #include // for size_t #include // for ceil #include +#ifdef OPENMC_MPI +#include "mpi.h" +#endif +#include "xtensor/xbuilder.hpp" #include "xtensor/xeval.hpp" #include "xtensor/xmath.hpp" #include "xtensor/xsort.hpp" @@ -13,6 +18,8 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" +#include "openmc/message_passing.h" +#include "openmc/search.h" #include "openmc/xml_interface.h" namespace openmc { @@ -521,6 +528,76 @@ void RegularMesh::to_hdf5(hid_t group) close_group(mesh_group); } +xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, + int n_energy, const double* energies, bool* outside) +{ + // Determine shape of array for counts + int m = xt::prod(shape_)(); + std::vector shape; + if (n_energy > 0) { + shape = {m, n_energy}; + } else { + shape = {m}; + } + + // Create array of zeros + xt::xarray cnt {shape, 0.0}; + bool outside_ = false; + + for (int64_t i = 0; i < n; ++i) { + // determine scoring bin for entropy mesh + int mesh_bin = get_bin({bank[i].xyz}); + + // if outside mesh, skip particle + if (mesh_bin < 0) { + outside_ = true; + continue; + } + + if (n_energy > 0) { + double E = bank[i].E; + int e_bin; + if (E >= energies[0] && E <= energies[n_energy - 1]) { + // determine energy bin + int e_bin = lower_bound_index(energies, energies + n_energy, E); + + // Add to appropriate bin + cnt(mesh_bin, e_bin) += bank[i].wgt; + } + } else { + // Add to appropriate bin + cnt(mesh_bin) += bank[i].wgt; + } + } + + // Create copy of count data + int total = cnt.size(); + double* cnt_reduced = new double[total]; + +#ifdef OPENMC_MPI + // collect values from all processors + MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, + mpi::intracomm); + if (outside) { + MPI_Reduce(&outside_, outside, 1, MPI_BOOL, MPI_LOR, 0, mpi::intracomm); + } + + // Check if there were sites outside the mesh for any processor + MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & + mpi_intracomm, mpi_err) + end if +#else + std::copy(cnt.data(), cnt.data() + total, cnt_reduced); + if (outside) *outside = outside_; +#endif + + // Adapt reduced values in array back into an xarray + auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape); + xt::xarray counts = arr; + + return counts; +} + //============================================================================== // C API functions //============================================================================== @@ -660,4 +737,50 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, return 0; } +//============================================================================== +// Fortran compatibility +//============================================================================== + +extern "C" { + int32_t mesh_id(RegularMesh* m) { return m->id_; } + + double mesh_volume_frac(RegularMesh* m) { return m->volume_frac_; } + + int mesh_n_dimension(RegularMesh* m) { return m->n_dimension_; } + + double* mesh_lower_left(RegularMesh* m) { return m->lower_left_.data(); } + + int mesh_get_bin(RegularMesh* m, const double* xyz) + { + return m->get_bin({xyz}); + } + + void meshes_to_hdf5(hid_t group) + { + // Write number of meshes + hid_t meshes_group = create_group(group, "meshes"); + int32_t n_meshes = meshes.size(); + write_attribute(meshes_group, "n_meshes", n_meshes); + + if (n_meshes > 0) { + // Write IDs of meshes + std::vector ids; + for (const auto& m : meshes) { + m->to_hdf5(meshes_group); + ids.push_back(m->id_); + } + write_attribute(meshes_group, "ids", ids); + } + + close_group(meshes_group); + } + + void free_memory_mesh() + { + meshes.clear(); + mesh_map.clear(); + } +} + + } // namespace openmc diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 52c220976c..8a8d16c084 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -11,7 +11,6 @@ module mesh_header implicit none private - public :: free_memory_mesh public :: openmc_extend_meshes public :: openmc_get_mesh_index public :: openmc_mesh_get_id @@ -27,22 +26,23 @@ module mesh_header !=============================================================================== type, public :: RegularMesh - integer :: id = -1 ! user-specified id - integer :: type = MESH_REGULAR ! rectangular, hexagonal + type(C_PTR) :: ptr integer(C_INT) :: n_dimension ! rank of mesh - real(8) :: volume_frac ! volume fraction of each cell integer(C_INT), allocatable :: dimension(:) ! number of cells in each direction real(C_DOUBLE), allocatable :: lower_left(:) ! lower-left corner of mesh real(C_DOUBLE), allocatable :: upper_right(:) ! upper-right corner of mesh real(C_DOUBLE), allocatable :: width(:) ! width of each mesh cell contains + procedure :: id => regular_id + procedure :: volume_frac => regular_volume_frac + + procedure :: from_xml => regular_from_xml 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 procedure :: intersects => regular_intersects - procedure :: to_hdf5 => regular_to_hdf5 end type RegularMesh integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes @@ -121,6 +121,36 @@ module mesh_header contains + function regular_id(this) result(id) + class(RegularMesh), intent(in) :: this + integer(C_INT32_T) :: id + + interface + 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 + end interface + + 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 + + interface + 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 + end interface + + volume_frac = mesh_volume_frac(this % ptr) + end function + subroutine regular_from_xml(this, node) class(RegularMesh), intent(inout) :: this type(XMLNode), intent(in) :: node @@ -598,35 +628,4 @@ contains end function mesh_intersects_3d -!=============================================================================== -! TO_HDF5 writes the mesh data to an HDF5 group -!=============================================================================== - - subroutine regular_to_hdf5(this, group) - class(RegularMesh), intent(in) :: this - integer(HID_T), intent(in) :: group - - integer(HID_T) :: mesh_group - - mesh_group = create_group(group, "mesh " // trim(to_str(this % id))) - - call write_dataset(mesh_group, "type", "regular") - call write_dataset(mesh_group, "dimension", this % dimension) - call write_dataset(mesh_group, "lower_left", this % lower_left) - call write_dataset(mesh_group, "upper_right", this % upper_right) - call write_dataset(mesh_group, "width", this % width) - - call close_group(mesh_group) - end subroutine regular_to_hdf5 - -!=============================================================================== -! FREE_MEMORY_MESH deallocates global arrays defined in this module -!=============================================================================== - - subroutine free_memory_mesh() - n_meshes = 0 - if (allocated(meshes)) deallocate(meshes) - call mesh_dict % clear() - end subroutine free_memory_mesh - end module mesh_header diff --git a/src/message_passing.cpp b/src/message_passing.cpp index aac677c6ef..ec2fc2d688 100644 --- a/src/message_passing.cpp +++ b/src/message_passing.cpp @@ -5,6 +5,7 @@ namespace mpi { int rank {0}; int n_procs {1}; +bool master {false}; #ifdef OPENMC_MPI MPI_Comm intracomm; diff --git a/src/output.F90 b/src/output.F90 index b7b418a39b..e286956670 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -33,6 +33,14 @@ module output integer :: ou = OUTPUT_UNIT integer :: eu = ERROR_UNIT + interface + function entropy(i) result(h) bind(C, name='entropy_c') + import C_INT, C_DOUBLE + integer(C_INT), value :: i + real(C_DOUBLE) :: h + end function + end interface + contains !=============================================================================== @@ -335,7 +343,7 @@ contains ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy % data(i) + entropy(i) if (n > 1) then write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5," +/-",F8.5)', ADVANCE='NO') & @@ -369,7 +377,7 @@ contains ! write out entropy info if (entropy_on) write(UNIT=OUTPUT_UNIT, FMT='(3X, F8.5)', ADVANCE='NO') & - entropy % data(i) + entropy(i) ! write out accumulated k-effective if after first active batch if (n > 1) then diff --git a/src/physics.F90 b/src/physics.F90 index 992e487675..f7d80261ad 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1187,6 +1187,14 @@ contains real(8) :: weight ! weight adjustment for ufs method type(Nuclide), pointer :: nuc + interface + function ufs_get_weight(p) result(weight) bind(C) + import Particle, C_DOUBLE + type(Particle), intent(in) :: p + real(C_DOUBLE) :: WEIGHT + end function + end interface + ! Get pointers nuc => nuclides(i_nuclide) @@ -1196,20 +1204,7 @@ contains ! the expected number of fission sites produced if (ufs) then - associate (m => meshes(index_ufs_mesh)) - ! Determine indices on ufs mesh for current location - call m % get_bin(p % coord(1) % xyz, mesh_bin) - if (mesh_bin == NO_BIN_FOUND) then - call particle_write_restart(p) - call fatal_error("Source site outside UFS mesh!") - end if - - if (source_frac(1, mesh_bin) /= ZERO) then - weight = m % volume_frac / source_frac(1, mesh_bin) - else - weight = ONE - end if - end associate + weight = ufs_get_weight(p) else weight = ONE end if diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 40df6e1cb5..84695153cb 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -174,27 +174,21 @@ contains real(8) :: phi ! fission neutron azimuthal angle real(8) :: weight ! weight adjustment for ufs method + interface + function ufs_get_weight(p) result(weight) bind(C) + import Particle, C_DOUBLE + type(Particle), intent(in) :: p + real(C_DOUBLE) :: WEIGHT + end function + end interface + ! TODO: Heat generation from fission ! If uniform fission source weighting is turned on, we increase of decrease ! the expected number of fission sites produced if (ufs) then - associate (m => meshes(index_ufs_mesh)) - ! Determine indices on ufs mesh for current location - call m % get_bin(p % coord(1) % xyz, mesh_bin) - - if (mesh_bin == NO_BIN_FOUND) then - call particle_write_restart(p) - call fatal_error("Source site outside UFS mesh!") - end if - - if (source_frac(1, mesh_bin) /= ZERO) then - weight = m % volume_frac / source_frac(1, mesh_bin) - else - weight = ONE - end if - end associate + weight = ufs_get_weight(p) else weight = ONE end if diff --git a/src/settings.cpp b/src/settings.cpp index 86f3c1fef4..7b00b9eca7 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -537,10 +537,6 @@ read_settings_xml() m.width_ = (m.upper_right_ - m.lower_left_) / m.shape_; } - // TODO: Allocate entropy_P - // Allocate space for storing number of fission sites in each mesh cell - //allocate(entropy_p(1, product(m % dimension))) - // Turn on Shannon entropy calculation settings::entropy_on = true; } @@ -574,10 +570,6 @@ read_settings_xml() } if (index_ufs_mesh >= 0) { - // Allocate array to store source fraction for UFS - // TODO: Allocate source_frac - //allocate(source_frac(1, product(meshes(index_ufs_mesh) % dimension))) - // Turn on uniform fission source weighting settings::ufs_on = true; } diff --git a/src/simulation.F90 b/src/simulation.F90 index d334b46b8c..8b7cfdb27a 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -10,7 +10,7 @@ module simulation use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd use cmfd_header, only: cmfd_on use constants, only: ZERO - use eigenvalue, only: count_source_for_ufs, calculate_average_keff, & + use eigenvalue, only: calculate_average_keff, & calculate_generation_keff, shannon_entropy, & synchronize_bank, keff_generation, k_sum #ifdef _OPENMP @@ -234,12 +234,17 @@ contains subroutine initialize_generation() + interface + subroutine ufs_count_sites() bind(C) + end subroutine + end interface + if (run_mode == MODE_EIGENVALUE) then ! Reset number of fission bank sites n_bank = 0 ! Count source sites if using uniform fission source weighting - if (ufs) call count_source_for_ufs() + if (ufs) call ufs_count_sites() ! Store current value of tracklength k keff_generation = global_tallies(RESULT_VALUE, K_TRACKLENGTH) @@ -256,6 +261,9 @@ contains interface subroutine fill_source_bank_fixedsource() bind(C) end subroutine + + subroutine shannon_entropy() bind(C) + end subroutine end interface ! Update global tallies with the omp private accumulation variables diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 0e38b76f25..9b7c68cd46 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -47,13 +47,6 @@ module simulation_header real(8) :: k_col_tra = ZERO ! sum over batches of k_collision * k_tracklength real(8) :: k_abs_tra = ZERO ! sum over batches of k_absorption * k_tracklength - ! Shannon entropy - type(VectorReal) :: entropy ! shannon entropy at each generation - real(8), allocatable :: entropy_p(:,:) ! % of source sites in each cell - - ! Uniform fission source weighting - real(8), allocatable :: source_frac(:,:) - ! ============================================================================ ! PARALLEL PROCESSING VARIABLES @@ -87,8 +80,6 @@ contains !=============================================================================== subroutine free_memory_simulation() - if (allocated(entropy_p)) deallocate(entropy_p) - if (allocated(source_frac)) deallocate(source_frac) if (allocated(work_index)) deallocate(work_index) call k_generation % clear() diff --git a/src/state_point.F90 b/src/state_point.F90 index bc6474de22..b7f4e13776 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -20,7 +20,6 @@ module state_point use endf, only: reaction_name use error, only: fatal_error, warning, write_message use hdf5_interface - use mesh_header, only: RegularMesh, meshes, n_meshes use message_passing use mgxs_interface use nuclide_header, only: nuclides @@ -79,6 +78,17 @@ contains character(MAX_WORD_LEN, kind=C_CHAR) :: temp_name logical :: parallel + interface + subroutine meshes_to_hdf5(group) bind(C) + import HID_T + integer(HID_T), value :: group + end subroutine + subroutine entropy_to_hdf5(group) bind(C) + import HID_T + integer(HID_T), value :: group + end subroutine + end interface + err = 0 ! Set the filename @@ -163,9 +173,7 @@ contains call write_dataset(file_id, "generations_per_batch", gen_per_batch) k = k_generation % size() call write_dataset(file_id, "k_generation", k_generation % data(1:k)) - if (entropy_on) then - call write_dataset(file_id, "entropy", entropy % data(1:k)) - end if + call entropy_to_hdf5() call write_dataset(file_id, "k_col_abs", k_col_abs) call write_dataset(file_id, "k_col_tra", k_col_tra) call write_dataset(file_id, "k_abs_tra", k_abs_tra) @@ -192,26 +200,8 @@ contains tallies_group = create_group(file_id, "tallies") - ! Write number of meshes - meshes_group = create_group(tallies_group, "meshes") - call write_attribute(meshes_group, "n_meshes", n_meshes) - - if (n_meshes > 0) then - ! Write IDs of meshes - allocate(id_array(n_meshes)) - do i = 1, n_meshes - id_array(i) = meshes(i) % id - end do - call write_attribute(meshes_group, "ids", id_array) - deallocate(id_array) - - ! Write information for meshes - MESH_LOOP: do i = 1, n_meshes - call meshes(i) % to_hdf5(meshes_group) - end do MESH_LOOP - end if - - call close_group(meshes_group) + ! Write meshes + call meshes_to_hdf5(tallies_group) ! Write information for derivatives. if (size(tally_derivs) > 0) then @@ -641,6 +631,11 @@ contains logical :: source_present character(MAX_WORD_LEN) :: word + interface + subroutine entropy_from_hdf5() bind(C) + end subroutine + end interface + ! Write message call write_message("Loading state point " // trim(path_state_point) & // "...", 5) @@ -722,10 +717,7 @@ contains call k_generation % resize(n) call read_dataset(k_generation % data(1:n), file_id, "k_generation") - if (entropy_on) then - call entropy % resize(n) - call read_dataset(entropy % data(1:n), file_id, "entropy") - end if + call entropy_from_hdf5() call read_dataset(k_col_abs, file_id, "k_col_abs") call read_dataset(k_col_tra, file_id, "k_col_tra") call read_dataset(k_abs_tra, file_id, "k_abs_tra") diff --git a/src/summary.F90 b/src/summary.F90 index 750a9899f4..0fafd7e5d0 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -6,7 +6,6 @@ module summary use geometry_header use hdf5_interface use material_header, only: Material, n_materials, openmc_material_get_volume - use mesh_header, only: RegularMesh use message_passing use mgxs_interface use nuclide_header diff --git a/src/tallies/tally.F90 b/src/tallies/tally.F90 index 9967606d17..a618f351e6 100644 --- a/src/tallies/tally.F90 +++ b/src/tallies/tally.F90 @@ -8,7 +8,6 @@ module tally use error, only: fatal_error use geometry_header use math, only: t_percentile - use mesh_header, only: RegularMesh, meshes use message_passing use mgxs_interface use nuclide_header From 18f9f91833f8267796f63d7dde826c177300c0d0 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 09:57:26 -0500 Subject: [PATCH 07/21] Read meshes from tallies file --- include/openmc/mesh.h | 8 ++++++++ include/openmc/output.h | 2 ++ include/openmc/settings.h | 5 ++--- src/input_xml.F90 | 31 ++++++++++--------------------- src/mesh.cpp | 15 +++++++++++++++ src/settings.cpp | 14 ++------------ 6 files changed, 39 insertions(+), 36 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 51a0a627f5..28f9c3eac9 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -51,6 +51,14 @@ private: bool intersects_3d(Position r0, Position r1); }; +//============================================================================== +// Non-member functions +//============================================================================== + +//! Read meshes from either settings/tallies +//! \param[in] root XML node +extern "C" void read_meshes(pugi::xml_node* root); + //============================================================================== // Global variables //============================================================================== diff --git a/include/openmc/output.h b/include/openmc/output.h index cabd39edeb..1ba5d65991 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -22,5 +22,7 @@ void header(const char* msg, int level); extern "C" void print_overlap_check(); +extern "C" void title(); + } // namespace openmc #endif // OPENMC_OUTPUT_H diff --git a/include/openmc/settings.h b/include/openmc/settings.h index be8661f4dd..e5684c7ae6 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -89,13 +89,12 @@ extern "C" double weight_survive; //!< Survival weight after Russian roul } // namespace settings -//============================================================================== //! Read settings from XML file //! \param[in] root XML node for -//============================================================================== - extern "C" void read_settings_xml(); +extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr); + } // namespace openmc #endif // OPENMC_SETTINGS_H diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 61b3208043..bfd604f9ee 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -65,17 +65,17 @@ module input_xml subroutine read_surfaces(node_ptr) bind(C) import C_PTR - type(C_PTR) :: node_ptr + type(C_PTR), value :: node_ptr end subroutine read_surfaces subroutine read_cells(node_ptr) bind(C) import C_PTR - type(C_PTR) :: node_ptr + type(C_PTR), value :: node_ptr end subroutine read_cells subroutine read_lattices(node_ptr) bind(C) import C_PTR - type(C_PTR) :: node_ptr + type(C_PTR), value :: node_ptr end subroutine read_lattices subroutine read_settings_xml() bind(C) @@ -83,9 +83,14 @@ module input_xml subroutine read_materials(node_ptr) bind(C) import C_PTR - type(C_PTR) :: node_ptr + type(C_PTR), value :: node_ptr end subroutine read_materials + subroutine read_meshes(node_ptr) bind(C) + import C_PTR + type(C_PTR), value :: node_ptr + end subroutine + function find_root_universe() bind(C) result(root) import C_INT32_T integer(C_INT32_T) :: root @@ -1202,9 +1207,6 @@ contains ! ========================================================================== ! DETERMINE SIZE OF ARRAYS AND ALLOCATE - ! Get pointer list to XML - call get_node_list(root, "mesh", node_mesh_list) - ! Get pointer list to XML call get_node_list(root, "filter", node_filt_list) @@ -1220,20 +1222,7 @@ contains ! READ MESH DATA ! Check for user meshes and allocate - n = size(node_mesh_list) - if (n > 0) then - err = openmc_extend_meshes(n, i_start, i_end) - end if - - do i = 1, n - m => meshes(i_start + i - 1) - - ! Instantiate mesh from XML node - call m % from_xml(node_mesh_list(i)) - - ! Add mesh to dictionary - call mesh_dict % set(m % id, i_start + i - 1) - end do + call read_meshes() ! We only need the mesh info for plotting if (run_mode == MODE_PLOTTING) then diff --git a/src/mesh.cpp b/src/mesh.cpp index 47f03ad3de..113968242e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -737,6 +737,21 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, return 0; } +//============================================================================== +// Non-member functions +//============================================================================== + +void read_meshes(pugi::xml_node* root) +{ + for (auto node : root->children("mesh")) { + // Read mesh and add to vector + meshes.emplace_back(new RegularMesh{node}); + + // Map ID to position in vector + mesh_map[meshes.back()->id_] = meshes.size() - 1; + } +} + //============================================================================== // Fortran compatibility //============================================================================== diff --git a/src/settings.cpp b/src/settings.cpp index 7b00b9eca7..4883d07a30 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -170,11 +170,7 @@ void get_run_parameters(pugi::xml_node node_base) } } -extern "C" void title(); -extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr); - -extern "C" void -read_settings_xml() +void read_settings_xml() { using namespace settings; using namespace pugi; @@ -488,13 +484,7 @@ read_settings_xml() } // Read meshes - for (auto node : root.children("mesh")) { - // Read mesh and add to vector - meshes.emplace_back(new RegularMesh{node}); - - // Map ID to position in vector - mesh_map[meshes.back()->id_] = meshes.size() - 1; - } + read_meshes(&root); // Shannon Entropy mesh if (check_for_node(root, "entropy_mesh")) { From 666553da2775084bb0598910abd27399e042d039 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 10:38:53 -0500 Subject: [PATCH 08/21] Add bindings on RegularMesh Fortran type --- src/mesh.cpp | 25 +- src/mesh_header.F90 | 581 ++++++++++---------------------------------- src/plot.F90 | 85 +++---- src/plot_header.F90 | 3 +- 4 files changed, 196 insertions(+), 498 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 113968242e..9c955fae76 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -757,19 +757,42 @@ void read_meshes(pugi::xml_node* root) //============================================================================== extern "C" { + RegularMesh* mesh_ptr(int i) { return meshes[i-1].get(); } + int32_t mesh_id(RegularMesh* m) { return m->id_; } double mesh_volume_frac(RegularMesh* m) { return m->volume_frac_; } int mesh_n_dimension(RegularMesh* m) { return m->n_dimension_; } - double* mesh_lower_left(RegularMesh* m) { return m->lower_left_.data(); } + int mesh_dimension(RegularMesh* m, int i) { return m->shape_(i - 1); } + + double mesh_lower_left(RegularMesh* m, int i) { return m->lower_left_(i - 1); } + + double mesh_upper_right(RegularMesh* m, int i) { return m->upper_right_(i - 1); } + + double mesh_width(RegularMesh* m, int i) { return m->width_(i - 1); } int mesh_get_bin(RegularMesh* m, const double* xyz) { return m->get_bin({xyz}); } + int mesh_get_bin_from_indices(RegularMesh* m, const int* ijk) + { + return m->get_bin_from_indices(ijk); + } + + void mesh_get_indices(RegularMesh* m, const double* xyz, int* ijk, bool* in_mesh) + { + m->get_indices({xyz}, ijk, in_mesh); + } + + void mesh_get_indices_from_bin(RegularMesh* m, int bin, int* ijk) + { + m->get_indices_from_bin(bin, ijk); + } + void meshes_to_hdf5(hid_t group) { // Write number of meshes diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 8a8d16c084..591263000b 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -27,28 +27,22 @@ module mesh_header type, public :: RegularMesh type(C_PTR) :: ptr - integer(C_INT) :: n_dimension ! rank of mesh - integer(C_INT), allocatable :: dimension(:) ! number of cells in each direction - real(C_DOUBLE), allocatable :: lower_left(:) ! lower-left corner of mesh - real(C_DOUBLE), allocatable :: upper_right(:) ! upper-right corner of mesh - real(C_DOUBLE), allocatable :: width(:) ! width of each mesh cell contains procedure :: id => regular_id procedure :: volume_frac => regular_volume_frac + procedure :: n_dimension => regular_n_dimension + procedure :: lower_left => regular_lower_left + procedure :: upper_right => regular_upper_right + procedure :: width => regular_width - - procedure :: from_xml => regular_from_xml 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 - procedure :: intersects => regular_intersects end type RegularMesh integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes - type(RegularMesh), public, allocatable, target :: meshes(:) - ! Dictionary that maps user IDs to indices in 'meshes' type(DictIntInt), public :: mesh_dict @@ -117,167 +111,143 @@ module mesh_header 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 + + function mesh_get_bin(ptr, xyz) result(bin) bind(C) + import C_PTR, C_DOUBLE, C_INT + type(C_PTR), value :: ptr + real(C_DOUBLE), intent(in) :: xyz(*) + integer(C_INT) :: bin + end function + + function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), intent(in) :: ijk(*) + integer(C_INT) :: bin + end function + + subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) + import C_PTR, C_DOUBLE, C_INT, C_BOOL + type(C_PTR), value :: ptr + real(C_DOUBLE), intent(in) :: xyz(*) + integer(C_INT), intent(out) :: ijk(*) + logical(C_BOOL), intent(out) :: in_mesh + end subroutine + + subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: ptr + integer(C_INT), 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 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 - - interface - 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 - end interface - 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 - - interface - 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 - end interface - volume_frac = mesh_volume_frac(this % ptr) end function - subroutine regular_from_xml(this, node) - class(RegularMesh), intent(inout) :: this - type(XMLNode), intent(in) :: node + function regular_n_dimension(this) result(n) + class(RegularMesh), intent(in) :: this + integer(C_INT) :: n + n = mesh_n_dimension(this % ptr) + end function - integer :: n - character(MAX_LINE_LEN) :: temp_str + 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 - ! Copy mesh id - if (check_for_node(node, "id")) then - call get_node_value(node, "id", this % id) + 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 - ! Check to make sure 'id' hasn't been used - if (mesh_dict % has(this % id)) then - call fatal_error("Two or more meshes use the same unique ID: " & - // to_str(this % id)) - end if - end if + 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 - ! Read mesh type - if (check_for_node(node, "type")) then - call get_node_value(node, "type", temp_str) - select case (to_lower(temp_str)) - case ('rect', 'rectangle', 'rectangular') - call warning("Mesh type '" // trim(temp_str) // "' is deprecated. & - &Please use 'regular' instead.") - this % type = MESH_REGULAR - case ('regular') - this % type = MESH_REGULAR - case default - call fatal_error("Invalid mesh type: " // trim(temp_str)) - end select - else - this % type = MESH_REGULAR - end if - - ! Determine number of dimensions for mesh - if (check_for_node(node, "dimension")) then - n = node_word_count(node, "dimension") - if (n /= 1 .and. n /= 2 .and. n /= 3) then - call fatal_error("Mesh must be one, two, or three dimensions.") - end if - this % n_dimension = n - - ! Allocate attribute arrays - allocate(this % dimension(n)) - - ! Check that dimensions are all greater than zero - call get_node_array(node, "dimension", this % dimension) - if (any(this % dimension <= 0)) then - call fatal_error("All entries on the element for a tally & - &mesh must be positive.") - end if - end if - - ! Check for lower-left coordinates - if (check_for_node(node, "lower_left")) then - n = node_word_count(node, "lower_left") - allocate(this % lower_left(n)) - - ! Read mesh lower-left corner location - call get_node_array(node, "lower_left", this % lower_left) - else - call fatal_error("Must specify on a mesh.") - end if - - if (check_for_node(node, "width")) then - ! Make sure both upper-right or width were specified - if (check_for_node(node, "upper_right")) then - call fatal_error("Cannot specify both and on a & - &mesh.") - end if - - n = node_word_count(node, "width") - allocate(this % width(n)) - allocate(this % upper_right(n)) - - ! Check to ensure width has same dimensions - if (n /= size(this % lower_left)) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - - ! Check for negative widths - call get_node_array(node, "width", this % width) - if (any(this % width < ZERO)) then - call fatal_error("Cannot have a negative on a tally mesh.") - end if - - ! Set width and upper right coordinate - this % upper_right = this % lower_left + this % dimension * this % width - - elseif (check_for_node(node, "upper_right")) then - n = node_word_count(node, "upper_right") - allocate(this % upper_right(n)) - allocate(this % width(n)) - - ! Check to ensure width has same dimensions - if (n /= size(this % lower_left)) then - call fatal_error("Number of entries on must be the & - &same as the number of entries on .") - end if - - ! Check that upper-right is above lower-left - call get_node_array(node, "upper_right", this % upper_right) - if (any(this % upper_right < this % lower_left)) then - call fatal_error("The coordinates must be greater than & - &the coordinates on a tally mesh.") - end if - - ! Set width and upper right coordinate - this % width = (this % upper_right - this % lower_left) / this % dimension - else - call fatal_error("Must specify either and on a & - &mesh.") - end if - - if (allocated(this % dimension)) then - if (size(this % dimension) /= size(this % lower_left)) then - call fatal_error("Number of entries on must be the same & - &as the number of entries on .") - end if - - ! Set volume fraction - this % volume_frac = ONE/real(product(this % dimension),8) - end if - - end subroutine regular_from_xml + 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 @@ -288,38 +258,8 @@ contains real(8), intent(in) :: xyz(:) ! coordinates integer, intent(out) :: bin ! tally bin - integer :: n ! size of mesh - integer :: d ! mesh dimension index - integer :: ijk(3) ! indices in mesh - logical :: in_mesh ! was given coordinate in mesh at all? - - ! Get number of dimensions - n = this % n_dimension - - ! Loop over the dimensions of the mesh - do d = 1, n - - ! Check for cases where particle is outside of mesh - if (xyz(d) < this % lower_left(d)) then - bin = NO_BIN_FOUND - return - elseif (xyz(d) > this % upper_right(d)) then - bin = NO_BIN_FOUND - return - end if - end do - - ! Determine indices - call this % get_indices(xyz, ijk, in_mesh) - - ! Convert indices to bin - if (in_mesh) then - bin = this % get_bin_from_indices(ijk) - else - bin = NO_BIN_FOUND - end if - - end subroutine regular_get_bin + bin = mesh_get_bin(this % ptr, xyz) + end subroutine !=============================================================================== ! GET_MESH_INDICES determines the indices of a particle in a structured mesh @@ -331,18 +271,10 @@ contains integer, intent(out) :: ijk(:) ! indices in mesh logical, intent(out) :: in_mesh ! were given coords in mesh? - ! Find particle in mesh - ijk(:this % n_dimension) = ceiling((xyz(:this % n_dimension) - & - this % lower_left)/this % width) - - ! Determine if particle is in mesh - if (any(ijk(:this % n_dimension) < 1) .or. & - any(ijk(:this % n_dimension) > this % dimension)) then - in_mesh = .false. - else - in_mesh = .true. - end if + logical(C_BOOL) :: in_mesh_ + call mesh_get_indices(this % ptr, xyz, ijk, in_mesh_) + in_mesh = in_mesh_ end subroutine regular_get_indices !=============================================================================== @@ -355,15 +287,7 @@ contains integer, intent(in) :: ijk(:) integer :: bin - if (this % n_dimension == 1) then - bin = ijk(1) - elseif (this % n_dimension == 2) then - bin = (ijk(2) - 1) * this % dimension(1) + ijk(1) - elseif (this % n_dimension == 3) then - bin = ((ijk(3) - 1) * this % dimension(2) + (ijk(2) - 1)) & - * this % dimension(1) + ijk(1) - end if - + bin = mesh_get_bin_from_indices(this % ptr, ijk) end function regular_get_bin_from_indices !=============================================================================== @@ -376,256 +300,7 @@ contains integer, intent(in) :: bin integer, intent(out) :: ijk(:) - if (this % n_dimension == 1) then - ijk(1) = bin - else if (this % n_dimension == 2) then - ijk(1) = mod(bin - 1, this % dimension(1)) + 1 - ijk(2) = (bin - 1)/this % dimension(1) + 1 - else if (this % n_dimension == 3) then - ijk(1) = mod(bin - 1, this % dimension(1)) + 1 - ijk(2) = mod(bin - 1, this % dimension(1) * this % dimension(2)) & - / this % dimension(1) + 1 - ijk(3) = (bin - 1)/(this % dimension(1) * this % dimension(2)) + 1 - end if - + call mesh_get_indices_from_bin(this % ptr, bin, ijk) end subroutine regular_get_indices_from_bin -!=============================================================================== -! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the -! outer boundary of the given mesh. This is important for determining whether a -! track will score to a mesh tally. -!=============================================================================== - - pure function regular_intersects(this, xyz0, xyz1) result(intersects) - class(RegularMesh), intent(in) :: this - real(8), intent(in) :: xyz0(:) - real(8), intent(in) :: xyz1(:) - logical :: intersects - - select case(this % n_dimension) - case (1) - intersects = mesh_intersects_1d(this, xyz0, xyz1) - case (2) - intersects = mesh_intersects_2d(this, xyz0, xyz1) - case (3) - intersects = mesh_intersects_3d(this, xyz0, xyz1) - end select - end function regular_intersects - - pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(:) - real(8), intent(in) :: xyz1(:) - logical :: intersects - - real(8) :: x0 ! track start point - real(8) :: x1 ! track end point - real(8) :: xm0 ! lower-left coordinates of mesh - real(8) :: xm1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - - ! Copy coordinates of ending point - x1 = xyz1(1) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - intersects = .true. - return - end if - - ! Check if line intersects right surface - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - intersects = .true. - return - end if - - end function mesh_intersects_1d - - pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(:) - real(8), intent(in) :: xyz1(:) - logical :: intersects - - real(8) :: x0, y0 ! track start point - real(8) :: x1, y1 ! track end point - real(8) :: xi, yi ! track intersection point with mesh - real(8) :: xm0, ym0 ! lower-left coordinates of mesh - real(8) :: xm1, ym1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - y0 = xyz0(2) - - ! Copy coordinates of ending point - x1 = xyz1(1) - y1 = xyz1(2) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - ym0 = m % lower_left(2) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - ym1 = m % upper_right(2) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface -- calculate the intersection point - ! y - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects back surface -- calculate the intersection point - ! x - if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then - xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects right surface -- calculate the intersection - ! point y - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects front surface -- calculate the intersection point - ! x - if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then - xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1) then - intersects = .true. - return - end if - end if - - end function mesh_intersects_2d - - pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(:) - real(8), intent(in) :: xyz1(:) - logical :: intersects - - real(8) :: x0, y0, z0 ! track start point - real(8) :: x1, y1, z1 ! track end point - real(8) :: xi, yi, zi ! track intersection point with mesh - real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh - real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - y0 = xyz0(2) - z0 = xyz0(3) - - ! Copy coordinates of ending point - x1 = xyz1(1) - y1 = xyz1(2) - z1 = xyz1(3) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - ym0 = m % lower_left(2) - zm0 = m % lower_left(3) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - ym1 = m % upper_right(2) - zm1 = m % upper_right(3) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface -- calculate the intersection point - ! (y,z) - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects back surface -- calculate the intersection point - ! (x,z) - if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then - xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects bottom surface -- calculate the intersection - ! point (x,y) - if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then - xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) - if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects right surface -- calculate the intersection point - ! (y,z) - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects front surface -- calculate the intersection point - ! (x,z) - if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then - xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects top surface -- calculate the intersection point - ! (x,y) - if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then - xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) - if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - end function mesh_intersects_3d - end module mesh_header diff --git a/src/plot.F90 b/src/plot.F90 index 517c27fa17..b8f3db0796 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -9,6 +9,7 @@ module plot use hdf5_interface use output, only: time_stamp use material_header, only: materials + use mesh_header, only: meshes use particle_header use plot_header use progress_header, only: ProgressBar @@ -184,7 +185,7 @@ contains !$omp end parallel do ! Draw tally mesh boundaries on the image if requested - if (associated(pl % meshlines_mesh)) call draw_mesh_lines(pl, data) + if (pl % index_meshlines_mesh >= 0) call draw_mesh_lines(pl, data) ! Write out the ppm to a file call output_ppm(pl, data) @@ -214,6 +215,7 @@ contains real(8) :: xyz_ur_plot(3) ! upper right xyz of plot image real(8) :: xyz_ll(3) ! lower left xyz real(8) :: xyz_ur(3) ! upper right xyz + type(RegularMesh) :: m rgb(:) = pl % meshlines_color % rgb @@ -239,57 +241,56 @@ contains width = xyz_ur_plot - xyz_ll_plot - associate (m => pl % meshlines_mesh) - call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) - call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) + m = meshes(pl % index_meshlines_mesh) + call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) + call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) - ! sweep through all meshbins on this plane and draw borders - do i = ijk_ll(outer), ijk_ur(outer) - do j = ijk_ll(inner), ijk_ur(inner) - ! check if we're in the mesh for this ijk - if (i > 0 .and. i <= m % dimension(outer) .and. & - j > 0 .and. j <= m % dimension(inner)) then + ! sweep through all meshbins on this plane and draw borders + do i = ijk_ll(outer), ijk_ur(outer) + do j = ijk_ll(inner), ijk_ur(inner) + ! check if we're in the mesh for this ijk + if (i > 0 .and. i <= m % dimension(outer) .and. & + j > 0 .and. j <= m % dimension(inner)) then - ! get xyz's of lower left and upper right of this mesh cell - xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) - xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1) - xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i - xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j + ! get xyz's of lower left and upper right of this mesh cell + xyz_ll(outer) = m % lower_left(outer) + m % width(outer) * (i - 1) + xyz_ll(inner) = m % lower_left(inner) + m % width(inner) * (j - 1) + xyz_ur(outer) = m % lower_left(outer) + m % width(outer) * i + xyz_ur(inner) = m % lower_left(inner) + m % width(inner) * j - ! map the xyz ranges to pixel ranges + ! map the xyz ranges to pixel ranges - frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer) - outrange(1) = int(frac * real(pl % pixels(1), 8)) - frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) - outrange(2) = int(frac * real(pl % pixels(1), 8)) + frac = (xyz_ll(outer) - xyz_ll_plot(outer)) / width(outer) + outrange(1) = int(frac * real(pl % pixels(1), 8)) + frac = (xyz_ur(outer) - xyz_ll_plot(outer)) / width(outer) + outrange(2) = int(frac * real(pl % pixels(1), 8)) - frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8)) - frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) - inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8)) + frac = (xyz_ur(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(1) = int((ONE - frac) * real(pl % pixels(2), 8)) + frac = (xyz_ll(inner) - xyz_ll_plot(inner)) / width(inner) + inrange(2) = int((ONE - frac) * real(pl % pixels(2), 8)) - ! draw lines - do out_ = outrange(1), outrange(2) - do plus = 0, pl % meshlines_width - data(:, out_ + 1, inrange(1) + plus + 1) = rgb - data(:, out_ + 1, inrange(2) + plus + 1) = rgb - data(:, out_ + 1, inrange(1) - plus + 1) = rgb - data(:, out_ + 1, inrange(2) - plus + 1) = rgb - end do + ! draw lines + do out_ = outrange(1), outrange(2) + do plus = 0, pl % meshlines_width + data(:, out_ + 1, inrange(1) + plus + 1) = rgb + data(:, out_ + 1, inrange(2) + plus + 1) = rgb + data(:, out_ + 1, inrange(1) - plus + 1) = rgb + data(:, out_ + 1, inrange(2) - plus + 1) = rgb end do - do in_ = inrange(1), inrange(2) - do plus = 0, pl % meshlines_width - data(:, outrange(1) + plus + 1, in_ + 1) = rgb - data(:, outrange(2) + plus + 1, in_ + 1) = rgb - data(:, outrange(1) - plus + 1, in_ + 1) = rgb - data(:, outrange(2) - plus + 1, in_ + 1) = rgb - end do + end do + do in_ = inrange(1), inrange(2) + do plus = 0, pl % meshlines_width + data(:, outrange(1) + plus + 1, in_ + 1) = rgb + data(:, outrange(2) + plus + 1, in_ + 1) = rgb + data(:, outrange(1) - plus + 1, in_ + 1) = rgb + data(:, outrange(2) - plus + 1, in_ + 1) = rgb end do + end do - end if - end do + end if end do - end associate + end do end subroutine draw_mesh_lines diff --git a/src/plot_header.F90 b/src/plot_header.F90 index 1091bd73a8..761481b332 100644 --- a/src/plot_header.F90 +++ b/src/plot_header.F90 @@ -4,7 +4,6 @@ module plot_header use constants use dict_header, only: DictIntInt - use mesh_header, only: RegularMesh implicit none @@ -31,7 +30,7 @@ module plot_header integer :: pixels(3) ! pixel width/height of plot slice integer :: meshlines_width ! pixel width of meshlines integer :: level ! universe depth to plot the cells of - type(RegularMesh), pointer :: meshlines_mesh => null() ! mesh to plot + integer :: index_meshlines_mesh = -1 ! index of mesh to plot type(ObjectColor) :: meshlines_color ! Color for meshlines type(ObjectColor) :: not_found ! color for positions where no cell found type(ObjectColor), allocatable :: colors(:) ! colors of cells/mats From 21a79eefd910ceb2ca3aecb0da4894261b909d9f Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 10:51:05 -0500 Subject: [PATCH 09/21] Remove some uses of mesh_dict --- src/input_xml.F90 | 15 ++++++--------- src/tallies/tally_filter_mesh.F90 | 22 +++++++++++----------- src/tallies/tally_filter_meshsurface.F90 | 22 ++++++++++------------ 3 files changed, 27 insertions(+), 32 deletions(-) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index bfd604f9ee..be12870ad5 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2064,6 +2064,7 @@ contains integer :: i, j integer :: n_cols, col_id, n_comp, n_masks, n_meshlines integer :: meshid + integer(C_INT) :: err, idx integer, allocatable :: iarray(:) logical :: file_exists ! does plots.xml file exist? character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml @@ -2386,7 +2387,7 @@ contains // trim(to_str(pl % id))) end if - pl % meshlines_mesh => meshes(index_ufs_mesh) + pl % index_meshlines_mesh = index_ufs_mesh case ('cmfd') @@ -2404,7 +2405,7 @@ contains // trim(to_str(pl % id))) end if - pl % meshlines_mesh => meshes(index_entropy_mesh) + pl % index_meshlines_mesh = index_entropy_mesh case ('tally') @@ -2417,17 +2418,13 @@ contains end if ! Check if the specified tally mesh exists - if (mesh_dict % has(meshid)) then - pl % meshlines_mesh => meshes(mesh_dict % get(meshid)) - if (meshes(meshid) % type /= MESH_REGULAR) then - call fatal_error("Non-rectangular mesh specified in & - &meshlines for plot " // trim(to_str(pl % id))) - end if - else + err = openmc_get_mesh_index(meshid, idx) + if (err /= 0) then call fatal_error("Could not find mesh " & // trim(to_str(meshid)) // " specified in meshlines for & &plot " // trim(to_str(pl % id))) end if + pl % index_meshlines_mesh = idx case default call fatal_error("Invalid type for meshlines on plot " & diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 486569da63..28eafcaaf5 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -5,7 +5,7 @@ module tally_filter_mesh use constants use dict_header, only: EMPTY use error - use mesh_header, only: RegularMesh, meshes, n_meshes, mesh_dict + use mesh_header use hdf5_interface use particle_header, only: Particle use string, only: to_str @@ -38,10 +38,11 @@ contains class(MeshFilter), intent(inout) :: this type(XMLNode), intent(in) :: node - integer :: i_mesh + integer :: i integer :: id integer :: n - integer :: val + integer(C_INT) :: err + type(RegularMesh) :: m n = node_word_count(node, "bins") @@ -52,19 +53,18 @@ contains call get_node_value(node, "bins", id) ! Get pointer to mesh - val = mesh_dict % get(id) - if (val /= EMPTY) then - i_mesh = val - else + err = openmc_get_mesh_index(id, this % mesh) + if (err /= 0) then call fatal_error("Could not find mesh " // trim(to_str(id)) & // " specified on filter.") end if ! Determine number of bins - this % n_bins = product(meshes(i_mesh) % dimension) - - ! Store the index of the mesh - this % mesh = i_mesh + m = meshes(this % mesh) + this % n_bins = 1 + do i = 1, m % n_dimension() + this % n_bins = this % n_bins * m % dimension(i) + end do end subroutine from_xml subroutine get_all_bins_mesh(this, p, estimator, match) diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index 801d4252ca..c7363a2db2 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -5,7 +5,7 @@ module tally_filter_meshsurface use constants use dict_header, only: EMPTY use error - use mesh_header, only: RegularMesh, meshes, n_meshes, mesh_dict + use mesh_header use hdf5_interface use particle_header, only: Particle use string, only: to_str @@ -38,11 +38,11 @@ contains class(MeshSurfaceFilter), intent(inout) :: this type(XMLNode), intent(in) :: node - integer :: i_mesh + integer :: i integer :: id integer :: n integer :: n_dim - integer :: val + type(RegularMesh) :: m n = node_word_count(node, "bins") @@ -53,20 +53,18 @@ contains call get_node_value(node, "bins", id) ! Get pointer to mesh - val = mesh_dict % get(id) - if (val /= EMPTY) then - i_mesh = val - else + err = openmc_get_mesh_index(id, this % mesh) + if (err /= 0) then call fatal_error("Could not find mesh " // trim(to_str(id)) & // " specified on filter.") end if ! Determine number of bins - n_dim = meshes(i_mesh) % n_dimension - this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension) - - ! Store the index of the mesh - this % mesh = i_mesh + m = meshes(this % mesh) + this % n_bins = 4 * m % n_dimension() + do i = 1, m % n_dimension() + this % n_bins = this % n_bins * m % dimension(i) + end do end subroutine from_xml subroutine get_all_bins(this, p, estimator, match) From 04e48224e6ce8272acd05080fded3326ad774159 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 30 Aug 2018 14:22:55 -0500 Subject: [PATCH 10/21] Almost complete mesh conversion --- CMakeLists.txt | 1 + include/openmc/eigenvalue.h | 4 + include/openmc/hdf5_interface.h | 8 + include/openmc/particle.h | 2 +- src/cmfd_data.F90 | 7 +- src/cmfd_header.F90 | 3 +- src/cmfd_input.F90 | 110 +------- src/eigenvalue.F90 | 1 - src/eigenvalue.cpp | 5 + src/input_xml.F90 | 12 +- src/mesh.cpp | 6 +- src/mesh_header.F90 | 42 ++- src/particle.cpp | 2 +- src/physics.F90 | 1 - src/physics_mg.F90 | 1 - src/plot.F90 | 6 +- src/settings.cpp | 1 + src/simulation.F90 | 5 +- src/simulation_header.F90 | 10 +- src/state_point.F90 | 2 +- src/tallies/tally_filter_mesh.F90 | 277 ++++++++++---------- src/tallies/tally_filter_meshsurface.F90 | 309 ++++++++++++----------- src/tallies/trigger.F90 | 11 +- 23 files changed, 380 insertions(+), 446 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 78d1caf4d8..713bf7c077 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -385,6 +385,7 @@ add_library(libopenmc SHARED src/distribution_energy.cpp src/distribution_multi.cpp src/distribution_spatial.cpp + src/eigenvalue.cpp src/endf.cpp src/initialize.cpp src/finalize.cpp diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index c2e76e26f6..cc879461ed 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -2,6 +2,9 @@ #define OPENMC_EIGENVALUE_H #include // for int64_t +#include + +#include "xtensor/xtensor.hpp" #include "openmc/particle.h" @@ -15,6 +18,7 @@ extern std::vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS extern "C" int64_t n_bank; +#pragma omp threadprivate(n_bank) //============================================================================== // Non-member functions diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 944fe75cb5..b53703fc94 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -317,6 +317,14 @@ write_attribute(hid_t obj_id, const char* name, const std::array& buffer) write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); } +template inline void +write_attribute(hid_t obj_id, const char* name, const std::vector& buffer) +{ + hsize_t dims[] {buffer.size()}; + write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); +} + + //============================================================================== // Templates/overloads for write_dataset //============================================================================== diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 846d2e5a6a..298b2878a0 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -161,7 +161,7 @@ extern "C" { {mark_as_lost(message.str());} //! create a particle restart HDF5 file - void write_restart(); + void write_restart() const; }; diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index ce4825426b..ec07a4d04c 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -80,7 +80,7 @@ contains integer :: i_mesh ! flattend index for mesh logical :: energy_filters! energy filters present real(8) :: flux ! temp variable for flux - type(RegularMesh), pointer :: m ! pointer for mesh object + type(RegularMesh) :: m ! pointer for mesh object ! Extract spatial and energy indices from object nx = cmfd % indices(1) @@ -99,7 +99,7 @@ contains select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m => meshes(filt % mesh) + m = meshes(filt % mesh) end select ! Set mesh widths @@ -354,9 +354,6 @@ contains ! Normalize openmc source distribution cmfd % openmc_src = cmfd % openmc_src/sum(cmfd % openmc_src)*cmfd%norm - ! Nullify all pointers - if (associated(m)) nullify(m) - end subroutine compute_xs !=============================================================================== diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index 2e6162b49e..eafba7b8ad 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -95,7 +95,8 @@ module cmfd_header ! Main object type(cmfd_type), public :: cmfd - type(RegularMesh), public, pointer :: cmfd_mesh => null() + integer, public :: index_cmfd_mesh + type(RegularMesh), public :: cmfd_mesh ! Pointers for different tallies type(TallyContainer), public, pointer :: cmfd_tallies(:) => null() diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index de73bc03f6..3208ad9d06 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -3,7 +3,7 @@ module cmfd_input use, intrinsic :: ISO_C_BINDING use cmfd_header - use mesh_header, only: mesh_dict + use mesh_header use mgxs_interface, only: energy_bins, num_energy_groups use tally use tally_header @@ -241,7 +241,7 @@ contains use constants, only: MAX_LINE_LEN use error, only: fatal_error, warning - use mesh_header, only: RegularMesh, openmc_extend_meshes + use mesh_header use string use tally, only: openmc_tally_allocate use tally_header, only: openmc_extend_tallies @@ -266,112 +266,22 @@ contains integer :: iarray3(3) ! temp integer array real(8) :: rarray3(3) ! temp double array real(C_DOUBLE), allocatable :: energies(:) - type(RegularMesh), pointer :: m type(XMLNode) :: node_mesh - err = openmc_extend_meshes(1, i_start) + ! Read CMFD mesh + call read_meshes(root % ptr) - ! Allocate mesh - cmfd_mesh => meshes(i_start) - m => meshes(i_start) + ! Get index of cmfd mesh and set ID + i_start = n_meshes() - 1 + err = openmc_mesh_set_id(i_start, i_start) - ! Set mesh id - m % id = i_start - - ! Set mesh type to rectangular - m % type = MESH_REGULAR + ! Save reference to CMFD mesh + index_cmfd_mesh = i_start + cmfd_mesh = meshes(i_start) ! Get pointer to mesh XML node node_mesh = root % child("mesh") - ! Determine number of dimensions for mesh - n = node_word_count(node_mesh, "dimension") - if (n /= 2 .and. n /= 3) then - call fatal_error("Mesh must be two or three dimensions.") - end if - m % n_dimension = n - - ! Allocate attribute arrays - allocate(m % dimension(n)) - allocate(m % lower_left(n)) - allocate(m % width(n)) - allocate(m % upper_right(n)) - - ! Check that dimensions are all greater than zero - call get_node_array(node_mesh, "dimension", iarray3(1:n)) - if (any(iarray3(1:n) <= 0)) then - call fatal_error("All entries on the element for a tally mesh& - & must be positive.") - end if - - ! Read dimensions in each direction - m % dimension = iarray3(1:n) - - ! Read mesh lower-left corner location - if (m % n_dimension /= node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same as & - &the number of entries on .") - end if - call get_node_array(node_mesh, "lower_left", m % lower_left) - - ! Make sure both upper-right or width were specified - if (check_for_node(node_mesh, "upper_right") .and. & - check_for_node(node_mesh, "width")) then - call fatal_error("Cannot specify both and on a & - &tally mesh.") - end if - - ! Make sure either upper-right or width was specified - if (.not.check_for_node(node_mesh, "upper_right") .and. & - .not.check_for_node(node_mesh, "width")) then - call fatal_error("Must specify either and on a & - &tally mesh.") - end if - - if (check_for_node(node_mesh, "width")) then - ! Check to ensure width has same dimensions - if (node_word_count(node_mesh, "width") /= & - node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same as the & - &number of entries on .") - end if - - ! Check for negative widths - call get_node_array(node_mesh, "width", rarray3(1:n)) - if (any(rarray3(1:n) < ZERO)) then - call fatal_error("Cannot have a negative on a tally mesh.") - end if - - ! Set width and upper right coordinate - m % width = rarray3(1:n) - m % upper_right = m % lower_left + m % dimension * m % width - - elseif (check_for_node(node_mesh, "upper_right")) then - ! Check to ensure width has same dimensions - if (node_word_count(node_mesh, "upper_right") /= & - node_word_count(node_mesh, "lower_left")) then - call fatal_error("Number of entries on must be the same & - &as the number of entries on .") - end if - - ! Check that upper-right is above lower-left - call get_node_array(node_mesh, "upper_right", rarray3(1:n)) - if (any(rarray3(1:n) < m % lower_left)) then - call fatal_error("The coordinates must be greater than & - &the coordinates on a tally mesh.") - end if - - ! Set upper right coordinate and width - m % upper_right = rarray3(1:n) - m % width = (m % upper_right - m % lower_left) / real(m % dimension, 8) - end if - - ! Set volume fraction - m % volume_frac = ONE/real(product(m % dimension),8) - - ! Add mesh to dictionary - call mesh_dict % set(m % id, i_start) - ! Determine number of filters energy_filters = check_for_node(node_mesh, "energy") n = merge(5, 3, energy_filters) diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 38284c2798..14fe59110e 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -7,7 +7,6 @@ module eigenvalue use error, only: fatal_error, warning use math, only: t_percentile use mesh, only: count_bank_sites - use mesh_header, only: RegularMesh, meshes use message_passing use random_lcg, only: prn, set_particle_seed, advance_prn_seed use settings diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 9ce5b9ef80..ef4159e782 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -147,5 +147,10 @@ extern "C" double entropy_c(int i) return entropy.at(i - 1); } +extern "C" double entropy_clear() +{ + entropy.clear(); +} + } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index be12870ad5..c3fedd881e 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4,7 +4,7 @@ module input_xml use algorithm, only: find use cmfd_input, only: configure_cmfd - use cmfd_header, only: cmfd_mesh + use cmfd_header, only: index_cmfd_mesh use constants use dict_header, only: DictIntInt, DictCharInt, DictEntryCI use endf, only: reaction_name @@ -86,11 +86,6 @@ module input_xml type(C_PTR), value :: node_ptr end subroutine read_materials - subroutine read_meshes(node_ptr) bind(C) - import C_PTR - type(C_PTR), value :: node_ptr - end subroutine - function find_root_universe() bind(C) result(root) import C_INT32_T integer(C_INT32_T) :: root @@ -232,7 +227,6 @@ contains if (run_mode == MODE_EIGENVALUE) then ! Preallocate space for keff and entropy by generation call k_generation % reserve(n_max_batches*gen_per_batch) - call entropy % reserve(n_max_batches*gen_per_batch) end if ! Particle tracks @@ -1222,7 +1216,7 @@ contains ! READ MESH DATA ! Check for user meshes and allocate - call read_meshes() + call read_meshes(root % ptr) ! We only need the mesh info for plotting if (run_mode == MODE_PLOTTING) then @@ -2396,7 +2390,7 @@ contains &meshlines on plot " // trim(to_str(pl % id))) end if - pl % meshlines_mesh => cmfd_mesh + pl % index_meshlines_mesh = index_cmfd_mesh case ('entropy') diff --git a/src/mesh.cpp b/src/mesh.cpp index 9c955fae76..a4526e5b60 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -532,10 +532,10 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, int n_energy, const double* energies, bool* outside) { // Determine shape of array for counts - int m = xt::prod(shape_)(); + std::size_t m = xt::prod(shape_)(); std::vector shape; if (n_energy > 0) { - shape = {m, n_energy}; + shape = {m, static_cast(n_energy)}; } else { shape = {m}; } @@ -757,6 +757,8 @@ void read_meshes(pugi::xml_node* root) //============================================================================== extern "C" { + int n_meshes() { return meshes.size(); } + RegularMesh* mesh_ptr(int i) { return meshes[i-1].get(); } int32_t mesh_id(RegularMesh* m) { return m->id_; } diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 591263000b..2a6abee7d8 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -2,35 +2,20 @@ module mesh_header use, intrinsic :: ISO_C_BINDING - use constants - use dict_header, only: DictIntInt - use error - use hdf5_interface - use string, only: to_str, to_lower - use xml_interface - implicit none - private - public :: openmc_extend_meshes - public :: openmc_get_mesh_index - public :: openmc_mesh_get_id - public :: openmc_mesh_get_dimension - public :: openmc_mesh_get_params - public :: openmc_mesh_set_id - public :: openmc_mesh_set_dimension - public :: openmc_mesh_set_params !=============================================================================== ! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by ! congruent squares or cubes !=============================================================================== - type, public :: RegularMesh + 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 @@ -41,11 +26,6 @@ module mesh_header procedure :: get_indices_from_bin => regular_get_indices_from_bin end type RegularMesh - integer(C_INT32_T), public, bind(C) :: n_meshes = 0 ! # of structured meshes - - ! Dictionary that maps user IDs to indices in 'meshes' - type(DictIntInt), public :: mesh_dict - interface function openmc_extend_meshes(n, index_start, index_end) result(err) bind(C) import C_INT32_T, C_INT @@ -158,21 +138,21 @@ module mesh_header real(C_DOUBLE) :: w end function - function mesh_get_bin(ptr, xyz) result(bin) bind(C) + pure function mesh_get_bin(ptr, xyz) result(bin) bind(C) import C_PTR, C_DOUBLE, C_INT type(C_PTR), value :: ptr real(C_DOUBLE), intent(in) :: xyz(*) integer(C_INT) :: bin end function - function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) + pure function mesh_get_bin_from_indices(ptr, ijk) result(bin) bind(C) import C_PTR, C_INT type(C_PTR), value :: ptr integer(C_INT), intent(in) :: ijk(*) integer(C_INT) :: bin end function - subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) + pure subroutine mesh_get_indices(ptr, xyz, ijk, in_mesh) bind(C) import C_PTR, C_DOUBLE, C_INT, C_BOOL type(C_PTR), value :: ptr real(C_DOUBLE), intent(in) :: xyz(*) @@ -180,7 +160,7 @@ module mesh_header logical(C_BOOL), intent(out) :: in_mesh end subroutine - subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) + pure subroutine mesh_get_indices_from_bin(ptr, bin, ijk) bind(C) import C_PTR, C_INT type(C_PTR), value :: ptr integer(C_INT), value :: bin @@ -192,6 +172,16 @@ module mesh_header integer(C_INT), value :: i type(C_PTR) :: ptr end function + + subroutine read_meshes(node_ptr) bind(C) + import C_PTR + type(C_PTR), value :: node_ptr + end subroutine + + function n_meshes() result(n) bind(C) + import C_INT + integer(C_INT) :: n + end function end interface contains diff --git a/src/particle.cpp b/src/particle.cpp index 729a860eb6..ca57c2faab 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -141,7 +141,7 @@ Particle::mark_as_lost(const char* message) } void -Particle::write_restart() +Particle::write_restart() const { // Dont write another restart file if in particle restart mode if (settings::run_mode == RUN_MODE_PARTICLE) return; diff --git a/src/physics.F90 b/src/physics.F90 index f7d80261ad..f1a0d7c57b 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -6,7 +6,6 @@ module physics use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math - use mesh_header, only: meshes use message_passing use nuclide_header use particle_header diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index 84695153cb..bcafc43142 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -7,7 +7,6 @@ module physics_mg use error, only: fatal_error, warning, write_message use material_header, only: Material, materials use math, only: rotate_angle - use mesh_header, only: meshes use mgxs_interface use message_passing use nuclide_header, only: material_xs diff --git a/src/plot.F90 b/src/plot.F90 index b8f3db0796..bfc8bc2251 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -9,7 +9,7 @@ module plot use hdf5_interface use output, only: time_stamp use material_header, only: materials - use mesh_header, only: meshes + use mesh_header, only: meshes, RegularMesh use particle_header use plot_header use progress_header, only: ProgressBar @@ -242,8 +242,8 @@ contains width = xyz_ur_plot - xyz_ll_plot m = meshes(pl % index_meshlines_mesh) - call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) - call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) + call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension()), in_mesh) + call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension()), in_mesh) ! sweep through all meshbins on this plane and draw borders do i = ijk_ll(outer), ijk_ur(outer) diff --git a/src/settings.cpp b/src/settings.cpp index 4883d07a30..db6aee7f54 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -16,6 +16,7 @@ #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/mesh.h" +#include "openmc/output.h" #include "openmc/random_lcg.h" #include "openmc/source.h" #include "openmc/string_utils.h" diff --git a/src/simulation.F90 b/src/simulation.F90 index 8b7cfdb27a..3ca5d9623a 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -10,8 +10,7 @@ module simulation use cmfd_execute, only: cmfd_init_batch, cmfd_tally_init, execute_cmfd use cmfd_header, only: cmfd_on use constants, only: ZERO - use eigenvalue, only: calculate_average_keff, & - calculate_generation_keff, shannon_entropy, & + use eigenvalue, only: calculate_average_keff, calculate_generation_keff, & synchronize_bank, keff_generation, k_sum #ifdef _OPENMP use eigenvalue, only: join_bank_from_threads @@ -479,7 +478,7 @@ contains ! will potentially populate k_generation and entropy) current_batch = 0 call k_generation % clear() - call entropy % clear() + call entropy_clear() need_depletion_rx = .false. ! If this is a restart run, load the state point data and binary source diff --git a/src/simulation_header.F90 b/src/simulation_header.F90 index 9b7c68cd46..ad5513c618 100644 --- a/src/simulation_header.F90 +++ b/src/simulation_header.F90 @@ -64,6 +64,12 @@ module simulation_header !$omp threadprivate(trace, thread_id, current_work) + interface + subroutine entropy_clear() bind(C) + end subroutine + end interface + + contains !=============================================================================== @@ -80,12 +86,12 @@ contains !=============================================================================== subroutine free_memory_simulation() + if (allocated(work_index)) deallocate(work_index) call k_generation % clear() call k_generation % shrink_to_fit() - call entropy % clear() - call entropy % shrink_to_fit() + call entropy_clear() end subroutine free_memory_simulation end module simulation_header diff --git a/src/state_point.F90 b/src/state_point.F90 index b7f4e13776..d20456d7b5 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -173,7 +173,7 @@ contains call write_dataset(file_id, "generations_per_batch", gen_per_batch) k = k_generation % size() call write_dataset(file_id, "k_generation", k_generation % data(1:k)) - call entropy_to_hdf5() + call entropy_to_hdf5(file_id) call write_dataset(file_id, "k_col_abs", k_col_abs) call write_dataset(file_id, "k_col_tra", k_col_tra) call write_dataset(file_id, "k_abs_tra", k_abs_tra) diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 28eafcaaf5..3e8eaa1c6d 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -92,155 +92,154 @@ contains real(8) :: distance ! distance traveled in mesh cell logical :: start_in_mesh ! starting coordinates inside mesh? logical :: end_in_mesh ! ending coordinates inside mesh? - type(RegularMesh), pointer :: m + type(RegularMesh) :: m ! Get a pointer to the mesh. - m => meshes(this % mesh) - n = m % n_dimension + m = meshes(this % mesh) if (estimator /= ESTIMATOR_TRACKLENGTH) then ! If this is an analog or collision tally, then there can only be one ! valid mesh bin. call m % get_bin(p % coord(1) % xyz, bin) - if (bin /= NO_BIN_FOUND) then + if (bin >= 0) then call match % bins % push_back(bin) call match % weights % push_back(ONE) end if return end if - ! A track can span multiple mesh bins so we need to handle a lot of - ! intersection logic for tracklength tallies. + ! ! A track can span multiple mesh bins so we need to handle a lot of + ! ! intersection logic for tracklength tallies. - ! ======================================================================== - ! Determine if the track intersects the tally mesh. + ! ! ======================================================================== + ! ! Determine if the track intersects the tally mesh. - ! Copy the starting and ending coordinates of the particle. Offset these - ! just a bit for the purposes of determining if there was an intersection - ! in case the mesh surfaces coincide with lattice/geometric surfaces which - ! might produce finite-precision errors. - xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw - xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw + ! ! Copy the starting and ending coordinates of the particle. Offset these + ! ! just a bit for the purposes of determining if there was an intersection + ! ! in case the mesh surfaces coincide with lattice/geometric surfaces which + ! ! might produce finite-precision errors. + ! xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw + ! xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - ! Determine indices for starting and ending location. - call m % get_indices(xyz0, ijk0(:n), start_in_mesh) - call m % get_indices(xyz1, ijk1(:n), end_in_mesh) + ! ! Determine indices for starting and ending location. + ! call m % get_indices(xyz0, ijk0(:n), start_in_mesh) + ! call m % get_indices(xyz1, ijk1(:n), end_in_mesh) - ! If this is the first iteration of the filter loop, check if the track - ! intersects any part of the mesh. - if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - if (.not. m % intersects(xyz0, xyz1)) return - end if + ! ! If this is the first iteration of the filter loop, check if the track + ! ! intersects any part of the mesh. + ! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + ! if (.not. m % intersects(xyz0, xyz1)) return + ! end if - ! ======================================================================== - ! Figure out which mesh cell to tally. + ! ! ======================================================================== + ! ! Figure out which mesh cell to tally. - ! Copy the un-modified coordinates the particle direction. - xyz0 = p % last_xyz - xyz1 = p % coord(1) % xyz - uvw = p % coord(1) % uvw + ! ! Copy the un-modified coordinates the particle direction. + ! xyz0 = p % last_xyz + ! xyz1 = p % coord(1) % xyz + ! uvw = p % coord(1) % uvw - ! Compute the length of the entire track. - total_distance = sqrt(sum((xyz1 - xyz0)**2)) + ! ! Compute the length of the entire track. + ! total_distance = sqrt(sum((xyz1 - xyz0)**2)) - ! We are looking for the first valid mesh bin. Check to see if the - ! particle starts inside the mesh. - if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then - ! The particle does not start in the mesh. Note that we nudged the - ! start and end coordinates by a TINY_BIT each so we will have - ! difficulty resolving tracks that are less than 2*TINY_BIT in length. - ! If the track is that short, it is also insignificant so we can - ! safely ignore it in the tallies. - if (total_distance < 2*TINY_BIT) return + ! ! We are looking for the first valid mesh bin. Check to see if the + ! ! particle starts inside the mesh. + ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then + ! ! The particle does not start in the mesh. Note that we nudged the + ! ! start and end coordinates by a TINY_BIT each so we will have + ! ! difficulty resolving tracks that are less than 2*TINY_BIT in length. + ! ! If the track is that short, it is also insignificant so we can + ! ! safely ignore it in the tallies. + ! if (total_distance < 2*TINY_BIT) return - ! The particle does not start in the mesh so keep iterating the ijk0 - ! indices to cross the nearest mesh surface until we've found a valid - ! bin. MAX_SEARCH_ITER prevents an infinite loop. - search_iter = 0 - do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) - if (search_iter == MAX_SEARCH_ITER) then - call warning("Failed to find a mesh intersection on a tally mesh & - &filter.") - return - end if + ! ! The particle does not start in the mesh so keep iterating the ijk0 + ! ! indices to cross the nearest mesh surface until we've found a valid + ! ! bin. MAX_SEARCH_ITER prevents an infinite loop. + ! search_iter = 0 + ! do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) + ! if (search_iter == MAX_SEARCH_ITER) then + ! call warning("Failed to find a mesh intersection on a tally mesh & + ! &filter.") + ! return + ! end if - do j = 1, n - if (abs(uvw(j)) < FP_PRECISION) then - d(j) = INFINITY - else if (uvw(j) > 0) then - xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - else - xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - end if - end do - j = minloc(d(:n), 1) - if (uvw(j) > ZERO) then - ijk0(j) = ijk0(j) + 1 - else - ijk0(j) = ijk0(j) - 1 - end if + ! do j = 1, n + ! if (abs(uvw(j)) < FP_PRECISION) then + ! d(j) = INFINITY + ! else if (uvw(j) > 0) then + ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! else + ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! end if + ! end do + ! j = minloc(d(:n), 1) + ! if (uvw(j) > ZERO) then + ! ijk0(j) = ijk0(j) + 1 + ! else + ! ijk0(j) = ijk0(j) - 1 + ! end if - search_iter = search_iter + 1 - end do - distance = d(j) - xyz0 = xyz0 + distance * uvw - end if + ! search_iter = search_iter + 1 + ! end do + ! distance = d(j) + ! xyz0 = xyz0 + distance * uvw + ! end if - do - ! ======================================================================== - ! Compute the length of the track segment in the appropiate mesh cell and - ! return. + ! do + ! ! ======================================================================== + ! ! Compute the length of the track segment in the appropiate mesh cell and + ! ! return. - if (all(ijk0(:n) == ijk1(:n))) then - ! The track ends in this cell. Use the particle end location rather - ! than the mesh surface. - distance = sqrt(sum((xyz1 - xyz0)**2)) - else - ! The track exits this cell. Determine the distance to the closest mesh - ! surface. - do j = 1, n - if (abs(uvw(j)) < FP_PRECISION) then - d(j) = INFINITY - else if (uvw(j) > 0) then - xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - else - xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - d(j) = (xyz_cross - xyz0(j)) / uvw(j) - end if - end do - j = minloc(d(:n), 1) - distance = d(j) - end if + ! if (all(ijk0(:n) == ijk1(:n))) then + ! ! The track ends in this cell. Use the particle end location rather + ! ! than the mesh surface. + ! distance = sqrt(sum((xyz1 - xyz0)**2)) + ! else + ! ! The track exits this cell. Determine the distance to the closest mesh + ! ! surface. + ! do j = 1, n + ! if (abs(uvw(j)) < FP_PRECISION) then + ! d(j) = INFINITY + ! else if (uvw(j) > 0) then + ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! else + ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) + ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) + ! end if + ! end do + ! j = minloc(d(:n), 1) + ! distance = d(j) + ! end if - ! Assign the next tally bin and the score. - bin = m % get_bin_from_indices(ijk0(:n)) - call match % bins % push_back(bin) - call match % weights % push_back(distance / total_distance) + ! ! Assign the next tally bin and the score. + ! bin = m % get_bin_from_indices(ijk0(:n)) + ! call match % bins % push_back(bin) + ! call match % weights % push_back(distance / total_distance) - ! Find the next mesh cell that the particle enters. + ! ! Find the next mesh cell that the particle enters. - ! If the particle track ends in that bin, then we are done. - if (all(ijk0(:n) == ijk1(:n))) exit + ! ! If the particle track ends in that bin, then we are done. + ! if (all(ijk0(:n) == ijk1(:n))) exit - ! Translate the starting coordintes by the distance to that face. This - ! should be the xyz that we computed the distance to in the last - ! iteration of the filter loop. - xyz0 = xyz0 + distance * uvw + ! ! Translate the starting coordintes by the distance to that face. This + ! ! should be the xyz that we computed the distance to in the last + ! ! iteration of the filter loop. + ! xyz0 = xyz0 + distance * uvw - ! Increment the indices into the next mesh cell. - if (uvw(j) > ZERO) then - ijk0(j) = ijk0(j) + 1 - else - ijk0(j) = ijk0(j) - 1 - end if + ! ! Increment the indices into the next mesh cell. + ! if (uvw(j) > ZERO) then + ! ijk0(j) = ijk0(j) + 1 + ! else + ! ijk0(j) = ijk0(j) - 1 + ! end if - ! If the next indices are invalid, then the track has left the mesh and - ! we are done. - if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit - end do + ! ! If the next indices are invalid, then the track has left the mesh and + ! ! we are done. + ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit + ! end do end subroutine get_all_bins_mesh @@ -248,9 +247,12 @@ contains class(MeshFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group + type(RegularMesh) :: m + + m = meshes(this % mesh) call write_dataset(filter_group, "type", "mesh") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", meshes(this % mesh) % id) + call write_dataset(filter_group, "bins", m % id()) end subroutine to_statepoint_mesh function text_label_mesh(this, bin) result(label) @@ -259,20 +261,20 @@ contains character(MAX_LINE_LEN) :: label integer, allocatable :: ijk(:) + type(RegularMesh) :: m - associate (m => meshes(this % mesh)) - allocate(ijk(m % n_dimension)) - call m % get_indices_from_bin(bin, ijk) - if (m % n_dimension == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if - end associate + m = meshes(this % mesh) + allocate(ijk(m % n_dimension())) + call m % get_indices_from_bin(bin, ijk) + if (m % n_dimension() == 1) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" + elseif (m % n_dimension() == 2) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" + elseif (m % n_dimension() == 3) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + end if end function text_label_mesh !=============================================================================== @@ -304,18 +306,25 @@ contains integer(C_INT32_T), value, intent(in) :: index_mesh integer(C_INT) :: err + type(RegularMesh) :: m + integer :: i + err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MeshFilter) - if (index_mesh >= 1 .and. index_mesh <= n_meshes) then + if (index_mesh >= 0 .and. index_mesh < n_meshes()) then f % mesh = index_mesh - f % n_bins = product(meshes(index_mesh) % dimension) + f % n_bins = 1 + m = meshes(index_mesh) + do i = 1, m % n_dimension() + f % n_bins = f % n_bins * m % dimension(i) + end do else err = E_OUT_OF_BOUNDS call set_errmsg("Index in 'meshes' array is out of bounds.") end if - class default + class default err = E_INVALID_TYPE call set_errmsg("Tried to set mesh on a non-mesh filter.") end select diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index c7363a2db2..d1500b8b39 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -42,6 +42,7 @@ contains integer :: id integer :: n integer :: n_dim + integer(C_INT) :: err type(RegularMesh) :: m n = node_word_count(node, "bins") @@ -90,133 +91,130 @@ contains real(8) :: distance ! actual distance traveled logical :: start_in_mesh ! particle's starting xyz in mesh? logical :: end_in_mesh ! particle's ending xyz in mesh? + type(RegularMesh) :: m - ! Copy starting and ending location of particle - xyz0 = p % last_xyz_current - xyz1 = p % coord(1) % xyz + ! ! Copy starting and ending location of particle + ! xyz0 = p % last_xyz_current + ! xyz1 = p % coord(1) % xyz - associate (m => meshes(this % mesh)) - n_dim = m % n_dimension + ! m = meshes(this % mesh) + ! n_dim = m % n_dimension() - ! Determine indices for starting and ending location - call m % get_indices(xyz0, ijk0, start_in_mesh) - call m % get_indices(xyz1, ijk1, end_in_mesh) + ! ! Determine indices for starting and ending location + ! call m % get_indices(xyz0, ijk0, start_in_mesh) + ! call m % get_indices(xyz1, ijk1, end_in_mesh) - ! Check to see if start or end is in mesh -- if not, check if track still - ! intersects with mesh - if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - if (.not. m % intersects(xyz0, xyz1)) return - end if + ! ! Check to see if start or end is in mesh -- if not, check if track still + ! ! intersects with mesh + ! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then + ! !if (.not. m % intersects(xyz0, xyz1)) return + ! end if - ! Calculate number of surface crossings - n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) - if (n_cross == 0) return + ! ! Calculate number of surface crossings + ! n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) + ! if (n_cross == 0) return - ! Copy particle's direction - uvw = p % coord(1) % uvw + ! ! Copy particle's direction + ! uvw = p % coord(1) % uvw - ! Bounding coordinates - do d1 = 1, n_dim - if (uvw(d1) > 0) then - xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) - else - xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) - end if - end do + ! ! Bounding coordinates + ! do d1 = 1, n_dim + ! if (uvw(d1) > 0) then + ! xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) + ! else + ! xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) + ! end if + ! end do - do j = 1, n_cross - ! Set the distances to infinity - d = INFINITY + ! do j = 1, n_cross + ! ! Set the distances to infinity + ! d = INFINITY - ! Calculate distance to each bounding surface. We need to treat - ! special case where the cosine of the angle is zero since this would - ! result in a divide-by-zero. - do d1 = 1, n_dim - if (uvw(d1) == 0) then - d(d1) = INFINITY - else - d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) - end if - end do + ! ! Calculate distance to each bounding surface. We need to treat + ! ! special case where the cosine of the angle is zero since this would + ! ! result in a divide-by-zero. + ! do d1 = 1, n_dim + ! if (uvw(d1) == 0) then + ! d(d1) = INFINITY + ! else + ! d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) + ! end if + ! end do - ! Determine the closest bounding surface of the mesh cell by - ! calculating the minimum distance. Then use the minimum distance and - ! direction of the particle to determine which surface was crossed. - distance = minval(d) + ! ! Determine the closest bounding surface of the mesh cell by + ! ! calculating the minimum distance. Then use the minimum distance and + ! ! direction of the particle to determine which surface was crossed. + ! distance = minval(d) - ! Loop over the dimensions - do d1 = 1, n_dim + ! ! Loop over the dimensions + ! do d1 = 1, n_dim - ! Check whether distance is the shortest distance - if (distance == d(d1)) then + ! ! Check whether distance is the shortest distance + ! if (distance == d(d1)) then - ! Check whether particle is moving in positive d1 direction - if (uvw(d1) > 0) then + ! ! Check whether particle is moving in positive d1 direction + ! if (uvw(d1) > 0) then - ! Outward current on d1 max surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 1 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! Outward current on d1 max surface + ! if (all(ijk0(:n_dim) >= 1) .and. & + ! all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 1 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - ! Advance position - ijk0(d1) = ijk0(d1) + 1 - xyz_cross(d1) = xyz_cross(d1) + m % width(d1) + ! ! Advance position + ! ijk0(d1) = ijk0(d1) + 1 + ! xyz_cross(d1) = xyz_cross(d1) + m % width(d1) - ! If the particle crossed the surface, tally the inward current on - ! d1 min surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 2 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! If the particle crossed the surface, tally the inward current on + ! ! d1 min surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 2 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - else - ! The particle is moving in the negative d1 direction + ! else + ! ! The particle is moving in the negative d1 direction - ! Outward current on d1 min surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - 3 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! Outward current on d1 min surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 - 3 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if - ! Advance position - ijk0(d1) = ijk0(d1) - 1 - xyz_cross(d1) = xyz_cross(d1) - m % width(d1) + ! ! Advance position + ! ijk0(d1) = ijk0(d1) - 1 + ! xyz_cross(d1) = xyz_cross(d1) - m % width(d1) - ! If the particle crossed the surface, tally the inward current on - ! d1 max surface - if (all(ijk0(:n_dim) >= 1) .and. & - all(ijk0(:n_dim) <= m % dimension)) then - i_surf = d1 * 4 - i_mesh = m % get_bin_from_indices(ijk0) - i_bin = 4*n_dim*(i_mesh - 1) + i_surf + ! ! If the particle crossed the surface, tally the inward current on + ! ! d1 max surface + ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then + ! i_surf = d1 * 4 + ! i_mesh = m % get_bin_from_indices(ijk0) + ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - call match % bins % push_back(i_bin) - call match % weights % push_back(ONE) - end if - end if - end if - end do + ! call match % bins % push_back(i_bin) + ! call match % weights % push_back(ONE) + ! end if + ! end if + ! end if + ! end do - ! Calculate new coordinates - xyz0 = xyz0 + distance * uvw - end do - end associate + ! ! Calculate new coordinates + ! xyz0 = xyz0 + distance * uvw + ! end do end subroutine get_all_bins @@ -224,9 +222,12 @@ contains class(MeshSurfaceFilter), intent(in) :: this integer(HID_T), intent(in) :: filter_group + type(RegularMesh) :: m + + m = meshes(this % mesh) call write_dataset(filter_group, "type", "meshsurface") call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", meshes(this % mesh) % id) + call write_dataset(filter_group, "bins", m % id()) end subroutine to_statepoint function text_label(this, bin) result(label) @@ -238,55 +239,55 @@ contains integer :: i_surf integer :: n_dim integer, allocatable :: ijk(:) + type(RegularMesh) :: m - associate (m => meshes(this % mesh)) - n_dim = m % n_dimension - allocate(ijk(n_dim)) + m = meshes(this % mesh) + n_dim = m % n_dimension() + allocate(ijk(n_dim)) - ! Get flattend mesh index and surface index - i_mesh = (bin - 1) / (4*n_dim) + 1 - i_surf = mod(bin - 1, 4*n_dim) + 1 + ! Get flattend mesh index and surface index + i_mesh = (bin - 1) / (4*n_dim) + 1 + i_surf = mod(bin - 1, 4*n_dim) + 1 - ! Get mesh index part of label - call m % get_indices_from_bin(i_mesh, ijk) - if (m % n_dimension == 1) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" - elseif (m % n_dimension == 2) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ")" - elseif (m % n_dimension == 3) then - label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & - trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" - end if + ! Get mesh index part of label + call m % get_indices_from_bin(i_mesh, ijk) + if (m % n_dimension() == 1) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" + elseif (m % n_dimension() == 2) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ")" + elseif (m % n_dimension() == 3) then + label = "Mesh Index (" // trim(to_str(ijk(1))) // ", " // & + trim(to_str(ijk(2))) // ", " // trim(to_str(ijk(3))) // ")" + end if - ! Get surface part of label - select case (i_surf) - case (OUT_LEFT) - label = trim(label) // " Outgoing, x-min" - case (IN_LEFT) - label = trim(label) // " Incoming, x-min" - case (OUT_RIGHT) - label = trim(label) // " Outgoing, x-max" - case (IN_RIGHT) - label = trim(label) // " Incoming, x-max" - case (OUT_BACK) - label = trim(label) // " Outgoing, y-min" - case (IN_BACK) - label = trim(label) // " Incoming, y-min" - case (OUT_FRONT) - label = trim(label) // " Outgoing, y-max" - case (IN_FRONT) - label = trim(label) // " Incoming, y-max" - case (OUT_BOTTOM) - label = trim(label) // " Outgoing, z-min" - case (IN_BOTTOM) - label = trim(label) // " Incoming, z-min" - case (OUT_TOP) - label = trim(label) // " Outgoing, z-max" - case (IN_TOP) - label = trim(label) // " Incoming, z-max" - end select - end associate + ! Get surface part of label + select case (i_surf) + case (OUT_LEFT) + label = trim(label) // " Outgoing, x-min" + case (IN_LEFT) + label = trim(label) // " Incoming, x-min" + case (OUT_RIGHT) + label = trim(label) // " Outgoing, x-max" + case (IN_RIGHT) + label = trim(label) // " Incoming, x-max" + case (OUT_BACK) + label = trim(label) // " Outgoing, y-min" + case (IN_BACK) + label = trim(label) // " Incoming, y-min" + case (OUT_FRONT) + label = trim(label) // " Outgoing, y-max" + case (IN_FRONT) + label = trim(label) // " Incoming, y-max" + case (OUT_BOTTOM) + label = trim(label) // " Outgoing, z-min" + case (IN_BOTTOM) + label = trim(label) // " Incoming, z-min" + case (OUT_TOP) + label = trim(label) // " Outgoing, z-max" + case (IN_TOP) + label = trim(label) // " Incoming, z-max" + end select end function text_label !=============================================================================== @@ -318,16 +319,22 @@ contains integer(C_INT32_T), value, intent(in) :: index_mesh integer(C_INT) :: err + integer :: i integer :: n_dim + type(RegularMesh) :: m err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MeshSurfaceFilter) - if (index_mesh >= 1 .and. index_mesh <= n_meshes) then + if (index_mesh >= 0 .and. index_mesh < n_meshes()) then f % mesh = index_mesh - n_dim = meshes(index_mesh) % n_dimension - f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension) + m = meshes(index_mesh) + n_dim = m % n_dimension() + f % n_bins = 4*n_dim + do i = 1, n_dim + f % n_bins = f % n_bins * m % dimension(i) + end do else err = E_OUT_OF_BOUNDS call set_errmsg("Index in 'meshes' array is out of bounds.") diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 2cb3bf8194..27aa5ec50e 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -257,14 +257,14 @@ contains logical :: print_ebin ! should incoming energy bin be displayed? real(8) :: rel_err = ZERO ! temporary relative error of result real(8) :: std_dev = ZERO ! temporary standard deviration of result - type(RegularMesh), pointer :: m ! surface current mesh + type(RegularMesh) :: m ! surface current mesh ! Get pointer to mesh i_filter_mesh = t % filter(t % find_filter(FILTER_MESH)) i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE)) select type(filt => filters(i_filter_mesh) % obj) type is (MeshFilter) - m => meshes(filt % mesh) + m = meshes(filt % mesh) end select ! initialize bins array @@ -285,8 +285,11 @@ contains end if ! Get the dimensions and number of cells in the mesh - n_dim = m % n_dimension - n_cells = product(m % dimension) + n_dim = m % n_dimension() + n_cells = 1 + do j = 1, n_dim + n_cells = n_cells * m % dimension(j) + end do ! Loop over all the mesh cells do i = 1, n_cells From a10ef739baf5e9a5d8ad390e9250d3b85f01b962 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 3 Sep 2018 14:01:59 -0500 Subject: [PATCH 11/21] Bindings for bins_crossed --- include/openmc/mesh.h | 7 ++ src/input_xml.F90 | 8 +- src/mesh.cpp | 51 +++++++--- src/mesh_header.F90 | 2 +- src/output.F90 | 9 ++ src/simulation.F90 | 8 ++ src/stl_vector.F90 | 26 +++++ src/tallies/tally_filter_header.F90 | 4 +- src/tallies/tally_filter_mesh.F90 | 148 +++------------------------- 9 files changed, 107 insertions(+), 156 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 28f9c3eac9..47258edd82 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -1,3 +1,6 @@ +//! \file mesh.h +//! \brief Mesh types used for tallies, Shannon entropy, CMFD, etc. + #ifndef OPENMC_MESH_H #define OPENMC_MESH_H @@ -59,6 +62,10 @@ private: //! \param[in] root XML node extern "C" void read_meshes(pugi::xml_node* root); +//! Write mesh data to an HDF5 group +//! \param[in] group HDF5 group +extern "C" void meshes_to_hdf5(hid_t group); + //============================================================================== // Global variables //============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c3fedd881e..4766650bf9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -65,17 +65,17 @@ module input_xml subroutine read_surfaces(node_ptr) bind(C) import C_PTR - type(C_PTR), value :: node_ptr + type(C_PTR) :: node_ptr end subroutine read_surfaces subroutine read_cells(node_ptr) bind(C) import C_PTR - type(C_PTR), value :: node_ptr + type(C_PTR) :: node_ptr end subroutine read_cells subroutine read_lattices(node_ptr) bind(C) import C_PTR - type(C_PTR), value :: node_ptr + type(C_PTR) :: node_ptr end subroutine read_lattices subroutine read_settings_xml() bind(C) @@ -83,7 +83,7 @@ module input_xml subroutine read_materials(node_ptr) bind(C) import C_PTR - type(C_PTR), value :: node_ptr + type(C_PTR) :: node_ptr end subroutine read_materials function find_root_universe() bind(C) result(root) diff --git a/src/mesh.cpp b/src/mesh.cpp index a4526e5b60..3749c55859 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -752,14 +752,38 @@ void read_meshes(pugi::xml_node* root) } } +void meshes_to_hdf5(hid_t group) +{ + // Write number of meshes + hid_t meshes_group = create_group(group, "meshes"); + int32_t n_meshes = meshes.size(); + write_attribute(meshes_group, "n_meshes", n_meshes); + + if (n_meshes > 0) { + // Write IDs of meshes + std::vector ids; + for (const auto& m : meshes) { + m->to_hdf5(meshes_group); + ids.push_back(m->id_); + } + write_attribute(meshes_group, "ids", ids); + } + + close_group(meshes_group); +} + //============================================================================== // Fortran compatibility //============================================================================== extern "C" { + // Declaration of Fortran procedures + void vector_int_push_back(void* ptr, int value); + void vector_real_push_back(void* ptr, double value); + int n_meshes() { return meshes.size(); } - RegularMesh* mesh_ptr(int i) { return meshes[i-1].get(); } + RegularMesh* mesh_ptr(int i) { return meshes.at(i).get(); } int32_t mesh_id(RegularMesh* m) { return m->id_; } @@ -795,24 +819,19 @@ extern "C" { m->get_indices_from_bin(bin, ijk); } - void meshes_to_hdf5(hid_t group) + void mesh_bins_crossed(RegularMesh* m, const Particle* p, void* match_bins, + void* match_weights) { - // Write number of meshes - hid_t meshes_group = create_group(group, "meshes"); - int32_t n_meshes = meshes.size(); - write_attribute(meshes_group, "n_meshes", n_meshes); + // Get bins crossed + std::vector bins; + std::vector lengths; + m->bins_crossed(p, bins, lengths); - if (n_meshes > 0) { - // Write IDs of meshes - std::vector ids; - for (const auto& m : meshes) { - m->to_hdf5(meshes_group); - ids.push_back(m->id_); - } - write_attribute(meshes_group, "ids", ids); + // Call bindings for VectorInt and VectorReal on Fortran side + for (int i = 0; i < bins.size(); ++i) { + vector_int_push_back(match_bins, bins[i]); + vector_real_push_back(match_weights, lengths[i]); } - - close_group(meshes_group); } void free_memory_mesh() diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 2a6abee7d8..d7e26fc6b8 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -175,7 +175,7 @@ module mesh_header subroutine read_meshes(node_ptr) bind(C) import C_PTR - type(C_PTR), value :: node_ptr + type(C_PTR) :: node_ptr end subroutine function n_meshes() result(n) bind(C) diff --git a/src/output.F90 b/src/output.F90 index e286956670..4363731f09 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -660,6 +660,10 @@ contains if (n_tallies == 0) return allocate(matches(n_filters)) + do i = 1, n_filters + allocate(matches(i) % bins) + allocate(matches(i) % weights) + end do ! Initialize names for scores score_names(abs(SCORE_FLUX)) = "Flux" @@ -855,6 +859,11 @@ contains close(UNIT=unit_tally) + do i = 1, n_filters + deallocate(matches(i) % bins) + deallocate(matches(i) % weights) + end do + end subroutine write_tallies !=============================================================================== diff --git a/src/simulation.F90 b/src/simulation.F90 index 3ca5d9623a..176dfb921b 100644 --- a/src/simulation.F90 +++ b/src/simulation.F90 @@ -472,6 +472,10 @@ contains ! Allocate array for matching filter bins allocate(filter_matches(n_filters)) + do i = 1, n_filters + allocate(filter_matches(i) % bins) + allocate(filter_matches(i) % weights) + end do !$omp end parallel ! Reset global variables -- this is done before loading state point (as that @@ -557,6 +561,10 @@ contains deallocate(materials(i) % mat_nuclide_index) end do !$omp parallel + do i = 1, size(filter_matches) + deallocate(filter_matches(i) % bins) + deallocate(filter_matches(i) % weights) + end do deallocate(micro_xs, micro_photon_xs, filter_matches) !$omp end parallel diff --git a/src/stl_vector.F90 b/src/stl_vector.F90 index 06f487dc1e..1562907814 100644 --- a/src/stl_vector.F90 +++ b/src/stl_vector.F90 @@ -35,6 +35,8 @@ module stl_vector ! ! size -- Returns the number of elements in the vector. + use, intrinsic :: ISO_C_BINDING + implicit none private @@ -521,4 +523,28 @@ contains size = this%size_ end function size_char +!=============================================================================== +! Procedures to be called from C++ +!=============================================================================== + + subroutine vector_int_push_back(ptr, val) bind(C) + type(C_PTR), value :: ptr + integer(C_INT), value :: val + + type(VectorInt), pointer :: vec + + call C_F_POINTER(ptr, vec) + call vec % push_back(val) + end subroutine + + subroutine vector_real_push_back(ptr, val) bind(C) + type(C_PTR), value :: ptr + real(C_DOUBLE), value :: val + + type(VectorReal), pointer :: vec + + call C_F_POINTER(ptr, vec) + call vec % push_back(val) + end subroutine + end module stl_vector diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index e57423ddfa..d6aef4abd9 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -28,8 +28,8 @@ module tally_filter_header type, public :: TallyFilterMatch ! Index of the bin and weight being used in the current filter combination integer :: i_bin - type(VectorInt) :: bins - type(VectorReal) :: weights + type(VectorInt), pointer :: bins + type(VectorReal), pointer :: weights ! Indicates whether all valid bins for this filter have been found logical :: bins_present = .false. diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 3e8eaa1c6d..bc770b199a 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -93,6 +93,17 @@ contains logical :: start_in_mesh ! starting coordinates inside mesh? logical :: end_in_mesh ! ending coordinates inside mesh? type(RegularMesh) :: m + type(C_PTR) :: ptr_bins, ptr_weights + + interface + subroutine mesh_bins_crossed(m, p, bins, weights) bind(C) + import C_PTR, Particle + type(C_PTR), value :: m + type(Particle), intent(in) :: p + type(C_PTR), value :: bins + type(C_PTR), value :: weights + end subroutine + end interface ! Get a pointer to the mesh. m = meshes(this % mesh) @@ -106,141 +117,12 @@ contains call match % weights % push_back(ONE) end if return + else + ptr_bins = C_LOC(match % bins) + ptr_weights = C_LOC(match % weights) + call mesh_bins_crossed(m % ptr, p, ptr_bins, ptr_weights) end if - ! ! A track can span multiple mesh bins so we need to handle a lot of - ! ! intersection logic for tracklength tallies. - - ! ! ======================================================================== - ! ! Determine if the track intersects the tally mesh. - - ! ! Copy the starting and ending coordinates of the particle. Offset these - ! ! just a bit for the purposes of determining if there was an intersection - ! ! in case the mesh surfaces coincide with lattice/geometric surfaces which - ! ! might produce finite-precision errors. - ! xyz0 = p % last_xyz + TINY_BIT * p % coord(1) % uvw - ! xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw - - ! ! Determine indices for starting and ending location. - ! call m % get_indices(xyz0, ijk0(:n), start_in_mesh) - ! call m % get_indices(xyz1, ijk1(:n), end_in_mesh) - - ! ! If this is the first iteration of the filter loop, check if the track - ! ! intersects any part of the mesh. - ! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - ! if (.not. m % intersects(xyz0, xyz1)) return - ! end if - - ! ! ======================================================================== - ! ! Figure out which mesh cell to tally. - - ! ! Copy the un-modified coordinates the particle direction. - ! xyz0 = p % last_xyz - ! xyz1 = p % coord(1) % xyz - ! uvw = p % coord(1) % uvw - - ! ! Compute the length of the entire track. - ! total_distance = sqrt(sum((xyz1 - xyz0)**2)) - - ! ! We are looking for the first valid mesh bin. Check to see if the - ! ! particle starts inside the mesh. - ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) then - ! ! The particle does not start in the mesh. Note that we nudged the - ! ! start and end coordinates by a TINY_BIT each so we will have - ! ! difficulty resolving tracks that are less than 2*TINY_BIT in length. - ! ! If the track is that short, it is also insignificant so we can - ! ! safely ignore it in the tallies. - ! if (total_distance < 2*TINY_BIT) return - - ! ! The particle does not start in the mesh so keep iterating the ijk0 - ! ! indices to cross the nearest mesh surface until we've found a valid - ! ! bin. MAX_SEARCH_ITER prevents an infinite loop. - ! search_iter = 0 - ! do while (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) - ! if (search_iter == MAX_SEARCH_ITER) then - ! call warning("Failed to find a mesh intersection on a tally mesh & - ! &filter.") - ! return - ! end if - - ! do j = 1, n - ! if (abs(uvw(j)) < FP_PRECISION) then - ! d(j) = INFINITY - ! else if (uvw(j) > 0) then - ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) - ! else - ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) - ! end if - ! end do - ! j = minloc(d(:n), 1) - ! if (uvw(j) > ZERO) then - ! ijk0(j) = ijk0(j) + 1 - ! else - ! ijk0(j) = ijk0(j) - 1 - ! end if - - ! search_iter = search_iter + 1 - ! end do - ! distance = d(j) - ! xyz0 = xyz0 + distance * uvw - ! end if - - ! do - ! ! ======================================================================== - ! ! Compute the length of the track segment in the appropiate mesh cell and - ! ! return. - - ! if (all(ijk0(:n) == ijk1(:n))) then - ! ! The track ends in this cell. Use the particle end location rather - ! ! than the mesh surface. - ! distance = sqrt(sum((xyz1 - xyz0)**2)) - ! else - ! ! The track exits this cell. Determine the distance to the closest mesh - ! ! surface. - ! do j = 1, n - ! if (abs(uvw(j)) < FP_PRECISION) then - ! d(j) = INFINITY - ! else if (uvw(j) > 0) then - ! xyz_cross = m % lower_left(j) + ijk0(j) * m % width(j) - ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) - ! else - ! xyz_cross = m % lower_left(j) + (ijk0(j) - 1) * m % width(j) - ! d(j) = (xyz_cross - xyz0(j)) / uvw(j) - ! end if - ! end do - ! j = minloc(d(:n), 1) - ! distance = d(j) - ! end if - - ! ! Assign the next tally bin and the score. - ! bin = m % get_bin_from_indices(ijk0(:n)) - ! call match % bins % push_back(bin) - ! call match % weights % push_back(distance / total_distance) - - ! ! Find the next mesh cell that the particle enters. - - ! ! If the particle track ends in that bin, then we are done. - ! if (all(ijk0(:n) == ijk1(:n))) exit - - ! ! Translate the starting coordintes by the distance to that face. This - ! ! should be the xyz that we computed the distance to in the last - ! ! iteration of the filter loop. - ! xyz0 = xyz0 + distance * uvw - - ! ! Increment the indices into the next mesh cell. - ! if (uvw(j) > ZERO) then - ! ijk0(j) = ijk0(j) + 1 - ! else - ! ijk0(j) = ijk0(j) - 1 - ! end if - - ! ! If the next indices are invalid, then the track has left the mesh and - ! ! we are done. - ! if (any(ijk0(:n) < 1) .or. any(ijk0(:n) > m % dimension)) exit - ! end do - end subroutine get_all_bins_mesh subroutine to_statepoint_mesh(this, filter_group) From 8af84064ca2f14a783b122153dde1b105ed687b9 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Sep 2018 07:01:27 -0500 Subject: [PATCH 12/21] Implement surface_bins_crossed for mesh --- include/openmc/mesh.h | 1 + src/mesh.cpp | 141 ++++++++++++++++++++- src/tallies/tally_filter_meshsurface.F90 | 154 +++-------------------- 3 files changed, 155 insertions(+), 141 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 47258edd82..c817564915 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -30,6 +30,7 @@ public: // Methods void bins_crossed(const Particle* p, std::vector& bins, std::vector& lengths); + void surface_bins_crossed(const Particle* p, std::vector& bins); int get_bin(Position r); int get_bin_from_indices(const int* ijk); void get_indices(Position r, int* ijk, bool* in_mesh); diff --git a/src/mesh.cpp b/src/mesh.cpp index 3749c55859..88aa6af126 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1,6 +1,6 @@ #include "openmc/mesh.h" -#include // for copy +#include // for copy, min #include // for size_t #include // for ceil #include @@ -515,6 +515,130 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } } +void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins) +{ + // ======================================================================== + // Determine if the track intersects the tally mesh. + + // Copy the starting and ending coordinates of the particle. + Position r0 {p->last_xyz_current}; + Position r1 {p->coord[0].xyz}; + Direction u {p->coord[0].uvw}; + + // Determine indices for starting and ending location. + int n = n_dimension_; + xt::xtensor ijk0 = xt::empty({n}); + bool start_in_mesh; + get_indices(r0, ijk0.data(), &start_in_mesh); + xt::xtensor ijk1 = xt::empty({n}); + bool end_in_mesh; + get_indices(r1, ijk1.data(), &end_in_mesh); + + // If this is the first iteration of the filter loop, check if the track + // intersects any part of the mesh. + if ((!start_in_mesh) && (!end_in_mesh)) { + if (!intersects(r0, r1)) return; + } + + // ======================================================================== + // Figure out which mesh cell to tally. + + // Calculate number of surface crossings + int n_cross = xt::sum(xt::abs(ijk1 - ijk0))(); + if (n_cross == 0) return; + + // Bounding coordinates + Position xyz_cross; + for (int i = 0; i < n; ++i) { + if (u[i] > 0.0) { + xyz_cross[i] = lower_left_[i] + ijk0[i] * width_[i]; + } else { + xyz_cross[i] = lower_left_[i] + (ijk0[i] - 1) * width_[i]; + } + } + + for (int j = 0; j < n_cross; ++j) { + // Set the distances to infinity + Position d {INFTY, INFTY, INFTY}; + + // Determine closest bounding surface. We need to treat + // special case where the cosine of the angle is zero since this would + // result in a divide-by-zero. + double distance = INFTY; + for (int i = 0; i < n; ++i) { + if (u[i] == 0) { + d[i] = INFINITY; + } else { + d[i] = (xyz_cross[i] - r0[i])/u[i]; + } + distance = std::min(distance, d[i]); + } + + // Loop over the dimensions + for (int i = 0; i < n; ++i) { + // Check whether distance is the shortest distance + if (distance == d[i]) { + + // Check whether particle is moving in positive i direction + if (u[i] > 0) { + + // Outward current on i max surface + if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { + int i_surf = 4*i + 3; + int i_mesh = get_bin_from_indices(ijk0.data()); + int i_bin = 4*n*(i_mesh - 1) + i_surf; + + bins.push_back(i_bin); + } + + // Advance position + ++ijk0[i]; + xyz_cross[i] += width_[i]; + + // If the particle crossed the surface, tally the inward current on + // i min surface + if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { + int i_surf = 4*i + 2; + int i_mesh = get_bin_from_indices(ijk0.data()); + int i_bin = 4*n*(i_mesh - 1) + i_surf; + + bins.push_back(i_bin); + } + + } else { + // The particle is moving in the negative i direction + + // Outward current on i min surface + if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_) ){ + int i_surf = 4*i + 1; + int i_mesh = get_bin_from_indices(ijk0.data()); + int i_bin = 4*n*(i_mesh - 1) + i_surf; + + bins.push_back(i_bin); + } + + // Advance position + --ijk0[i]; + xyz_cross[i] -= width_[i]; + + // If the particle crossed the surface, tally the inward current on + // i max surface + if (xt::all(ijk0 >= 1) && xt::all(ijk0 <= shape_)) { + int i_surf = 4*i + 4; + int i_mesh = get_bin_from_indices(ijk0.data()); + int i_bin = 4*n*(i_mesh - 1) + i_surf; + + bins.push_back(i_bin); + } + } + } + } + + // Calculate new coordinates + r0 += distance * u; + } +} + void RegularMesh::to_hdf5(hid_t group) { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); @@ -585,7 +709,6 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, // Check if there were sites outside the mesh for any processor MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & mpi_intracomm, mpi_err) - end if #else std::copy(cnt.data(), cnt.data() + total, cnt_reduced); if (outside) *outside = outside_; @@ -834,6 +957,20 @@ extern "C" { } } + void mesh_surface_bins_crossed(RegularMesh* m, const Particle* p, + void* match_bins, void* match_weights) + { + // Get surface bins crossed + std::vector bins; + m->surface_bins_crossed(p, bins); + + // Call bindings for VectorInt and VectorReal + for (auto b : bins) { + vector_int_push_back(match_bins, b); + vector_real_push_back(match_weights, 1.0); + } + } + void free_memory_mesh() { meshes.clear(); diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index d1500b8b39..f08bc35efb 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -3,7 +3,6 @@ module tally_filter_meshsurface use, intrinsic :: ISO_C_BINDING use constants - use dict_header, only: EMPTY use error use mesh_header use hdf5_interface @@ -41,7 +40,6 @@ contains integer :: i integer :: id integer :: n - integer :: n_dim integer(C_INT) :: err type(RegularMesh) :: m @@ -74,147 +72,25 @@ contains integer, intent(in) :: estimator type(TallyFilterMatch), intent(inout) :: match - integer :: j ! loop indices - integer :: n_dim ! num dimensions of the mesh - integer :: d1 ! dimension index - integer :: ijk0(3) ! indices of starting coordinates - integer :: ijk1(3) ! indices of ending coordinates - integer :: n_cross ! number of surface crossings - integer :: i_mesh ! flattened mesh bin index - integer :: i_surf ! surface index (1--12) - integer :: i_bin ! actual index for filter - real(8) :: uvw(3) ! cosine of angle of particle - real(8) :: xyz0(3) ! starting/intermediate coordinates - real(8) :: xyz1(3) ! ending coordinates of particle - real(8) :: xyz_cross(3) ! coordinates of bounding surfaces - real(8) :: d(3) ! distance to each bounding surface - real(8) :: distance ! actual distance traveled - logical :: start_in_mesh ! particle's starting xyz in mesh? - logical :: end_in_mesh ! particle's ending xyz in mesh? type(RegularMesh) :: m + type(C_PTR) :: ptr_bins, ptr_weights - ! ! Copy starting and ending location of particle - ! xyz0 = p % last_xyz_current - ! xyz1 = p % coord(1) % xyz + interface + subroutine mesh_surface_bins_crossed(m, p, bins, weights) bind(C) + import C_PTR, Particle + type(C_PTR), value :: m + type(Particle), intent(in) :: p + type(C_PTR), value :: bins + type(C_PTR), value :: weights + end subroutine + end interface - ! m = meshes(this % mesh) - ! n_dim = m % n_dimension() + ! Get a pointer to the mesh. + m = meshes(this % mesh) - ! ! Determine indices for starting and ending location - ! call m % get_indices(xyz0, ijk0, start_in_mesh) - ! call m % get_indices(xyz1, ijk1, end_in_mesh) - - ! ! Check to see if start or end is in mesh -- if not, check if track still - ! ! intersects with mesh - ! if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - ! !if (.not. m % intersects(xyz0, xyz1)) return - ! end if - - ! ! Calculate number of surface crossings - ! n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim))) - ! if (n_cross == 0) return - - ! ! Copy particle's direction - ! uvw = p % coord(1) % uvw - - ! ! Bounding coordinates - ! do d1 = 1, n_dim - ! if (uvw(d1) > 0) then - ! xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1) - ! else - ! xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1) - ! end if - ! end do - - ! do j = 1, n_cross - ! ! Set the distances to infinity - ! d = INFINITY - - ! ! Calculate distance to each bounding surface. We need to treat - ! ! special case where the cosine of the angle is zero since this would - ! ! result in a divide-by-zero. - ! do d1 = 1, n_dim - ! if (uvw(d1) == 0) then - ! d(d1) = INFINITY - ! else - ! d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1) - ! end if - ! end do - - ! ! Determine the closest bounding surface of the mesh cell by - ! ! calculating the minimum distance. Then use the minimum distance and - ! ! direction of the particle to determine which surface was crossed. - ! distance = minval(d) - - ! ! Loop over the dimensions - ! do d1 = 1, n_dim - - ! ! Check whether distance is the shortest distance - ! if (distance == d(d1)) then - - ! ! Check whether particle is moving in positive d1 direction - ! if (uvw(d1) > 0) then - - ! ! Outward current on d1 max surface - ! if (all(ijk0(:n_dim) >= 1) .and. & - ! all(ijk0(:n_dim) <= m % dimension)) then - ! i_surf = d1 * 4 - 1 - ! i_mesh = m % get_bin_from_indices(ijk0) - ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - - ! call match % bins % push_back(i_bin) - ! call match % weights % push_back(ONE) - ! end if - - ! ! Advance position - ! ijk0(d1) = ijk0(d1) + 1 - ! xyz_cross(d1) = xyz_cross(d1) + m % width(d1) - - ! ! If the particle crossed the surface, tally the inward current on - ! ! d1 min surface - ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then - ! i_surf = d1 * 4 - 2 - ! i_mesh = m % get_bin_from_indices(ijk0) - ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - - ! call match % bins % push_back(i_bin) - ! call match % weights % push_back(ONE) - ! end if - - ! else - ! ! The particle is moving in the negative d1 direction - - ! ! Outward current on d1 min surface - ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then - ! i_surf = d1 * 4 - 3 - ! i_mesh = m % get_bin_from_indices(ijk0) - ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - - ! call match % bins % push_back(i_bin) - ! call match % weights % push_back(ONE) - ! end if - - ! ! Advance position - ! ijk0(d1) = ijk0(d1) - 1 - ! xyz_cross(d1) = xyz_cross(d1) - m % width(d1) - - ! ! If the particle crossed the surface, tally the inward current on - ! ! d1 max surface - ! if (all(ijk0(:n_dim) >= 1)) then ! .and. all(ijk0(:n_dim) <= m % dimension)) then - ! i_surf = d1 * 4 - ! i_mesh = m % get_bin_from_indices(ijk0) - ! i_bin = 4*n_dim*(i_mesh - 1) + i_surf - - ! call match % bins % push_back(i_bin) - ! call match % weights % push_back(ONE) - ! end if - ! end if - ! end if - ! end do - - ! ! Calculate new coordinates - ! xyz0 = xyz0 + distance * uvw - ! end do + ptr_bins = C_LOC(match % bins) + ptr_weights = C_LOC(match % weights) + call mesh_surface_bins_crossed(m % ptr, p, ptr_bins, ptr_weights) end subroutine get_all_bins From 1b42fa20c6e2ec5cefb955f677c927c19677e0d5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 4 Sep 2018 07:26:57 -0500 Subject: [PATCH 13/21] Make sure to use std::fabs --- src/mesh.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 88aa6af126..87d5f1e289 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -432,9 +432,9 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } for (j = 0; j < n; ++j) { - if (std::abs(u[j]) < FP_PRECISION) { + if (std::fabs(u[j]) < FP_PRECISION) { d(j) = INFTY; - } else if (u[j] > 0) { + } else if (u[j] > 0.0) { double xyz_cross = lower_left_[j] + ijk0(j) * width_[j]; d(j) = (xyz_cross - r0[j]) / u[j]; } else { @@ -473,7 +473,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, // surface. xt::xtensor d = xt::zeros({n}); for (int j = 0; j < n; ++j) { - if (abs(u[j]) < FP_PRECISION) { + if (std::fabs(u[j]) < FP_PRECISION) { d(j) = INFTY; } else if (u[j] > 0) { double xyz_cross = lower_left_[j] + ijk0(j) * width_[j]; From 829b4c4f913e7aeb069ec35de16c046687d4f83c Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 06:35:38 -0500 Subject: [PATCH 14/21] Use C++ based count_sites for CMFD. Fix several bugs --- CMakeLists.txt | 1 + src/cmfd_execute.F90 | 18 ++++++++++++++---- src/cmfd_execute.cpp | 34 ++++++++++++++++++++++++++++++++++ src/cmfd_header.F90 | 8 +++++--- src/eigenvalue.F90 | 1 - src/eigenvalue.cpp | 1 - src/mesh.cpp | 6 +++--- src/message_passing.cpp | 2 +- src/settings.cpp | 2 +- 9 files changed, 59 insertions(+), 14 deletions(-) create mode 100644 src/cmfd_execute.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 713bf7c077..f4ff9ae56d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -380,6 +380,7 @@ add_library(libopenmc SHARED src/tallies/trigger.F90 src/tallies/trigger_header.F90 src/cell.cpp + src/cmfd_execute.cpp src/distribution.cpp src/distribution_angle.cpp src/distribution_energy.cpp diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 74264a036e..6709f782ac 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -223,7 +223,7 @@ contains integer :: nx ! maximum number of cells in x direction integer :: ny ! maximum number of cells in y direction integer :: nz ! maximum number of cells in z direction - integer :: ng ! maximum number of energy groups + integer(C_INT) :: ng ! maximum number of energy groups integer :: i ! iteration counter integer :: g ! index for group integer :: ijk(3) ! spatial bin location @@ -231,12 +231,22 @@ contains integer :: mesh_bin ! mesh bin of soruce particle integer :: n_groups ! number of energy groups real(8) :: norm ! normalization factor - logical :: outside ! any source sites outside mesh + logical(C_BOOL) :: outside ! any source sites outside mesh logical :: in_mesh ! source site is inside mesh #ifdef OPENMC_MPI integer :: mpi_err #endif + interface + subroutine cmfd_populate_sourcecounts(ng, energies, source_counts, outside) bind(C) + import C_INT, C_DOUBLE, C_BOOL + integer(C_INT), value :: ng + real(C_DOUBLE), intent(in) :: energies + real(C_DOUBLE), intent(out) :: source_counts + logical(C_BOOL), intent(out) :: outside + end subroutine + end interface + ! Get maximum of spatial and group indices nx = cmfd % indices(1) ny = cmfd % indices(2) @@ -260,8 +270,8 @@ contains cmfd%weightfactors = ONE ! Count bank sites in mesh and reverse due to egrid structure - call count_bank_sites(cmfd_mesh, source_bank, cmfd%sourcecounts, & - cmfd % egrid, sites_outside=outside, size_bank=work) + call cmfd_populate_sourcecounts(ng + 1, cmfd % egrid(1), & + cmfd % sourcecounts(1,1), outside) ! Check for sites outside of the mesh if (master .and. outside) then diff --git a/src/cmfd_execute.cpp b/src/cmfd_execute.cpp new file mode 100644 index 0000000000..7774cfe1d2 --- /dev/null +++ b/src/cmfd_execute.cpp @@ -0,0 +1,34 @@ +#include // for copy +#include +#include + +#include "xtensor/xarray.hpp" +#include "xtensor/xio.hpp" + +#include "openmc/capi.h" +#include "openmc/mesh.h" + +namespace openmc { + +extern "C" int index_cmfd_mesh; + +extern "C" void +cmfd_populate_sourcecounts(int n_energy, const double* energies, + double* source_counts, bool* outside) +{ + // Get pointer to source bank + Bank* source_bank; + int64_t n; + openmc_source_bank(&source_bank, &n); + + // Get source counts in each mesh bin / energy bin + auto& m = meshes.at(index_cmfd_mesh); + xt::xarray counts = m->count_sites(openmc_work, source_bank, n_energy, energies, outside); + + std::cout << counts << "\n"; + + // Copy data from the xarray into the source counts array + std::copy(counts.begin(), counts.end(), source_counts); +} + +} // namespace openmc diff --git a/src/cmfd_header.F90 b/src/cmfd_header.F90 index eafba7b8ad..e5e5cc423e 100644 --- a/src/cmfd_header.F90 +++ b/src/cmfd_header.F90 @@ -1,5 +1,7 @@ module cmfd_header + use, intrinsic :: ISO_C_BINDING + use constants, only: CMFD_NOACCEL, ZERO, ONE use mesh_header, only: RegularMesh use set_header, only: SetInt @@ -24,7 +26,7 @@ module cmfd_header integer :: mat_dim = CMFD_NOACCEL ! Energy grid - real(8), allocatable :: egrid(:) + real(C_DOUBLE), allocatable :: egrid(:) ! Cross sections real(8), allocatable :: totalxs(:,:,:,:) @@ -53,7 +55,7 @@ module cmfd_header real(8), allocatable :: openmc_src(:,:,:,:) ! Source sites in each mesh box - real(8), allocatable :: sourcecounts(:,:) + real(C_DOUBLE), allocatable :: sourcecounts(:,:) ! Weight adjustment factors real(8), allocatable :: weightfactors(:,:,:,:) @@ -95,7 +97,7 @@ module cmfd_header ! Main object type(cmfd_type), public :: cmfd - integer, public :: index_cmfd_mesh + integer(C_INT), public, bind(C) :: index_cmfd_mesh type(RegularMesh), public :: cmfd_mesh ! Pointers for different tallies diff --git a/src/eigenvalue.F90 b/src/eigenvalue.F90 index 14fe59110e..bd2a20b413 100644 --- a/src/eigenvalue.F90 +++ b/src/eigenvalue.F90 @@ -6,7 +6,6 @@ module eigenvalue use constants, only: ZERO use error, only: fatal_error, warning use math, only: t_percentile - use mesh, only: count_bank_sites use message_passing use random_lcg, only: prn, set_particle_seed, advance_prn_seed use settings diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index ef4159e782..fa995e6a03 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -80,7 +80,6 @@ void ufs_count_sites() int64_t n; openmc_source_bank(&source_bank, &n); - // count number of source sites in each ufs mesh cell bool sites_outside; source_frac = m->count_sites(openmc_work, source_bank, 0, nullptr, diff --git a/src/mesh.cpp b/src/mesh.cpp index 87d5f1e289..58f166f625 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -659,7 +659,7 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, std::size_t m = xt::prod(shape_)(); std::vector shape; if (n_energy > 0) { - shape = {m, static_cast(n_energy)}; + shape = {m, static_cast(n_energy - 1)}; } else { shape = {m}; } @@ -670,7 +670,8 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, for (int64_t i = 0; i < n; ++i) { // determine scoring bin for entropy mesh - int mesh_bin = get_bin({bank[i].xyz}); + // TODO: off-by-one + int mesh_bin = get_bin({bank[i].xyz}) - 1; // if outside mesh, skip particle if (mesh_bin < 0) { @@ -680,7 +681,6 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, if (n_energy > 0) { double E = bank[i].E; - int e_bin; if (E >= energies[0] && E <= energies[n_energy - 1]) { // determine energy bin int e_bin = lower_bound_index(energies, energies + n_energy, E); diff --git a/src/message_passing.cpp b/src/message_passing.cpp index ec2fc2d688..dc287e3b3a 100644 --- a/src/message_passing.cpp +++ b/src/message_passing.cpp @@ -5,7 +5,7 @@ namespace mpi { int rank {0}; int n_procs {1}; -bool master {false}; +bool master {true}; #ifdef OPENMC_MPI MPI_Comm intracomm; diff --git a/src/settings.cpp b/src/settings.cpp index db6aee7f54..7da8b31669 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -520,7 +520,7 @@ void read_settings_xml() // If the user did not specify how many mesh cells are to be used in // each direction, we automatically determine an appropriate number of // cells - int n = std::ceil(std::pow(settings::n_particles / 20, 1.0/3.0)); + int n = std::ceil(std::pow(settings::n_particles / 20.0, 1.0/3.0)); m.shape_ = {n, n, n}; m.n_dimension_ = 3; From 49cb5da76906344ad2989236808e5fedeb22559a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 07:02:55 -0500 Subject: [PATCH 15/21] Fix bug in UFS --- src/eigenvalue.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index fa995e6a03..71ce1c295a 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -110,7 +110,7 @@ void ufs_count_sites() double ufs_get_weight(const Particle* p) { - auto& m = meshes[settings::index_entropy_mesh]; + auto& m = meshes[settings::index_ufs_mesh]; // Determine indices on ufs mesh for current location // TODO: off by one From 7a6d36c279799581834fb740be9f666f748483a6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 09:08:16 -0500 Subject: [PATCH 16/21] Cosmetic changes --- src/mesh.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 58f166f625..c173cabe11 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -394,7 +394,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, // If this is the first iteration of the filter loop, check if the track // intersects any part of the mesh. - if ((!start_in_mesh) && (!end_in_mesh)) { + if (!start_in_mesh && !end_in_mesh) { if (!intersects(r0, r1)) return; } @@ -402,8 +402,8 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, // Figure out which mesh cell to tally. // Copy the un-modified coordinates the particle direction. - r0 = Position{p->last_xyz}; - r1 = Position{p->coord[0].xyz}; + r0 = last_r; + r1 = r; // Compute the length of the entire track. double total_distance = (r1 - r0).norm(); @@ -536,7 +536,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins // If this is the first iteration of the filter loop, check if the track // intersects any part of the mesh. - if ((!start_in_mesh) && (!end_in_mesh)) { + if (!start_in_mesh && !end_in_mesh) { if (!intersects(r0, r1)) return; } From f1fb400f925fbb9429b9f46f3004c2583cf43f69 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 11:23:20 -0500 Subject: [PATCH 17/21] Fix last mesh bugs (all tests pass now) --- CMakeLists.txt | 1 - openmc/capi/mesh.py | 6 +- src/cmfd_execute.F90 | 1 - src/mesh.F90 | 106 ------------------------------ src/mesh.cpp | 7 +- src/tallies/tally_filter_mesh.F90 | 18 ----- 6 files changed, 8 insertions(+), 131 deletions(-) delete mode 100644 src/mesh.F90 diff --git a/CMakeLists.txt b/CMakeLists.txt index f4ff9ae56d..d6821c4143 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -312,7 +312,6 @@ add_library(libopenmc SHARED src/material_header.F90 src/math.F90 src/matrix_header.F90 - src/mesh.F90 src/mesh_header.F90 src/message_passing.F90 src/mgxs_data.F90 diff --git a/openmc/capi/mesh.py b/openmc/capi/mesh.py index 091c5194b9..70a4345b09 100644 --- a/openmc/capi/mesh.py +++ b/openmc/capi/mesh.py @@ -41,6 +41,8 @@ _dll.openmc_mesh_set_params.errcheck = _error_handler _dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_get_mesh_index.restype = c_int _dll.openmc_get_mesh_index.errcheck = _error_handler +_dll.n_meshes.argtypes = [] +_dll.n_meshes.restype = c_int class Mesh(_FortranObjectWithID): @@ -172,10 +174,10 @@ class _MeshMapping(Mapping): def __iter__(self): for i in range(len(self)): - yield Mesh(index=i + 1).id + yield Mesh(index=i).id def __len__(self): - return c_int32.in_dll(_dll, 'n_meshes').value + return _dll.n_meshes() def __repr__(self): return repr(dict(self)) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index 6709f782ac..18f5ebec19 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -214,7 +214,6 @@ contains use bank_header, only: source_bank use constants, only: ZERO, ONE use error, only: warning, fatal_error - use mesh, only: count_bank_sites use message_passing use string, only: to_str diff --git a/src/mesh.F90 b/src/mesh.F90 deleted file mode 100644 index 790dc7d889..0000000000 --- a/src/mesh.F90 +++ /dev/null @@ -1,106 +0,0 @@ -module mesh - - use algorithm, only: binary_search - use bank_header, only: bank - use constants - use mesh_header - use message_passing - - implicit none - -contains - -!=============================================================================== -! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a -! given mesh as well as an optional energy group structure. This can be used for -! a variety of purposes (Shannon entropy, CMFD, uniform fission source -! weighting) -!=============================================================================== - - subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, & - sites_outside) - - type(RegularMesh), intent(in) :: m ! mesh to count sites - type(Bank), intent(in) :: bank_array(:) ! fission or source bank - real(8), intent(out) :: cnt(:,:) ! weight of sites in each - ! cell and energy group - real(8), intent(in), optional :: energies(:) ! energy grid to search - integer(8), intent(in), optional :: size_bank ! # of bank sites (on each proc) - logical, intent(inout), optional :: sites_outside ! were there sites outside mesh? - real(8), allocatable :: cnt_(:,:) - - integer :: i ! loop index for local fission sites - integer :: n_sites ! size of bank array - integer :: n ! number of energy groups / size - integer :: mesh_bin ! mesh bin - integer :: e_bin ! energy bin -#ifdef OPENMC_MPI - integer :: mpi_err ! MPI error code -#endif - logical :: outside ! was any site outside mesh? - - ! initialize variables - allocate(cnt_(size(cnt,1), size(cnt,2))) - cnt_ = ZERO - outside = .false. - - ! Set size of bank - if (present(size_bank)) then - n_sites = int(size_bank,4) - else - n_sites = size(bank_array) - end if - - ! Determine number of energies in group structure - if (present(energies)) then - n = size(energies) - 1 - else - n = 1 - end if - - ! loop over fission sites and count how many are in each mesh box - FISSION_SITES: do i = 1, n_sites - ! determine scoring bin for entropy mesh - call m % get_bin(bank_array(i) % xyz, mesh_bin) - - ! if outside mesh, skip particle - if (mesh_bin == NO_BIN_FOUND) then - outside = .true. - cycle - end if - - ! determine energy bin - if (present(energies)) then - if (bank_array(i) % E < energies(1)) then - e_bin = 1 - elseif (bank_array(i) % E > energies(n + 1)) then - e_bin = n - else - e_bin = binary_search(energies, n + 1, bank_array(i) % E) - end if - else - e_bin = 1 - end if - - ! add to appropriate mesh box - cnt_(e_bin, mesh_bin) = cnt_(e_bin, mesh_bin) + bank_array(i) % wgt - end do FISSION_SITES - -#ifdef OPENMC_MPI - ! collect values from all processors - n = size(cnt_) - call MPI_REDUCE(cnt_, cnt, n, MPI_REAL8, MPI_SUM, 0, mpi_intracomm, mpi_err) - - ! Check if there were sites outside the mesh for any processor - if (present(sites_outside)) then - call MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & - mpi_intracomm, mpi_err) - end if -#else - sites_outside = outside - cnt = cnt_ -#endif - - end subroutine count_bank_sites - -end module mesh diff --git a/src/mesh.cpp b/src/mesh.cpp index c173cabe11..f2bd7fe791 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -443,7 +443,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } } - int j = xt::argmin(d)(0); + j = xt::argmin(d)(0); if (u[j] > 0.0) { ++ijk0(j); } else { @@ -729,11 +729,11 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, extern "C" int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end) { - if (!index_start) *index_start = meshes.size(); + if (index_start) *index_start = meshes.size(); for (int i = 0; i < n; ++i) { meshes.emplace_back(new RegularMesh{}); } - if (!index_end) *index_end = meshes.size() - 1; + if (index_end) *index_end = meshes.size() - 1; return 0; } @@ -824,6 +824,7 @@ openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, *ll = m->lower_left_.data(); *ur = m->upper_right_.data(); + *width = m->width_.data(); *n = m->n_dimension_; return 0; } diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index bc770b199a..dda0bc0ce2 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -73,25 +73,7 @@ contains integer, intent(in) :: estimator type(TallyFilterMatch), intent(inout) :: match - integer, parameter :: MAX_SEARCH_ITER = 100 ! Maximum number of times we can - ! can loop while trying to find - ! the first intersection. - - integer :: j ! loop index for direction - integer :: n - integer :: ijk0(3) ! indices of starting coordinates - integer :: ijk1(3) ! indices of ending coordinates - integer :: search_iter ! loop count for intersection search integer :: bin - real(8) :: uvw(3) ! cosine of angle of particle - real(8) :: xyz0(3) ! starting/intermediate coordinates - real(8) :: xyz1(3) ! ending coordinates of particle - real(8) :: xyz_cross ! coordinates of next boundary - real(8) :: d(3) ! distance to each bounding surface - real(8) :: total_distance ! distance of entire particle track - real(8) :: distance ! distance traveled in mesh cell - logical :: start_in_mesh ! starting coordinates inside mesh? - logical :: end_in_mesh ! ending coordinates inside mesh? type(RegularMesh) :: m type(C_PTR) :: ptr_bins, ptr_weights From 50dcae177225b4ddd86eef7515caf3758da244a2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 12:27:18 -0500 Subject: [PATCH 18/21] Add Doxygen comments, remove unused variables --- include/openmc/capi.h | 1 - include/openmc/mesh.h | 62 ++++++++++++++++++++++++++++++++++++---- src/cmfd_input.F90 | 2 -- src/input_xml.F90 | 10 ------- src/multipole_header.F90 | 4 +-- src/nuclide_header.F90 | 1 - src/physics.F90 | 1 - src/physics_mg.F90 | 1 - src/settings.cpp | 13 ++++----- src/state_point.F90 | 2 +- 10 files changed, 64 insertions(+), 33 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index c155e7bd13..c48b212bd7 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -140,7 +140,6 @@ extern "C" { extern int32_t n_filters; extern int32_t n_lattices; extern int32_t n_materials; - extern int32_t n_meshes; extern int n_nuclides; extern int32_t n_plots; extern int32_t n_realizations; diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index c817564915..dca028287d 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -28,26 +28,76 @@ public: RegularMesh(pugi::xml_node node); // Methods + + //! Determine which bins were crossed by a particle + //! + //! \param[in] p Particle to check + //! \param[out] bins Bins that were crossed + //! \param[out] lengths Fraction of tracklength in each bin void bins_crossed(const Particle* p, std::vector& bins, std::vector& lengths); + + //! Determine which surface bins were crossed by a particle + //! + //! \param[in] p Particle to check + //! \param[out] bins Surface bins that were crossed void surface_bins_crossed(const Particle* p, std::vector& bins); + + //! Get bin at a given position in space + //! + //! \param[in] r Position to get bin for + //! \return Mesh bin int get_bin(Position r); + + //! Get bin given mesh indices + //! + //! \param[in] Array of mesh indices + //! \return Mesh bin int get_bin_from_indices(const int* ijk); + + //! Get mesh indices given a position + //! + //! \param[in] r Position to get indices for + //! \param[out] ijk Array of mesh indices + //! \param[out] in_mesh Whether position is in mesh void get_indices(Position r, int* ijk, bool* in_mesh); + + //! Get mesh indices corresponding to a mesh bin + //! + //! \param[in] bin Mesh bin + //! \param[out] ijk Mesh indices void get_indices_from_bin(int bin, int* ijk); + + //! Check if a line connected by two points intersects the mesh + //! + //! \param[in] r0 Starting position + //! \param[in] r1 Ending position + //! \return Whether line connecting r0 and r1 intersects mesh bool intersects(Position r0, Position r1); + + //! Write mesh data to an HDF5 group + //! + //! \param[in] group HDF5 group void to_hdf5(hid_t group); + //! Count number of bank sites in each mesh bin / energy bin + //! + //! \param[in] n Number of bank sites + //! \param[in] bank Array of bank sites + //! \param[in] n_energy Number of energies + //! \param[in] energies Array of energies + //! \param[out] Whether any bank sites are outside the mesh + //! \return Array indicating number of sites in each mesh/energy bin xt::xarray count_sites(int64_t n, const Bank* bank, int n_energy, const double* energies, bool* outside); int id_ {-1}; //!< User-specified ID - int n_dimension_; - double volume_frac_; - xt::xarray shape_; - xt::xarray lower_left_; - xt::xarray upper_right_; - xt::xarray width_; + int n_dimension_; //!< Number of dimensions + double volume_frac_; //!< Volume fraction of each mesh element + xt::xarray shape_; //!< Number of mesh elements in each dimension + xt::xarray lower_left_; //!< Lower-left coordinates of mesh + xt::xarray upper_right_; //!< Upper-right coordinates of mesh + xt::xarray width_; //!< Width of each mesh element private: bool intersects_1d(Position r0, Position r1); diff --git a/src/cmfd_input.F90 b/src/cmfd_input.F90 index 3208ad9d06..d1e0f6a9a7 100644 --- a/src/cmfd_input.F90 +++ b/src/cmfd_input.F90 @@ -263,8 +263,6 @@ contains integer :: i_filt ! index in filters array integer :: filt_id integer :: tally_id - integer :: iarray3(3) ! temp integer array - real(8) :: rarray3(3) ! temp double array real(C_DOUBLE), allocatable :: energies(:) type(XMLNode) :: node_mesh diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4766650bf9..e756a6aa59 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -203,22 +203,14 @@ contains subroutine read_settings_xml_f(root_ptr) bind(C) type(C_PTR), value :: root_ptr - character(MAX_LINE_LEN) :: temp_str integer :: i integer :: n - integer :: temp_int - integer :: temp_int_array3(3) - integer(C_INT32_T) :: i_start, i_end - integer(C_INT) :: err integer, allocatable :: temp_int_array(:) integer :: n_tracks type(XMLNode) :: root - type(XMLNode) :: node_entropy - type(XMLNode) :: node_ufs type(XMLNode) :: node_sp type(XMLNode) :: node_res_scat type(XMLNode) :: node_vol - type(XMLNode), allocatable :: node_mesh_list(:) type(XMLNode), allocatable :: node_vol_list(:) ! Get proper XMLNode type given pointer @@ -1164,13 +1156,11 @@ contains character(MAX_WORD_LEN), allocatable :: sarray(:) type(DictCharInt) :: trigger_scores type(TallyFilterContainer), pointer :: f - type(RegularMesh), pointer :: m type(XMLDocument) :: doc type(XMLNode) :: root type(XMLNode) :: node_tal type(XMLNode) :: node_filt type(XMLNode) :: node_trigger - type(XMLNode), allocatable :: node_mesh_list(:) type(XMLNode), allocatable :: node_tal_list(:) type(XMLNode), allocatable :: node_filt_list(:) type(XMLNode), allocatable :: node_trigger_list(:) diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index 05eb62e1aa..64610e4716 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -1,7 +1,6 @@ module multipole_header use constants - use dict_header, only: DictIntInt use error, only: fatal_error use hdf5_interface @@ -74,12 +73,11 @@ contains character(len=*), intent(in) :: filename character(len=10) :: version - integer :: i, n_poles, n_residues, n_windows + integer :: n_poles, n_residues, n_windows integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3) integer(HID_T) :: file_id integer(HID_T) :: group_id integer(HID_T) :: dset - type(DictIntInt) :: l_val_dict ! Open file for reading and move into the /isotope group file_id = file_open(filename, 'r', parallel=.true.) diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 481419fdc3..47a2d700a1 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -600,7 +600,6 @@ contains integer :: i, j, k, l integer :: t - integer :: m integer :: n integer :: n_grid integer :: i_fission diff --git a/src/physics.F90 b/src/physics.F90 index f1a0d7c57b..29a8064951 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -1181,7 +1181,6 @@ contains integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born integer :: i ! loop index integer :: nu ! actual number of neutrons produced - integer :: mesh_bin ! mesh bin for source site real(8) :: nu_t ! total nu real(8) :: weight ! weight adjustment for ufs method type(Nuclide), pointer :: nuc diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index bcafc43142..1d60b371ee 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -167,7 +167,6 @@ contains integer :: dg ! delayed group integer :: gout ! group out integer :: nu ! actual number of neutrons produced - integer :: mesh_bin ! mesh bin for source site real(8) :: nu_t ! total nu real(8) :: mu ! fission neutron angular cosine real(8) :: phi ! fission neutron azimuthal angle diff --git a/src/settings.cpp b/src/settings.cpp index 7da8b31669..007c515b0d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -3,7 +3,6 @@ #include // for ceil, pow #include // for numeric_limits #include -#include // for out_of_range #include #include @@ -490,13 +489,13 @@ void read_settings_xml() // Shannon Entropy mesh if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); - try { - index_entropy_mesh = mesh_map.at(temp); - } catch (const std::out_of_range& e) { + if (mesh_map.find(temp) == mesh_map.end()) { std::stringstream msg; msg << "Mesh " << temp << " specified for Shannon entropy does not exist."; fatal_error(msg); } + index_entropy_mesh = mesh_map.at(temp); + } else if (check_for_node(root, "entropy")) { warning("Specifying a Shannon entropy mesh via the element " "is deprecated. Please create a mesh using and then reference " @@ -535,14 +534,14 @@ void read_settings_xml() // Uniform fission source weighting mesh if (check_for_node(root, "ufs_mesh")) { auto temp = std::stoi(get_node_value(root, "ufs_mesh")); - try { - index_ufs_mesh = mesh_map.at(temp); - } catch (const std::out_of_range& e) { + if (mesh_map.find(temp) == mesh_map.end()) { std::stringstream msg; msg << "Mesh " << temp << " specified for uniform fission site method " "does not exist."; fatal_error(msg); } + index_ufs_mesh = mesh_map.at(temp); + } else if (check_for_node(root, "uniform_fs")) { warning("Specifying a UFS mesh via the element " "is deprecated. Please create a mesh using and then reference " diff --git a/src/state_point.F90 b/src/state_point.F90 index d20456d7b5..526066c529 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -67,7 +67,7 @@ contains integer :: i_xs integer, allocatable :: id_array(:) integer(HID_T) :: file_id - integer(HID_T) :: cmfd_group, tallies_group, tally_group, meshes_group, & + integer(HID_T) :: cmfd_group, tallies_group, tally_group, & filters_group, filter_group, derivs_group, & deriv_group, runtime_group integer(C_INT) :: ignored_err From 6174461daa08b36972b56d7af81ec1ba81300eea Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 13:50:19 -0500 Subject: [PATCH 19/21] Fix MPI-related bugs --- src/eigenvalue.cpp | 4 ++-- src/mesh.cpp | 8 +++----- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 71ce1c295a..d4f629e06b 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -92,8 +92,8 @@ void ufs_count_sites() #ifdef OPENMC_MPI // Send source fraction to all processors - int n = xt::prod(m->shape_)(); - MPI_Bcast(source_frac.data(), n, MPI_DOUBLE, 0, mpi::intracomm) + int n_bins = xt::prod(m->shape_)(); + MPI_Bcast(source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); #endif // Normalize to total weight to get fraction of source in each cell diff --git a/src/mesh.cpp b/src/mesh.cpp index f2bd7fe791..b922408de3 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -702,13 +702,11 @@ xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, // collect values from all processors MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0, mpi::intracomm); - if (outside) { - MPI_Reduce(&outside_, outside, 1, MPI_BOOL, MPI_LOR, 0, mpi::intracomm); - } // Check if there were sites outside the mesh for any processor - MPI_REDUCE(outside, sites_outside, 1, MPI_LOGICAL, MPI_LOR, 0, & - mpi_intracomm, mpi_err) + if (outside) { + MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm); + } #else std::copy(cnt.data(), cnt.data() + total, cnt_reduced); if (outside) *outside = outside_; From 7807c169c4ad5a1269ed2d464571d819b8d054fa Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Wed, 5 Sep 2018 15:06:49 -0500 Subject: [PATCH 20/21] Try using coveralls-python instead of python-coveralls --- tools/ci/travis-install.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tools/ci/travis-install.sh b/tools/ci/travis-install.sh index 5ad2381698..88c59d051b 100755 --- a/tools/ci/travis-install.sh +++ b/tools/ci/travis-install.sh @@ -22,4 +22,4 @@ python tools/ci/travis-install.py pip install -e .[test] # For uploading to coveralls -pip install python-coveralls +pip install coveralls From 4c47d868b4e5d481dc391aba26674199a99dbc7a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 6 Sep 2018 10:19:41 -0500 Subject: [PATCH 21/21] Make sure RegularMesh methods are const. Address other comments by @smharper --- include/openmc/mesh.h | 24 ++++++++++++------------ src/mesh.cpp | 39 +++++++++++++++++---------------------- 2 files changed, 29 insertions(+), 34 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index dca028287d..3d31fb4891 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -35,50 +35,50 @@ public: //! \param[out] bins Bins that were crossed //! \param[out] lengths Fraction of tracklength in each bin void bins_crossed(const Particle* p, std::vector& bins, - std::vector& lengths); + std::vector& lengths) const; //! Determine which surface bins were crossed by a particle //! //! \param[in] p Particle to check //! \param[out] bins Surface bins that were crossed - void surface_bins_crossed(const Particle* p, std::vector& bins); + void surface_bins_crossed(const Particle* p, std::vector& bins) const; //! Get bin at a given position in space //! //! \param[in] r Position to get bin for //! \return Mesh bin - int get_bin(Position r); + int get_bin(Position r) const; //! Get bin given mesh indices //! //! \param[in] Array of mesh indices //! \return Mesh bin - int get_bin_from_indices(const int* ijk); + int get_bin_from_indices(const int* ijk) const; //! Get mesh indices given a position //! //! \param[in] r Position to get indices for //! \param[out] ijk Array of mesh indices //! \param[out] in_mesh Whether position is in mesh - void get_indices(Position r, int* ijk, bool* in_mesh); + void get_indices(Position r, int* ijk, bool* in_mesh) const; //! Get mesh indices corresponding to a mesh bin //! //! \param[in] bin Mesh bin //! \param[out] ijk Mesh indices - void get_indices_from_bin(int bin, int* ijk); + void get_indices_from_bin(int bin, int* ijk) const; //! Check if a line connected by two points intersects the mesh //! //! \param[in] r0 Starting position //! \param[in] r1 Ending position //! \return Whether line connecting r0 and r1 intersects mesh - bool intersects(Position r0, Position r1); + bool intersects(Position r0, Position r1) const; //! Write mesh data to an HDF5 group //! //! \param[in] group HDF5 group - void to_hdf5(hid_t group); + void to_hdf5(hid_t group) const; //! Count number of bank sites in each mesh bin / energy bin //! @@ -89,7 +89,7 @@ public: //! \param[out] Whether any bank sites are outside the mesh //! \return Array indicating number of sites in each mesh/energy bin xt::xarray count_sites(int64_t n, const Bank* bank, - int n_energy, const double* energies, bool* outside); + int n_energy, const double* energies, bool* outside) const; int id_ {-1}; //!< User-specified ID int n_dimension_; //!< Number of dimensions @@ -100,9 +100,9 @@ public: xt::xarray width_; //!< Width of each mesh element private: - bool intersects_1d(Position r0, Position r1); - bool intersects_2d(Position r0, Position r1); - bool intersects_3d(Position r0, Position r1); + bool intersects_1d(Position r0, Position r1) const; + bool intersects_2d(Position r0, Position r1) const; + bool intersects_3d(Position r0, Position r1) const; }; //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index b922408de3..c1bd830a6f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -138,7 +138,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) } } -int RegularMesh::get_bin(Position r) +int RegularMesh::get_bin(Position r) const { // Loop over the dimensions of the mesh for (int i = 0; i < n_dimension_; ++i) { @@ -160,7 +160,7 @@ int RegularMesh::get_bin(Position r) return get_bin_from_indices(ijk); } -int RegularMesh::get_bin_from_indices(const int* ijk) +int RegularMesh::get_bin_from_indices(const int* ijk) const { if (n_dimension_ == 1) { return ijk[0]; @@ -171,7 +171,7 @@ int RegularMesh::get_bin_from_indices(const int* ijk) } } -void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) +void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) const { // Find particle in mesh *in_mesh = true; @@ -183,7 +183,7 @@ void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) } } -void RegularMesh::get_indices_from_bin(int bin, int* ijk) +void RegularMesh::get_indices_from_bin(int bin, int* ijk) const { if (n_dimension_ == 1) { ijk[0] = bin; @@ -192,12 +192,12 @@ void RegularMesh::get_indices_from_bin(int bin, int* ijk) ijk[1] = (bin - 1) / shape_[0] + 1; } else if (n_dimension_ == 3) { ijk[0] = (bin - 1) % shape_[0] + 1; - ijk[1] = (bin - 1) % (shape_[0] * shape_[1]) / shape_[0] + 1; + ijk[1] = ((bin - 1) % (shape_[0] * shape_[1])) / shape_[0] + 1; ijk[2] = (bin - 1) / (shape_[0] * shape_[1]) + 1; } } -bool RegularMesh::intersects(Position r0, Position r1) +bool RegularMesh::intersects(Position r0, Position r1) const { switch(n_dimension_) { case 1: @@ -209,7 +209,7 @@ bool RegularMesh::intersects(Position r0, Position r1) } } -bool RegularMesh::intersects_1d(Position r0, Position r1) +bool RegularMesh::intersects_1d(Position r0, Position r1) const { // Copy coordinates of mesh lower_left and upper_right double left = lower_left_[0]; @@ -225,7 +225,7 @@ bool RegularMesh::intersects_1d(Position r0, Position r1) } } -bool RegularMesh::intersects_2d(Position r0, Position r1) +bool RegularMesh::intersects_2d(Position r0, Position r1) const { // Copy coordinates of starting point double x0 = r0.x; @@ -280,7 +280,7 @@ bool RegularMesh::intersects_2d(Position r0, Position r1) return false; } -bool RegularMesh::intersects_3d(Position r0, Position r1) +bool RegularMesh::intersects_3d(Position r0, Position r1) const { // Copy coordinates of starting point double x0 = r0.x; @@ -365,7 +365,7 @@ bool RegularMesh::intersects_3d(Position r0, Position r1) } void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, - std::vector& lengths) + std::vector& lengths) const { constexpr int MAX_SEARCH_ITER = 100; @@ -392,8 +392,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, bool end_in_mesh; get_indices(r1, ijk1.data(), &end_in_mesh); - // If this is the first iteration of the filter loop, check if the track - // intersects any part of the mesh. + // Check if the track intersects any part of the mesh. if (!start_in_mesh && !end_in_mesh) { if (!intersects(r0, r1)) return; } @@ -459,8 +458,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, while (true) { // ======================================================================== - // Compute the length of the track segment in the appropiate mesh cell and - // return. + // Compute the length of the track segment in the each mesh cell and return double distance; int j; @@ -492,9 +490,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, bins.push_back(bin); lengths.push_back(distance / total_distance); - // Find the next mesh cell that the particle enters. - - // If the particle track ends in that bin, { we are done. + // If the particle track ends in that bin, then we are done. if (ijk0 == ijk1) break; // Translate the starting coordintes by the distance to that face. This @@ -515,7 +511,7 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } } -void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins) +void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins) const { // ======================================================================== // Determine if the track intersects the tally mesh. @@ -534,8 +530,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins bool end_in_mesh; get_indices(r1, ijk1.data(), &end_in_mesh); - // If this is the first iteration of the filter loop, check if the track - // intersects any part of the mesh. + // Check if the track intersects any part of the mesh. if (!start_in_mesh && !end_in_mesh) { if (!intersects(r0, r1)) return; } @@ -639,7 +634,7 @@ void RegularMesh::surface_bins_crossed(const Particle* p, std::vector& bins } } -void RegularMesh::to_hdf5(hid_t group) +void RegularMesh::to_hdf5(hid_t group) const { hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); @@ -653,7 +648,7 @@ void RegularMesh::to_hdf5(hid_t group) } xt::xarray RegularMesh::count_sites(int64_t n, const Bank* bank, - int n_energy, const double* energies, bool* outside) + int n_energy, const double* energies, bool* outside) const { // Determine shape of array for counts std::size_t m = xt::prod(shape_)();