From 5a73710110a354fac42993a4ea2bbc6edd8b2e6f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 7 Oct 2018 15:12:04 -0400 Subject: [PATCH] Move mesh and meshsurface filter internals to C++ --- include/openmc/tallies/tally_filter_mesh.h | 26 +++- .../openmc/tallies/tally_filter_meshsurface.h | 58 ++++++++- src/cmfd_data.F90 | 10 +- src/tallies/tally_filter_header.F90 | 1 + src/tallies/tally_filter_mesh.F90 | 94 +++------------ src/tallies/tally_filter_meshsurface.F90 | 111 ------------------ src/tallies/trigger.F90 | 2 +- 7 files changed, 107 insertions(+), 195 deletions(-) diff --git a/include/openmc/tallies/tally_filter_mesh.h b/include/openmc/tallies/tally_filter_mesh.h index d48e8e35b..01ffbae6f 100644 --- a/include/openmc/tallies/tally_filter_mesh.h +++ b/include/openmc/tallies/tally_filter_mesh.h @@ -53,7 +53,31 @@ public: } } - virtual std::string text_label(int bin) const {}; + virtual void + to_statepoint(hid_t filter_group) const override + { + write_dataset(filter_group, "type", "mesh"); + write_dataset(filter_group, "n_bins", n_bins_); + write_dataset(filter_group, "bins", meshes[mesh_]->id_); + } + + virtual std::string + text_label(int bin) const override + { + auto& mesh = *meshes[mesh_]; + int n_dim = mesh.n_dimension_; + + int ijk[n_dim]; + mesh.get_indices_from_bin(bin, ijk); + + std::stringstream out; + out << "Mesh Index (" << ijk[0]; + if (n_dim > 1) out << ", " << ijk[1]; + if (n_dim > 2) out << ", " << ijk[2]; + out << ")"; + + return out.str(); + } int32_t mesh_; }; diff --git a/include/openmc/tallies/tally_filter_meshsurface.h b/include/openmc/tallies/tally_filter_meshsurface.h index 5bf06e095..7b374b0a6 100644 --- a/include/openmc/tallies/tally_filter_meshsurface.h +++ b/include/openmc/tallies/tally_filter_meshsurface.h @@ -4,6 +4,7 @@ #include #include +#include "openmc/constants.h" #include "openmc/error.h" #include "openmc/mesh.h" #include "openmc/tallies/tally_filter.h" @@ -45,7 +46,62 @@ public: for (auto b : match.bins) match.weights.push_back(1.0); } - virtual std::string text_label(int bin) const {}; + virtual void + to_statepoint(hid_t filter_group) const override + { + write_dataset(filter_group, "type", "meshsurface"); + write_dataset(filter_group, "n_bins", n_bins_); + write_dataset(filter_group, "bins", meshes[mesh_]->id_); + } + + virtual std::string + text_label(int bin) const override + { + auto& mesh = *meshes[mesh_]; + int n_dim = mesh.n_dimension_; + + // Get flattend mesh index and surface index. + int i_mesh = (bin - 1) / (4 * n_dim) + 1; + int i_surf = ((bin - 1) % (4 * n_dim)) + 1; + + // Get mesh index part of label. + int ijk[n_dim]; + mesh.get_indices_from_bin(i_mesh, ijk); + std::stringstream out; + out << "Mesh Index (" << ijk[0]; + if (n_dim > 1) out << ", " << ijk[1]; + if (n_dim > 2) out << ", " << ijk[2]; + out << ")"; + + // Get surface part of label. + if (i_surf == OUT_LEFT) { + out << " Outgoing, x-min"; + } else if (i_surf == IN_LEFT) { + out << " Incoming, x-min"; + } else if (i_surf == OUT_RIGHT) { + out << " Outgoing, x-max"; + } else if (i_surf == IN_RIGHT) { + out << " Incoming, x-max"; + } else if (i_surf == OUT_BACK) { + out << " Outgoing, y-min"; + } else if (i_surf == IN_BACK) { + out << " Incoming, y-min"; + } else if (i_surf == OUT_FRONT) { + out << " Outgoing, y-max"; + } else if (i_surf == IN_FRONT) { + out << " Incoming, y-max"; + } else if (i_surf == OUT_BOTTOM) { + out << " Outgoing, z-min"; + } else if (i_surf == IN_BOTTOM) { + out << " Incoming, z-min"; + } else if (i_surf == OUT_TOP) { + out << " Outgoing, z-max"; + } else if (i_surf == IN_TOP) { + out << " Incoming, z-max"; + } + + return out.str(); + } int32_t mesh_; }; diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 07276d900..a12ca217e 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -82,6 +82,14 @@ contains real(8) :: flux ! temp variable for flux type(RegularMesh) :: m ! pointer for mesh object + interface + function mesh_filter_get_mesh(filt) result(index_mesh) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: index_mesh + end function mesh_filter_get_mesh + end interface + ! Extract spatial and energy indices from object nx = cmfd % indices(1) ny = cmfd % indices(2) @@ -99,7 +107,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 diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index 162157489..93e1a82d6 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -379,6 +379,7 @@ contains class(CppTallyFilter), intent(inout) :: this type(XMLNode), intent(in) :: node call this % from_xml_cpp_inner(node) + this % n_bins = this % n_bins_cpp() end subroutine from_xml_cpp_default subroutine get_all_bins_cpp_default(this, p, estimator, match) diff --git a/src/tallies/tally_filter_mesh.F90 b/src/tallies/tally_filter_mesh.F90 index 816a58e53..5d2307aef 100644 --- a/src/tallies/tally_filter_mesh.F90 +++ b/src/tallies/tally_filter_mesh.F90 @@ -3,14 +3,9 @@ module tally_filter_mesh use, intrinsic :: ISO_C_BINDING use constants - use dict_header, only: EMPTY use error use mesh_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none private @@ -24,77 +19,24 @@ module tally_filter_mesh !=============================================================================== type, public, extends(CppTallyFilter) :: MeshFilter - integer :: mesh contains - procedure :: from_xml - procedure :: to_statepoint => to_statepoint_mesh - procedure :: text_label => text_label_mesh + procedure :: mesh => get_mesh end type MeshFilter contains - subroutine from_xml(this, node) - class(MeshFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: id - integer :: n - integer(C_INT) :: err - - call this % from_xml_cpp_inner(node) - - n = node_word_count(node, "bins") - - if (n /= 1) call fatal_error("Only one mesh can be & - &specified per mesh filter.") - - ! Determine id of mesh - call get_node_value(node, "bins", id) - - ! Get pointer to mesh - 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 = this % n_bins_cpp() - end subroutine from_xml - - subroutine to_statepoint_mesh(this, filter_group) - 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", m % id()) - end subroutine to_statepoint_mesh - - function text_label_mesh(this, bin) result(label) - class(MeshFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer, allocatable :: ijk(:) - type(RegularMesh) :: m - - 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 + function get_mesh(this) result(mesh) + class(MeshFilter) :: this + integer :: mesh + interface + function mesh_filter_get_mesh(filt) result(index_mesh) bind(C) + import C_PTR, C_INT + type(C_PTR), value :: filt + integer(C_INT) :: index_mesh + end function mesh_filter_get_mesh + end interface + mesh = mesh_filter_get_mesh(this % ptr) + end function get_mesh !=============================================================================== ! C API FUNCTIONS @@ -106,19 +48,11 @@ contains integer(C_INT32_T), intent(out) :: index_mesh integer(C_INT) :: err - interface - function mesh_filter_get_mesh(filt) result(index_mesh) bind(C) - import C_PTR, C_INT - type(C_PTR), value :: filt - integer(C_INT) :: index_mesh - end function mesh_filter_get_mesh - end interface - err = verify_filter(index) if (err == 0) then select type (f => filters(index) % obj) type is (MeshFilter) - index_mesh = mesh_filter_get_mesh(f % ptr) + index_mesh = f % mesh() class default err = E_INVALID_TYPE call set_errmsg("Tried to set mesh on a non-mesh filter.") diff --git a/src/tallies/tally_filter_meshsurface.F90 b/src/tallies/tally_filter_meshsurface.F90 index d2ff898b0..22453a40f 100644 --- a/src/tallies/tally_filter_meshsurface.F90 +++ b/src/tallies/tally_filter_meshsurface.F90 @@ -5,11 +5,7 @@ module tally_filter_meshsurface use constants use error use mesh_header - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none private @@ -23,116 +19,10 @@ module tally_filter_meshsurface !=============================================================================== type, public, extends(CppTallyFilter) :: MeshSurfaceFilter - integer :: mesh - contains - procedure :: from_xml - procedure :: to_statepoint - procedure :: text_label end type MeshSurfaceFilter contains - subroutine from_xml(this, node) - class(MeshSurfaceFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: id - integer :: n - integer(C_INT) :: err - - call this % from_xml_cpp_inner(node) - - n = node_word_count(node, "bins") - - if (n /= 1) call fatal_error("Only one mesh can be & - &specified per meshsurface filter.") - - ! Determine id of mesh - call get_node_value(node, "bins", id) - - ! Get pointer to mesh - 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 = this % n_bins_cpp() - end subroutine from_xml - - subroutine to_statepoint(this, filter_group) - 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", m % id()) - end subroutine to_statepoint - - function text_label(this, bin) result(label) - class(MeshSurfaceFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: i_mesh - integer :: i_surf - integer :: n_dim - integer, allocatable :: ijk(:) - type(RegularMesh) :: m - - 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 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 function text_label - !=============================================================================== ! C API FUNCTIONS !=============================================================================== @@ -155,7 +45,6 @@ contains if (err == 0) then select type (f => filters(index) % obj) type is (MeshSurfaceFilter) - index_mesh = f % mesh index_mesh = meshsurface_filter_get_mesh(f % ptr) class default err = E_INVALID_TYPE diff --git a/src/tallies/trigger.F90 b/src/tallies/trigger.F90 index 4450a292a..fbecaecd9 100644 --- a/src/tallies/trigger.F90 +++ b/src/tallies/trigger.F90 @@ -264,7 +264,7 @@ contains 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