crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move distribcell offsets to C++

Browse source
This commit is contained in:
Sterling Harper 2018-05-21 16:32:44 -04:00
parent a686ad33b2
commit 9da5e36daf
11 changed files with 265 additions and 185 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/cell.cpp
View file

@ -485,6 +485,8 @@ extern "C" {
i_surf = out.second;
}
int32_t cell_offset(Cell *c, int map) {return c->offset[map];}
void cell_to_hdf5(Cell *c, hid_t group) {c->to_hdf5(group);}
void extend_cells_c(int32_t n)

2
src/cell.h
View file

@ -71,6 +71,8 @@ public:
std::vector<std::int32_t> rpn;
bool simple; //!< Does the region contain only intersections?
std::vector<int32_t> offset; //!< Distribcell offset table
Cell() {};
explicit Cell(pugi::xml_node cell_node);

92
src/geometry.F90
View file

@ -26,7 +26,7 @@ module geometry
end function cell_contains_c
function count_universe_instances(search_univ, target_univ_id) bind(C) &
result(count)
result(count)
import C_INT32_T, C_INT
integer(C_INT32_T), intent(in), value :: search_univ
integer(C_INT32_T), intent(in), value :: target_univ_id
@ -174,19 +174,19 @@ contains
offset = 0
do k = 1, p % n_coord
if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then
offset = offset + cells(p % coord(k) % cell) % &
offset(distribcell_index)
offset = offset + cells(p % coord(k) % cell) &
% offset(distribcell_index-1)
elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then
if (lattices(p % coord(k + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])) then
offset = offset + lattices(p % coord(k + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z)
offset = offset + lattices(p % coord(k + 1) % lattice) % obj &
% offset(distribcell_index - 1, &
[p % coord(k + 1) % lattice_x - 1, &
p % coord(k + 1) % lattice_y - 1, &
p % coord(k + 1) % lattice_z - 1])
end if
end if
end do
@ -534,82 +534,6 @@ contains
end subroutine neighbor_lists
!===============================================================================
! CALC_OFFSETS calculates and stores the offsets in all fill cells. This
! routine is called once upon initialization.
!===============================================================================
subroutine calc_offsets(univ_id, map, univ)
integer, intent(in) :: univ_id ! target universe ID
integer, intent(in) :: map ! map index in vector of maps
type(Universe), intent(in) :: univ ! universe searching in
integer :: i ! index over cells
integer :: j, k, m ! indices in lattice
integer :: offset ! total offset for a given cell
integer :: cell_index ! index in cells array
type(Cell), pointer :: c ! pointer to current cell
class(Lattice), pointer :: lat ! pointer to current lattice
offset = 0
do i = 1, size(univ % cells)
cell_index = univ % cells(i)
! get pointer to cell
c => cells(cell_index)
! ====================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
if (c % type() == FILL_MATERIAL) then
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
elseif (c % type() == FILL_UNIVERSE) then
c % offset(map) = offset
offset = offset + count_universe_instances(c % fill - 1, univ_id)
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
lat => lattices(c % fill) % obj
select type (lat)
type is (RectLattice)
do m = 1, lat % n_cells(3)
do k = 1, lat % n_cells(2)
do j = 1, lat % n_cells(1)
lat % offset(map, j, k, m) = offset
offset = offset + count_universe_instances(&
lat % get([j-1, k-1, m-1]), univ_id)
end do
end do
end do
type is (HexLattice)
do m = 1, lat % n_axial
do k = 1, 2*lat % n_rings - 1
do j = 1, 2*lat % n_rings - 1
if (lat % are_valid_indices([j, k, m])) then
lat % offset(map, j, k, m) = offset
offset = offset + count_universe_instances(&
lat % get([j-1, k-1, m-1]), univ_id)
end if
end do
end do
end do
end select
end if
end do
end subroutine calc_offsets
!===============================================================================
! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in
! the geometry

68
src/geometry_aux.cpp
View file

@ -1,5 +1,4 @@
//! \file geometry_aux.cpp
//! Auxilary functions for geometry initialization and general data handling.
#include "geometry_aux.h"
#include <sstream>
#include <unordered_set>
@ -14,11 +13,9 @@
namespace openmc {
//==============================================================================
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
extern "C" void
void
adjust_indices_c()
{
// Adjust material/fill idices.
@ -65,15 +62,9 @@ adjust_indices_c()
}
}
//==============================================================================
//! Figure out which Universe is the root universe.
//!
//! This function looks for a universe that is not listed in a Cell::fill or in
//! a Lattice.
//! @return The index of the root universe.
//==============================================================================
extern "C" int32_t
int32_t
find_root_universe()
{
// Find all the universes listed as a cell fill.
@ -115,15 +106,24 @@ find_root_universe()
}
//==============================================================================
//! Recursively search through the geometry and count cell instances.
//!
//! This function will update the Cell::n_instances value for each cell in the
//! geometry.
//! @param univ_indx The index of the universe to begin searching from (probably
//! the root universe).
void
allocate_offset_tables(int n_maps)
{
for (Cell *c : cells_c) {
if (c->type != FILL_MATERIAL) {
c->offset.resize(n_maps, C_NONE);
}
}
for (Lattice *lat : lattices_c) {
lat->allocate_offset_table(n_maps);
}
}
//==============================================================================
extern "C" void
void
count_cell_instances(int32_t univ_indx)
{
for (int32_t cell_indx : universes_c[univ_indx]->cells) {
@ -144,14 +144,9 @@ count_cell_instances(int32_t univ_indx)
}
}
//==============================================================================
//! Recursively search through universes and count the number of instances of
//! the target universe in the geometry tree.
//! @param search_univ The index of the universe to begin searching from.
//! @param target_univ_id The ID of the universe to be counted.
//==============================================================================
extern "C" int
int
count_universe_instances(int32_t search_univ, int32_t target_univ_id)
{
// If this is the target, it can't contain itself.
@ -179,4 +174,27 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id)
return count;
}
//==============================================================================
void
fill_offset_tables(int32_t target_univ_id, int map)
{
for (Universe *univ : universes_c) {
int32_t offset {0}; // TODO: is this a bug? It matches F90 implementation.
for (int32_t cell_indx : univ->cells) {
Cell &c = *cells_c[cell_indx];
if (c.type == FILL_UNIVERSE) {
c.offset[map] = offset;
int32_t search_univ = c.fill - 1; // TODO: off-by-one
offset += count_universe_instances(search_univ, target_univ_id);
} else if (c.type == FILL_LATTICE) {
Lattice &lat = *lattices_c[c.fill - 1]; // TODO: off-by-one
offset = lat.fill_offset_table(offset, target_univ_id, map);
}
}
}
}
} // namespace openmc

66
src/geometry_aux.h Normal file
View file

@ -0,0 +1,66 @@
//! \file geometry_aux.h
//! Auxilary functions for geometry initialization and general data handling.
#ifndef GEOMETRY_AUX_H
#define GEOMETRY_AUX_H
#include <cstdint>
namespace openmc {
//==============================================================================
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
extern "C" void adjust_indices_c();
//==============================================================================
//! Figure out which Universe is the root universe.
//!
//! This function looks for a universe that is not listed in a Cell::fill or in
//! a Lattice.
//! @return The index of the root universe.
//==============================================================================
extern "C" int32_t find_root_universe();
//==============================================================================
//! Allocate storage in Lattice and Cell objects for distribcell offset tables.
//==============================================================================
extern "C" void allocate_offset_tables(int n_maps);
//==============================================================================
//! Recursively search through the geometry and count cell instances.
//!
//! This function will update the Cell::n_instances value for each cell in the
//! geometry.
//! @param univ_indx The index of the universe to begin searching from (probably
//! the root universe).
//==============================================================================
extern "C" void count_cell_instances(int32_t univ_indx);
//==============================================================================
//! Recursively search through universes and count universe instances.
//! @param search_univ The index of the universe to begin searching from.
//! @param target_univ_id The ID of the universe to be counted.
//! @return The number of instances of target_univ_id in the geometry tree under
//! search_univ.
//==============================================================================
extern "C" int
count_universe_instances(int32_t search_univ, int32_t target_univ_id);
//==============================================================================
//! Populate Cell and Lattice distribcell offset tables.
//! @param target_univ_id The ID of the universe to be counted.
//! @param map The index of the distribcell map that defines the offsets for the
//! target universe.
//==============================================================================
extern "C" void fill_offset_tables(int32_t target_univ_id, int map);
} // namespace openmc
#endif // GEOMETRY_AUX_H

49
src/geometry_header.F90
View file

@ -68,7 +68,7 @@ module geometry_header
end subroutine cell_set_universe_c
function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') &
result(n_instances)
result(n_instances)
import C_PTR, C_INT32_T
implicit none
type(C_PTR), intent(in), value :: cell_ptr
@ -94,9 +94,16 @@ module geometry_header
integer(C_INT32_T), intent(out) :: i_surf
end subroutine cell_distance_c
function cell_offset_c(cell_ptr, map) bind(C, name="cell_offset") &
result(offset)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), intent(in), value :: cell_ptr
integer(C_INT), intent(in), value :: map
integer(C_INT32_T) :: offset
end function cell_offset_c
subroutine cell_to_hdf5_c(cell_ptr, group) bind(C, name='cell_to_hdf5')
import HID_T, C_PTR
implicit none
type(C_PTR), intent(in), value :: cell_ptr
integer(HID_T), intent(in), value :: group
end subroutine cell_to_hdf5_c
@ -104,14 +111,12 @@ module geometry_header
function lattice_pointer_c(lat_ind) bind(C, name='lattice_pointer') &
result(ptr)
import C_PTR, C_INT32_T
implicit none
integer(C_INT32_T), intent(in), value :: lat_ind
type(C_PTR) :: ptr
end function lattice_pointer_c
function lattice_id_c(lat_ptr) bind(C, name='lattice_id') result(id)
import C_PTR, C_INT32_T
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT32_T) :: id
end function lattice_id_c
@ -119,7 +124,6 @@ module geometry_header
function lattice_are_valid_indices_c(lat_ptr, i_xyz) &
bind(C, name='lattice_are_valid_indices') result (is_valid)
import C_PTR, C_INT, C_BOOL
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT), intent(in) :: i_xyz(3)
logical(C_BOOL) :: is_valid
@ -128,7 +132,6 @@ module geometry_header
subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) &
bind(C, name='lattice_distance')
import C_PTR, C_INT, C_DOUBLE
implicit none
type(C_PTR), intent(in), value :: lat_ptr
real(C_DOUBLE), intent(in) :: xyz(3)
real(C_DOUBLE), intent(in) :: uvw(3)
@ -140,7 +143,6 @@ module geometry_header
subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) &
bind(C, name='lattice_get_indices')
import C_PTR, C_INT, C_DOUBLE
implicit none
type(C_PTR), intent(in), value :: lat_ptr
real(C_DOUBLE), intent(in) :: xyz(3)
integer(C_INT), intent(out) :: i_xyz(3)
@ -149,7 +151,6 @@ module geometry_header
subroutine lattice_get_local_xyz_c(lat_ptr, global_xyz, i_xyz, local_xyz) &
bind(C, name='lattice_get_local_xyz')
import C_PTR, C_INT, C_DOUBLE
implicit none
type(C_PTR), intent(in), value :: lat_ptr
real(C_DOUBLE), intent(in) :: global_xyz(3)
integer(C_INT), intent(in) :: i_xyz(3)
@ -158,11 +159,19 @@ module geometry_header
subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5')
import HID_T, C_PTR
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(HID_T), intent(in), value :: group
end subroutine lattice_to_hdf5_c
function lattice_offset_c(lat_ptr, map, i_xyz) &
bind(C, name='lattice_offset') result(offset)
import C_PTR, C_INT, C_INT32_T
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT), intent(in), value :: map
integer(C_INT), intent(in) :: i_xyz(3)
integer(C_INT32_T) :: offset
end function lattice_offset_c
function lattice_outer_c(lat_ptr) bind(C, name='lattice_outer') &
result(outer)
import C_PTR, C_INT32_T
@ -173,7 +182,6 @@ module geometry_header
function lattice_universe_c(lat_ptr, i_xyz) &
bind(C, name='lattice_universe') result(univ)
import C_PTR, C_INT32_t, C_INT
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT), intent(in) :: i_xyz(3)
integer(C_INT32_T) :: univ
@ -181,7 +189,6 @@ module geometry_header
subroutine extend_cells_c(n) bind(C)
import C_INT32_t
implicit none
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_cells_c
end interface
@ -207,7 +214,6 @@ module geometry_header
type(C_PTR) :: ptr
logical :: is_3d ! Lattice has cells on z axis
integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets
contains
procedure :: id => lattice_id
@ -216,6 +222,7 @@ module geometry_header
procedure :: get => lattice_get
procedure :: get_indices => lattice_get_indices
procedure :: get_local_xyz => lattice_get_local_xyz
procedure :: offset => lattice_offset
procedure :: outer => lattice_outer
procedure :: to_hdf5 => lattice_to_hdf5
end type Lattice
@ -256,8 +263,6 @@ module geometry_header
integer, allocatable :: material(:) ! Material within cell. Multiple
! materials for distribcell
! instances. 0 signifies a universe
integer, allocatable :: offset(:) ! Distribcell offset for tally
! counter
integer, allocatable :: region(:) ! Definition of spatial region as
! Boolean expression of half-spaces
integer :: distribcell_index ! Index corresponding to this cell in
@ -282,6 +287,7 @@ module geometry_header
procedure :: n_instances => cell_n_instances
procedure :: simple => cell_simple
procedure :: distance => cell_distance
procedure :: offset => cell_offset
procedure :: to_hdf5 => cell_to_hdf5
end type Cell
@ -350,6 +356,14 @@ contains
call lattice_get_local_xyz_c(this % ptr, global_xyz, i_xyz, local_xyz)
end function lattice_get_local_xyz
function lattice_offset(this, map, i_xyz) result(offset)
class(Lattice), intent(in) :: this
integer(C_INT), intent(in) :: map
integer(C_INT), intent(in) :: i_xyz(3)
integer(C_INT32_T) :: offset
offset = lattice_offset_c(this % ptr, map, i_xyz)
end function lattice_offset
function lattice_outer(this) result(outer)
class(Lattice), intent(in) :: this
integer(C_INT32_T) :: outer
@ -422,6 +436,13 @@ contains
call cell_distance_c(this % ptr, xyz, uvw, on_surface, min_dist, i_surf)
end subroutine cell_distance
function cell_offset(this, map) result(offset)
class(Cell), intent(in) :: this
integer(C_INT), intent(in) :: map
integer(C_INT32_T) :: offset
offset = cell_offset_c(this % ptr, map)
end function cell_offset
subroutine cell_to_hdf5(this, group)
class(Cell), intent(in) :: this
integer(HID_T), intent(in) :: group

44
src/input_xml.F90
View file

@ -11,7 +11,7 @@ module input_xml
use distribution_univariate
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message, openmc_err_msg
use geometry, only: calc_offsets, maximum_levels, neighbor_lists
use geometry, only: maximum_levels, neighbor_lists
use geometry_header
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
@ -50,10 +50,21 @@ module input_xml
subroutine adjust_indices_c() bind(C)
end subroutine adjust_indices_c
subroutine allocate_offset_tables(n_maps) bind(C)
import C_INT
integer(C_INT), intent(in), value :: n_maps
end subroutine allocate_offset_tables
subroutine fill_offset_tables(target_univ_id, map) bind(C)
import C_INT32_T, C_INT
integer(C_INT32_T), intent(in), value :: target_univ_id
integer(C_INT), intent(in), value :: map
end subroutine fill_offset_tables
subroutine count_cell_instances(univ_indx) bind(C)
import C_INT32_T
integer(C_INT32_T), intent(in), value :: univ_indx
end subroutine
end subroutine count_cell_instances
subroutine read_surfaces(node_ptr) bind(C)
import C_PTR
@ -3884,12 +3895,11 @@ contains
! Allocate offset maps at each level in the geometry
call allocate_offsets(univ_list)
call allocate_offset_tables(n_maps)
! Calculate offsets for each target distribcell
do i = 1, n_maps
do j = 1, n_universes
call calc_offsets(univ_list(i), i, universes(j))
end do
call fill_offset_tables(univ_list(i), i-1)
end do
end subroutine prepare_distribcell
@ -3952,30 +3962,6 @@ contains
end do
end do
! Allocate the offset tables for lattices
do i = 1, n_lattices
associate(lat => lattices(i) % obj)
select type(lat)
type is (RectLattice)
allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), &
lat % n_cells(3)))
type is (HexLattice)
allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, &
2 * lat % n_rings - 1, lat % n_axial))
end select
lat % offset(:, :, :, :) = 0
end associate
end do
! Allocate offset table for fill cells
do i = 1, n_cells
if (cells(i) % type() /= FILL_MATERIAL) then
allocate(cells(i) % offset(n_maps))
end if
end do
! Free up memory
call cell_list % clear()

48
src/lattice.cpp
View file

@ -7,6 +7,7 @@
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "geometry_aux.h"
#include "hdf5_interface.h"
#include "xml_interface.h"
@ -98,6 +99,26 @@ Lattice::adjust_indices()
//==============================================================================
void
Lattice::allocate_offset_table(int n_maps)
{
offsets.resize(n_maps * universes.size(), C_NONE);
}
//==============================================================================
int32_t
Lattice::fill_offset_table(int32_t offset, int32_t target_univ_id, int map)
{
for (auto it = begin(); it != end(); ++it) {
offsets[map * universes.size() + it.indx] = offset;
offset += count_universe_instances(*it, target_univ_id);
}
return offset;
}
//==============================================================================
void
Lattice::to_hdf5(hid_t lattices_group) const
{
@ -316,6 +337,17 @@ RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const
//==============================================================================
int32_t&
RectLattice::offset(int map, const int i_xyz[3])
{
int nx = n_cells[0];
int ny = n_cells[1];
int nz = n_cells[2];
return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]];
}
//==============================================================================
void
RectLattice::to_hdf5_inner(hid_t lat_group) const
{
@ -763,6 +795,17 @@ HexLattice::is_valid_index(int indx) const
//==============================================================================
int32_t&
HexLattice::offset(int map, const int i_xyz[3])
{
int nx = 2*n_rings - 1;
int ny = 2*n_rings - 1;
int nz = n_axial;
return offsets[nx*ny*nz*map + nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]];
}
//==============================================================================
void
HexLattice::to_hdf5_inner(hid_t lat_group) const
{
@ -874,6 +917,11 @@ extern "C" {
local_xyz[2] = xyz[2];
}
int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3])
{
return lat->offset(map, i_xyz);
}
int32_t lattice_outer(Lattice *lat) {return lat->outer;}
void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}

40
src/lattice.h
View file

@ -41,11 +41,11 @@ class LatticeIter;
class Lattice
{
public:
int32_t id; //! Universe ID number
std::string name; //! User-defined name
std::vector<int32_t> universes; //! Universes filling each lattice tile
int32_t outer{NO_OUTER_UNIVERSE}; //! Universe tiled outside the lattice
//std::vector<int32_t> offset; //! Distribcell offsets
int32_t id; //!< Universe ID number
std::string name; //!< User-defined name
std::vector<int32_t> universes; //!< Universes filling each lattice tile
int32_t outer{NO_OUTER_UNIVERSE}; //!< Universe tiled outside the lattice
std::vector<int32_t> offsets; //!< Distribcell offset table
explicit Lattice(pugi::xml_node lat_node);
@ -59,6 +59,11 @@ public:
//! Convert internal universe values from IDs to indices using universe_dict.
void adjust_indices();
//! Allocate offset table for distribcell.
void allocate_offset_table(int n_maps);
int32_t fill_offset_table(int32_t offset, int32_t target_univ_id, int map);
//! Check lattice indices.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return true if the given indices fit within the lattice bounds. False
@ -96,6 +101,13 @@ public:
return (indx > 0) && (indx < universes.size());
}
//! Get the distribcell offset for a lattice tile.
//! @param The map index for the target cell.
//! @param i_xyz[3] The indices for a lattice tile.
//! @return Distribcell offset i.e. the largest instance number for the target
//! cell found in the geometry tree under this lattice tile.
virtual int32_t& offset(int map, const int i_xyz[3]) = 0;
//! Write all information needed to reconstruct the lattice to an HDF5 group.
//! @param group_id An HDF5 group id.
void to_hdf5(hid_t group_id) const;
@ -166,12 +178,14 @@ public:
std::array<double, 3>
get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;
int32_t& offset(int map, const int i_xyz[3]);
void to_hdf5_inner(hid_t group_id) const;
protected:
std::array<int, 3> n_cells; //! Number of cells along each axis
std::array<double, 3> lower_left; //! Global lower-left corner of the lattice
std::array<double, 3> pitch; //! Lattice tile width along each axis
std::array<int, 3> n_cells; //!< Number of cells along each axis
std::array<double, 3> lower_left; //!< Global lower-left corner of the lattice
std::array<double, 3> pitch; //!< Lattice tile width along each axis
};
//==============================================================================
@ -200,13 +214,15 @@ public:
bool is_valid_index(int indx) const;
int32_t& offset(int map, const int i_xyz[3]);
void to_hdf5_inner(hid_t group_id) const;
protected:
int n_rings; //! Number of radial tile positions
int n_axial; //! Number of axial tile positions
std::array<double, 3> center; //! Global center of lattice
std::array<double, 2> pitch; //! Lattice tile width and height
int n_rings; //!< Number of radial tile positions
int n_axial; //!< Number of axial tile positions
std::array<double, 3> center; //!< Global center of lattice
std::array<double, 2> pitch; //!< Lattice tile width and height
};
} // namespace openmc

5
src/summary.F90
View file

@ -226,11 +226,6 @@ contains
case (FILL_UNIVERSE)
call write_dataset(cell_group, "fill_type", "universe")
call write_dataset(cell_group, "fill", universes(c%fill)%id)
if (allocated(c%offset)) then
if (size(c%offset) > 0) then
call write_dataset(cell_group, "offset", c%offset)
end if
end if
if (allocated(c%translation)) then
call write_dataset(cell_group, "translation", c%translation)

34
src/tallies/tally_filter_distribcell.F90
View file

@ -59,19 +59,19 @@ contains
offset = 0
do i = 1, p % n_coord
if (cells(p % coord(i) % cell) % type() == FILL_UNIVERSE) then
offset = offset + cells(p % coord(i) % cell) % &
offset(distribcell_index)
offset = offset + cells(p % coord(i) % cell) &
% offset(distribcell_index-1)
elseif (cells(p % coord(i) % cell) % type() == FILL_LATTICE) then
if (lattices(p % coord(i + 1) % lattice) % obj &
% are_valid_indices([&
p % coord(i + 1) % lattice_x, &
p % coord(i + 1) % lattice_y, &
p % coord(i + 1) % lattice_z])) then
offset = offset + lattices(p % coord(i + 1) % lattice) % obj % &
offset(distribcell_index, &
p % coord(i + 1) % lattice_x, &
p % coord(i + 1) % lattice_y, &
p % coord(i + 1) % lattice_z)
offset = offset + lattices(p % coord(i + 1) % lattice) % obj &
% offset(distribcell_index - 1, &
[p % coord(i + 1) % lattice_x - 1, &
p % coord(i + 1) % lattice_y - 1, &
p % coord(i + 1) % lattice_z - 1])
end if
end if
if (this % cell == p % coord(i) % cell) then
@ -211,12 +211,12 @@ contains
! Two cases, lattice or fill cell
if (c % type() == FILL_UNIVERSE) then
temp_offset = c % offset(map)
temp_offset = c % offset(map-1)
! Get the offset of the first lattice location
else
lat => lattices(c % fill) % obj
temp_offset = lat % offset(map, 1, 1, 1)
temp_offset = lat % offset(map-1, [0, 0, 0])
end if
! If the final offset is in the range of offset - temp_offset+offset
@ -248,7 +248,7 @@ contains
if (c % type() == FILL_UNIVERSE) then
! Enter this cell to update the current offset
offset = c % offset(map) + offset
offset = c % offset(map-1) + offset
next_univ => universes(c % fill)
call find_offset(i_cell, next_univ, target_offset, offset, path)
@ -282,10 +282,11 @@ contains
do l = 1, n_y
do k = 1, n_x
if (target_offset >= lat % offset(map, k, l, m) + offset) then
if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) &
+ offset) then
if (k == n_x .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
lat_offset = lat % offset(map-1, [k-1, l-1, m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
if (lat % is_3d) then
@ -306,7 +307,7 @@ contains
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
if (lat % is_3d) then
@ -352,10 +353,11 @@ contains
cycle
end if
if (target_offset >= lat % offset(map, k, l, m) + offset) then
if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) &
+ offset) then
if (k == lat % n_rings .and. l == n_y .and. m == n_z) then
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
lat_offset = lat % offset(map-1, [k-1, l-1, m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
if (lat % is_3d) then
@ -378,7 +380,7 @@ contains
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
if (lat % is_3d) then