Bindings for bins_crossed

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
 //==============================================================================

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)

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<int> 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<int> bins;
+    std::vector<double> lengths;
+    m->bins_crossed(p, bins, lengths);
 
-    if (n_meshes > 0) {
-      // Write IDs of meshes
-      std::vector<int> 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()

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)

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
 
 !===============================================================================

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
 

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

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.

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)