crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-25 20:45:35 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Fix rotation matrix in transport()

Browse source
This commit is contained in:
Paul Romano 2019-02-14 15:30:00 -06:00
parent 968e9a95d4
commit 47c19353e4
3 changed files with 5 additions and 43 deletions
Download patch file Download diff file Expand all files Collapse all files

5
src/geometry_header.F90
View file

@ -82,11 +82,6 @@ module geometry_header
type Cell
type(C_PTR) :: ptr
! Rotation matrix and translation vector
real(8), allocatable :: rotation(:)
real(8), allocatable :: rotation_matrix(:,:)
contains
procedure :: id => cell_id

35
src/input_xml.F90
View file

@ -217,41 +217,6 @@ contains
// to_str(c % id()))
end if
! Rotation matrix
if (check_for_node(node_cell, "rotation")) then
! Rotations can only be applied to cells that are being filled with
! another universe
if (c % fill() == C_NONE) then
call fatal_error("Cannot apply a rotation to cell " // trim(to_str(&
&c % id())) // " because it is not filled with another universe")
end if
! Read number of rotation parameters
n = node_word_count(node_cell, "rotation")
if (n /= 3) then
call fatal_error("Incorrect number of rotation parameters on cell " &
// to_str(c % id()))
end if
! Copy rotation angles in x,y,z directions
allocate(c % rotation(3))
call get_node_array(node_cell, "rotation", c % rotation)
phi = -c % rotation(1) * PI/180.0_8
theta = -c % rotation(2) * PI/180.0_8
psi = -c % rotation(3) * PI/180.0_8
! Calculate rotation matrix based on angles given
allocate(c % rotation_matrix(3,3))
c % rotation_matrix = reshape((/ &
cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta), &
-cos(phi)*sin(psi) + sin(phi)*sin(theta)*cos(psi), &
cos(phi)*cos(psi) + sin(phi)*sin(theta)*sin(psi), &
sin(phi)*cos(theta), &
sin(phi)*sin(psi) + cos(phi)*sin(theta)*cos(psi), &
-sin(phi)*cos(psi) + cos(phi)*sin(theta)*sin(psi), &
cos(phi)*cos(theta) /), (/ 3,3 /))
end if
! Add cell to dictionary
call cell_dict % set(c % id(), i)

8
src/particle.cpp
View file

@ -357,9 +357,11 @@ Particle::transport()
for (int j = 0; j < n_coord - 1; ++j) {
if (coord[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
// FIXME
// coord[j + 1].uvw = matmul(cells(coord[j].cell + 1) % &
// rotation_matrix, coord[j].uvw)
const auto& m {model::cells[coord[j].cell]->rotation_};
Direction u {coord[j].uvw};
coord[j + 1].uvw[0] = m[3]*u.x + m[4]*u.y + m[5]*u.z;
coord[j + 1].uvw[1] = m[6]*u.x + m[7]*u.y + m[8]*u.z;
coord[j + 1].uvw[2] = m[9]*u.x + m[10]*u.y + m[11]*u.z;
} else {
// Otherwise, copy this level's direction
std::copy(coord[j].uvw, coord[j].uvw + 3, coord[j + 1].uvw);