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

Merge branch 'multipole2' into diff_tally5

Browse source
This commit is contained in:
Sterling Harper 2016-05-25 18:42:51 -04:00
commit 6feaa9be00
1 changed files with 30 additions and 73 deletions
Download patch file Download diff file Expand all files Collapse all files

103
src/multipole.F90
View file

@ -27,11 +27,9 @@ contains
integer(HID_T) :: file_id
integer(HID_T) :: group_id
integer :: is_fissionable
real(8) :: insert_pts(4) ! New points in the energy grid
integer :: cut1, cut2 ! Old indices just outside MP region
integer :: new_n_grid ! Number of points in new E grid
real(8), allocatable :: new_energy(:) ! New energy grid
real(8) :: f1, f2 ! Interpolation near cut1 & cut2
real(8), allocatable :: new_xs(:) ! New cross sections
integer :: i
integer :: IE ! Reaction threshold
@ -90,84 +88,49 @@ contains
! interpolated to these points. The other two points are used to zero the
! cross sections inside the multipole region.
! Define the four new inserted points.
insert_pts(:) = [multipole % start_E / 1e6_8, &
multipole % start_E / 1e6_8 + 1e-12_8, &
multipole % end_E / 1e6_8 - 1e-12_8, &
multipole % end_E / 1e6_8]
! Find the points just outside the multipole region.
cut1 = binary_search(nuc % energy, nuc % n_grid, insert_pts(1))
cut2 = binary_search(nuc % energy, nuc % n_grid, insert_pts(4)) + 1
if (nuc % energy(cut1) == insert_pts(1)) cut1 = cut1 - 1
if (nuc % energy(cut2) == insert_pts(4)) cut2 = cut2 + 1
cut1 = binary_search(nuc % energy, nuc % n_grid, insert_pts(1)) + 2
cut2 = binary_search(nuc % energy, nuc % n_grid, insert_pts(4)) - 1
! Generate the new energy grid.
new_n_grid = nuc % n_grid - (cut2 - cut1 - 1) + 4
new_n_grid = nuc % n_grid - (cut2 - cut1 - 1)
allocate(new_energy(new_n_grid))
new_energy(1:cut1) = nuc % energy(1:cut1)
new_energy(cut1+1:cut1+4) = insert_pts(:)
new_energy(cut1+5:new_n_grid) = nuc % energy(cut2:nuc % n_grid)
! Compute interpolation factors for the new energy points.
f1 = (insert_pts(1) - nuc % energy(cut1)) &
/ (nuc % energy(cut1+1) - nuc % energy(cut1))
f2 = (insert_pts(4) - nuc % energy(cut2-1)) &
/ (nuc % energy(cut2) - nuc % energy(cut2-1))
new_energy(cut1+1:new_n_grid) = nuc % energy(cut2:nuc % n_grid)
! Adjust the total cross section.
allocate(new_xs(new_n_grid))
new_xs(1:cut1) = nuc % total(1:cut1)
new_xs(cut1+1) = (ONE - f1) * nuc % total(cut1) &
+ f1 * nuc % total(cut1+1)
new_xs(cut1+2:cut1+3) = ZERO
new_xs(cut1+4) = (ONE - f2) * nuc % total(cut2-1) &
+ f2 * nuc % total(cut2)
new_xs(cut1+5:new_n_grid) = nuc % total(cut2:nuc % n_grid)
new_xs(1:cut1-1) = nuc % total(1:cut1-1)
new_xs(cut1:cut1+1) = ZERO
new_xs(cut1+2:new_n_grid) = nuc % total(cut2+1:nuc % n_grid)
call move_alloc(new_xs, nuc % total)
! Adjust the elastic cross section.
allocate(new_xs(new_n_grid))
new_xs(1:cut1) = nuc % elastic(1:cut1)
new_xs(cut1+1) = (ONE - f1) * nuc % elastic(cut1) &
+ f1 * nuc % elastic(cut1+1)
new_xs(cut1+2:cut1+3) = ZERO
new_xs(cut1+4) = (ONE - f2) * nuc % elastic(cut2-1) &
+ f2 * nuc % elastic(cut2)
new_xs(cut1+5:new_n_grid) = nuc % elastic(cut2:nuc % n_grid)
new_xs(1:cut1-1) = nuc % elastic(1:cut1-1)
new_xs(cut1:cut1+1) = ZERO
new_xs(cut1+2:new_n_grid) = nuc % elastic(cut2+1:nuc % n_grid)
call move_alloc(new_xs, nuc % elastic)
! Adjust the fission cross section.
allocate(new_xs(new_n_grid))
new_xs(1:cut1) = nuc % fission(1:cut1)
new_xs(cut1+1) = (ONE - f1) * nuc % fission(cut1) &
+ f1 * nuc % fission(cut1+1)
new_xs(cut1+2:cut1+3) = ZERO
new_xs(cut1+4) = (ONE - f2) * nuc % fission(cut2-1) &
+ f2 * nuc % fission(cut2)
new_xs(cut1+5:new_n_grid) = nuc % fission(cut2:nuc % n_grid)
new_xs(1:cut1-1) = nuc % fission(1:cut1-1)
new_xs(cut1:cut1+1) = ZERO
new_xs(cut1+2:new_n_grid) = nuc % fission(cut2+1:nuc % n_grid)
call move_alloc(new_xs, nuc % fission)
! Adjust the nu-fission cross section.
allocate(new_xs(new_n_grid))
new_xs(1:cut1) = nuc % nu_fission(1:cut1)
new_xs(cut1+1) = (ONE - f1) * nuc % nu_fission(cut1) &
+ f1 * nuc % nu_fission(cut1+1)
new_xs(cut1+2:cut1+3) = ZERO
new_xs(cut1+4) = (ONE - f2) * nuc % nu_fission(cut2-1) &
+ f2 * nuc % nu_fission(cut2)
new_xs(cut1+5:new_n_grid) = nuc % nu_fission(cut2:nuc % n_grid)
new_xs(1:cut1-1) = nuc % nu_fission(1:cut1-1)
new_xs(cut1:cut1+1) = ZERO
new_xs(cut1+2:new_n_grid) = nuc % nu_fission(cut2+1:nuc % n_grid)
call move_alloc(new_xs, nuc % nu_fission)
! Adjust the absorption cross section.
allocate(new_xs(new_n_grid))
new_xs(1:cut1) = nuc % absorption(1:cut1)
new_xs(cut1+1) = (ONE - f1) * nuc % absorption(cut1) &
+ f1 * nuc % absorption(cut1+1)
new_xs(cut1+2:cut1+3) = ZERO
new_xs(cut1+4) = (ONE - f2) * nuc % absorption(cut2-1) &
+ f2 * nuc % absorption(cut2)
new_xs(cut1+5:new_n_grid) = nuc % absorption(cut2:nuc % n_grid)
new_xs(1:cut1-1) = nuc % absorption(1:cut1-1)
new_xs(cut1:cut1+1) = ZERO
new_xs(cut1+2:new_n_grid) = nuc % absorption(cut2+1:nuc % n_grid)
call move_alloc(new_xs, nuc % absorption)
! Adjust other cross sections.
@ -178,30 +141,24 @@ contains
if (rxn % threshold >= cut2) then
! The threshold is above the multipole range. All we need to do
! is adjust the threshold index to match the new grid.
rxn % threshold = rxn % threshold - (cut2 - cut1 - 1) + 4
rxn % threshold = rxn % threshold - (cut2 - cut1 - 1)
else if (rxn % threshold <= cut1) then
! The threhold is below the multipole range. Remove the multipole
! region just like we did with the other reactions.
! The new grid removed (cut2 - cut1 - 1) points and added 4.
allocate(new_xs(size(rxn % sigma) - (cut2 - cut1 - 1) + 4))
new_xs(1:cut1-IE+1) = rxn % sigma(1:cut1-IE+1)
new_xs(cut1-IE+2) = (ONE - f1) * rxn % sigma(cut1-IE+1) &
+ f1 * rxn % sigma(cut1-IE+2)
new_xs(cut1-IE+3:cut1-IE+4) = ZERO
new_xs(cut1-IE+5) = (ONE - f2) * rxn % sigma(cut2-IE) &
+ f2 * rxn % sigma(cut2-IE+1)
new_xs(cut1-IE+6:size(new_xs)) = &
rxn % sigma(cut2-IE+1:size(rxn % sigma))
! The new grid removed (cut2 - cut1 - 1) points.
allocate(new_xs(size(rxn % sigma) - (cut2 - cut1 - 1)))
new_xs(1:cut1-IE) = rxn % sigma(1:cut1-IE)
new_xs(cut1-IE+1:cut1-IE+2) = ZERO
new_xs(cut1-IE+3:size(new_xs)) = &
rxn % sigma(cut2-IE+2:size(rxn % sigma))
call move_alloc(new_xs, rxn % sigma)
else
! The threshold lies within the multipole range. Remove the first
! cut2-IE points and add an interpolated point
allocate(new_xs(size(rxn % sigma) - (cut2-IE) + 1))
new_xs(1) = (ONE - f2) * rxn % sigma(cut2-IE) &
+ f2 * rxn % sigma(cut2-IE+1)
new_xs(2:size(new_xs)) = rxn % sigma(cut2-IE+1:size(rxn % sigma))
! cut2-IE+1 points.
allocate(new_xs(size(rxn % sigma) - (cut2-IE+1)))
new_xs(1:size(new_xs)) = rxn % sigma(cut2-IE+2:size(rxn % sigma))
call move_alloc(new_xs, rxn % sigma)
rxn % threshold = cut1 + 4
rxn % threshold = cut2 + 1
end if
end associate
end do