From d15090099a781c77bd5679d82fb9e75fdb596390 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 28 May 2019 15:51:51 -0400 Subject: [PATCH] Implement RectilinearMesh::surface_bins_crossed --- openmc/filter.py | 5 +- src/mesh.cpp | 139 ++++++++++++++++------------- src/tallies/filter_meshsurface.cpp | 19 ++-- 3 files changed, 95 insertions(+), 68 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 77fe84c3ae..7ad901b59d 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -773,13 +773,12 @@ class MeshSurfaceFilter(MeshFilter): @MeshFilter.mesh.setter def mesh(self, mesh): - cv.check_type('filter mesh', mesh, openmc.Mesh) + cv.check_type('filter mesh', mesh, openmc.MeshBase) self._mesh = mesh # Take the product of mesh indices and current names - n_dim = len(mesh.dimension) self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in - product(mesh.indices, _CURRENT_NAMES[:4*n_dim])] + product(mesh.indices, _CURRENT_NAMES[:4*3])] def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. diff --git a/src/mesh.cpp b/src/mesh.cpp index ea2b353949..0e3198c63f 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -911,12 +911,68 @@ void RectilinearMesh::surface_bins_crossed(const Particle* p, bool end_in_mesh; get_indices(r1, ijk1, &end_in_mesh); - // Check if the track intersects any part of the mesh. + // If the starting coordinates do not lie in the mesh, compute the coords and + // mesh indices of the first intersection, and add the bin for this first + // intersection. Return if the particle does not intersect the mesh at all. if (!start_in_mesh) { - Position r0_copy = r0; - int ijk0_copy[3]; - for (int i = 0; i < 3; ++i) ijk0_copy[i] = ijk0[i]; - if (!intersects(r0_copy, r1, ijk0_copy)) return; + // Compute the incoming intersection coordinates and indices. + if (!intersects(r0, r1, ijk0)) return; + + // Determine which surface the particle entered. + double min_dist = INFTY; + int i_surf; + for (int i = 0; i < 3; ++i) { + if (u[i] > 0.0 && ijk0[i] == 1) { + double d = std::abs(r0[i] - grid_[i][0]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 2; + } + } else if (u[i] < 0.0 && ijk0[i] == shape_[i]) { + double d = std::abs(r0[i] - grid_[i][shape_[i]]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 4; + } + } // u[i] == 0 intentionally skipped + } + + // Add the incoming current bin. + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); + } + + // If the ending coordinates do not lie in the mesh, compute the coords and + // mesh indices of the last intersection, and add the bin for this last + // intersection. + if (!end_in_mesh) { + // Compute the outgoing intersection coordinates and indices. + intersects(r1, r0, ijk1); + + // Determine which surface the particle exited. + double min_dist = INFTY; + int i_surf; + for (int i = 0; i < 3; ++i) { + if (u[i] > 0.0 && ijk1[i] == shape_[i]) { + double d = std::abs(r1[i] - grid_[i][shape_[i]]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 3; + } + } else if (u[i] < 0.0 && ijk1[i] == 1) { + double d = std::abs(r1[i] - grid_[i][0]); + if (d < min_dist) { + min_dist = d; + i_surf = 4*i + 1; + } + } // u[i] == 0 intentionally skipped + } + + // Add the outgoing current bin. + int i_mesh = get_bin_from_indices(ijk1); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); } // ======================================================================== @@ -959,80 +1015,43 @@ void RectilinearMesh::surface_bins_crossed(const Particle* p, // Check whether distance is the shortest distance if (distance == d[i]) { - // Check whether the current indices are within the mesh bounds - bool in_mesh = true; - for (int j = 0; j < 3; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - // Check whether particle is moving in positive i direction if (u[i] > 0) { // Outward current on i max surface - if (in_mesh) { - int i_surf = 4*i + 3; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4*3*i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } + int i_surf = 4*i + 3; + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); // Advance position ++ijk0[i]; xyz_cross[i] = grid_[i][ijk0[i]]; - in_mesh = true; - for (int j = 0; j < 3; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - // If the particle crossed the surface, tally the inward current on - // i min surface - if (in_mesh) { - int i_surf = 4*i + 2; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4*3*i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } + // Inward current on i min surface + i_surf = 4*i + 2; + i_mesh = get_bin_from_indices(ijk0); + i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); } else { // The particle is moving in the negative i direction // Outward current on i min surface - if (in_mesh) { - int i_surf = 4*i + 1; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4*3*i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } + int i_surf = 4*i + 1; + int i_mesh = get_bin_from_indices(ijk0); + int i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); // Advance position --ijk0[i]; xyz_cross[i] = grid_[i][ijk0[i] - 1]; - in_mesh = true; - for (int j = 0; j < 3; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - // If the particle crossed the surface, tally the inward current on - // i max surface - if (in_mesh) { - int i_surf = 4*i + 4; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4*3*i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } + // Inward current on i min surface + i_surf = 4*i + 4; + i_mesh = get_bin_from_indices(ijk0); + i_bin = 4*3*i_mesh + i_surf - 1; + bins.push_back(i_bin); } } } diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index 8f4ae4e356..aa64520139 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -18,8 +18,12 @@ MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator, std::string MeshSurfaceFilter::text_label(int bin) const { - auto& mesh = *dynamic_cast(model::meshes[mesh_].get()); - int n_dim = mesh.n_dimension_; + int n_dim; + if (auto* mesh = dynamic_cast(model::meshes[mesh_].get())) { + n_dim = mesh->n_dimension_; + } else if (auto* mesh = dynamic_cast(model::meshes[mesh_].get())) { + n_dim = 3; + } // Get flattend mesh index and surface index. int i_mesh = bin / (4 * n_dim); @@ -75,9 +79,14 @@ void MeshSurfaceFilter::set_mesh(int32_t mesh) { mesh_ = mesh; - auto &m = *dynamic_cast(model::meshes[mesh_].get()); - n_bins_ = 4 * m.n_dimension_; - for (auto dim : m.shape_) n_bins_ *= dim; + auto* m_ptr = model::meshes[mesh_].get(); + if (auto* m = dynamic_cast(m_ptr)) { + n_bins_ = 4 * m->n_dimension_; + for (auto dim : m->shape_) n_bins_ *= dim; + } else if (auto* m = dynamic_cast(m_ptr)) { + n_bins_ = 4 * 3; + for (auto dim : m->shape_) n_bins_ *= dim; + } } //==============================================================================