mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 05:35:49 -04:00
Don't include 0-index mesh cells for mesh surface filter
This commit is contained in:
Paul Romano
parent
794ac9255a
commit
616d492d01
4 changed files with 130 additions and 349 deletions
|
|
@ -2927,6 +2927,7 @@ contains
|
|||
err = openmc_tally_set_filters(i_start + i - 1, n_filter, temp_filter)
|
||||
deallocate(temp_filter)
|
||||
else
|
||||
t % type = TALLY_MESH_CURRENT
|
||||
t % score_bins(j) = SCORE_CURRENT
|
||||
end if
|
||||
|
||||
|
|
|
|||
|
|
@ -771,12 +771,6 @@ contains
|
|||
end associate
|
||||
end if
|
||||
|
||||
! Handle surface current tallies separately
|
||||
if (t % type == TALLY_MESH_CURRENT) then
|
||||
call write_surface_current(t, unit_tally)
|
||||
cycle
|
||||
end if
|
||||
|
||||
! WARNING: Admittedly, the logic for moving for printing results is
|
||||
! extremely confusing and took quite a bit of time to get correct. The
|
||||
! logic is structured this way since it is not practical to have a do
|
||||
|
|
|
|||
|
|
@ -3064,7 +3064,7 @@ contains
|
|||
! tally total or partial currents between two cells
|
||||
!===============================================================================
|
||||
|
||||
subroutine score_surface_tally(p)
|
||||
subroutine score_surface_tally(p)
|
||||
|
||||
type(Particle), intent(in) :: p
|
||||
|
||||
|
|
@ -3199,277 +3199,121 @@ contains
|
|||
|
||||
integer :: i
|
||||
integer :: i_tally
|
||||
integer :: j, k ! loop indices
|
||||
integer :: n_dim ! num dimensions of the mesh
|
||||
integer :: d1 ! dimension index
|
||||
integer :: d2 ! dimension index
|
||||
integer :: d3 ! dimension index
|
||||
integer :: ijk0(3) ! indices of starting coordinates
|
||||
integer :: ijk1(3) ! indices of ending coordinates
|
||||
integer :: n_cross ! number of surface crossings
|
||||
integer :: filter_index ! index of scoring bin
|
||||
integer :: i_filter_mesh ! index of mesh filter in filters array
|
||||
integer :: i_filter_surf ! index of surface filter in filters
|
||||
integer :: i_filter_energy ! index of energy filter in filters
|
||||
real(8) :: uvw(3) ! cosine of angle of particle
|
||||
real(8) :: xyz0(3) ! starting/intermediate coordinates
|
||||
real(8) :: xyz1(3) ! ending coordinates of particle
|
||||
real(8) :: xyz_cross(3) ! coordinates of bounding surfaces
|
||||
real(8) :: d(3) ! distance to each bounding surface
|
||||
real(8) :: distance ! actual distance traveled
|
||||
integer :: matching_bin ! next valid filter bin
|
||||
logical :: start_in_mesh ! particle's starting xyz in mesh?
|
||||
logical :: end_in_mesh ! particle's ending xyz in mesh?
|
||||
logical :: cross_surface ! whether the particle crosses a surface
|
||||
logical :: energy_filter ! energy filter present
|
||||
type(RegularMesh), pointer :: m
|
||||
integer :: i_filt
|
||||
integer :: i_bin
|
||||
integer :: q ! loop index for scoring bins
|
||||
integer :: k ! working index for expand and score
|
||||
integer :: score_bin ! scoring bin, e.g. SCORE_FLUX
|
||||
integer :: score_index ! scoring bin index
|
||||
integer :: j ! loop index for scoring bins
|
||||
integer :: filter_index ! single index for single bin
|
||||
real(8) :: flux ! collision estimate of flux
|
||||
real(8) :: filter_weight ! combined weight of all filters
|
||||
real(8) :: score ! analog tally score
|
||||
logical :: finished ! found all valid bin combinations
|
||||
|
||||
! No collision, so no weight change when survival biasing
|
||||
flux = p % wgt
|
||||
|
||||
TALLY_LOOP: do i = 1, active_current_tallies % size()
|
||||
! Copy starting and ending location of particle
|
||||
xyz0 = p % last_xyz_current
|
||||
xyz1 = p % coord(1) % xyz
|
||||
|
||||
! Get pointer to tally
|
||||
! Get index of tally and pointer to tally
|
||||
i_tally = active_current_tallies % data(i)
|
||||
associate (t => tallies(i_tally) % obj)
|
||||
|
||||
! Check for energy filter
|
||||
energy_filter = (t % find_filter(FILTER_ENERGYIN) > 0)
|
||||
|
||||
! Get index for mesh, surface, and energy filters
|
||||
i_filter_mesh = t % filter(t % find_filter(FILTER_MESH))
|
||||
i_filter_surf = t % filter(t % find_filter(FILTER_SURFACE))
|
||||
if (energy_filter) then
|
||||
i_filter_energy = t % filter(t % find_filter(FILTER_ENERGYIN))
|
||||
end if
|
||||
|
||||
! Reset the matching bins arrays
|
||||
call filter_matches(i_filter_mesh) % bins % resize(1)
|
||||
call filter_matches(i_filter_surf) % bins % resize(1)
|
||||
if (energy_filter) then
|
||||
call filter_matches(i_filter_energy) % bins % resize(1)
|
||||
end if
|
||||
|
||||
! Get pointer to mesh
|
||||
select type(filt => filters(i_filter_mesh) % obj)
|
||||
type is (MeshFilter)
|
||||
m => meshes(filt % mesh)
|
||||
end select
|
||||
|
||||
n_dim = m % n_dimension
|
||||
|
||||
! Determine indices for starting and ending location
|
||||
call m % get_indices(xyz0, ijk0, start_in_mesh)
|
||||
call m % get_indices(xyz1, ijk1, end_in_mesh)
|
||||
|
||||
! Check to see if start or end is in mesh -- if not, check if track still
|
||||
! intersects with mesh
|
||||
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
|
||||
if (.not. m % intersects(xyz0, xyz1)) cycle
|
||||
end if
|
||||
|
||||
! Calculate number of surface crossings
|
||||
n_cross = sum(abs(ijk1(:n_dim) - ijk0(:n_dim)))
|
||||
if (n_cross == 0) then
|
||||
cycle
|
||||
end if
|
||||
|
||||
! Copy particle's direction
|
||||
uvw = p % coord(1) % uvw
|
||||
|
||||
! Determine incoming energy bin. We need to tell the energy filter this
|
||||
! is a tracklength tally so it uses the pre-collision energy.
|
||||
if (energy_filter) then
|
||||
call filter_matches(i_filter_energy) % bins % clear()
|
||||
call filter_matches(i_filter_energy) % weights % clear()
|
||||
call filters(i_filter_energy) % obj % get_all_bins(p, &
|
||||
ESTIMATOR_TRACKLENGTH, filter_matches(i_filter_energy))
|
||||
if (filter_matches(i_filter_energy) % bins % size() == 0) cycle
|
||||
matching_bin = filter_matches(i_filter_energy) % bins % data(1)
|
||||
filter_matches(i_filter_energy) % bins % data(1) = matching_bin
|
||||
end if
|
||||
|
||||
! Bounding coordinates
|
||||
do d1 = 1, n_dim
|
||||
if (uvw(d1) > 0) then
|
||||
xyz_cross(d1) = m % lower_left(d1) + ijk0(d1) * m % width(d1)
|
||||
else
|
||||
xyz_cross(d1) = m % lower_left(d1) + (ijk0(d1) - 1) * m % width(d1)
|
||||
! Find all valid bins in each filter if they have not already been found
|
||||
! for a previous tally.
|
||||
do j = 1, size(t % filter)
|
||||
i_filt = t % filter(j)
|
||||
if (.not. filter_matches(i_filt) % bins_present) then
|
||||
call filter_matches(i_filt) % bins % clear()
|
||||
call filter_matches(i_filt) % weights % clear()
|
||||
call filters(i_filt) % obj % get_all_bins(p, t % estimator, &
|
||||
filter_matches(i_filt))
|
||||
filter_matches(i_filt) % bins_present = .true.
|
||||
end if
|
||||
! If there are no valid bins for this filter, then there is nothing to
|
||||
! score and we can move on to the next tally.
|
||||
if (filter_matches(i_filt) % bins % size() == 0) cycle TALLY_LOOP
|
||||
|
||||
! Set the index of the bin used in the first filter combination
|
||||
filter_matches(i_filt) % i_bin = 1
|
||||
end do
|
||||
|
||||
do j = 1, n_cross
|
||||
! Reset scoring bin index
|
||||
filter_matches(i_filter_surf) % bins % data(1) = 0
|
||||
! ========================================================================
|
||||
! Loop until we've covered all valid bins on each of the filters.
|
||||
|
||||
! Set the distances to infinity
|
||||
d = INFINITY
|
||||
FILTER_LOOP: do
|
||||
|
||||
! Calculate distance to each bounding surface. We need to treat
|
||||
! special case where the cosine of the angle is zero since this would
|
||||
! result in a divide-by-zero.
|
||||
do d1 = 1, n_dim
|
||||
if (uvw(d1) == 0) then
|
||||
d(d1) = INFINITY
|
||||
! Reset scoring index and weight
|
||||
filter_index = 1
|
||||
filter_weight = ONE
|
||||
|
||||
! Determine scoring index and weight for this filter combination
|
||||
do j = 1, size(t % filter)
|
||||
i_filt = t % filter(j)
|
||||
i_bin = filter_matches(i_filt) % i_bin
|
||||
filter_index = filter_index + (filter_matches(i_filt) % bins % &
|
||||
data(i_bin) - 1) * t % stride(j)
|
||||
filter_weight = filter_weight * filter_matches(i_filt) % weights % &
|
||||
data(i_bin)
|
||||
end do
|
||||
|
||||
! Determine score
|
||||
score = flux * filter_weight
|
||||
|
||||
! Currently only one score type
|
||||
k = 0
|
||||
SCORE_LOOP: do q = 1, t % n_user_score_bins
|
||||
k = k + 1
|
||||
|
||||
! determine what type of score bin
|
||||
score_bin = t % score_bins(q)
|
||||
|
||||
! determine scoring bin index, no offset from nuclide bins
|
||||
score_index = q
|
||||
|
||||
call expand_and_score(p, t, score_index, filter_index, score_bin, &
|
||||
score, k)
|
||||
end do SCORE_LOOP
|
||||
|
||||
! ======================================================================
|
||||
! Filter logic
|
||||
|
||||
! Increment the filter bins, starting with the last filter to find the
|
||||
! next valid bin combination
|
||||
finished = .true.
|
||||
do j = size(t % filter), 1, -1
|
||||
i_filt = t % filter(j)
|
||||
if (filter_matches(i_filt) % i_bin < filter_matches(i_filt) % &
|
||||
bins % size()) then
|
||||
filter_matches(i_filt) % i_bin = filter_matches(i_filt) % i_bin + 1
|
||||
finished = .false.
|
||||
exit
|
||||
else
|
||||
d(d1) = (xyz_cross(d1) - xyz0(d1))/uvw(d1)
|
||||
filter_matches(i_filt) % i_bin = 1
|
||||
end if
|
||||
end do
|
||||
|
||||
! Determine the closest bounding surface of the mesh cell by
|
||||
! calculating the minimum distance. Then use the minimum distance and
|
||||
! direction of the particle to determine which surface was crossed.
|
||||
distance = minval(d)
|
||||
! Once we have finished all valid bins for each of the filters, exit
|
||||
! the loop.
|
||||
if (finished) exit FILTER_LOOP
|
||||
|
||||
! Loop over the dimensions
|
||||
do d1 = 1, n_dim
|
||||
|
||||
! Get the other dimensions.
|
||||
if (d1 == 1) then
|
||||
d2 = mod(d1, 3) + 1
|
||||
d3 = mod(d1 + 1, 3) + 1
|
||||
else
|
||||
d2 = mod(d1 + 1, 3) + 1
|
||||
d3 = mod(d1, 3) + 1
|
||||
end if
|
||||
|
||||
! Check whether distance is the shortest distance
|
||||
if (distance == d(d1)) then
|
||||
|
||||
! Check whether particle is moving in positive d1 direction
|
||||
if (uvw(d1) > 0) then
|
||||
|
||||
! Outward current on d1 max surface
|
||||
if (all(ijk0(:n_dim) >= 1) .and. &
|
||||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 1
|
||||
filter_matches(i_filter_mesh) % bins % data(1) = &
|
||||
m % get_bin_from_indices(ijk0)
|
||||
filter_index = 1
|
||||
do k = 1, size(t % filter)
|
||||
filter_index = filter_index + (filter_matches(t % &
|
||||
filter(k)) % bins % data(1) - 1) * t % stride(k)
|
||||
end do
|
||||
!$omp atomic
|
||||
t % results(RESULT_VALUE, 1, filter_index) = &
|
||||
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
|
||||
end if
|
||||
|
||||
! Inward current on d1 min surface
|
||||
cross_surface = .false.
|
||||
select case(n_dim)
|
||||
|
||||
case (1)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (2)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. &
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (3)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. &
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
|
||||
ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
end select
|
||||
|
||||
! If the particle crossed the surface, tally the current
|
||||
if (cross_surface) then
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 2
|
||||
filter_matches(i_filter_mesh) % bins % data(1) = &
|
||||
m % get_bin_from_indices(ijk0)
|
||||
filter_index = 1
|
||||
do k = 1, size(t % filter)
|
||||
filter_index = filter_index + (filter_matches(t % &
|
||||
filter(k)) % bins % data(1) - 1) * t % stride(k)
|
||||
end do
|
||||
!$omp atomic
|
||||
t % results(RESULT_VALUE, 1, filter_index) = &
|
||||
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
end if
|
||||
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
|
||||
|
||||
! The particle is moving in the negative d1 direction
|
||||
else
|
||||
|
||||
! Outward current on d1 min surface
|
||||
if (all(ijk0(:n_dim) >= 1) .and. &
|
||||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 3
|
||||
filter_matches(i_filter_mesh) % bins % data(1) = &
|
||||
m % get_bin_from_indices(ijk0)
|
||||
filter_index = 1
|
||||
do k = 1, size(t % filter)
|
||||
filter_index = filter_index + (filter_matches(t % &
|
||||
filter(k)) % bins % data(1) - 1) * t % stride(k)
|
||||
end do
|
||||
!$omp atomic
|
||||
t % results(RESULT_VALUE, 1, filter_index) = &
|
||||
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
|
||||
end if
|
||||
|
||||
! Inward current on d1 max surface
|
||||
cross_surface = .false.
|
||||
select case(n_dim)
|
||||
|
||||
case (1)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (2)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (3)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
|
||||
ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
end select
|
||||
|
||||
! If the particle crossed the surface, tally the current
|
||||
if (cross_surface) then
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
filter_matches(i_filter_surf) % bins % data(1) = d1 * 4
|
||||
filter_matches(i_filter_mesh) % bins % data(1) = &
|
||||
m % get_bin_from_indices(ijk0)
|
||||
filter_index = 1
|
||||
do k = 1, size(t % filter)
|
||||
filter_index = filter_index + (filter_matches(t % &
|
||||
filter(k)) % bins % data(1) - 1) * t % stride(k)
|
||||
end do
|
||||
!$omp atomic
|
||||
t % results(RESULT_VALUE, 1, filter_index) = &
|
||||
t % results(RESULT_VALUE, 1, filter_index) + p % wgt
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
end if
|
||||
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
||||
! Calculate new coordinates
|
||||
xyz0 = xyz0 + distance * uvw
|
||||
end do
|
||||
end do FILTER_LOOP
|
||||
|
||||
end associate
|
||||
|
||||
! If the user has specified that we can assume all tallies are spatially
|
||||
! separate, this implies that once a tally has been scored to, we needn't
|
||||
! check the others. This cuts down on overhead when there are many
|
||||
! tallies specified
|
||||
|
||||
if (assume_separate) exit TALLY_LOOP
|
||||
|
||||
end do TALLY_LOOP
|
||||
|
||||
! Reset filter matches flag
|
||||
filter_matches(:) % bins_present = .false.
|
||||
|
||||
end subroutine score_surface_current
|
||||
|
||||
!===============================================================================
|
||||
|
|
|
|||
|
|
@ -64,7 +64,7 @@ contains
|
|||
|
||||
! Determine number of bins
|
||||
n_dim = meshes(i_mesh) % n_dimension
|
||||
this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension + 1)
|
||||
this % n_bins = 4*n_dim*product(meshes(i_mesh) % dimension)
|
||||
|
||||
! Store the index of the mesh
|
||||
this % mesh = i_mesh
|
||||
|
|
@ -79,8 +79,6 @@ contains
|
|||
integer :: j ! loop indices
|
||||
integer :: n_dim ! num dimensions of the mesh
|
||||
integer :: d1 ! dimension index
|
||||
integer :: d2 ! dimension index
|
||||
integer :: d3 ! dimension index
|
||||
integer :: ijk0(3) ! indices of starting coordinates
|
||||
integer :: ijk1(3) ! indices of ending coordinates
|
||||
integer :: n_cross ! number of surface crossings
|
||||
|
|
@ -95,7 +93,6 @@ contains
|
|||
real(8) :: distance ! actual distance traveled
|
||||
logical :: start_in_mesh ! particle's starting xyz in mesh?
|
||||
logical :: end_in_mesh ! particle's ending xyz in mesh?
|
||||
logical :: cross_surface ! whether the particle crosses a surface
|
||||
|
||||
! Copy starting and ending location of particle
|
||||
xyz0 = p % last_xyz_current
|
||||
|
|
@ -153,15 +150,6 @@ contains
|
|||
! Loop over the dimensions
|
||||
do d1 = 1, n_dim
|
||||
|
||||
! Get the other dimensions.
|
||||
if (d1 == 1) then
|
||||
d2 = mod(d1, 3) + 1
|
||||
d3 = mod(d1 + 1, 3) + 1
|
||||
else
|
||||
d2 = mod(d1 + 1, 3) + 1
|
||||
d3 = mod(d1, 3) + 1
|
||||
end if
|
||||
|
||||
! Check whether distance is the shortest distance
|
||||
if (distance == d(d1)) then
|
||||
|
||||
|
|
@ -173,103 +161,57 @@ contains
|
|||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
i_surf = d1 * 4 - 1
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*i_mesh + i_surf + 1
|
||||
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
end if
|
||||
|
||||
! Inward current on d1 min surface
|
||||
cross_surface = .false.
|
||||
select case(n_dim)
|
||||
|
||||
case (1)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (2)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. &
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (3)
|
||||
if (ijk0(d1) >= 0 .and. ijk0(d1) < m % dimension(d1) .and. &
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
|
||||
ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
end select
|
||||
|
||||
! If the particle crossed the surface, tally the current
|
||||
if (cross_surface) then
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
i_surf = d1 * 4 - 2
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*i_mesh + i_surf + 1
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
end if
|
||||
|
||||
! Advance position
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
xyz_cross(d1) = xyz_cross(d1) + m % width(d1)
|
||||
|
||||
! The particle is moving in the negative d1 direction
|
||||
! If the particle crossed the surface, tally the inward current on
|
||||
! d1 min surface
|
||||
if (all(ijk0(:n_dim) >= 1) .and. &
|
||||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
i_surf = d1 * 4 - 2
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
end if
|
||||
|
||||
else
|
||||
! The particle is moving in the negative d1 direction
|
||||
|
||||
! Outward current on d1 min surface
|
||||
if (all(ijk0(:n_dim) >= 1) .and. &
|
||||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
i_surf = d1 * 4 - 3
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*i_mesh + i_surf + 1
|
||||
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
end if
|
||||
|
||||
! Inward current on d1 max surface
|
||||
cross_surface = .false.
|
||||
select case(n_dim)
|
||||
|
||||
case (1)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (2)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
|
||||
case (3)
|
||||
if (ijk0(d1) > 1 .and. ijk0(d1) <= m % dimension(d1) + 1 .and.&
|
||||
ijk0(d2) >= 1 .and. ijk0(d2) <= m % dimension(d2) .and. &
|
||||
ijk0(d3) >= 1 .and. ijk0(d3) <= m % dimension(d3)) then
|
||||
cross_surface = .true.
|
||||
end if
|
||||
end select
|
||||
|
||||
! If the particle crossed the surface, tally the current
|
||||
if (cross_surface) then
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
i_surf = d1 * 4
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*i_mesh + i_surf + 1
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
|
||||
ijk0(d1) = ijk0(d1) + 1
|
||||
end if
|
||||
|
||||
! Advance position
|
||||
ijk0(d1) = ijk0(d1) - 1
|
||||
xyz_cross(d1) = xyz_cross(d1) - m % width(d1)
|
||||
|
||||
! If the particle crossed the surface, tally the inward current on
|
||||
! d1 max surface
|
||||
if (all(ijk0(:n_dim) >= 1) .and. &
|
||||
all(ijk0(:n_dim) <= m % dimension)) then
|
||||
i_surf = d1 * 4
|
||||
i_mesh = m % get_bin_from_indices(ijk0)
|
||||
i_bin = 4*n_dim*(i_mesh - 1) + i_surf
|
||||
|
||||
call match % bins % push_back(i_bin)
|
||||
call match % weights % push_back(ONE)
|
||||
end if
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
|
|
@ -305,7 +247,7 @@ contains
|
|||
allocate(ijk(n_dim))
|
||||
|
||||
! Get flattend mesh index and surface index
|
||||
i_mesh = (bin - 1) / (4*n_dim)
|
||||
i_mesh = (bin - 1) / (4*n_dim) + 1
|
||||
i_surf = mod(bin - 1, 4*n_dim) + 1
|
||||
|
||||
! Get mesh index part of label
|
||||
|
|
@ -370,7 +312,7 @@ contains
|
|||
if (index_mesh >= 1 .and. index_mesh <= n_meshes) then
|
||||
f % mesh = index_mesh
|
||||
n_dim = meshes(index_mesh) % n_dimension
|
||||
f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension + 1)
|
||||
f % n_bins = 4*n_dim*product(meshes(index_mesh) % dimension)
|
||||
else
|
||||
err = E_OUT_OF_BOUNDS
|
||||
call set_errmsg("Index in 'meshes' array is out of bounds.")
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue