Move most active_tally lists to C++

include/openmc/tallies/tally.h

@@ -30,6 +30,8 @@ public:
 
   int32_t n_filter_bins() const {return n_filter_bins_;}
 
+  int type_ {TALLY_VOLUME}; //!< volume, surface current
+
   //! Event type that contributes to this tally
   int estimator_ {ESTIMATOR_TRACKLENGTH};
 
@@ -64,12 +66,12 @@ extern "C" double total_weight;
 namespace model {
   extern std::vector<std::unique_ptr<Tally>> tallies;
 
+  extern std::vector<int> active_tallies;
   extern std::vector<int> active_analog_tallies;
   extern std::vector<int> active_tracklength_tallies;
-  extern std::vector<int> active_meshsurf_tallies;
   extern std::vector<int> active_collision_tallies;
-  extern std::vector<int> active_tallies;
-  extern std::vector<int> active_surface_tallies;
+  //extern std::vector<int> active_meshsurf_tallies;
+  //extern std::vector<int> active_surface_tallies;
 }
 
 // Threadprivate variables

src/input_xml.F90

@@ -1394,12 +1394,12 @@ contains
                      &in the same tally as normal surface currents")
               end if
 
-              t % type = TALLY_SURFACE
+              call t % set_type(TALLY_SURFACE)
               t % score_bins(j) = SCORE_CURRENT
 
             else if (has_meshsurface) then
               t % score_bins(j) = SCORE_CURRENT
-              t % type = TALLY_MESH_SURFACE
+              call t % set_type(TALLY_MESH_SURFACE)
 
               ! Check to make sure that current is the only desired response
               ! for this tally

src/state_point.F90

@@ -346,7 +346,7 @@ contains
       call write_dataset(file_id, "global_tallies", global_tallies)
 
       ! Write tallies
-      if (active_tallies % size() > 0) then
+      if (active_tallies_size() > 0) then
         ! Indicate that tallies are on
         call write_attribute(file_id, "tallies_present", 1)
 

src/tallies/tally.F90

@@ -2066,9 +2066,9 @@ contains
     ! A loop over all tallies is necessary because we need to simultaneously
     ! determine different filter bins for the same tally in order to score to it
 
-    TALLY_LOOP: do i = 1, active_analog_tallies % size()
+    TALLY_LOOP: do i = 1, active_analog_tallies_size()
       ! Get index of tally and pointer to tally
-      i_tally = active_analog_tallies % data(i)
+      i_tally = active_analog_tallies_data(i)
       associate (t => tallies(i_tally) % obj)
 
       ! Find all valid bins in each filter if they have not already been found
@@ -2212,9 +2212,9 @@ contains
     ! A loop over all tallies is necessary because we need to simultaneously
     ! determine different filter bins for the same tally in order to score to it
 
-    TALLY_LOOP: do i = 1, active_analog_tallies % size()
+    TALLY_LOOP: do i = 1, active_analog_tallies_size()
       ! Get index of tally and pointer to tally
-      i_tally = active_analog_tallies % data(i)
+      i_tally = active_analog_tallies_data(i)
       associate (t => tallies(i_tally) % obj)
 
       ! Find all valid bins in each filter if they have not already been found
@@ -2577,9 +2577,9 @@ contains
     ! A loop over all tallies is necessary because we need to simultaneously
     ! determine different filter bins for the same tally in order to score to it
 
-    TALLY_LOOP: do i = 1, active_tracklength_tallies % size()
+    TALLY_LOOP: do i = 1, active_tracklength_tallies_size()
       ! Get index of tally and pointer to tally
-      i_tally = active_tracklength_tallies % data(i)
+      i_tally = active_tracklength_tallies_data(i)
       associate (t => tallies(i_tally) % obj)
 
       ! Find all valid bins in each filter if they have not already been found
@@ -2734,9 +2734,9 @@ contains
     ! A loop over all tallies is necessary because we need to simultaneously
     ! determine different filter bins for the same tally in order to score to it
 
-    TALLY_LOOP: do i = 1, active_collision_tallies % size()
+    TALLY_LOOP: do i = 1, active_collision_tallies_size()
       ! Get index of tally and pointer to tally
-      i_tally = active_collision_tallies % data(i)
+      i_tally = active_collision_tallies_data(i)
       associate (t => tallies(i_tally) % obj)
 
       ! Find all valid bins in each filter if they have not already been found
@@ -3774,8 +3774,8 @@ contains
     end if
 
     ! Accumulate results for each tally
-    do i = 1, active_tallies % size()
-      call tallies(active_tallies % data(i)) % obj % accumulate()
+    do i = 1, active_tallies_size()
+      call tallies(active_tallies_data(i)) % obj % accumulate()
     end do
 
   end subroutine accumulate_tallies
@@ -3790,10 +3790,13 @@ contains
     integer(C_INT) :: err
     logical(C_BOOL) :: active
 
-    call active_tallies % clear()
-    call active_analog_tallies % clear()
-    call active_collision_tallies % clear()
-    call active_tracklength_tallies % clear()
+    interface
+      subroutine setup_active_tallies_c() bind(C)
+      end subroutine
+    end interface
+
+    call setup_active_tallies_c()
+
     call active_surface_tallies % clear()
     call active_meshsurf_tallies % clear()
 
@@ -3801,21 +3804,10 @@ contains
       associate (t => tallies(i) % obj)
         err = openmc_tally_get_active(i, active)
         if (active) then
-          ! Add tally to active tallies
-          call active_tallies % push_back(i)
-
           ! Check what type of tally this is and add it to the appropriate list
-          if (t % type == TALLY_VOLUME) then
-            if (t % estimator() == ESTIMATOR_ANALOG) then
-              call active_analog_tallies % push_back(i)
-            elseif (t % estimator() == ESTIMATOR_TRACKLENGTH) then
-              call active_tracklength_tallies % push_back(i)
-            elseif (t % estimator() == ESTIMATOR_COLLISION) then
-              call active_collision_tallies % push_back(i)
-            end if
-          elseif (t % type == TALLY_MESH_SURFACE) then
+          if (t % type() == TALLY_MESH_SURFACE) then
             call active_meshsurf_tallies % push_back(i)
-          elseif (t % type == TALLY_SURFACE) then
+          elseif (t % type() == TALLY_SURFACE) then
             call active_surface_tallies % push_back(i)
           end if
 

src/tallies/tally.cpp

@@ -27,12 +27,12 @@ namespace openmc {
 namespace model {
   std::vector<std::unique_ptr<Tally>> tallies;
 
+  std::vector<int> active_tallies;
   std::vector<int> active_analog_tallies;
   std::vector<int> active_tracklength_tallies;
-  std::vector<int> active_meshsurf_tallies;
   std::vector<int> active_collision_tallies;
-  std::vector<int> active_tallies;
-  std::vector<int> active_surface_tallies;
+  //std::vector<int> active_meshsurf_tallies;
+  //std::vector<int> active_surface_tallies;
 }
 
 double global_tally_absorption;
@@ -177,6 +177,48 @@ void reduce_tally_results()
 }
 #endif
 
+extern "C" void
+setup_active_tallies_c()
+{
+  model::active_tallies.clear();
+  model::active_analog_tallies.clear();
+  model::active_tracklength_tallies.clear();
+  model::active_collision_tallies.clear();
+  //model::active_meshsurf_tallies.clear();
+  //model::active_surface_tallies.clear();
+
+  //TODO: off-by-one all through here
+  for (auto i = 0; i < model::tallies.size(); ++i) {
+    const auto& tally {*model::tallies[i]};
+
+    if (tally.active_) {
+      model::active_tallies.push_back(i + 1);
+      switch (tally.type_) {
+
+      case TALLY_VOLUME:
+        switch (tally.estimator_) {
+          case ESTIMATOR_ANALOG:
+            model::active_analog_tallies.push_back(i + 1);
+            break;
+          case ESTIMATOR_TRACKLENGTH:
+            model::active_tracklength_tallies.push_back(i + 1);
+            break;
+          case ESTIMATOR_COLLISION:
+            model::active_collision_tallies.push_back(i + 1);
+        }
+        break;
+
+      //case TALLY_MESH_SURFACE:
+      //  model::active_meshsurf_tallies.push_back(i + 1);
+      //  break;
+
+      //case TALLY_SURFACE:
+      //  model::active_surface_tallies.push_back(i + 1);
+      }
+    }
+  }
+}
+
 extern "C" void
 free_memory_tally_c()
 {
@@ -189,12 +231,51 @@ free_memory_tally_c()
   model::tally_filters.clear();
 
   model::tallies.clear();
+
+  model::active_tallies.clear();
+  model::active_analog_tallies.clear();
+  model::active_tracklength_tallies.clear();
+  model::active_collision_tallies.clear();
+  //model::active_meshsurf_tallies.clear();
+  //model::active_surface_tallies.clear();
 }
 
 //==============================================================================
 // C-API functions
 //==============================================================================
 
+extern "C" int
+openmc_tally_get_type(int32_t index, int32_t* type)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  //TODO: off-by-one
+  *type = model::tallies[index-1]->type_;
+}
+
+extern "C" int
+openmc_tally_set_type(int32_t index, const char* type)
+{
+  if (index < 1 || index > model::tallies.size()) {
+    set_errmsg("Index in tallies array is out of bounds.");
+    return OPENMC_E_OUT_OF_BOUNDS;
+  }
+  if (strcmp(type, "volume") == 0) {
+    model::tallies[index-1]->type_ = TALLY_VOLUME;
+  } else if (strcmp(type, "mesh-surface") == 0) {
+    model::tallies[index-1]->type_ = TALLY_MESH_SURFACE;
+  } else if (strcmp(type, "surface") == 0) {
+    model::tallies[index-1]->type_ = TALLY_SURFACE;
+  } else {
+    std::stringstream errmsg;
+    errmsg << "Unknown tally type: " << type;
+    set_errmsg(errmsg);
+    return OPENMC_E_INVALID_ARGUMENT;
+  }
+}
+
 extern "C" int
 openmc_tally_get_active(int32_t index, bool* active)
 {
@@ -266,6 +347,34 @@ extern "C" {
       model::tallies.push_back(std::make_unique<Tally>());
   }
 
+  int active_tallies_data(int i)
+  {return model::active_tallies[i-1];}
+
+  int active_tallies_size()
+  {return model::active_tallies.size();}
+
+  int active_analog_tallies_data(int i)
+  {return model::active_analog_tallies[i-1];}
+
+  int active_analog_tallies_size()
+  {return model::active_analog_tallies.size();}
+
+  int active_tracklength_tallies_data(int i)
+  {return model::active_tracklength_tallies[i-1];}
+
+  int active_tracklength_tallies_size()
+  {return model::active_tracklength_tallies.size();}
+
+  int active_collision_tallies_data(int i)
+  {return model::active_collision_tallies[i-1];}
+
+  int active_collision_tallies_size()
+  {return model::active_collision_tallies.size();}
+
+  int tally_get_type_c(Tally* tally) {return tally->type_;}
+
+  void tally_set_type_c(Tally* tally, int type) {tally->type_ = type;}
+
   int tally_get_estimator_c(Tally* tally) {return tally->estimator_;}
 
   void tally_set_estimator_c(Tally* tally, int e) {tally->estimator_ = e;}

src/tallies/tally_header.F90

@@ -16,29 +16,6 @@ module tally_header
   use trigger_header,      only: TriggerObject
 
   implicit none
-  private
-  public :: allocate_tally_results
-  public :: free_memory_tally
-  public :: openmc_extend_tallies
-  public :: openmc_get_tally_index
-  public :: openmc_get_tally_next_id
-  public :: openmc_global_tallies
-  public :: openmc_tally_get_active
-  public :: openmc_tally_get_estimator
-  public :: openmc_tally_get_id
-  public :: openmc_tally_get_n_realizations
-  public :: openmc_tally_get_nuclides
-  public :: openmc_tally_get_scores
-  public :: openmc_tally_get_type
-  public :: openmc_tally_reset
-  public :: openmc_tally_results
-  public :: openmc_tally_set_active
-  public :: openmc_tally_set_estimator
-  public :: openmc_tally_set_filters
-  public :: openmc_tally_set_id
-  public :: openmc_tally_set_nuclides
-  public :: openmc_tally_set_scores
-  public :: openmc_tally_set_type
 
   interface
     function openmc_tally_set_filters(index, n, filter_indices) result(err) bind(C)
@@ -62,6 +39,50 @@ module tally_header
       logical(C_BOOL), intent(out) :: active
       integer(C_INT) :: err
     end function
+
+    function active_tallies_size() result(size) bind(C)
+      import C_INT
+      integer(C_INT) :: size
+    end function
+
+    function active_tallies_data(i) result(tally) bind(C)
+      import C_INT
+      integer(C_INT), value :: i
+      integer(C_INT) :: tally
+    end function
+
+    function active_analog_tallies_size() result(size) bind(C)
+      import C_INT
+      integer(C_INT) :: size
+    end function
+
+    function active_analog_tallies_data(i) result(tally) bind(C)
+      import C_INT
+      integer(C_INT), value :: i
+      integer(C_INT) :: tally
+    end function
+
+    function active_tracklength_tallies_size() result(size) bind(C)
+      import C_INT
+      integer(C_INT) :: size
+    end function
+
+    function active_tracklength_tallies_data(i) result(tally) bind(C)
+      import C_INT
+      integer(C_INT), value :: i
+      integer(C_INT) :: tally
+    end function
+
+    function active_collision_tallies_size() result(size) bind(C)
+      import C_INT
+      integer(C_INT) :: size
+    end function
+
+    function active_collision_tallies_data(i) result(tally) bind(C)
+      import C_INT
+      integer(C_INT), value :: i
+      integer(C_INT) :: tally
+    end function
   end interface
 
 !===============================================================================
@@ -77,10 +98,7 @@ module tally_header
 
     integer :: id                   ! user-defined identifier
     character(len=104) :: name = "" ! user-defined name
-    integer :: type = TALLY_VOLUME  ! volume, surface current
-    !integer :: estimator = ESTIMATOR_TRACKLENGTH ! collision, track-length
     real(8) :: volume               ! volume of region
-    !logical :: active = .false.
     logical :: depletion_rx = .false. ! has depletion reactions, e.g. (n,2n)
 
     ! Individual nuclides to tally
@@ -112,6 +130,8 @@ module tally_header
     procedure :: allocate_results => tally_allocate_results
     procedure :: read_results_hdf5 => tally_read_results_hdf5
     procedure :: write_results_hdf5 => tally_write_results_hdf5
+    procedure :: type => tally_get_type
+    procedure :: set_type => tally_set_type
     procedure :: estimator => tally_get_estimator
     procedure :: set_estimator => tally_set_estimator
     procedure :: n_filters => tally_get_n_filters
@@ -163,11 +183,7 @@ module tally_header
 !$omp&              global_tally_tracklength, global_tally_leakage)
 
   ! Active tally lists
-  type(VectorInt), public :: active_analog_tallies
-  type(VectorInt), public :: active_tracklength_tallies
   type(VectorInt), public :: active_meshsurf_tallies
-  type(VectorInt), public :: active_collision_tallies
-  type(VectorInt), public :: active_tallies
   type(VectorInt), public :: active_surface_tallies
 
   ! Normalization for statistics
@@ -298,6 +314,32 @@ contains
 
   end subroutine tally_allocate_results
 
+  function tally_get_type(this) result(t)
+    class(TallyObject) :: this
+    integer(C_INT) :: t
+    interface
+      function tally_get_type_c(tally) result(t) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: tally
+        integer(C_INT) :: t
+      end function
+    end interface
+    t = tally_get_type_c(this % ptr)
+  end function
+
+  subroutine tally_set_type(this, t)
+    class(TallyObject) :: this
+    integer(C_INT) :: t
+    interface
+      subroutine tally_set_type_c(tally, t) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: tally
+        integer(C_INT), value :: t
+      end subroutine
+    end interface
+    call tally_set_type_c(this % ptr, t)
+  end subroutine
+
   function tally_get_estimator(this) result(e)
     class(TallyObject) :: this
     integer(C_INT) :: e
@@ -513,12 +555,8 @@ contains
     if (allocated(global_tallies)) deallocate(global_tallies)
 
     ! Deallocate tally node lists
-    call active_analog_tallies % clear()
-    call active_tracklength_tallies % clear()
     call active_meshsurf_tallies % clear()
-    call active_collision_tallies % clear()
     call active_surface_tallies % clear()
-    call active_tallies % clear()
   end subroutine free_memory_tally
 
 !===============================================================================
@@ -701,22 +739,6 @@ contains
   end function openmc_tally_get_scores
 
 
-  function openmc_tally_get_type(index, type) result(err) bind(C)
-    ! Return the type of a tally
-    integer(C_INT32_T), value    :: index
-    integer(C_INT32_T), intent(out) :: type
-    integer(C_INT) :: err
-
-    if (index >= 1 .and. index <= size(tallies)) then
-      type = tallies(index) % obj % type
-      err = 0
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg('Index in tallies array is out of bounds.')
-    end if
-  end function openmc_tally_get_type
-
-
   function openmc_tally_reset(index) result(err) bind(C)
     ! Reset tally results and number of realizations
     integer(C_INT32_T), intent(in), value :: index
@@ -1048,37 +1070,6 @@ contains
   end function openmc_tally_set_scores
 
 
-  function openmc_tally_set_type(index, type) result(err) bind(C)
-    ! Update the type of a tally that is already allocated
-    integer(C_INT32_T), value, intent(in) :: index
-    character(kind=C_CHAR), intent(in) :: type(*)
-    integer(C_INT) :: err
-
-    character(:), allocatable :: type_
-
-    ! Convert C string to Fortran string
-    type_ = to_f_string(type)
-
-    err = 0
-    if (index >= 1 .and. index <= size(tallies)) then
-      select case (type_)
-      case ('volume')
-        tallies(index) % obj % type = TALLY_VOLUME
-      case ('mesh-surface')
-        tallies(index) % obj % type = TALLY_MESH_SURFACE
-      case ('surface')
-        tallies(index) % obj % type = TALLY_SURFACE
-      case default
-        err = E_INVALID_ARGUMENT
-        call set_errmsg("Unknown tally type: " // trim(type_))
-      end select
-    else
-      err = E_OUT_OF_BOUNDS
-      call set_errmsg("Index in tally array is out of bounds.")
-    end if
-  end function openmc_tally_set_type
-
-
   subroutine openmc_get_tally_next_id(id) bind(C)
     ! Returns an ID number that has not been used by any other tallies.
     integer(C_INT32_T), intent(out) :: id

src/tallies/trigger.F90

@@ -166,7 +166,7 @@ contains
           trigger % variance = ZERO
 
           ! Mesh current tally triggers require special treatment
-          if (t % type == TALLY_MESH_SURFACE) then
+          if (t % type() == TALLY_MESH_SURFACE) then
             call compute_tally_current(t, trigger)
 
           else

src/tracking.F90

@@ -86,7 +86,7 @@ contains
     endif
 
     ! Every particle starts with no accumulated flux derivative.
-    if (active_tallies % size() > 0) call zero_flux_derivs()
+    if (active_tallies_size() > 0) call zero_flux_derivs()
 
     EVENT_LOOP: do
       ! Set the random number stream
@@ -172,7 +172,7 @@ contains
       end do
 
       ! Score track-length tallies
-      if (active_tracklength_tallies % size() > 0) then
+      if (active_tracklength_tallies_size() > 0) then
         call score_tracklength_tally(p, distance)
       end if
 
@@ -183,7 +183,7 @@ contains
       end if
 
       ! Score flux derivative accumulators for differential tallies.
-      if (active_tallies % size() > 0) call score_track_derivative(p, distance)
+      if (active_tallies_size() > 0) call score_track_derivative(p, distance)
 
       if (d_collision > d_boundary) then
         ! ====================================================================
@@ -241,8 +241,8 @@ contains
         ! Score collision estimator tallies -- this is done after a collision
         ! has occurred rather than before because we need information on the
         ! outgoing energy for any tallies with an outgoing energy filter
-        if (active_collision_tallies % size() > 0) call score_collision_tally(p)
-        if (active_analog_tallies % size() > 0) call score_analog_tally(p)
+        if (active_collision_tallies_size() > 0) call score_collision_tally(p)
+        if (active_analog_tallies_size() > 0) call score_analog_tally(p)
 
         ! Reset banked weight during collision
         p % n_bank   = 0
@@ -273,7 +273,7 @@ contains
         end do
 
         ! Score flux derivative accumulators for differential tallies.
-        if (active_tallies % size() > 0) call score_collision_derivative(p)
+        if (active_tallies_size() > 0) call score_collision_derivative(p)
       end if
 
       ! If particle has too many events, display warning and kill it

tests/unit_tests/test_capi.py

@@ -162,6 +162,7 @@ def test_tally_mapping(capi_init):
 
 def test_tally(capi_init):
     t = openmc.capi.tallies[1]
+    assert t.type == 'volume'
     assert len(t.filters) == 2
     assert isinstance(t.filters[0], openmc.capi.MaterialFilter)
     assert isinstance(t.filters[1], openmc.capi.EnergyFilter)