diff --git a/src/constants.F90 b/src/constants.F90 index 9c52572db6..77d2dcb228 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -305,7 +305,7 @@ module constants EVENT_ABSORB = 2 ! Tally score type - integer, parameter :: N_SCORE_TYPES = 25 + integer, parameter :: N_SCORE_TYPES = 26 integer, parameter :: & SCORE_FLUX = -1, & ! flux SCORE_TOTAL = -2, & ! total reaction rate @@ -331,7 +331,8 @@ module constants SCORE_FISS_Q_PROMPT = -22, & ! prompt fission Q-value SCORE_FISS_Q_RECOV = -23, & ! recoverable fission Q-value SCORE_DECAY_RATE = -24, & ! delayed neutron precursor decay rate - SCORE_CELL_TO_CELL_TYPE = -25 ! cell to cell partial current + SCORE_CELL_TO_CELL = -25, & ! cell to cell partial current + SCORE_CELL_TO_CELL_NORMAL_PROJECTION = -26 ! cell to cell partial current projected on surface normal ! Maximum scattering order supported integer, parameter :: MAX_ANG_ORDER = 10 diff --git a/src/endf.F90 b/src/endf.F90 index 0ba2db388b..f672776a94 100644 --- a/src/endf.F90 +++ b/src/endf.F90 @@ -66,6 +66,10 @@ contains string = "fission-q-prompt" case (SCORE_FISS_Q_RECOV) string = "fission-q-recoverable" + case (SCORE_CELL_TO_CELL) + string = "cell-to-cell-partial-current" + case (SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + string = "cell-to-cell-partial-current-projection-on-normal" ! Normal ENDF-based reactions case (TOTAL_XS) diff --git a/src/geometry.F90 b/src/geometry.F90 index 7283748a0c..03a7de7888 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -12,7 +12,7 @@ module geometry use stl_vector, only: VectorInt use string, only: to_str use tally, only: score_surface_current, & - &score_cell_to_cell + &score_partial_current implicit none @@ -399,8 +399,6 @@ contains logical :: rotational ! if rotational periodic BC applied logical :: found ! particle found in universe? class(Surface), pointer :: surf - - !print *, "cross surface called" i_surface = abs(p % surface) surf => surfaces(i_surface)%obj @@ -613,23 +611,33 @@ contains ! cells on the positive side call find_cell(p, found, surf%neighbor_pos) - - ! /CHANGE/ Score cell_to_cell -! print *, "Positive side, calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) - + + ! Save cosine of angle between surface normal and particle direction + p % normal_proj = dot_product(p % coord(1) % uvw, & + & surf % normal(p % coord(1) % xyz)) / & + & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& + &, surf % normal(p % coord(1) % xyz))) + + ! Score cell to cell partial currents + call score_partial_current(p) + if (found) return elseif (p % surface < 0 .and. allocated(surf%neighbor_neg)) then ! If coming from positive side of surface, search all the neighboring ! cells on the negative side + ! Save cosine of angle between surface normal and particle direction + p % normal_proj = dot_product(p % coord(1) % uvw, & + & surf % normal(p % coord(1) % xyz)) / & + & sqrt(dot_product(surf % normal(p % coord(1) % xyz)& + &, surf % normal(p % coord(1) % xyz))) + call find_cell(p, found, surf%neighbor_neg) - - ! /CHANGE/ Score cell_to_cell -! print *, "Negative side, calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) - + + ! Score cell to cell partial currents + call score_partial_current(p) + if (found) return end if @@ -662,10 +670,9 @@ contains return end if end if - - ! /CHANGE/ Score cell_to_cell -! print *, "Searched all cells then calling score_cell_to_cell" - call score_cell_to_cell(p, last_cell) + + ! Score cell to cell partial currents + call score_partial_current(p) end subroutine cross_surface diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 7ded1060cb..e74824e1d4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3004,7 +3004,7 @@ contains end if n_words = node_word_count(node_filt, "bins") case ("mesh", "universe", "material", "cell", "distribcell", & - "cellborn", "surface", "delayedgroup") + "cellborn", "cellto", "cellfrom", "surface", "delayedgroup") if (.not. check_for_node(node_filt, "bins")) then call fatal_error("Bins not set in filter " // trim(to_str(filter_id))) end if @@ -3035,7 +3035,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellfrom') ! Allocate and declare the filter type allocate(CellFromFilter :: f % obj) @@ -3046,7 +3046,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellto') ! Allocate and declare the filter type allocate(CellToFilter :: f % obj) @@ -3057,7 +3057,7 @@ contains allocate(filt % cells(n_words)) call get_node_array(node_filt, "bins", filt % cells) end select - + case ('cellborn') ! Allocate and declare the filter type allocate(CellbornFilter :: f % obj) @@ -3414,17 +3414,10 @@ contains t % find_filter(FILTER_DISTRIBCELL) = j type is (CellFilter) t % find_filter(FILTER_CELL) = j - type is (CellFromFilter) t % find_filter(FILTER_CELL_TO_CELL) = j - - t % estimator = ESTIMATOR_ANALOG - type is (CellToFilter) t % find_filter(FILTER_CELL_TO_CELL) = j - - t % estimator = ESTIMATOR_ANALOG - type is (CellbornFilter) t % find_filter(FILTER_CELLBORN) = j type is (MaterialFilter) @@ -3927,13 +3920,14 @@ contains t % find_filter(FILTER_SURFACE) = n_filter t % filter(n_filter) = i_filt - case ('events') - t % score_bins(j) = SCORE_EVENTS - case ('partial_current') t % type = TALLY_CELL_TO_CELL - t % score_bins(j) = SCORE_CELL_TO_CELL_TYPE - + t % score_bins(j) = SCORE_CELL_TO_CELL + case ('projected_partial_current') + t % type = TALLY_CELL_TO_CELL + t % score_bins(j) = SCORE_CELL_TO_CELL_NORMAL_PROJECTION + case ('events') + t % score_bins(j) = SCORE_EVENTS case ('elastic', '(n,elastic)') t % score_bins(j) = ELASTIC case ('(n,2nd)') diff --git a/src/output.F90 b/src/output.F90 index 10d4f0511c..46e8827564 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -733,6 +733,9 @@ contains score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" + score_names(abs(SCORE_CELL_TO_CELL)) = "Partial current" + score_names(abs(SCORE_CELL_TO_CELL_NORMAL_PROJECTION)) = "Partial & + ¤t projected on normal of surface" ! Create filename for tally output filename = trim(path_output) // "tallies.out" diff --git a/src/particle_header.F90 b/src/particle_header.F90 index c6d478540d..bf22a2c6af 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -15,9 +15,7 @@ module particle_header !=============================================================================== type LocalCoord - - integer :: last_cell = NONE - + ! Indices in various arrays for this level integer :: cell = NONE integer :: universe = NONE @@ -60,13 +58,11 @@ module particle_header real(8) :: wgt ! particle weight real(8) :: mu ! angle of scatter logical :: alive ! is particle alive? - - ! Previous cell before crossing surface - integer :: last_cell - ! Particle coordinates - integer :: last_n_coord ! number of current coordinates -! type(LocalCoord) :: last_coord(MAX_COORD) ! coordinates for all levels - + + ! Crossing surface data + integer :: last_cell ! last cell the particle was in before crossing + real(8) :: normal_proj ! cos of angle to surface normal when crossing + ! Pre-collision physical data real(8) :: last_xyz_current(3) ! coordinates of the last collision or ! reflective/periodic surface crossing @@ -136,6 +132,7 @@ contains ! clear attributes this % surface = NONE + this % last_cell = NONE this % cell_born = NONE this % material = NONE this % last_material = NONE @@ -143,6 +140,7 @@ contains this % wgt = ONE this % last_wgt = ONE this % absorb_wgt = ZERO + this % normal_proj = 1 this % n_bank = 0 this % wgt_bank = ZERO this % sqrtkT = ERROR_REAL @@ -181,8 +179,6 @@ contains elemental subroutine reset_coord(this) class(LocalCoord), intent(inout) :: this - this % last_cell = NONE - this % cell = NONE this % universe = NONE this % lattice = NONE diff --git a/src/tally.F90 b/src/tally.F90 index 3d18d162f0..ccfd828f15 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -116,12 +116,8 @@ contains !######################################################################### ! Determine appropirate scoring value. - -! print *, "Score bin", score_bin select case(score_bin) - - ! use default case for cell_to_from tallies, with t % estimator == ESTIMATOR_ANALOG case (SCORE_FLUX, SCORE_FLUX_YN) if (t % estimator == ESTIMATOR_ANALOG) then @@ -1180,16 +1176,13 @@ contains end if end if - case(SCORE_CELL_TO_CELL_TYPE) - ! Not using default because score_bin and event_MT thing - ! score_bin is -25, and event_MT seems to match the cell id - ! print*, p % event_MT, score_bin -! print *, "Scoring", p % last_wgt, flux + case(SCORE_CELL_TO_CELL) score = p % last_wgt * flux - + + case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + score = p % last_wgt * flux * p % normal_proj + case default - !print *, "Found right score type" - !print * , t % estimator if (t % estimator == ESTIMATOR_ANALOG) then ! Any other score is assumed to be a MT number. Thus, we just need ! to check if it matches the MT number of the event @@ -2091,14 +2084,13 @@ contains case (SCORE_EVENTS) ! Simply count number of scoring events score = ONE - - case(SCORE_CELL_TO_CELL_TYPE) - ! Not using default because score_bin and event_MT thing - ! score_bin is -25, and event_MT seems to match the cell id - ! print*, p % event_MT, score_bin -! print *, "Scoring", p % last_wgt, flux + + case(SCORE_CELL_TO_CELL) score = p % last_wgt * flux + case(SCORE_CELL_TO_CELL_NORMAL_PROJECTION) + score = p % last_wgt * flux * p % normal_proj + end select !######################################################################### @@ -2130,8 +2122,6 @@ contains integer :: num_nm ! Number of N,M orders in harmonic integer :: n ! Moment loop index real(8) :: uvw(3) - -! print *, "Expand and score, score bin", score_bin, "score", score, "filter_index", filter_index select case(score_bin) case (SCORE_SCATTER_N, SCORE_NU_SCATTER_N) @@ -2215,7 +2205,6 @@ contains case default -! print *, "Adding to total in expand-and-score", score !$omp atomic t % results(RESULT_VALUE, score_index, filter_index) = & t % results(RESULT_VALUE, score_index, filter_index) + score @@ -3160,18 +3149,17 @@ contains end subroutine score_collision_tally !=============================================================================== -! SCORE_CELL_TO_CELL tallies partial currents from cell_to to cell_from +! score_partial_current tallies partial currents from cell_to to cell_from !=============================================================================== - subroutine score_cell_to_cell(p, last_cell) + subroutine score_partial_current(p) type(Particle), intent(in) :: p - integer, intent(in) :: last_cell - + ! if particle has leaked - ! DF could preserve leakage aiight - !if p % alive = .false. then - + ! DF could preserve leakage + ! if p % alive = .false. then + !end if integer :: i integer :: i_tally @@ -3183,9 +3171,8 @@ contains real(8) :: flux ! collision estimate of flux real(8) :: filter_weight ! combined weight of all filters logical :: finished ! found all valid bin combinations - real(8) :: dummy_real8 type(TallyObject), pointer :: t - + ! Determine collision estimate of flux if (survival_biasing) then ! We need to account for the fact that some weight was already absorbed @@ -3193,13 +3180,15 @@ contains else flux = p % last_wgt end if - + ! A loop over all tallies is necessary because we need to simultaneously ! determine different filter bins for the same tally in order to score to it + ! However, we can restrict the list to cell_to_cell tallies since we know we + ! are crossing a surface - TALLY_LOOP: do i = 1, active_collision_tallies % size() + TALLY_LOOP: do i = 1, active_cell_to_cell_tallies % size() ! Get index of tally and pointer to tally - i_tally = active_collision_tallies % data(i) + i_tally = active_cell_to_cell_tallies % data(i) t => tallies(i_tally) ! Find all valid bins in each filter if they have not already been found @@ -3245,11 +3234,10 @@ contains filter_weight = filter_weight * filter_matches(i_filt) % weights % & data(i_bin) end do - - + call score_general(p, t, 0, filter_index, & - 0, dummy_real8, filter_weight) - + 0, 0d0, flux * filter_weight) + ! ====================================================================== ! Filter logic @@ -3282,15 +3270,15 @@ contains if (assume_separate) exit TALLY_LOOP end do TALLY_LOOP - + ! Reset filter matches flag filter_matches(:) % bins_present = .false. ! Reset tally map positioning position = 0 - - end subroutine score_cell_to_cell - + + end subroutine score_partial_current + !=============================================================================== ! SCORE_SURFACE_CURRENT tallies surface crossings in a mesh tally by manually ! determining which mesh surfaces were crossed @@ -4584,9 +4572,8 @@ contains call active_current_tallies % push_back(i_user_tallies + i) elseif (user_tallies(i) % type == TALLY_CELL_TO_CELL) then call active_cell_to_cell_tallies % push_back(i_user_tallies + i) - call active_current_tallies % push_back(i_user_tallies + i) end if - + end do call active_tallies % shrink_to_fit() @@ -4594,6 +4581,7 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() + call active_cell_to_cell_tallies % shrink_to_fit() end subroutine setup_active_usertallies @@ -4617,6 +4605,9 @@ contains else if (active_current_tallies % size() > 0) then call fatal_error("Active current tallies should not exist before CMFD & &tallies!") + else if (active_cell_to_cell_tallies % size() > 0) then + call fatal_error("Active cell to cell tallies should not exist before & + &CMFD tallies!") end if do i = 1, n_cmfd_tallies @@ -4627,6 +4618,7 @@ contains if (cmfd_tallies(i) % type == TALLY_VOLUME) then if (cmfd_tallies(i) % estimator == ESTIMATOR_ANALOG) then call active_analog_tallies % push_back(i_cmfd_tallies + i) + call active_cell_to_cell_tallies % push_back(i_cmfd_tallies + 1) elseif (cmfd_tallies(i) % estimator == ESTIMATOR_TRACKLENGTH) then call active_tracklength_tallies % push_back(i_cmfd_tallies + i) end if @@ -4640,6 +4632,7 @@ contains call active_tracklength_tallies % shrink_to_fit() call active_collision_tallies % shrink_to_fit() call active_current_tallies % shrink_to_fit() + call active_cell_to_cell_tallies % shrink_to_fit() end subroutine setup_active_cmfdtallies diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 341e44deb5..b70528bf4e 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -70,20 +70,7 @@ module tally_filter procedure :: text_label => text_label_cell procedure :: initialize => initialize_cell end type CellFilter - -!=============================================================================== -! CELLFROMTOFILTER specifies which geometric cells particles exit. -!=============================================================================== - type, extends(TallyFilter) :: CellToCellFilter - integer, allocatable :: cells(:) - type(DictIntInt) :: map - contains - procedure :: get_next_bin => get_next_bin_cell_to_cell - procedure :: to_statepoint => to_statepoint_cell_to_cell - procedure :: text_label => text_label_cell_to_cell - procedure :: initialize => initialize_cell_to_cell - end type CellToCellFilter - + !=============================================================================== ! CELLFROMFILTER specifies which geometric cells particles exit. !=============================================================================== @@ -96,7 +83,7 @@ module tally_filter procedure :: text_label => text_label_cell_from procedure :: initialize => initialize_cell_from end type CellFromFilter - + !=============================================================================== ! CELLTOFILTER specifies which geometric cells particles enter. !=============================================================================== @@ -757,99 +744,10 @@ contains end function text_label_cell !=============================================================================== -! Cell_to_cell_Filter methods +! CellFromFilter methods !=============================================================================== - subroutine get_next_bin_cell_to_cell(this, p, estimator, current_bin, & + subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & next_bin, weight) - - class(CellToCellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - integer, value, intent(in) :: current_bin - integer, intent(out) :: next_bin - real(8), intent(out) :: weight - - integer :: i, start - - ! Find the coordinate level of the last bin we found. - if (current_bin == NO_BIN_FOUND) then - start = 1 - else - do i = 1, p % n_coord - if (p % coord(i) % cell == this % cells(current_bin)) then - start = i + 1 - exit - end if - end do - end if - - ! Starting one coordinate level deeper, find the next bin. !!!!!!!!!!!!!!!!!! OLD !!!!!!!!!!!!!!!!!!!!!!!!!! - next_bin = NO_BIN_FOUND - weight = ERROR_REAL - do i = start, p % n_coord - if (this % map % has_key(p % coord(i) % cell) .and. & - this % map % has_key(p % last_cell)) then - next_bin = this % map % get_key(p % coord(i) % cell) - weight = ONE - exit - end if - end do - end subroutine get_next_bin_cell_to_cell - - subroutine to_statepoint_cell_to_cell(this, filter_group) - class(CellToCellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - integer :: i - integer, allocatable :: cell_ids(:) - - call write_dataset(filter_group, "type", "cell") - call write_dataset(filter_group, "n_bins", this % n_bins) - - allocate(cell_ids(size(this % cells))) - do i = 1, size(this % cells) - cell_ids(i) = cells(this % cells(i)) % id - end do - call write_dataset(filter_group, "bins", cell_ids) - end subroutine to_statepoint_cell_to_cell - - subroutine initialize_cell_to_cell(this) - class(CellToCellFilter), intent(inout) :: this - - integer :: i, id - - ! Convert ids to indices. - do i = 1, this % n_bins - id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) - if (cell_dict % has_key(id)) then - this % cells(i) = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end do - - ! Generate mapping from cell indices to filter bins. - do i = 1, this % n_bins - call this % map % add_key(this % cells(i), i) - end do - end subroutine initialize_cell_to_cell - - function text_label_cell_to_cell(this, bin) result(label) - class(CellToCellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - label = "Cell to from" // to_str(cells(this % cells(bin)) % id) - end function text_label_cell_to_cell - -!=============================================================================== -! Cell_From_Filter methods -!=============================================================================== - subroutine get_next_bin_cell_from(this, p, estimator, current_bin, & - next_bin, weight) - class(CellFromFilter), intent(in) :: this type(Particle), intent(in) :: p integer, intent(in) :: estimator @@ -859,45 +757,24 @@ contains integer :: i, start - ! Find the coordinate level of the last bin we found. -! if (current_bin == NO_BIN_FOUND) then -! start = 1 -! else -! do i = 1, p % last_n_coord !shouldnt be done this way, modify later!! -! ! print *, p % last_cell, this % cells(current_bin) -! if (p % last_cell == this % cells(current_bin)) then -! start = i + 1 -! exit -! end if -! end do -! end if - - ! Starting one coordinate level deeper, find the next bin. + ! Particle can only have one cell_from + ! If current_bin is already a bin, then no need to look for another bin next_bin = NO_BIN_FOUND weight = ERROR_REAL - weight = ONE if (current_bin == NO_BIN_FOUND) then -! print *, "Looking in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) if (this % map % has_key(p % last_cell)) then -! print *, "Found in map for cell at", p % last_cell," :",this % map % get_key(p % last_cell) next_bin = this % map % get_key(p % last_cell) + weight = ONE end if else next_bin = NO_BIN_FOUND end if - -! do i = start, p % last_n_coord -! if (this % map % has_key(p % last_cell)) then -! next_bin = this % map % get_key(p % coord(i) % last_cell) !need to use coord(i) to dig one level deeper -! weight = ONE -! exit -! end if -! end do + end subroutine get_next_bin_cell_from subroutine to_statepoint_cell_from(this, filter_group) class(CellFromFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group + integer(HID_T), intent(in) :: filter_group integer :: i integer, allocatable :: cell_ids(:) @@ -920,7 +797,6 @@ contains ! Convert ids to indices. do i = 1, this % n_bins id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) if (cell_dict % has_key(id)) then this % cells(i) = cell_dict % get_key(id) else @@ -931,31 +807,29 @@ contains ! Generate mapping from cell indices to filter bins. do i = 1, this % n_bins -! print *, "Filling bin/cell map : cell",this % cells(i), "bin", i call this % map % add_key(this % cells(i), i) end do end subroutine initialize_cell_from - + function text_label_cell_from(this, bin) result(label) class(CellFromFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label - label = "Cell to from" // to_str(cells(this % cells(bin)) % id) + label = "Cell from " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_from - + !=============================================================================== -! Cell_To_Filter methods +! CellToFilter methods !=============================================================================== - subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & + subroutine get_next_bin_cell_to(this, p, estimator, current_bin, & next_bin, weight) - class(CellToFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - integer, value, intent(in) :: current_bin - integer, intent(out) :: next_bin - real(8), intent(out) :: weight + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight integer :: i, start @@ -985,7 +859,7 @@ contains subroutine to_statepoint_cell_to(this, filter_group) class(CellToFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group + integer(HID_T), intent(in) :: filter_group integer :: i integer, allocatable :: cell_ids(:) @@ -1008,7 +882,7 @@ contains ! Convert ids to indices. do i = 1, this % n_bins id = this % cells(i) -! print *,this % n_bins, this % cells(i), cell_dict % get_key(id) + if (cell_dict % has_key(id)) then this % cells(i) = cell_dict % get_key(id) else @@ -1022,16 +896,15 @@ contains call this % map % add_key(this % cells(i), i) end do end subroutine initialize_cell_to - + function text_label_cell_to(this, bin) result(label) class(CellToFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label - label = "Cell to to" // to_str(cells(this % cells(bin)) % id) + label = "Cell to " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_to - - + !=============================================================================== ! DistribcellFilter methods !=============================================================================== diff --git a/src/tracking.F90 b/src/tracking.F90 index e9fbc86623..4a87b4ffa8 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -149,11 +149,10 @@ contains ! PARTICLE CROSSES SURFACE if (next_level > 0) p % n_coord = next_level - + ! Saving previous cell data last_cell = p % coord(p % n_coord) % cell - last_n_coord = p % n_coord - + p % coord(p % n_coord) % cell = NONE if (any(lattice_translation /= 0)) then ! Particle crosses lattice boundary @@ -163,12 +162,10 @@ contains else ! Particle crosses surface p % surface = surface_crossed - + ! /CHANGE/ Saving stuff on particle p % last_cell = last_cell - p % last_n_coord = last_n_coord - p % coord(next_level) % last_cell = last_cell - + call cross_surface(p, last_cell) p % event = EVENT_SURFACE end if