From c27c9bad6c761bc3fb959e2955de66e78f56a8c8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 25 May 2016 18:41:09 -0400 Subject: [PATCH] Use simpler scheme for removing non-WMP XS --- src/multipole.F90 | 103 ++++++++++++++-------------------------------- 1 file changed, 30 insertions(+), 73 deletions(-) diff --git a/src/multipole.F90 b/src/multipole.F90 index 59c7fed69..0ded8846a 100644 --- a/src/multipole.F90 +++ b/src/multipole.F90 @@ -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