Move FILL_MATERIAL chunk of find_cell to C++

This commit is contained in:
Sterling Harper 2018-08-20 12:24:08 -04:00
parent ef8bf2b820
commit cb5db790d7
3 changed files with 76 additions and 62 deletions

View file

@ -69,9 +69,7 @@ contains
type(Particle), intent(inout) :: p
logical, intent(inout) :: found
integer, optional :: search_cells(:)
integer :: j, k ! coordinate level index
integer :: offset ! instance # of a distributed cell
integer :: distribcell_index
integer :: j ! coordinate level index
integer :: i_xyz(3) ! indices in lattice
integer :: i_cell ! index in cells array
@ -85,61 +83,7 @@ contains
if (found) then
associate(c => cells(i_cell))
CELL_TYPE: if (c % type() == FILL_MATERIAL) then
! ======================================================================
! AT LOWEST UNIVERSE, TERMINATE SEARCH
! Save previous material and temperature
p % last_material = p % material
p % last_sqrtkT = p % sqrtkT
! Get distributed offset
if (c % material_size() > 1 .or. c % sqrtkT_size() > 1) then
! Distributed instances of this cell have different
! materials/temperatures. Determine which instance this is for
! assigning the matching material/temperature.
distribcell_index = c % distribcell_index()
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-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 - 1, &
[p % coord(k + 1) % lattice_x, &
p % coord(k + 1) % lattice_y, &
p % coord(k + 1) % lattice_z])
end if
end if
end do
! Keep track of which instance of the cell the particle is in
p % cell_instance = offset + 1
else
p % cell_instance = 1
end if
! Save the material
if (c % material_size() > 1) then
p % material = c % material(offset + 1)
else
p % material = c % material(1)
end if
! Save the temperature
if (c % sqrtkT_size() > 1) then
p % sqrtkT = c % sqrtkT(offset)
else
p % sqrtkT = c % sqrtkT(0)
end if
elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE
CELL_TYPE: if (c % type() == FILL_UNIVERSE) then
! ======================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL

View file

@ -5,7 +5,7 @@
#include "cell.h"
#include "constants.h"
#include "error.h"
#include "particle.h"
#include "lattice.h"
#include "settings.h"
//TODO: remove this include
@ -16,6 +16,8 @@ namespace openmc {
std::vector<int64_t> overlap_check_count;
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p) {
int n_coord = p->n_coord;
@ -72,10 +74,11 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) {
// Find which cell of this universe the particle is in.
bool found = false;
int32_t i_cell;
for (int i = 0; i < n_search_cells; i++) {
int32_t i_cell = search_cells[i];
i_cell = search_cells[i];
// Make sure the search cell is in the same universe
// Make sure the search cell is in the same universe.
//TODO: off-by-one indexing
if (global_cells[i_cell]->universe - 1 != i_universe) continue;
@ -96,6 +99,58 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) {
}
}
if (found) {
Cell& c {*global_cells[i_cell]};
if (c.type == FILL_MATERIAL) {
// Find the distribcell instance number.
if (c.material.size() > 1 || c.sqrtkT.size() > 1) {
//=====================================================================
//! Found a material cell which means this is the lowest coord level.
//TODO: off-by-one indexing
int distribcell_index = c.distribcell_index - 1;
int offset = 0;
for (int i = 0; i < p->n_coord; i++) {
Cell& c_i {*global_cells[p->coord[i].cell-1]};
if (c_i.type == FILL_UNIVERSE) {
offset += c_i.offset[distribcell_index];
} else if (c_i.type == FILL_LATTICE) {
Lattice& lat {*lattices_c[p->coord[i+1].lattice-1]};
int i_xyz[3] {p->coord[i+1].lattice_x,
p->coord[i+1].lattice_y,
p->coord[i+1].lattice_z};
if (lat.are_valid_indices(i_xyz)) {
offset += lat.offset(distribcell_index, i_xyz);
}
}
}
p->cell_instance = offset + 1;
} else {
p->cell_instance = 1;
}
// Set the material and temperature.
p->last_material = p->material;
int32_t mat;
if (c.material.size() > 1) {
mat = c.material[p->cell_instance-1];
} else {
mat = c.material[0];
}
if (mat == MATERIAL_VOID) {
p->material = MATERIAL_VOID;
} else {
p->material = mat + 1;
}
p->last_sqrtkT = p->sqrtkT;
if (c.sqrtkT.size() > 1) {
p->sqrtkT = c.sqrtkT[p->cell_instance-1];
} else {
p->sqrtkT = c.sqrtkT[0];
}
}
}
return found;
}

View file

@ -4,13 +4,28 @@
#include <cstdint>
#include <vector>
#include "particle.h"
namespace openmc {
extern "C" int openmc_root_universe;
//TODO: free this memory
extern std::vector<int64_t> overlap_check_count;
//==============================================================================
//! Check for overlapping cells at the particle's position.
//==============================================================================
extern "C" bool
check_cell_overlap(Particle* p);
//==============================================================================
//! Locate the particle in the geometry tree and set its geometry data fields.
//==============================================================================
extern "C" bool
find_cell(Particle* p, int n_search_cells, int* search_cells);
} // namespace openmc
#endif // OPENMC_GEOMETRY_H