mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge branch 'nuclidescatt' into mgxs_integration
This commit is contained in:
Adam Nelson
commit
c779ca42c4
3 changed files with 242 additions and 233 deletions
|
|
@ -201,7 +201,7 @@ module mgxs_header
|
|||
! the xsdata object node itself.
|
||||
!===============================================================================
|
||||
|
||||
subroutine mgxs_init_file(this,node_xsdata,i_listing)
|
||||
subroutine mgxs_init_file(this, node_xsdata, i_listing)
|
||||
class(Mgxs), intent(inout) :: this ! Working Object
|
||||
type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml
|
||||
integer, intent(in) :: i_listing ! Index in listings array
|
||||
|
|
@ -236,7 +236,7 @@ module mgxs_header
|
|||
else if (temp_str == 'tabular') then
|
||||
this % scatt_type = ANGLE_TABULAR
|
||||
else
|
||||
call fatal_error("Invalid Scatt Type Option!")
|
||||
call fatal_error("Invalid scatt_type option!")
|
||||
end if
|
||||
else
|
||||
this % scatt_type = ANGLE_LEGENDRE
|
||||
|
|
@ -259,8 +259,8 @@ module mgxs_header
|
|||
|
||||
end subroutine mgxs_init_file
|
||||
|
||||
subroutine mgxsiso_init_file(this,node_xsdata,groups,get_kfiss,get_fiss, &
|
||||
max_order,i_listing)
|
||||
subroutine mgxsiso_init_file(this, node_xsdata, groups, get_kfiss, get_fiss, &
|
||||
max_order, i_listing)
|
||||
class(MgxsIso), intent(inout) :: this ! Working Object
|
||||
type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml
|
||||
integer, intent(in) :: groups ! Number of Energy groups
|
||||
|
|
@ -282,7 +282,7 @@ module mgxs_header
|
|||
integer :: legendre_mu_points, imu
|
||||
|
||||
! Call generic data gathering routine (will populate the metadata)
|
||||
call mgxs_init_file(this,node_xsdata,i_listing)
|
||||
call mgxs_init_file(this, node_xsdata, i_listing)
|
||||
|
||||
! Load the more specific data
|
||||
allocate(this % nu_fission(groups))
|
||||
|
|
@ -292,7 +292,7 @@ module mgxs_header
|
|||
! Chi was provided, that means they are giving chi and nu-fission
|
||||
! vectors
|
||||
! Get chi
|
||||
allocate(temp_arr(1 * groups))
|
||||
allocate(temp_arr(groups))
|
||||
call get_node_array(node_xsdata, "chi", temp_arr)
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
|
|
@ -329,7 +329,7 @@ module mgxs_header
|
|||
end do
|
||||
|
||||
! Now pull out information needed for chi
|
||||
this % chi = temp_2d
|
||||
this % chi(:, :) = temp_2d
|
||||
! Normalize chi so its CDF goes to 1
|
||||
do gin = 1, groups
|
||||
this % chi(:, gin) = this % chi(:, gin) / sum(this % chi(:, gin))
|
||||
|
|
@ -375,14 +375,14 @@ module mgxs_header
|
|||
if (arr_len == groups * groups) then
|
||||
allocate(temp_arr(arr_len))
|
||||
call get_node_array(node_xsdata, "multiplicity", temp_arr)
|
||||
temp_mult = reshape(temp_arr, (/groups, groups/))
|
||||
temp_mult(:, :) = reshape(temp_arr, (/groups, groups/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("Multiplicity length not same as number of groups&
|
||||
& squared!")
|
||||
end if
|
||||
else
|
||||
temp_mult = ONE
|
||||
temp_mult(:, :) = ONE
|
||||
end if
|
||||
|
||||
! Get scattering treatment information
|
||||
|
|
@ -426,7 +426,7 @@ module mgxs_header
|
|||
if (check_for_node(node_xsdata, "order")) then
|
||||
call get_node_value(node_xsdata, "order", order)
|
||||
else
|
||||
call fatal_error("Order Must Be Provided!")
|
||||
call fatal_error("Order must be provided!")
|
||||
end if
|
||||
|
||||
! Before retrieving the data, store the dimensionality of the data in
|
||||
|
|
@ -546,11 +546,11 @@ module mgxs_header
|
|||
if (check_for_node(node_xsdata, "total")) then
|
||||
call get_node_array(node_xsdata, "total", this % total)
|
||||
else
|
||||
this % total = this % absorption + this % scatter % scattxs
|
||||
this % total(:) = this % absorption(:) + this % scatter % scattxs(:)
|
||||
end if
|
||||
|
||||
! Deallocate temporaries for the next material
|
||||
deallocate(input_scatt,scatt_coeffs,temp_mult)
|
||||
deallocate(input_scatt, scatt_coeffs, temp_mult)
|
||||
|
||||
! Finally, check sigT to ensure it is not 0 since it is
|
||||
! often divided by in the tally routines
|
||||
|
|
@ -561,8 +561,8 @@ module mgxs_header
|
|||
|
||||
end subroutine mgxsiso_init_file
|
||||
|
||||
subroutine mgxsang_init_file(this,node_xsdata,groups,get_kfiss,get_fiss, &
|
||||
max_order,i_listing)
|
||||
subroutine mgxsang_init_file(this, node_xsdata, groups, get_kfiss, get_fiss, &
|
||||
max_order, i_listing)
|
||||
class(MgxsAngle), intent(inout) :: this ! Working Object
|
||||
type(Node), pointer, intent(in) :: node_xsdata ! Data from MGXS xml
|
||||
integer, intent(in) :: groups ! Number of Energy groups
|
||||
|
|
@ -584,18 +584,18 @@ module mgxs_header
|
|||
integer :: legendre_mu_points, imu, ipol, iazi
|
||||
|
||||
! Call generic data gathering routine (will populate the metadata)
|
||||
call mgxs_init_file(this,node_xsdata,i_listing)
|
||||
call mgxs_init_file(this, node_xsdata, i_listing)
|
||||
|
||||
if (check_for_node(node_xsdata, "num_polar")) then
|
||||
call get_node_value(node_xsdata, "num_polar", this % n_pol)
|
||||
else
|
||||
call fatal_error("num_polar Must Be Provided!")
|
||||
call fatal_error("num_polar must be provided!")
|
||||
end if
|
||||
|
||||
if (check_for_node(node_xsdata, "num_azimuthal")) then
|
||||
call get_node_value(node_xsdata, "num_azimuthal", this % n_azi)
|
||||
else
|
||||
call fatal_error("num_azimuthal Must Be Provided!")
|
||||
call fatal_error("num_azimuthal must be provided!")
|
||||
end if
|
||||
|
||||
! Load angle data, if present (else equally spaced)
|
||||
|
|
@ -663,8 +663,8 @@ module mgxs_header
|
|||
if (check_for_node(node_xsdata, "nu_fission")) then
|
||||
allocate(temp_arr(groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "nu_fission", temp_arr)
|
||||
this % nu_fission = reshape(temp_arr,(/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
this % nu_fission(:, :, :) = reshape(temp_arr, (/groups, &
|
||||
this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("If fissionable, must provide nu_fission!")
|
||||
|
|
@ -677,8 +677,8 @@ module mgxs_header
|
|||
allocate(temp_arr(groups * groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "nu_fission", temp_arr)
|
||||
allocate(temp_4d(groups, groups, this % n_azi,this % n_pol))
|
||||
temp_4d = reshape(temp_arr, (/groups, groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
temp_4d(:, :, :, :) = reshape(temp_arr, (/groups, groups, &
|
||||
this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("If fissionable, must provide nu_fission!")
|
||||
|
|
@ -716,8 +716,8 @@ module mgxs_header
|
|||
allocate(temp_arr(groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "fission", temp_arr)
|
||||
allocate(this % fission(groups, this % n_azi, this % n_pol))
|
||||
this % fission = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
this % fission(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("Fission data missing, required due to fission&
|
||||
|
|
@ -729,8 +729,8 @@ module mgxs_header
|
|||
allocate(temp_arr(groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "kappa_fission", temp_arr)
|
||||
allocate(this % k_fission(groups, this % n_azi, this % n_pol))
|
||||
this % k_fission = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
this % k_fission(:, :, :) = reshape(temp_arr, (/groups, &
|
||||
this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("kappa_fission data missing, required due to &
|
||||
|
|
@ -738,16 +738,16 @@ module mgxs_header
|
|||
end if
|
||||
end if
|
||||
else
|
||||
this % nu_fission = ZERO
|
||||
this % chi = ZERO
|
||||
this % nu_fission(:, :, :) = ZERO
|
||||
this % chi(:, :, :, :) = ZERO
|
||||
end if
|
||||
|
||||
if (check_for_node(node_xsdata, "absorption")) then
|
||||
allocate(temp_arr(groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "absorption", temp_arr)
|
||||
allocate(this % absorption(groups, this % n_azi, this % n_pol))
|
||||
this % absorption = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
this % absorption(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("Must provide absorption!")
|
||||
|
|
@ -760,15 +760,15 @@ module mgxs_header
|
|||
if (arr_len == groups * groups * this % n_azi * this % n_pol) then
|
||||
allocate(temp_arr(arr_len))
|
||||
call get_node_array(node_xsdata, "multiplicity", temp_arr)
|
||||
temp_mult = reshape(temp_arr, (/groups, groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
temp_mult(:, :, :, :) = reshape(temp_arr, (/groups, groups, &
|
||||
this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
call fatal_error("Multiplicity length not same as number of groups&
|
||||
& squared!")
|
||||
end if
|
||||
else
|
||||
temp_mult = ONE
|
||||
temp_mult(:, :, :, :) = ONE
|
||||
end if
|
||||
|
||||
! Get scattering treatment information
|
||||
|
|
@ -812,7 +812,7 @@ module mgxs_header
|
|||
if (check_for_node(node_xsdata, "order")) then
|
||||
call get_node_value(node_xsdata, "order", order)
|
||||
else
|
||||
call fatal_error("Order Must Be Provided!")
|
||||
call fatal_error("Order must be provided!")
|
||||
end if
|
||||
|
||||
! Before retrieving the data, store the dimensionality of the data in
|
||||
|
|
@ -836,8 +836,8 @@ module mgxs_header
|
|||
allocate(temp_arr(groups * groups * order_dim * this % n_azi * &
|
||||
this % n_pol))
|
||||
call get_node_array(node_xsdata, "scatter", temp_arr)
|
||||
input_scatt = reshape(temp_arr, (/groups, groups, order_dim, &
|
||||
this % n_azi, this % n_pol/))
|
||||
input_scatt(:, :, :, :, :) = reshape(temp_arr, (/groups, groups, &
|
||||
order_dim, this % n_azi, this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
|
||||
! Compare the number of orders given with the maximum order of the
|
||||
|
|
@ -951,8 +951,8 @@ module mgxs_header
|
|||
if (check_for_node(node_xsdata, "total")) then
|
||||
allocate(temp_arr(groups * this % n_azi * this % n_pol))
|
||||
call get_node_array(node_xsdata, "total", temp_arr)
|
||||
this % total = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
this % total(:, :, :) = reshape(temp_arr, (/groups, this % n_azi, &
|
||||
this % n_pol/))
|
||||
deallocate(temp_arr)
|
||||
else
|
||||
do ipol = 1, this % n_pol
|
||||
|
|
@ -1251,7 +1251,7 @@ module mgxs_header
|
|||
if (present(gout)) then
|
||||
xs = this % chi(gout, gin, iazi, ipol)
|
||||
else
|
||||
! Not sure youd want a 1 or a 0, but here you go!
|
||||
! Not sure you would want a 1 or a 0, but here you go!
|
||||
xs = sum(this % chi(:, gin, iazi, ipol))
|
||||
end if
|
||||
case('scatter')
|
||||
|
|
@ -1316,7 +1316,7 @@ module mgxs_header
|
|||
! objects
|
||||
!===============================================================================
|
||||
|
||||
subroutine mgxs_combine(this,mat,scatt_type,i_listing)
|
||||
subroutine mgxs_combine(this, mat, scatt_type, i_listing)
|
||||
class(Mgxs), intent(inout) :: this ! The Mgxs to initialize
|
||||
type(Material), pointer, intent(in) :: mat ! base material
|
||||
integer, intent(in) :: scatt_type ! How is data presented
|
||||
|
|
@ -1340,7 +1340,7 @@ module mgxs_header
|
|||
|
||||
end subroutine mgxs_combine
|
||||
|
||||
subroutine mgxsiso_combine(this,mat,nuclides,groups,max_order,scatt_type, &
|
||||
subroutine mgxsiso_combine(this, mat, nuclides, groups, max_order, scatt_type, &
|
||||
i_listing)
|
||||
class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize
|
||||
type(Material), pointer, intent(in) :: mat ! base material
|
||||
|
|
@ -1373,9 +1373,9 @@ module mgxs_header
|
|||
do i = 2, mat % n_nuclides
|
||||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsIso)
|
||||
if (order /= size(nuc % scatter % dist(1) % data,dim=1)) &
|
||||
call fatal_error("All Histogram Scattering Entries Must Be&
|
||||
& Same Length!")
|
||||
if (order /= size(nuc % scatter % dist(1) % data, dim=1)) &
|
||||
call fatal_error("All histogram scattering entries must be&
|
||||
& same length!")
|
||||
end select
|
||||
end do
|
||||
! Ok, got our order, store the dimensionality
|
||||
|
|
@ -1390,8 +1390,8 @@ module mgxs_header
|
|||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsIso)
|
||||
if (order /= size(nuc % scatter % dist(1) % data, dim=1)) &
|
||||
call fatal_error("All Tabular Scattering Entries Must Be&
|
||||
& Same Length!")
|
||||
call fatal_error("All tabular scattering entries must be&
|
||||
& same length!")
|
||||
end select
|
||||
end do
|
||||
! Ok, got our order, store the dimensionality
|
||||
|
|
@ -1406,7 +1406,7 @@ module mgxs_header
|
|||
do i = 1, mat % n_nuclides
|
||||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsIso)
|
||||
if (size(nuc % scatter % dist(1) % data,dim=1) > mat_max_order) &
|
||||
if (size(nuc % scatter % dist(1) % data, dim=1) > mat_max_order) &
|
||||
mat_max_order = size(nuc % scatter % dist(1) % data, dim=1)
|
||||
end select
|
||||
end do
|
||||
|
|
@ -1423,25 +1423,25 @@ module mgxs_header
|
|||
|
||||
! Allocate and initialize data needed for macro_xs(i_mat) object
|
||||
allocate(this % total(groups))
|
||||
this % total = ZERO
|
||||
this % total(:) = ZERO
|
||||
allocate(this % absorption(groups))
|
||||
this % absorption = ZERO
|
||||
this % absorption(:) = ZERO
|
||||
allocate(this % fission(groups))
|
||||
this % fission = ZERO
|
||||
this % fission(:) = ZERO
|
||||
allocate(this % k_fission(groups))
|
||||
this % k_fission = ZERO
|
||||
this % k_fission(:) = ZERO
|
||||
allocate(this % nu_fission(groups))
|
||||
this % nu_fission = ZERO
|
||||
this % nu_fission(:) = ZERO
|
||||
allocate(this % chi(groups,groups))
|
||||
this % chi = ZERO
|
||||
this % chi(:, :) = ZERO
|
||||
allocate(temp_mult(groups,groups))
|
||||
temp_mult = ZERO
|
||||
temp_mult(:, :) = ZERO
|
||||
allocate(mult_num(groups,groups))
|
||||
mult_num = ZERO
|
||||
mult_num(:, :) = ZERO
|
||||
allocate(mult_denom(groups,groups))
|
||||
mult_denom = ZERO
|
||||
mult_denom(:, :) = ZERO
|
||||
allocate(scatt_coeffs(order_dim,groups,groups))
|
||||
scatt_coeffs = ZERO
|
||||
scatt_coeffs(:, :, :) = ZERO
|
||||
|
||||
! Add contribution from each nuclide in material
|
||||
do i = 1, mat % n_nuclides
|
||||
|
|
@ -1452,19 +1452,19 @@ module mgxs_header
|
|||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsIso)
|
||||
! Add contributions to total, absorption, and fission data (if necessary)
|
||||
this % total = this % total + atom_density * nuc % total
|
||||
this % absorption = this % absorption + &
|
||||
atom_density * nuc % absorption
|
||||
this % total(:) = this % total(:) + atom_density * nuc % total(:)
|
||||
this % absorption(:) = this % absorption(:) + &
|
||||
atom_density * nuc % absorption(:)
|
||||
if (nuc % fissionable) then
|
||||
this % chi = this % chi + atom_density * nuc % chi
|
||||
this % nu_fission = this % nu_fission + atom_density * &
|
||||
nuc % nu_fission
|
||||
this % chi(:, :) = this % chi(:, :) + atom_density * nuc % chi(:, :)
|
||||
this % nu_fission(:) = this % nu_fission(:)+ atom_density * &
|
||||
nuc % nu_fission(:)
|
||||
if (allocated(nuc % fission)) then
|
||||
this % fission = this % fission + atom_density * nuc % fission
|
||||
this % fission(:) = this % fission(:) + atom_density * nuc % fission(:)
|
||||
end if
|
||||
if (allocated(nuc % k_fission)) then
|
||||
this % k_fission = this % k_fission + atom_density * &
|
||||
nuc % k_fission
|
||||
this % k_fission(:) = this % k_fission(:) + atom_density * &
|
||||
nuc % k_fission(:)
|
||||
end if
|
||||
end if
|
||||
|
||||
|
|
@ -1498,7 +1498,7 @@ module mgxs_header
|
|||
nuc % scatter % get_matrix(min(nuc_order_dim, order_dim))
|
||||
|
||||
type is (MgxsAngle)
|
||||
call fatal_error("Invalid Passing of MgxsAngle to MgxsIso Object")
|
||||
call fatal_error("Invalid passing of MgxsAngle to MgxsIso object")
|
||||
end select
|
||||
end do
|
||||
|
||||
|
|
@ -1531,7 +1531,7 @@ module mgxs_header
|
|||
|
||||
end subroutine mgxsiso_combine
|
||||
|
||||
subroutine mgxsang_combine(this,mat,nuclides,groups,max_order,scatt_type,&
|
||||
subroutine mgxsang_combine(this, mat, nuclides, groups, max_order, scatt_type, &
|
||||
i_listing)
|
||||
class(MgxsAngle), intent(inout) :: this ! The Mgxs to initialize
|
||||
type(Material), pointer, intent(in) :: mat ! base material
|
||||
|
|
@ -1551,7 +1551,7 @@ module mgxs_header
|
|||
real(8), allocatable :: mult_denom(:, :, :, :), scatt_coeffs(:, :, :, :, :)
|
||||
|
||||
! Set the meta-data
|
||||
call mgxs_combine(this,mat,scatt_type,i_listing)
|
||||
call mgxs_combine(this, mat, scatt_type, i_listing)
|
||||
|
||||
! Get the number of each polar and azi angles and make sure all the
|
||||
! NuclideAngle types have the same number of these angles
|
||||
|
|
@ -1564,12 +1564,12 @@ module mgxs_header
|
|||
n_pol = nuc % n_pol
|
||||
n_azi = nuc % n_azi
|
||||
allocate(this % polar(n_pol))
|
||||
this % polar = nuc % polar
|
||||
this % polar(:) = nuc % polar(:)
|
||||
allocate(this % azimuthal(n_azi))
|
||||
this % azimuthal = nuc % azimuthal
|
||||
this % azimuthal(:) = nuc % azimuthal(:)
|
||||
else
|
||||
if ((n_pol /= nuc % n_pol) .or. (n_azi /= nuc % n_azi)) then
|
||||
call fatal_error("All Angular Data Must Be Same Length!")
|
||||
call fatal_error("All angular data must be same length!")
|
||||
end if
|
||||
end if
|
||||
end select
|
||||
|
|
@ -1589,8 +1589,8 @@ module mgxs_header
|
|||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsAngle)
|
||||
if (order /= size(nuc % scatter(1,1) % obj % dist(1) % data, dim=1)) &
|
||||
call fatal_error("All Histogram Scattering Entries Must Be&
|
||||
& Same Length!")
|
||||
call fatal_error("All histogram scattering entries must be&
|
||||
& same length!")
|
||||
end select
|
||||
end do
|
||||
! Ok, got our order, store the dimensionality
|
||||
|
|
@ -1609,9 +1609,9 @@ module mgxs_header
|
|||
do i = 2, mat % n_nuclides
|
||||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsAngle)
|
||||
if (order /= size(nuc % scatter(1, 1) % obj % dist(1) % data,dim=1)) &
|
||||
call fatal_error("All Tabular Scattering Entries Must Be&
|
||||
& Same Length!")
|
||||
if (order /= size(nuc % scatter(1, 1) % obj % dist(1) % data, dim=1)) &
|
||||
call fatal_error("All tabular scattering entries must be&
|
||||
& same length!")
|
||||
end select
|
||||
end do
|
||||
! Ok, got our order, store the dimensionality
|
||||
|
|
@ -1653,25 +1653,25 @@ module mgxs_header
|
|||
|
||||
! Allocate and initialize data within macro_xs(i_mat) object
|
||||
allocate(this % total(groups, n_azi, n_pol))
|
||||
this % total = ZERO
|
||||
this % total(:, :, :) = ZERO
|
||||
allocate(this % absorption(groups, n_azi, n_pol))
|
||||
this % absorption = ZERO
|
||||
this % absorption(:, :, :) = ZERO
|
||||
allocate(this % fission(groups, n_azi, n_pol))
|
||||
this % fission = ZERO
|
||||
this % fission(:, :, :) = ZERO
|
||||
allocate(this % k_fission(groups, n_azi, n_pol))
|
||||
this % k_fission = ZERO
|
||||
this % k_fission(:, :, :) = ZERO
|
||||
allocate(this % nu_fission(groups, n_azi, n_pol))
|
||||
this % nu_fission = ZERO
|
||||
this % nu_fission(:, :, :) = ZERO
|
||||
allocate(this % chi(groups, groups, n_azi, n_pol))
|
||||
this % chi = ZERO
|
||||
this % chi(:, :, :, :) = ZERO
|
||||
allocate(temp_mult(groups, groups, n_azi, n_pol))
|
||||
temp_mult = ZERO
|
||||
temp_mult(:, :, :, :) = ZERO
|
||||
allocate(mult_num(groups, groups, n_azi, n_pol))
|
||||
mult_num = ZERO
|
||||
mult_num(:, :, :, :) = ZERO
|
||||
allocate(mult_denom(groups, groups, n_azi, n_pol))
|
||||
mult_denom = ZERO
|
||||
mult_denom(:, :, :, :) = ZERO
|
||||
allocate(scatt_coeffs(order_dim, groups, groups, n_azi, n_pol))
|
||||
scatt_coeffs = ZERO
|
||||
scatt_coeffs(:, :, :, :, :) = ZERO
|
||||
|
||||
! Add contribution from each nuclide in material
|
||||
do i = 1, mat % n_nuclides
|
||||
|
|
@ -1681,22 +1681,24 @@ module mgxs_header
|
|||
! Perform our operations which depend upon the type
|
||||
select type(nuc => nuclides(mat % nuclide(i)) % obj)
|
||||
type is (MgxsIso)
|
||||
call fatal_error("Invalid Passing of MgxsIso to MgxsAngle Object")
|
||||
call fatal_error("Invalid passing of MgxsIso to MgxsAngle object")
|
||||
type is (MgxsAngle)
|
||||
! Add contributions to total, absorption, and fission data (if necessary)
|
||||
this % total = this % total + atom_density * nuc % total
|
||||
this % absorption = this % absorption + &
|
||||
atom_density * nuc % absorption
|
||||
this % total(:, :, :) = this % total(:, :, :) + &
|
||||
atom_density * nuc % total(:, :, :)
|
||||
this % absorption(:, :, :) = this % absorption(:, :, :) + &
|
||||
atom_density * nuc % absorption(:, :, :)
|
||||
if (nuc % fissionable) then
|
||||
this % chi = this % chi + atom_density * nuc % chi
|
||||
this % nu_fission = this % nu_fission + atom_density * &
|
||||
nuc % nu_fission
|
||||
this % nu_fission(:, :, :) = this % nu_fission(:, :, :) + &
|
||||
atom_density * nuc % nu_fission(:, :, :)
|
||||
if (allocated(nuc % fission)) then
|
||||
this % fission = this % fission + atom_density * nuc % fission
|
||||
this % fission(:, :, :) = this % fission(:, :, :) + &
|
||||
atom_density * nuc % fission(:, :, :)
|
||||
end if
|
||||
if (allocated(nuc % k_fission)) then
|
||||
this % k_fission = this % k_fission + atom_density * &
|
||||
nuc % k_fission
|
||||
this % k_fission(:, :, :) = this % k_fission(:, :, :) + &
|
||||
atom_density * nuc % k_fission(:, :, :)
|
||||
end if
|
||||
end if
|
||||
|
||||
|
|
@ -1730,7 +1732,7 @@ module mgxs_header
|
|||
end do
|
||||
|
||||
! Get the complete scattering matrix
|
||||
nuc_order_dim = size(nuc % scatter(1,1) % obj % dist(1) % data,dim=1)
|
||||
nuc_order_dim = size(nuc % scatter(1, 1) % obj % dist(1) % data, dim=1)
|
||||
do ipol = 1, n_pol
|
||||
do iazi = 1, n_azi
|
||||
scatt_coeffs(1:min(nuc_order_dim, order_dim), :, :, iazi, ipol) = &
|
||||
|
|
|
|||
|
|
@ -14,7 +14,7 @@ module scattdata_header
|
|||
!===============================================================================
|
||||
|
||||
type :: Jagged2D
|
||||
real(8), allocatable :: data(:,:)
|
||||
real(8), allocatable :: data(:, :)
|
||||
end type Jagged2D
|
||||
|
||||
type :: Jagged1D
|
||||
|
|
@ -27,12 +27,15 @@ module scattdata_header
|
|||
!===============================================================================
|
||||
|
||||
type, abstract :: ScattData
|
||||
! normalized p0 matrix on its own for sampling energy
|
||||
! The data attribute of the energy, mult, and dist arrays
|
||||
! are not necessarily 1-indexed as they instead will be allocated
|
||||
! from a minimum outgoing group to an outgoing minimum group.
|
||||
! Normalized p0 matrix on its own for sampling energy
|
||||
type(Jagged1D), allocatable :: energy(:) ! (Gin % data(Gout))
|
||||
! nu-scatter multiplication (i.e. nu-scatt/scatt)
|
||||
! Nu-scatter multiplication (i.e. nu-scatt/scatt)
|
||||
type(Jagged1D), allocatable :: mult(:) ! (Gin % data(Gout))
|
||||
! Angular distribution
|
||||
type(Jagged2D), allocatable :: dist(:) ! (Gin % data(Order/Nmu x Gout)
|
||||
type(Jagged2D), allocatable :: dist(:) ! (Gin % data(Order/Nmu, Gout)
|
||||
integer, allocatable :: gmin(:) ! Minimum outgoing group
|
||||
integer, allocatable :: gmax(:) ! Maximum outgoing group
|
||||
real(8), allocatable :: scattxs(:) ! Isotropic Sigma_{s,g_{in}}
|
||||
|
|
@ -47,9 +50,9 @@ module scattdata_header
|
|||
abstract interface
|
||||
subroutine scattdata_init_(this, mult, coeffs)
|
||||
import ScattData
|
||||
class(ScattData), intent(inout) :: this ! Scattering Object to work with
|
||||
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
|
||||
class(ScattData), intent(inout) :: this ! Object to work with
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
end subroutine scattdata_init_
|
||||
|
||||
pure function scattdata_calc_f_(this, gin, gout, mu) result(f)
|
||||
|
|
@ -120,10 +123,10 @@ contains
|
|||
!===============================================================================
|
||||
|
||||
subroutine scattdata_init(this, order, energy, mult)
|
||||
class(ScattData), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: order ! Data Order
|
||||
real(8), intent(inout) :: energy(:,:) ! Energy Transfer Matrix
|
||||
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
|
||||
class(ScattData), intent(inout) :: this ! Object to work on
|
||||
integer, intent(in) :: order ! Data Order
|
||||
real(8), intent(inout) :: energy(:, :) ! Energy Transfer Matrix
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
|
||||
integer :: groups, gmin, gmax, gin
|
||||
real(8) :: norm
|
||||
|
|
@ -139,15 +142,15 @@ contains
|
|||
! Also set energy values when doing it
|
||||
do gin = 1, groups
|
||||
! Make sure energy is normalized (i.e., CDF is 1)
|
||||
norm = sum(energy(:,gin))
|
||||
if (norm /= ZERO) energy(:,gin) = energy(:,gin) / norm
|
||||
norm = sum(energy(:, gin))
|
||||
if (norm /= ZERO) energy(:, gin) = energy(:, gin) / norm
|
||||
! Find gmin by checking the P0 moment
|
||||
do gmin = 1, groups
|
||||
if (energy(gmin,gin) > ZERO) exit
|
||||
if (energy(gmin, gin) > ZERO) exit
|
||||
end do
|
||||
! Find gmax by checking the P0 moment
|
||||
do gmax = groups, 1, -1
|
||||
if (energy(gmax,gin) > ZERO) exit
|
||||
if (energy(gmax, gin) > ZERO) exit
|
||||
end do
|
||||
! Treat the case of all zeros
|
||||
if (gmin > gmax) then
|
||||
|
|
@ -156,10 +159,10 @@ contains
|
|||
! By not changing energy(gin) here we are leaving it as zero
|
||||
end if
|
||||
allocate(this % energy(gin) % data(gmin:gmax))
|
||||
this % energy(gin) % data(gmin:gmax) = energy(gmin:gmax,gin)
|
||||
this % energy(gin) % data(gmin:gmax) = energy(gmin:gmax, gin)
|
||||
allocate(this % mult(gin) % data(gmin:gmax))
|
||||
this % mult(gin) % data(gmin:gmax) = mult(gmin:gmax,gin)
|
||||
allocate(this % dist(gin) % data(order,gmin:gmax))
|
||||
this % mult(gin) % data(gmin:gmax) = mult(gmin:gmax, gin)
|
||||
allocate(this % dist(gin) % data(order, gmin:gmax))
|
||||
this % dist(gin) % data = ZERO
|
||||
this % gmin(gin) = gmin
|
||||
this % gmax(gin) = gmax
|
||||
|
|
@ -167,38 +170,38 @@ contains
|
|||
end subroutine scattdata_init
|
||||
|
||||
subroutine scattdatalegendre_init(this, mult, coeffs)
|
||||
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
|
||||
class(ScattDataLegendre), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
|
||||
real(8) :: dmu, mu, f, norm
|
||||
integer :: imu, Nmu, gout, gin, groups, order
|
||||
real(8), allocatable :: energy(:,:)
|
||||
real(8), allocatable :: matrix(:,:,:)
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
|
||||
groups = size(coeffs,dim=3)
|
||||
order = size(coeffs,dim=1)
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order,groups,groups))
|
||||
matrix = coeffs
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix (:, :, :)= coeffs
|
||||
|
||||
! Get scattxs value
|
||||
allocate(this % scattxs(groups))
|
||||
! Get this by summing the un-normalized P0 coefficient in matrix
|
||||
! over all outgoing groups
|
||||
this % scattxs = sum(matrix(1,:,:),dim=1)
|
||||
this % scattxs(:) = sum(matrix(1, :, :), dim=1)
|
||||
|
||||
allocate(energy(groups,groups))
|
||||
energy = ZERO
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = matrix(1,gout,gin)
|
||||
energy(gout,gin) = norm
|
||||
norm = matrix(1, gout, gin)
|
||||
energy(gout, gin) = norm
|
||||
if (norm /= ZERO) then
|
||||
matrix(:,gout,gin) = matrix(:,gout,gin) / norm
|
||||
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
|
@ -209,16 +212,16 @@ contains
|
|||
! Set dist values from matrix and initialize max_val
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
this % dist(gin) % data(:,gout) = matrix(:,gout,gin)
|
||||
this % dist(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
end do
|
||||
allocate(this % max_val(gin) % data(this % gmin(gin):this % gmax(gin)))
|
||||
this % max_val(gin) % data = ZERO
|
||||
this % max_val(gin) % data(:) = ZERO
|
||||
end do
|
||||
|
||||
! Step through the polynomial with fixed number of points to identify
|
||||
! the maximal value.
|
||||
Nmu = 1001
|
||||
dmu = TWO / real(Nmu - 1,8)
|
||||
dmu = TWO / real(Nmu - 1, 8)
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
do imu = 1, Nmu
|
||||
|
|
@ -228,7 +231,7 @@ contains
|
|||
else if (imu == Nmu) then
|
||||
mu = ONE
|
||||
else
|
||||
mu = -ONE + real(imu - 1,8) * dmu
|
||||
mu = -ONE + real(imu - 1, 8) * dmu
|
||||
end if
|
||||
! Calculate probability
|
||||
f = this % calc_f(gin,gout,mu)
|
||||
|
|
@ -245,37 +248,37 @@ contains
|
|||
|
||||
subroutine scattdatahistogram_init(this, mult, coeffs)
|
||||
class(ScattDataHistogram), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
|
||||
integer :: imu, gin, gout, groups, order
|
||||
real(8) :: norm
|
||||
real(8), allocatable :: energy(:,:)
|
||||
real(8), allocatable :: matrix(:,:,:)
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
|
||||
groups = size(coeffs,dim=3)
|
||||
order = size(coeffs,dim=1)
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order,groups,groups))
|
||||
matrix = coeffs
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix(:, :, :) = coeffs
|
||||
|
||||
! Get scattxs value
|
||||
allocate(this % scattxs(groups))
|
||||
! Get this by summing the un-normalized P0 coefficient in matrix
|
||||
! over all outgoing groups
|
||||
this % scattxs = sum(sum(matrix(:,:,:),dim=1),dim=1)
|
||||
this % scattxs(:) = sum(sum(matrix(:, :, :), dim=1), dim=1)
|
||||
|
||||
allocate(energy(groups,groups))
|
||||
energy = ZERO
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = sum(matrix(:,gout,gin))
|
||||
energy(gout,gin) = norm
|
||||
norm = sum(matrix(:, gout, gin))
|
||||
energy(gout, gin) = norm
|
||||
if (norm /= ZERO) then
|
||||
matrix(:,gout,gin) = matrix(:,gout,gin) / norm
|
||||
matrix(:, gout, gin) = matrix(:, gout, gin) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
|
@ -283,10 +286,10 @@ contains
|
|||
call scattdata_init(this, order, energy, mult)
|
||||
|
||||
allocate(this % mu(order))
|
||||
this % dmu = TWO / real(order,8)
|
||||
this % dmu = TWO / real(order, 8)
|
||||
this % mu(1) = -ONE
|
||||
do imu = 2, order
|
||||
this % mu(imu) = -ONE + real(imu - 1,8) * this % dmu
|
||||
this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
|
||||
end do
|
||||
|
||||
! Integrate this histogram so we can avoid rejection sampling while
|
||||
|
|
@ -297,21 +300,23 @@ contains
|
|||
this % gmin(gin):this % gmax(gin)))
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
! Store the histogram
|
||||
this % fmu(gin) % data(:,gout) = matrix(:,gout,gin)
|
||||
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
! Integrate the histogram
|
||||
this % dist(gin) % data(1,gout) = this % dmu * matrix(1,gout,gin)
|
||||
this % dist(gin) % data(1, gout) = &
|
||||
this % dmu * matrix(1, gout, gin)
|
||||
do imu = 2, order
|
||||
this % dist(gin) % data(imu,gout) = this % dmu * matrix(imu,gout,gin) + &
|
||||
this % dist(gin) % data(imu - 1,gout)
|
||||
this % dist(gin) % data(imu, gout) = &
|
||||
this % dmu * matrix(imu, gout, gin) + &
|
||||
this % dist(gin) % data(imu - 1, gout)
|
||||
end do
|
||||
|
||||
! Now make sure integral norms to zero
|
||||
norm = this % dist(gin) % data(order,gout)
|
||||
norm = this % dist(gin) % data(order, gout)
|
||||
if (norm > ZERO) then
|
||||
this % fmu(gin) % data(:,gout) = &
|
||||
this % fmu(gin) % data(:,gout) / norm
|
||||
this % dist(gin) % data(:,gout) = &
|
||||
this % dist(gin) % data(:,gout) / norm
|
||||
this % fmu(gin) % data(:, gout) = &
|
||||
this % fmu(gin) % data(:, gout) / norm
|
||||
this % dist(gin) % data(:, gout) = &
|
||||
this % dist(gin) % data(:, gout) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
|
@ -320,27 +325,27 @@ contains
|
|||
|
||||
subroutine scattdatatabular_init(this, mult, coeffs)
|
||||
class(ScattDataTabular), intent(inout) :: this ! Object to work on
|
||||
real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use
|
||||
real(8), intent(in) :: mult(:, :) ! Scatter Prod'n Matrix
|
||||
real(8), intent(in) :: coeffs(:, :, :) ! Coefficients to use
|
||||
|
||||
integer :: imu, gin, gout, groups, order
|
||||
real(8) :: norm
|
||||
real(8), allocatable :: energy(:,:)
|
||||
real(8), allocatable :: matrix(:,:,:)
|
||||
real(8), allocatable :: energy(:, :)
|
||||
real(8), allocatable :: matrix(:, :, :)
|
||||
|
||||
groups = size(coeffs,dim=3)
|
||||
order = size(coeffs,dim=1)
|
||||
groups = size(coeffs, dim=3)
|
||||
order = size(coeffs, dim=1)
|
||||
|
||||
! make a copy of coeffs that we can use to extract data and normalize
|
||||
allocate(matrix(order,groups,groups))
|
||||
matrix = coeffs
|
||||
allocate(matrix(order, groups, groups))
|
||||
matrix(:, :, :) = coeffs
|
||||
|
||||
! Build the angular distribution mu values
|
||||
allocate(this % mu(order))
|
||||
this % dmu = TWO / real(order - 1,8)
|
||||
this % dmu = TWO / real(order - 1, 8)
|
||||
this % mu(1) = -ONE
|
||||
do imu = 2, order - 1
|
||||
this % mu(imu) = -ONE + real(imu - 1,8) * this % dmu
|
||||
this % mu(imu) = -ONE + real(imu - 1, 8) * this % dmu
|
||||
end do
|
||||
this % mu(order) = ONE
|
||||
|
||||
|
|
@ -353,24 +358,24 @@ contains
|
|||
norm = ZERO
|
||||
do gout = 1, groups
|
||||
do imu = 2, order
|
||||
norm = norm + HALF * this % dmu * (matrix(imu - 1,gout,gin) + &
|
||||
matrix(imu,gout,gin))
|
||||
norm = norm + HALF * this % dmu * (matrix(imu - 1, gout, gin) + &
|
||||
matrix(imu, gout, gin))
|
||||
end do
|
||||
end do
|
||||
this % scattxs(gin) = norm
|
||||
end do
|
||||
|
||||
allocate(energy(groups,groups))
|
||||
energy = ZERO
|
||||
allocate(energy(groups, groups))
|
||||
energy(:, :) = ZERO
|
||||
! Build energy transfer probability matrix from data in matrix
|
||||
do gin = 1, groups
|
||||
do gout = 1, groups
|
||||
norm = ZERO
|
||||
do imu = 2, order
|
||||
norm = norm + HALF * this % dmu * &
|
||||
(matrix(imu - 1,gout,gin) + matrix(imu,gout,gin))
|
||||
(matrix(imu - 1, gout, gin) + matrix(imu, gout, gin))
|
||||
end do
|
||||
energy(gout,gin) = norm
|
||||
energy(gout, gin) = norm
|
||||
end do
|
||||
end do
|
||||
call scattdata_init(this, order, energy, mult)
|
||||
|
|
@ -383,12 +388,12 @@ contains
|
|||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
! Coeffs contain f(mu), put in f(mu) as that is where the
|
||||
! PDF lives
|
||||
this % fmu(gin) % data(:,gout) = matrix(:,gout,gin)
|
||||
this % fmu(gin) % data(:, gout) = matrix(:, gout, gin)
|
||||
|
||||
! Force positivity
|
||||
do imu = 1, order
|
||||
if (this % fmu(gin) % data(imu,gout) < ZERO) then
|
||||
this % fmu(gin) % data(imu,gout) = ZERO
|
||||
if (this % fmu(gin) % data(imu, gout) < ZERO) then
|
||||
this % fmu(gin) % data(imu, gout) = ZERO
|
||||
end if
|
||||
end do
|
||||
|
||||
|
|
@ -397,27 +402,27 @@ contains
|
|||
norm = ZERO
|
||||
do imu = 2, order
|
||||
norm = norm + HALF * this % dmu * &
|
||||
(this % fmu(gin) % data(imu - 1,gout) + &
|
||||
this % fmu(gin) % data(imu,gout))
|
||||
(this % fmu(gin) % data(imu - 1, gout) + &
|
||||
this % fmu(gin) % data(imu, gout))
|
||||
end do
|
||||
if (norm > ZERO) then
|
||||
this % fmu(gin) % data(:,gout) = &
|
||||
this % fmu(gin) % data(:,gout) / norm
|
||||
this % fmu(gin) % data(:, gout) = &
|
||||
this % fmu(gin) % data(:, gout) / norm
|
||||
end if
|
||||
|
||||
! Now create CDF from fmu with trapezoidal rule
|
||||
this % dist(gin) % data(1,gout) = ZERO
|
||||
this % dist(gin) % data(1, gout) = ZERO
|
||||
do imu = 2, order
|
||||
this % dist(gin) % data(imu,gout) = &
|
||||
this % dist(gin) % data(imu - 1,gout) + &
|
||||
HALF * this % dmu * (this % fmu(gin) % data(imu - 1,gout) + &
|
||||
this % fmu(gin) % data(imu,gout))
|
||||
this % dist(gin) % data(imu, gout) = &
|
||||
this % dist(gin) % data(imu - 1, gout) + &
|
||||
HALF * this % dmu * (this % fmu(gin) % data(imu - 1, gout) + &
|
||||
this % fmu(gin) % data(imu, gout))
|
||||
end do
|
||||
! Ensure we normalize to 1 still
|
||||
norm = this % dist(gin) % data(order,gout)
|
||||
norm = this % dist(gin) % data(order, gout)
|
||||
if (norm > ZERO) then
|
||||
this % dist(gin) % data(:,gout) = &
|
||||
this % dist(gin) % data(:,gout) / norm
|
||||
this % dist(gin) % data(:, gout) = &
|
||||
this % dist(gin) % data(:, gout) / norm
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
|
|
@ -438,7 +443,7 @@ contains
|
|||
if (gout < this % gmin(gin) .or. gout > this % gmax(gin)) then
|
||||
f = ZERO
|
||||
else
|
||||
f = evaluate_legendre(this % dist(gin) % data(:,gout),mu)
|
||||
f = evaluate_legendre(this % dist(gin) % data(:, gout), mu)
|
||||
end if
|
||||
|
||||
end function scattdatalegendre_calc_f
|
||||
|
|
@ -457,12 +462,12 @@ contains
|
|||
else
|
||||
! Find mu bin
|
||||
if (mu == ONE) then
|
||||
imu = size(this % fmu(gin) % data,dim=1)
|
||||
imu = size(this % fmu(gin) % data, dim=1)
|
||||
else
|
||||
imu = floor((mu + ONE)/ this % dmu + ONE)
|
||||
imu = floor((mu + ONE) / this % dmu + ONE)
|
||||
end if
|
||||
|
||||
f = this % fmu(gin) % data(imu,gout)
|
||||
f = this % fmu(gin) % data(imu, gout)
|
||||
end if
|
||||
|
||||
end function scattdatahistogram_calc_f
|
||||
|
|
@ -482,15 +487,15 @@ contains
|
|||
else
|
||||
! Find mu bin
|
||||
if (mu == ONE) then
|
||||
imu = size(this % fmu(gin) % data,dim=1) - 1
|
||||
imu = size(this % fmu(gin) % data, dim=1) - 1
|
||||
else
|
||||
imu = floor((mu + ONE)/ this % dmu + ONE)
|
||||
imu = floor((mu + ONE) / this % dmu + ONE)
|
||||
end if
|
||||
|
||||
! Now interpolate to find f(mu)
|
||||
r = (mu - this % mu(imu)) / (this % mu(imu + 1) - this % mu(imu))
|
||||
f = (ONE - r) * this % fmu(gin) % data(imu,gout) + &
|
||||
r * this % fmu(gin) % data(imu + 1,gout)
|
||||
f = (ONE - r) * this % fmu(gin) % data(imu, gout) + &
|
||||
r * this % fmu(gin) % data(imu + 1, gout)
|
||||
end if
|
||||
|
||||
end function scattdatatabular_calc_f
|
||||
|
|
@ -529,7 +534,7 @@ contains
|
|||
samples = 0
|
||||
do
|
||||
mu = TWO * prn() - ONE
|
||||
f = this % calc_f(gin,gout,mu)
|
||||
f = this % calc_f(gin, gout, mu)
|
||||
if (f > ZERO) then
|
||||
u = prn() * M
|
||||
if (u <= f) then
|
||||
|
|
@ -567,11 +572,11 @@ contains
|
|||
end do
|
||||
|
||||
xi = prn()
|
||||
if (xi < this % dist(gin) % data(1,gout)) then
|
||||
if (xi < this % dist(gin) % data(1, gout)) then
|
||||
imu = 1
|
||||
else
|
||||
imu = binary_search(this % dist(gin) % data(:,gout), &
|
||||
size(this % dist(gin) % data(:,gout)), xi)
|
||||
imu = binary_search(this % dist(gin) % data(:, gout), &
|
||||
size(this % dist(gin) % data(:, gout)), xi)
|
||||
end if
|
||||
|
||||
! Randomly select a mu in this bin.
|
||||
|
|
@ -604,12 +609,12 @@ contains
|
|||
end do
|
||||
|
||||
! determine outgoing cosine bin
|
||||
NP = size(this % dist(gin) % data(:,gout))
|
||||
NP = size(this % dist(gin) % data(:, gout))
|
||||
xi = prn()
|
||||
|
||||
c_k = this % dist(gin) % data(1,gout)
|
||||
c_k = this % dist(gin) % data(1, gout)
|
||||
do k = 1, NP - 1
|
||||
c_k1 = this % dist(gin) % data(k + 1,gout)
|
||||
c_k1 = this % dist(gin) % data(k + 1, gout)
|
||||
if (xi < c_k1) exit
|
||||
c_k = c_k1
|
||||
end do
|
||||
|
|
@ -617,18 +622,19 @@ contains
|
|||
! check to make sure k is <= NP - 1
|
||||
k = min(k, NP - 1)
|
||||
|
||||
p0 = this % fmu(gin) % data(k,gout)
|
||||
p0 = this % fmu(gin) % data(k, gout)
|
||||
mu0 = this % mu(k)
|
||||
! Linear-linear interpolation to find mu value w/in bin.
|
||||
p1 = this % fmu(gin) % data(k + 1,gout)
|
||||
p1 = this % fmu(gin) % data(k + 1, gout)
|
||||
mu1 = this % mu(k + 1)
|
||||
|
||||
frac = (p1 - p0)/(mu1 - mu0)
|
||||
frac = (p1 - p0) / (mu1 - mu0)
|
||||
|
||||
if (frac == ZERO) then
|
||||
mu = mu0 + (xi - c_k)/p0
|
||||
mu = mu0 + (xi - c_k) / p0
|
||||
else
|
||||
mu = mu0 + (sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac
|
||||
mu = mu0 + &
|
||||
(sqrt(max(ZERO, p0 * p0 + TWO * frac * (xi - c_k))) - p0) / frac
|
||||
end if
|
||||
|
||||
if (mu <= -ONE) then
|
||||
|
|
@ -649,23 +655,24 @@ contains
|
|||
pure function scattdata_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattData), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:,:,:) ! Resultant matrix just built
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
groups = size(this % energy)
|
||||
order = min(req_order,size(this % dist(1) % data(:,1)))
|
||||
! Set gin and gout for getting the order
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order,groups,groups))
|
||||
allocate(matrix(order, groups, groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix = ZERO
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:,gout,gin) = this % scattxs(gin) * &
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % dist(gin) % data(1:order,gout)
|
||||
this % dist(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdata_get_matrix
|
||||
|
|
@ -673,23 +680,23 @@ contains
|
|||
pure function scattdatahistogram_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattDataHistogram), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:,:,:) ! Resultant matrix just built
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
groups = size(this % energy)
|
||||
order = min(req_order,size(this % dist(1) % data(:,1)))
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order,groups,groups))
|
||||
allocate(matrix(order, groups, groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix = ZERO
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:,gout,gin) = this % scattxs(gin) * &
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % fmu(gin) % data(1:order,gout)
|
||||
this % fmu(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdatahistogram_get_matrix
|
||||
|
|
@ -697,23 +704,23 @@ contains
|
|||
pure function scattdatatabular_get_matrix(this, req_order) result(matrix)
|
||||
class(ScattDataTabular), intent(in) :: this ! Scattering Object to work with
|
||||
integer, intent(in) :: req_order ! Requested order of matrix
|
||||
real(8), allocatable :: matrix(:,:,:) ! Resultant matrix just built
|
||||
real(8), allocatable :: matrix(:, :, :) ! Resultant matrix just built
|
||||
|
||||
integer :: order, groups, gin, gout
|
||||
|
||||
groups = size(this % energy)
|
||||
order = min(req_order,size(this % dist(1) % data(:,1)))
|
||||
order = min(req_order, size(this % dist(1) % data, dim=1))
|
||||
|
||||
allocate(matrix(order,groups,groups))
|
||||
allocate(matrix(order, groups, groups))
|
||||
! Initialize to 0; this way the zero entries in the dense matrix dont
|
||||
! need to be explicitly set, requiring a significant increase in the
|
||||
! lines of code.
|
||||
matrix = ZERO
|
||||
matrix(:, :, :) = ZERO
|
||||
do gin = 1, groups
|
||||
do gout = this % gmin(gin), this % gmax(gin)
|
||||
matrix(:,gout,gin) = this % scattxs(gin) * &
|
||||
matrix(:, gout, gin) = this % scattxs(gin) * &
|
||||
this % energy(gin) % data(gout) * &
|
||||
this % fmu(gin) % data(1:order,gout)
|
||||
this % fmu(gin) % data(1:order, gout)
|
||||
end do
|
||||
end do
|
||||
end function scattdatatabular_get_matrix
|
||||
|
|
|
|||
|
|
@ -961,9 +961,9 @@ contains
|
|||
if (i_nuclide > 0) then
|
||||
score = score * atom_density * &
|
||||
nucxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
|
||||
UVW=p_uvw,MU=p % mu) / &
|
||||
UVW=p_uvw, MU=p % mu) / &
|
||||
matxs % get_xs('scatter*f_mu/mult', p % last_g, p % g, &
|
||||
UVW=p_uvw,MU=p % mu)
|
||||
UVW=p_uvw, MU=p % mu)
|
||||
end if
|
||||
|
||||
else
|
||||
|
|
@ -1080,11 +1080,11 @@ contains
|
|||
end if
|
||||
if (i_nuclide > 0) then
|
||||
score = score * atom_density * &
|
||||
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
|
||||
nucxs % get_xs('fission', p_g, UVW=p_uvw) / &
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
else
|
||||
score = score * &
|
||||
matxs % get_xs('fission', p_g, UVW=p_uvw) / &
|
||||
matxs % get_xs('fission', p_g, UVW=p_uvw) / &
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
end if
|
||||
else
|
||||
|
|
@ -1170,11 +1170,11 @@ contains
|
|||
if (i_nuclide > 0) then
|
||||
score = score * atom_density * &
|
||||
nucxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / &
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
else
|
||||
score = score * &
|
||||
matxs % get_xs('kappa_fission', p_g, UVW=p_uvw) / &
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
matxs % get_xs('absorption', p_g, UVW=p_uvw)
|
||||
end if
|
||||
else
|
||||
if (i_nuclide > 0) then
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue