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

Remove Fortran lattice % universes

Browse source
This commit is contained in:
Sterling Harper 2018-04-29 20:23:41 -04:00
parent 735c59650b
commit 98009a268d
10 changed files with 152 additions and 168 deletions
Download patch file Download diff file Expand all files Collapse all files

19
src/cell.cpp
View file

@ -8,6 +8,7 @@
#include "constants.h"
#include "error.h"
#include "hdf5_interface.h"
#include "lattice.h"
#include "surface.h"
#include "xml_interface.h"
@ -29,6 +30,19 @@ constexpr int32_t OP_UNION {std::numeric_limits<int32_t>::max() - 4};
extern "C" double FP_PRECISION;
//==============================================================================
// Global variables
//==============================================================================
// Braces force n_cells to be defined here, not just declared.
extern "C" {int32_t n_cells {0};}
std::vector<Cell*> cells_c;
std::map<int32_t, int32_t> cell_dict;
std::vector<Universe*> universes_c;
std::map<int32_t, int32_t> universe_dict;
//==============================================================================
//! Convert region specification string to integer tokens.
//!
@ -444,6 +458,11 @@ adjust_indices_c()
fatal_error(err_msg);
}
}
// Change all lattice universe values from IDs to indices.
for (Lattice *l : lattices_c) {
l->adjust_indices();
}
}
//==============================================================================

11
src/cell.h
View file

@ -16,16 +16,15 @@ namespace openmc {
// Global variables
//==============================================================================
// Braces force n_cells to be defined here, not just declared.
extern "C" {int32_t n_cells {0};}
extern "C" int32_t n_cells;
class Cell;
std::vector<Cell*> cells_c;
std::map<int32_t, int32_t> cell_dict;
extern std::vector<Cell*> cells_c;
extern std::map<int32_t, int32_t> cell_dict;
class Universe;
std::vector<Universe*> universes_c;
std::map<int32_t, int32_t> universe_dict;
extern std::vector<Universe*> universes_c;
extern std::map<int32_t, int32_t> universe_dict;
//==============================================================================
//! A geometry primitive that fills all space and contains cells.

3
src/constants.h
View file

@ -5,7 +5,8 @@
namespace openmc{
extern "C" double FP_COINCIDENT;
constexpr double INFTY{std::numeric_limits<double>::max()};
extern "C" double FP_PRECISION;
constexpr double INFTY {std::numeric_limits<double>::max()};
constexpr int C_NONE {-1};
} // namespace openmc

22
src/geometry.F90
View file

@ -254,7 +254,7 @@ contains
if (lat % are_valid_indices(i_xyz)) then
! Particle is inside the lattice.
p % coord(j + 1) % universe = &
lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
lat % get([i_xyz(1)-1, i_xyz(2)-1, i_xyz(3)-1]) + 1
else
! Particle is outside the lattice.
@ -331,7 +331,8 @@ contains
else OUTSIDE_LAT
! Find cell in next lattice element
p % coord(j) % universe = lat % universes(i_xyz(1), i_xyz(2), i_xyz(3))
p % coord(j) % universe = &
lat % get([i_xyz(1)-1, i_xyz(2)-1, i_xyz(3)-1]) + 1
call find_cell(p, found)
if (.not. found) then
@ -591,7 +592,7 @@ contains
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
offset = offset + &
count_target(next_univ, counts, found, univ_id, map)
end do
@ -612,7 +613,7 @@ contains
cycle
else
lat % offset(map, j, k, m) = offset
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
offset = offset + &
count_target(next_univ, counts, found, univ_id, map)
end if
@ -707,7 +708,7 @@ contains
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
! Found target - stop since target cannot contain itself
if (next_univ % id == univ_id) then
@ -735,7 +736,7 @@ contains
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
! Found target - stop since target cannot contain itself
if (next_univ % id == univ_id) then
@ -799,7 +800,7 @@ contains
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
call count_instance(universes(lat % universes(j, k, m)))
call count_instance(universes(lat % get([j-1, k-1, m-1])+1))
end do
end do
end do
@ -817,7 +818,8 @@ contains
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
call count_instance(universes(lat % universes(j, k, m)))
call count_instance(universes(lat % get([j-1, k-1, m-1]) &
+ 1))
end if
end do
end do
@ -874,7 +876,7 @@ contains
do j = 1, lat % n_cells(1)
do k = 1, lat % n_cells(2)
do m = 1, lat % n_cells(3)
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
levels_below = max(levels_below, maximum_levels(next_univ))
end do
end do
@ -893,7 +895,7 @@ contains
else if (j + k > 3*lat % n_rings - 1) then
cycle
else
next_univ => universes(lat % universes(j, k, m))
next_univ => universes(lat % get([j-1, k-1, m-1]) + 1)
levels_below = max(levels_below, maximum_levels(next_univ))
end if
end do

19
src/geometry_header.F90
View file

@ -148,6 +148,15 @@ module geometry_header
implicit none
integer(C_INT32_T), intent(in), value :: n
end subroutine extend_cells_c
function lattice_universe_c(lat_ptr, i_xyz) &
bind(C, name='lattice_universe') result(univ)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: lat_ptr
integer(C_INT), intent(in) :: i_xyz(3)
integer(C_INT32_T) :: univ
end function lattice_universe_c
end interface
!===============================================================================
@ -170,9 +179,7 @@ module geometry_header
type, abstract :: Lattice
type(C_PTR) :: ptr
!integer :: id ! Universe number for lattice
real(8), allocatable :: pitch(:) ! Pitch along each axis
integer, allocatable :: universes(:,:,:) ! Specified universes
integer :: outside ! Material to fill area outside
integer :: outer ! universe to tile outside the lat
logical :: is_3d ! Lattice has cells on z axis
@ -182,6 +189,7 @@ module geometry_header
procedure :: id => lattice_id
procedure :: are_valid_indices => lattice_are_valid_indices
procedure :: distance => lattice_distance
procedure :: get => lattice_get
procedure :: get_indices => lattice_get_indices
procedure :: get_local_xyz => lattice_get_local_xyz
procedure :: to_hdf5 => lattice_to_hdf5
@ -297,6 +305,13 @@ contains
call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans)
end subroutine lattice_distance
function lattice_get(this, i_xyz) result(univ)
class(Lattice), intent(in) :: this
integer(C_INT), intent(in) :: i_xyz(3)
integer(C_INT32_T) :: univ
univ = lattice_universe_c(this % ptr, i_xyz)
end function lattice_get
function lattice_get_indices(this, xyz) result(i_xyz)
class(Lattice), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3)

145
src/input_xml.F90
View file

@ -1331,7 +1331,6 @@ contains
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
n_z = lat % n_cells(3)
allocate(lat % universes(n_x, n_y, n_z))
! Check that number of universes matches size
n = node_word_count(node_lat, "universes")
@ -1340,21 +1339,16 @@ contains
&size of lattice " // trim(to_str(lat % id())) // ".")
end if
allocate(temp_int_array(n))
call get_node_array(node_lat, "universes", temp_int_array)
! Read universes
do m = 1, n_z
do m = 0, n_z-1
do k = 0, n_y - 1
do j = 1, n_x
lat % universes(j, n_y - k, m) = &
temp_int_array(j + n_x*k + n_x*n_y*(m-1))
if (find(fill_univ_ids, lat % universes(j, n_y - k, m)) == -1) &
call fill_univ_ids % push_back(lat % universes(j, n_y - k, m))
do j = 0, n_x - 1
univ_id = lat % get([j, k, m])
if (find(fill_univ_ids, univ_id) == -1) &
call fill_univ_ids % push_back(univ_id)
end do
end do
end do
deallocate(temp_int_array)
! Read outer universe for area outside lattice.
lat % outer = NO_OUTER_UNIVERSE
@ -1427,7 +1421,6 @@ contains
! Copy number of dimensions
n_rings = lat % n_rings
n_z = lat % n_axial
allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z))
! Check that number of universes matches size
n = node_word_count(node_lat, "universes")
@ -1440,91 +1433,15 @@ contains
call get_node_array(node_lat, "universes", temp_int_array)
! Read universes
! Universes in hexagonal lattices are stored in a manner that represents
! a skewed coordinate system: (x, alpha) rather than (x, y). There is
! no obvious, direct relationship between the order of universes in the
! input and the order that they will be stored in the skewed array so
! the following code walks a set of index values across the skewed array
! in a manner that matches the input order. Note that i_x = 0, i_a = 0
! corresponds to the center of the hexagonal lattice.
input_index = 1
do m = 1, n_z
! Initialize lattice indecies.
i_x = 1
i_a = n_rings - 1
! Map upper triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
do j = 1, k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
do m = 0, n_z-1
do k = 0, 2*n_rings-2
do j = 0, 2*n_rings-2
univ_id = lat % get([j, k, m])
if (find(fill_univ_ids, univ_id) == -1) &
call fill_univ_ids % push_back(univ_id)
end do
! Return lattice index to start of current row.
i_x = i_x - 2*k
i_a = i_a + k
end do
! Map middle square region of hexagonal lattice.
do k = 1, 2*n_rings - 1
if (mod(k, 2) == 1) then
! Walk index to lower-left neighbor of last row start.
i_x = i_x - 1
else
! Walk index to lower-right neighbor of last row start
i_x = i_x + 1
i_a = i_a - 1
end if
do j = 1, n_rings - mod(k-1, 2)
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - mod(k-1, 2))
i_a = i_a + n_rings - mod(k-1, 2)
end do
! Map lower triangular region of hexagonal lattice.
do k = 1, n_rings-1
! Walk index to lower-right neighbor of last row start.
i_x = i_x + 1
i_a = i_a - 1
do j = 1, n_rings - k
! Place universe in array.
lat % universes(i_x + n_rings, i_a + n_rings, m) = &
temp_int_array(input_index)
if (find(fill_univ_ids, temp_int_array(input_index)) == -1) &
call fill_univ_ids % push_back(temp_int_array(input_index))
! Walk index to closest non-adjacent right neighbor.
i_x = i_x + 2
i_a = i_a - 1
! Increment XML array index.
input_index = input_index + 1
end do
! Return lattice index to start of current row.
i_x = i_x - 2*(n_rings - k)
i_a = i_a + n_rings - k
end do
end do
deallocate(temp_int_array)
! Read outer universe for area outside lattice.
lat % outer = NO_OUTER_UNIVERSE
@ -4312,46 +4229,6 @@ contains
do i = 1, n_lattices
lat => lattices(i) % 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)
id = lat % universes(j,k,m)
if (universe_dict % has(id)) then
lat % universes(j,k,m) = universe_dict % get(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat % id())))
end if
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 (j + k < lat % n_rings + 1) then
cycle
else if (j + k > 3*lat % n_rings - 1) then
cycle
end if
id = lat % universes(j, k, m)
if (universe_dict % has(id)) then
lat % universes(j, k, m) = universe_dict % get(id)
else
call fatal_error("Invalid universe number " &
&// trim(to_str(id)) // " specified on lattice " &
&// trim(to_str(lat % id())))
end if
end do
end do
end do
end select
if (lat % outer /= NO_OUTER_UNIVERSE) then
if (universe_dict % has(lat % outer)) then

66
src/lattice.cpp
View file

@ -4,6 +4,8 @@
#include <sstream>
#include <vector>
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "hdf5_interface.h"
#include "xml_interface.h"
@ -131,6 +133,36 @@ RectLattice::RectLattice(pugi::xml_node lat_node)
//==============================================================================
int32_t&
RectLattice::operator[](const int i_xyz[3])
{
int nx = n_cells[0];
int ny = n_cells[1];
int nz = n_cells[2];
int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0];
return universes[indx];
}
//==============================================================================
void
RectLattice::adjust_indices()
{
int nx = n_cells[0];
int ny = n_cells[1];
int nz = n_cells[2];
for (int iz = 0; iz < nz; iz++) {
for (int iy = 0; iy < ny; iy++) {
for (int ix = 0; ix < nx; ix++) {
int indx = nx*ny*iz + nx*iy + ix;
universes[indx] = universe_dict[universes[indx]];
}
}
}
}
//==============================================================================
bool
RectLattice::are_valid_indices(const int i_xyz[3]) const
{
@ -386,6 +418,37 @@ HexLattice::HexLattice(pugi::xml_node lat_node)
//==============================================================================
int32_t&
HexLattice::operator[](const int i_xyz[3])
{
int indx = (2*n_rings-1)*(2*n_rings-1) * i_xyz[2]
+ (2*n_rings-1) * i_xyz[1]
+ i_xyz[0];
return universes[indx];
}
//==============================================================================
void
HexLattice::adjust_indices()
{
for (int iz = 0; iz < n_axial; iz++) {
for (int ia = 0; ia < 2*n_rings-1; ia++) {
for (int ix = 0; ix < 2*n_rings-1; ix++) {
int i_xyz[3] {ix+1, ia+1, iz+1};
if (are_valid_indices(i_xyz)) {
int indx = (2*n_rings-1)*(2*n_rings-1) * iz
+ (2*n_rings-1) * ia
+ ix;
universes[indx] = universe_dict[universes[indx]];
}
}
}
}
}
//==============================================================================
bool
HexLattice::are_valid_indices(const int i_xyz[3]) const
{
@ -655,6 +718,9 @@ extern "C" {
}
void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);}
int32_t lattice_universe(Lattice *lat, const int i_xyz[3])
{return (*lat)[i_xyz];}
}
} // namespace openmc

23
src/lattice.h
View file

@ -15,15 +15,7 @@
namespace openmc {
//==============================================================================
// Constants that should eventually be moved out of this file
//==============================================================================
extern "C" double FP_PRECISION;
constexpr double INFTY{std::numeric_limits<double>::max()};
constexpr int C_NONE {-1};
//==============================================================================
// Constants
// Module constants
//==============================================================================
constexpr int32_t NO_OUTER_UNIVERSE{-1};
@ -58,6 +50,11 @@ public:
virtual ~Lattice() {}
virtual int32_t& operator[](const int i_xyz[3]) = 0;
//! Convert internal universe values from IDs to indices using universe_dict.
virtual void adjust_indices() = 0;
//! 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
@ -104,6 +101,10 @@ public:
virtual ~RectLattice() {}
int32_t& operator[](const int i_xyz[3]);
void adjust_indices();
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>
@ -127,6 +128,10 @@ public:
virtual ~HexLattice() {}
int32_t& operator[](const int i_xyz[3]);
void adjust_indices();
bool are_valid_indices(const int i_xyz[3]) const;
std::pair<double, std::array<int, 3>>

4
src/summary.F90
View file

@ -288,7 +288,7 @@ contains
do k = 0, lat%n_cells(2) - 1
do m = 1, lat%n_cells(3)
lattice_universes(j, k+1, m) = &
universes(lat%universes(j, lat%n_cells(2) - k, m))%id
universes(lat%get([j-1, lat%n_cells(2)-k-1, m-1])+1)%id
end do
end do
end do
@ -319,7 +319,7 @@ contains
lattice_universes(j,k,m) = -1
cycle
end if
lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id
lattice_universes(j,k,m) = universes(lat%get([j-1,k-1,m-1])+1)%id
end do
end do
end do

8
src/tallies/tally_filter_distribcell.F90
View file

@ -289,7 +289,7 @@ contains
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
if (lat % is_3d) then
path = trim(path) // "(" // trim(to_str(k-1)) // &
"," // trim(to_str(l-1)) // "," // &
@ -310,7 +310,7 @@ contains
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
if (lat % is_3d) then
path = trim(path) // "(" // trim(to_str(old_k-1)) // &
"," // trim(to_str(old_l-1)) // "," // &
@ -359,7 +359,7 @@ contains
! This is last lattice cell, so target must be here
lat_offset = lat % offset(map, k, l, m)
offset = offset + lat_offset
next_univ => universes(lat % universes(k, l, m))
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
if (lat % is_3d) then
path = trim(path) // "(" // &
trim(to_str(k - lat % n_rings)) // "," // &
@ -382,7 +382,7 @@ contains
! Target is at this lattice position
lat_offset = lat % offset(map, old_k, old_l, old_m)
offset = offset + lat_offset
next_univ => universes(lat % universes(old_k, old_l, old_m))
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
if (lat % is_3d) then
path = trim(path) // "(" // &
trim(to_str(old_k - lat % n_rings)) // "," // &