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

Convert scattering routines to C++

Browse source
This commit is contained in:
Paul Romano 2018-11-21 11:13:55 -06:00
parent 9105cd195b
commit e8af7a2d16
14 changed files with 620 additions and 1100 deletions
Download patch file Download diff file Expand all files Collapse all files

11
include/openmc/endf.h
View file

@ -22,6 +22,17 @@ Interpolation int2interp(int i);
//! \return Whether corresponding reaction is a fission reaction
bool is_fission(int MT);
//! Determine if a given MT number is that of a disappearance reaction, i.e., a
//! reaction with no neutron in the exit channel
//! \param[in] MT ENDF MT value
//! \return Whether corresponding reaction is a disappearance reaction
bool is_disappearance(int MT);
//! Determine if a given MT number is that of an inelastic scattering reaction
//! \param[in] MT ENDF MT value
//! \return Whether corresponding reaction is an inelastic scattering reaction
bool is_inelastic_scatter(int MT);
//==============================================================================
//! Abstract one-dimensional function
//==============================================================================

6
include/openmc/material.h
View file

@ -43,5 +43,11 @@ public:
explicit Material(pugi::xml_node material_node);
};
//==============================================================================
// Fortran compatibility
//==============================================================================
extern "C" bool material_isotropic(int i_material, int i_nuc_mat);
} // namespace openmc
#endif // OPENMC_MATERIAL_H

26
include/openmc/nuclide.h
View file

@ -17,19 +17,35 @@
namespace openmc {
//==============================================================================
// Constants
//==============================================================================
constexpr double CACHE_INVALID {-1.0};
//===============================================================================
// Data for a nuclide
//===============================================================================
class Nuclide {
public:
// Types, aliases
using EmissionMode = ReactionProduct::EmissionMode;
struct EnergyGrid {
std::vector<int> grid_index;
std::vector<double> energy;
};
// Constructors
Nuclide(hid_t group, const double* temperature, int n);
// Methods
double nu(double E, EmissionMode mode, int group=0);
double nu(double E, EmissionMode mode, int group=0) const;
void calculate_elastic_xs(int i_nuclide) const;
//! Determines the microscopic 0K elastic cross section at a trial relative
//! energy used in resonance scattering
double elastic_xs_0K(double E) const;
// Data members
std::string name_; //! Name of nuclide, e.g. "U235"
@ -38,6 +54,7 @@ public:
int metastable_; //! Metastable state
double awr_; //! Atomic weight ratio
std::vector<double> kTs_; //! temperatures in eV (k*T)
std::vector<EnergyGrid> grid_; //! Energy grid at each temperature
bool fissionable_ {false}; //! Whether nuclide is fissionable
bool has_partial_fission_ {false}; //! has partial fission reactions?
@ -45,7 +62,14 @@ public:
int n_precursor_ {0}; //! Number of delayed neutron precursors
std::unique_ptr<Function1D> total_nu_; //! Total neutron yield
// Resonance scattering information
bool resonant_ {false};
std::vector<double> energy_0K_;
std::vector<double> elastic_0K_;
std::vector<double> xs_cdf_;
std::vector<std::unique_ptr<Reaction>> reactions_; //! Reactions
std::vector<int> index_inelastic_scatter_;
private:
void create_derived();

30
include/openmc/physics.h
View file

@ -2,6 +2,7 @@
#define OPENMC_PHYSICS_H
#include "openmc/bank.h"
#include "openmc/nuclide.h"
#include "openmc/particle.h"
#include "openmc/position.h"
#include "openmc/reaction.h"
@ -51,22 +52,33 @@ void sample_photon_product(int i_nuclide, double E, int* i_rx, int* i_product);
void absorption(Particle* p, int i_nuclide);
extern "C" void scatter(Particle*, int i_nuclide, int i_nuc_mat);
void scatter(Particle*, int i_nuclide, int i_nuc_mat);
// void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double* E,
// Direction* u, double* mu_lab, double* wgt);
//! Treats the elastic scattering of a neutron with a target.
void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
double* uvw, double* mu_lab, double* wgt);
// void sab_scatter(int i_nuclide, int i_sab, double* E, Direction* u, double* mu);
extern "C" void sab_scatter(int i_nuclide, int i_sab, double* E,
double* uvw, double* mu);
// void sample_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u,
// Direction v_neut, double* wgt, double xs_eff, double kT);
//! samples the target velocity. The constant cross section free gas model is
//! the default method. Methods for correctly accounting for the energy
//! dependence of cross sections in treating resonance elastic scattering such
//! as the DBRC, WCM, and a new, accelerated scheme are also implemented here.
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_neut, double xs_eff, double kT, double* wgt);
// void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u,
// double kT);
//! samples a target velocity based on the free gas scattering formulation, used
//! by most Monte Carlo codes, in which cross section is assumed to be constant
//! in energy. Excellent documentation for this method can be found in
//! FRA-TM-123.
Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT);
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site);
// void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p);
//! handles all reactions with a single secondary neutron (other than fission),
//! i.e. level scattering, (n,np), (n,na), etc.
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p);
void sample_secondary_photons(Particle* p, int i_nuclide);

6
include/openmc/position.h
View file

@ -25,6 +25,8 @@ struct Position {
Position& operator-=(double);
Position& operator*=(Position);
Position& operator*=(double);
Position& operator/=(Position);
Position& operator/=(double);
const double& operator[](int i) const {
switch (i) {
@ -76,6 +78,10 @@ inline Position operator*(Position a, Position b) { return a *= b; }
inline Position operator*(Position a, double b) { return a *= b; }
inline Position operator*(double a, Position b) { return b *= a; }
inline Position operator/(Position a, Position b) { return a /= b; }
inline Position operator/(Position a, double b) { return a /= b; }
inline Position operator/(double a, Position b) { return b /= a; }
inline bool operator==(Position a, Position b)
{return a.x == b.x && a.y == b.y && a.z == b.z;}

29
src/endf.cpp
View file

@ -45,6 +45,35 @@ bool is_fission(int mt)
return mt == 18 || mt == 19 || mt == 20 || mt == 21 || mt == 38;
}
bool is_disappearance(int mt)
{
if (mt >= N_DISAPPEAR && mt <= N_DA) {
return true;
} else if (mt >= N_P0 && mt <= N_AC) {
return true;
} else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE
|| mt == N_3HEA || mt == N_3P) {
return true;
} else {
return false;
}
}
bool is_inelastic_scatter(int mt)
{
if (mt < 100) {
if (is_fission(mt)) {
return false;
} else {
return mt >= MISC && mt != 27;
}
} else if (mt <= 200) {
return !is_disappearance(mt);
} else {
return false;
}
}
//==============================================================================
// Polynomial implementation
//==============================================================================

4
src/input_xml.F90
View file

@ -2113,9 +2113,6 @@ contains
call close_group(group_id)
call file_close(file_id)
! Assign resonant scattering data
if (res_scat_on) call nuclides(i_nuclide) % assign_0K_elastic_scattering()
! Determine if minimum/maximum energy for this nuclide is greater/less
! than the previous
if (size(nuclides(i_nuclide) % grid) >= 1) then
@ -2297,7 +2294,6 @@ contains
logical :: file_exists ! Does multipole library exist?
character(7) :: readable ! Is multipole library readable?
character(MAX_FILE_LEN) :: filename ! Path to multipole xs library
character(kind=C_CHAR), pointer :: string(:)
integer(HID_T) :: file_id
integer(HID_T) :: group_id

17
src/material_header.F90
View file

@ -1041,4 +1041,21 @@ contains
end subroutine bremsstrahlung_init
!===============================================================================
! Fortran compatibility
!===============================================================================
function material_isotropic(i_material, i_nuc_mat) result(iso) bind(C)
integer(C_INT), value :: i_material
integer(C_INT), value :: i_nuc_mat
logical(C_BOOL) :: iso
iso = .false.
associate (mat => materials(i_material))
if (mat % has_isotropic_nuclides) then
iso = mat % p0(i_nuc_mat)
end if
end associate
end function
end module material_header

111
src/nuclide.cpp
View file

@ -5,6 +5,7 @@
#include "openmc/error.h"
#include "openmc/hdf5_interface.h"
#include "openmc/message_passing.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/string_utils.h"
@ -143,22 +144,48 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n)
// Sort temperatures to read
std::sort(temps_to_read.begin(), temps_to_read.end());
// Determine exact kT values
hid_t energy_group = open_group(group, "energy");
for (const auto& T : temps_to_read) {
std::string dset {std::to_string(T) + "K"};
// Determine exact kT values
double kT;
read_dataset(kT_group, dset.c_str(), kT);
kTs_.push_back(kT);
// Read energy grid
grid_.emplace_back();
read_dataset(energy_group, dset.c_str(), grid_.back().energy);
}
close_group(kT_group);
// Check for 0K energy grid
if (object_exists(energy_group, "0K")) {
read_dataset(energy_group, "0K", energy_0K_);
}
close_group(energy_group);
// Read reactions
hid_t rxs_group = open_group(group, "reactions");
for (auto name : group_names(rxs_group)) {
if (starts_with(name, "reaction_")) {
hid_t rx_group = open_group(rxs_group, name.c_str());
reactions_.push_back(std::make_unique<Reaction>(rx_group, temps_to_read));
// Check for 0K elastic scattering
if (reactions_.back()->mt_ == ELASTIC) {
if (object_exists(rx_group, "0K")) {
hid_t temp_group = open_group(rx_group, "0K");
read_dataset(temp_group, "xs", elastic_0K_);
close_group(temp_group);
}
}
close_group(rx_group);
// Determine reaction indices for inelastic scattering reactions
if (is_inelastic_scatter(reactions_.back()->mt_)) {
index_inelastic_scatter_.push_back(reactions_.size() - 1);
}
}
}
close_group(rxs_group);
@ -211,9 +238,52 @@ void Nuclide::create_derived()
}
}
}
if (settings::res_scat_on) {
// Determine if this nuclide should be treated as a resonant scatterer
if (!settings::res_scat_nuclides.empty()) {
// If resonant nuclides were specified, check the list explicitly
for (const auto& name : settings::res_scat_nuclides) {
if (name_ == name) {
resonant_ = true;
// Make sure nuclide has 0K data
if (energy_0K_.empty()) {
fatal_error("Cannot treat " + name_ + " as a resonant scatterer "
"because 0 K elastic scattering data is not present.");
}
break;
}
}
} else {
// Otherwise, assume that any that have 0 K elastic scattering data are
// resonant
resonant_ = !energy_0K_.empty();
}
if (resonant_) {
// Build CDF for 0K elastic scattering
double xs_cdf_sum = 0.0;
xs_cdf_.resize(energy_0K_.size());
xs_cdf_[0] = 0.0;
const auto& E = energy_0K_;
auto& xs = elastic_0K_;
for (int i = 0; i < E.size() - 1; ++i) {
// Negative cross sections result in a CDF that is not monotonically
// increasing. Set all negative xs values to zero.
if (xs[i] < 0.0) xs[i] = 0.0;
// build xs cdf
xs_cdf_sum += (std::sqrt(E[i])*xs[i] + std::sqrt(E[i+1])*xs[i+1])
/ 2.0 * (E[i+1] - E[i]);
xs_cdf_[i] = xs_cdf_sum;
}
}
}
}
double Nuclide::nu(double E, EmissionMode mode, int group)
double Nuclide::nu(double E, EmissionMode mode, int group) const
{
if (!fissionable_) return 0.0;
@ -253,6 +323,43 @@ double Nuclide::nu(double E, EmissionMode mode, int group)
}
}
void Nuclide::calculate_elastic_xs(int i_nuclide) const
{
// Get temperature index, grid index, and interpolation factor
auto& micro = simulation::micro_xs[i_nuclide-1];
int i_temp = micro.index_temp - 1;
int i_grid = micro.index_grid - 1;
double f = micro.interp_factor;
if (i_temp >= 0) {
const auto& xs = reactions_[0]->xs_[i_temp].value;
micro.elastic = (1.0 - f)*xs[i_grid] + f*xs[i_grid + 1];
}
}
double Nuclide::elastic_xs_0K(double E) const
{
// Determine index on nuclide energy grid
int i_grid;
if (E < energy_0K_.front()) {
i_grid = 0;
} else if (E > energy_0K_.back()) {
i_grid = energy_0K_.size() - 2;
} else {
i_grid = lower_bound_index(energy_0K_.begin(), energy_0K_.end(), E);
}
// check for rare case where two energy points are the same
if (energy_0K_[i_grid] == energy_0K_[i_grid+1]) ++i_grid;
// calculate interpolation factor
double f = (E - energy_0K_[i_grid]) /
(energy_0K_[i_grid + 1] - energy_0K_[i_grid]);
// Calculate microscopic nuclide elastic cross section
return (1.0 - f)*elastic_0K_[i_grid] + f*elastic_0K_[i_grid + 1];
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================

162
src/nuclide_header.F90
View file

@ -71,15 +71,7 @@ module nuclide_header
! Microscopic cross sections
type(SumXS), allocatable :: xs(:)
! Resonance scattering info
logical :: resonant = .false. ! resonant scatterer?
real(8), allocatable :: energy_0K(:) ! energy grid for 0K xs
real(8), allocatable :: elastic_0K(:) ! Microscopic elastic cross section
real(8), allocatable :: xs_cdf(:) ! CDF of v_rel times cross section
! Fission information
logical :: has_partial_fission = .false. ! nuclide has partial fission reactions?
integer :: n_fission = 0 ! # of fission reactions
integer :: n_precursor = 0 ! # of delayed neutron precursors
integer, allocatable :: index_fission(:) ! indices in reactions
class(Function1D), allocatable :: total_nu
@ -95,7 +87,6 @@ module nuclide_header
! Reactions
type(Reaction), allocatable :: reactions(:)
integer, allocatable :: index_inelastic_scatter(:)
! Array that maps MT values to index in reactions; used at tally-time. Note
! that ENDF-102 does not assign any MT values above 891.
@ -108,7 +99,6 @@ module nuclide_header
type(C_PTR) :: ptr
contains
procedure :: assign_0K_elastic_scattering
procedure :: clear => nuclide_clear
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
@ -238,79 +228,6 @@ contains
ptr = C_LOC(micro_xs(1))
end function
!===============================================================================
! ASSIGN_0K_ELASTIC_SCATTERING
!===============================================================================
subroutine assign_0K_elastic_scattering(this)
class(Nuclide), intent(inout) :: this
integer :: i
real(8) :: xs_cdf_sum
interface
function res_scat_nuclides_empty() result(empty) bind(C)
import C_BOOL
logical(C_BOOL) :: empty
end function
function res_scat_nuclides_size() result(n) bind(C)
import C_INT
integer(C_INT) :: n
end function
function res_scat_nuclides_cmp(i, name) result(b) bind(C)
import C_INT, C_CHAR, C_BOOL
integer(C_INT), value :: i
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: b
end function
end interface
this % resonant = .false.
if (.not. res_scat_nuclides_empty()) then
! If resonant nuclides were specified, check the list explicitly
do i = 1, res_scat_nuclides_size()
if (res_scat_nuclides_cmp(i, to_c_string(this % name))) then
this % resonant = .true.
! Make sure nuclide has 0K data
if (.not. allocated(this % energy_0K)) then
call fatal_error("Cannot treat " // trim(this % name) // " as a &
&resonant scatterer because 0 K elastic scattering data is &
&not present.")
end if
exit
end if
end do
else
! Otherwise, assume that any that have 0 K elastic scattering data are
! resonant
this % resonant = allocated(this % energy_0K)
end if
if (this % resonant) then
! Build CDF for 0K elastic scattering
xs_cdf_sum = ZERO
allocate(this % xs_cdf(0:size(this % energy_0K)))
this % xs_cdf(0) = ZERO
associate (E => this % energy_0K, xs => this % elastic_0K)
do i = 1, size(E) - 1
! Negative cross sections result in a CDF that is not monotonically
! increasing. Set all negative xs values to zero.
if (xs(i) < ZERO) xs(i) = ZERO
! build xs cdf
xs_cdf_sum = xs_cdf_sum + (sqrt(E(i))*xs(i) + sqrt(E(i+1))*xs(i+1))&
/ TWO * (E(i+1) - E(i))
this % xs_cdf(i) = xs_cdf_sum
end do
end associate
end if
end subroutine assign_0K_elastic_scattering
!===============================================================================
! NUCLIDE_CLEAR resets and deallocates data in Nuclide
!===============================================================================
@ -341,7 +258,6 @@ contains
integer(HID_T) :: kT_group
integer(HID_T) :: rxs_group
integer(HID_T) :: rx_group
integer(HID_T) :: xs, temp_group
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
@ -355,7 +271,6 @@ contains
real(8) :: temp_actual
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
type(VectorInt) :: index_inelastic_scatter
interface
function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C)
@ -494,15 +409,6 @@ contains
call read_dataset(this % grid(i) % energy, energy_dset)
call close_dataset(energy_dset)
end do
! Check for 0K energy grid
if (object_exists(energy_group, '0K')) then
energy_dset = open_dataset(energy_group, '0K')
call get_shape(energy_dset, dims)
allocate(this % energy_0K(int(dims(1), 4)))
call read_dataset(this % energy_0K, energy_dset)
call close_dataset(energy_dset)
end if
call close_group(energy_group)
! Get MT values based on group names
@ -523,34 +429,10 @@ contains
! Set pointer for each reaction
call this % reactions(i) % init(this % ptr, i)
! Check for 0K elastic scattering
if (this % reactions(i) % MT == 2) then
if (object_exists(rx_group, '0K')) then
temp_group = open_group(rx_group, '0K')
xs = open_dataset(temp_group, 'xs')
call get_shape(xs, dims)
allocate(this % elastic_0K(int(dims(1), 4)))
call read_dataset(this % elastic_0K, xs)
call close_dataset(xs)
call close_group(temp_group)
end if
end if
! Add the reaction index to the scattering array if this is an inelastic
! scatter reaction
if (is_inelastic_scatter(MTs % data(i))) then
call index_inelastic_scatter % push_back(i)
end if
call close_group(rx_group)
end do
call close_group(rxs_group)
! Recast to a regular array to save space
allocate(this % index_inelastic_scatter(index_inelastic_scatter % size()))
this % index_inelastic_scatter = &
index_inelastic_scatter % data(1: index_inelastic_scatter % size())
! Read unresolved resonance probability tables if present
if (object_exists(group_id, 'urr')) then
this % urr_present = .true.
@ -726,7 +608,6 @@ contains
allocate(this % index_fission(1))
elseif (rx % MT == N_F) then
allocate(this % index_fission(PARTIAL_FISSION_MAX))
this % has_partial_fission = .true.
end if
end if
@ -746,7 +627,6 @@ contains
if (t == 1) then
i_fission = i_fission + 1
this % index_fission(i_fission) = i
this % n_fission = this % n_fission + 1
end if
end if ! fission
end do ! temperature
@ -1361,45 +1241,6 @@ contains
end if
end subroutine multipole_deriv_eval
!===============================================================================
! 0K_ELASTIC_XS determines the microscopic 0K elastic cross section
! for a given nuclide at the trial relative energy used in resonance scattering
!===============================================================================
pure function elastic_xs_0K(E, nuc) result(xs_out)
real(8), intent(in) :: E ! trial energy
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature
real(8) :: xs_out ! 0K xs at trial energy
integer :: i_grid ! index on nuclide energy grid
integer :: n_grid
real(8) :: f ! interp factor on nuclide energy grid
! Determine index on nuclide energy grid
n_grid = size(nuc % energy_0K)
if (E < nuc % energy_0K(1)) then
i_grid = 1
elseif (E > nuc % energy_0K(n_grid)) then
i_grid = n_grid - 1
else
i_grid = binary_search(nuc % energy_0K, n_grid, E)
end if
! check for rare case where two energy points are the same
if (nuc % energy_0K(i_grid) == nuc % energy_0K(i_grid+1)) then
i_grid = i_grid + 1
end if
! calculate interpolation factor
f = (E - nuc % energy_0K(i_grid)) &
& / (nuc % energy_0K(i_grid + 1) - nuc % energy_0K(i_grid))
! Calculate microscopic nuclide elastic cross section
xs_out = (ONE - f) * nuc % elastic_0K(i_grid) &
& + f * nuc % elastic_0K(i_grid + 1)
end function elastic_xs_0K
!===============================================================================
! CALCULATE_URR_XS determines cross sections in the unresolved resonance range
! from probability tables
@ -1678,9 +1519,6 @@ contains
call nuclide_dict % set(to_lower(name_), n)
n_nuclides = n
! Assign resonant scattering data
if (res_scat_on) call nuclides(n) % assign_0K_elastic_scattering()
! Initialize nuclide grid
call nuclides(n) % init_grid(energy_min(NEUTRON), &
energy_max(NEUTRON), n_log_bins)

567
src/physics.F90
View file

@ -1,9 +1,6 @@
module physics
use algorithm, only: binary_search
use bank_header
use constants
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message
use material_header, only: Material, materials
use math
@ -15,12 +12,10 @@ module physics
atomic_relaxation, pair_production, &
thick_target_bremsstrahlung
use physics_common
use random_lcg, only: prn, advance_prn_seed, prn_set_stream
use reaction_header, only: Reaction
use random_lcg, only: prn
use sab_header, only: sab_tables
use settings
use simulation_header
use string, only: to_str
use tally_header
implicit none
@ -254,229 +249,17 @@ contains
end function sample_element
!===============================================================================
! SCATTER
!===============================================================================
subroutine scatter(p, i_nuclide, i_nuc_mat) bind(C)
type(Particle), intent(inout) :: p
integer(C_INT), value :: i_nuclide
integer(C_INT), value :: i_nuc_mat
integer :: i
integer :: j
integer :: i_temp
integer :: i_grid
integer :: threshold
real(8) :: f
real(8) :: prob
real(8) :: cutoff
real(8) :: uvw_new(3) ! outgoing uvw for iso-in-lab scattering
real(8) :: uvw_old(3) ! incoming uvw for iso-in-lab scattering
real(8) :: phi ! azimuthal angle for iso-in-lab scattering
real(8) :: kT ! temperature in eV
logical :: sampled ! whether or not a reaction type has been sampled
type(Nuclide), pointer :: nuc
! copy incoming direction
uvw_old(:) = p % coord(1) % uvw
! Get pointer to nuclide and grid index/interpolation factor
nuc => nuclides(i_nuclide)
i_temp = micro_xs(i_nuclide) % index_temp
i_grid = micro_xs(i_nuclide) % index_grid
f = micro_xs(i_nuclide) % interp_factor
! For tallying purposes, this routine might be called directly. In that
! case, we need to sample a reaction via the cutoff variable
cutoff = prn() * (micro_xs(i_nuclide) % total - &
micro_xs(i_nuclide) % absorption)
sampled = .false.
! Calculate elastic cross section if it wasn't precalculated
if (micro_xs(i_nuclide) % elastic == CACHE_INVALID) then
call nuc % calculate_elastic_xs(micro_xs(i_nuclide))
end if
prob = micro_xs(i_nuclide) % elastic - micro_xs(i_nuclide) % thermal
if (prob > cutoff) then
! =======================================================================
! NON-S(A,B) ELASTIC SCATTERING
! Determine temperature
if (nuc % mp_present) then
kT = p % sqrtkT**2
else
kT = nuc % kTs(micro_xs(i_nuclide) % index_temp)
end if
! Perform collision physics for elastic scattering
call elastic_scatter(i_nuclide, nuc % reactions(1), kT, p % E, &
p % coord(1) % uvw, p % mu, p % wgt)
p % event_MT = ELASTIC
sampled = .true.
end if
prob = micro_xs(i_nuclide) % elastic
if (prob > cutoff .and. .not. sampled) then
! =======================================================================
! S(A,B) SCATTERING
call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab, p % E, &
p % coord(1) % uvw, p % mu)
p % event_MT = ELASTIC
sampled = .true.
end if
if (.not. sampled) then
! =======================================================================
! INELASTIC SCATTERING
j = 0
do while (prob < cutoff)
j = j + 1
i = nuc % index_inelastic_scatter(j)
! Check to make sure inelastic scattering reaction sampled
if (i > size(nuc % reactions)) then
call particle_write_restart(p)
call fatal_error("Did not sample any reaction for nuclide " &
&// trim(nuc % name))
end if
associate (rx => nuc % reactions(i))
! if energy is below threshold for this reaction, skip it
threshold = rx % xs_threshold(i_temp)
if (i_grid < threshold) cycle
! add to cumulative probability
prob = prob + ((ONE - f)*rx % xs(i_temp, i_grid - threshold + 1) &
+ f*(rx % xs(i_temp, i_grid - threshold + 2)))
end associate
end do
! Perform collision physics for inelastic scattering
call inelastic_scatter(nuc, nuc%reactions(i), p)
p % event_MT = nuc % reactions(i) % MT
end if
! Set event component
p % event = EVENT_SCATTER
! Sample new outgoing angle for isotropic-in-lab scattering
associate (mat => materials(p % material))
if (mat % has_isotropic_nuclides) then
if (materials(p % material) % p0(i_nuc_mat)) then
! Sample isotropic-in-lab outgoing direction
uvw_new(1) = TWO * prn() - ONE
phi = TWO * PI * prn()
uvw_new(2) = cos(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
uvw_new(3) = sin(phi) * sqrt(ONE - uvw_new(1)*uvw_new(1))
p % mu = dot_product(uvw_old, uvw_new)
! Change direction of particle
p % coord(1) % uvw = uvw_new
end if
end if
end associate
end subroutine scatter
!===============================================================================
! ELASTIC_SCATTER treats the elastic scattering of a neutron with a
! target.
!===============================================================================
subroutine elastic_scatter(i_nuclide, rxn, kT, E, uvw, mu_lab, wgt)
integer, intent(in) :: i_nuclide
type(Reaction), intent(in) :: rxn
real(8), intent(in) :: kT ! temperature in eV
real(8), intent(inout) :: E
real(8), intent(inout) :: uvw(3)
real(8), intent(out) :: mu_lab
real(8), intent(inout) :: wgt
real(8) :: awr ! atomic weight ratio of target
real(8) :: mu_cm ! cosine of polar angle in center-of-mass
real(8) :: vel ! magnitude of velocity
real(8) :: v_n(3) ! velocity of neutron
real(8) :: v_cm(3) ! velocity of center-of-mass
real(8) :: v_t(3) ! velocity of target nucleus
real(8) :: uvw_cm(3) ! directional cosines in center-of-mass
type(Nuclide), pointer :: nuc
! get pointer to nuclide
nuc => nuclides(i_nuclide)
vel = sqrt(E)
awr = nuc % awr
! Neutron velocity in LAB
v_n = vel * uvw
! Sample velocity of target nucleus
if (.not. micro_xs(i_nuclide) % use_ptable) then
call sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
micro_xs(i_nuclide) % elastic, kT)
else
v_t = ZERO
end if
! Velocity of center-of-mass
v_cm = (v_n + awr*v_t)/(awr + ONE)
! Transform to CM frame
v_n = v_n - v_cm
! Find speed of neutron in CM
vel = sqrt(dot_product(v_n, v_n))
! Sample scattering angle
mu_cm = rxn % sample_elastic_mu(E)
! Determine direction cosines in CM
uvw_cm = v_n/vel
! Rotate neutron velocity vector to new angle -- note that the speed of the
! neutron in CM does not change in elastic scattering. However, the speed
! will change when we convert back to LAB
v_n = vel * rotate_angle(uvw_cm, mu_cm)
! Transform back to LAB frame
v_n = v_n + v_cm
E = dot_product(v_n, v_n)
vel = sqrt(E)
! compute cosine of scattering angle in LAB frame by taking dot product of
! neutron's pre- and post-collision angle
mu_lab = dot_product(uvw, v_n) / vel
! Set energy and direction of particle in LAB frame
uvw = v_n / vel
! Because of floating-point roundoff, it may be possible for mu_lab to be
! outside of the range [-1,1). In these cases, we just set mu_lab to exactly
! -1 or 1
if (abs(mu_lab) > ONE) mu_lab = sign(ONE,mu_lab)
end subroutine elastic_scatter
!===============================================================================
! SAB_SCATTER performs thermal scattering of a particle with a bound scatterer
! according to a specified S(a,b) table.
!===============================================================================
subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu)
integer, intent(in) :: i_nuclide ! index in micro_xs
integer, intent(in) :: i_sab ! index in sab_tables
real(8), intent(inout) :: E ! incoming/outgoing energy
real(8), intent(inout) :: uvw(3) ! directional cosines
real(8), intent(out) :: mu ! scattering cosine
subroutine sab_scatter(i_nuclide, i_sab, E, uvw, mu) bind(C)
integer(C_INT), value :: i_nuclide ! index in micro_xs
integer(C_INT), value :: i_sab ! index in sab_tables
real(C_DOUBLE), intent(inout) :: E ! incoming/outgoing energy
real(C_DOUBLE), intent(inout) :: uvw(3) ! directional cosines
real(C_DOUBLE), intent(out) :: mu ! scattering cosine
real(C_DOUBLE) :: E_out
type(C_PTR) :: ptr
@ -490,340 +273,4 @@ contains
uvw = rotate_angle(uvw, mu)
end subroutine sab_scatter
!===============================================================================
! SAMPLE_TARGET_VELOCITY samples the target velocity. The constant cross section
! free gas model is the default method. Methods for correctly accounting
! for the energy dependence of cross sections in treating resonance elastic
! scattering such as the DBRC, WCM, and a new, accelerated scheme are also
! implemented here.
!===============================================================================
subroutine sample_target_velocity(nuc, v_target, E, uvw, v_neut, wgt, xs_eff, kT)
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature T
real(8), intent(out) :: v_target(3) ! target velocity
real(8), intent(in) :: E ! particle energy
real(8), intent(in) :: uvw(3) ! direction cosines
real(8), intent(in) :: v_neut(3) ! neutron velocity
real(8), intent(inout) :: wgt ! particle weight
real(8), intent(in) :: xs_eff ! effective elastic xs at temperature T
real(8), intent(in) :: kT ! equilibrium temperature of target in eV
real(8) :: awr ! target/neutron mass ratio
real(8) :: E_rel ! trial relative energy
real(8) :: xs_0K ! 0K xs at E_rel
real(8) :: wcf ! weight correction factor
real(8) :: E_red ! reduced energy (same as used by Cullen in SIGMA1)
real(8) :: E_low ! lowest practical relative energy
real(8) :: E_up ! highest practical relative energy
real(8) :: E_t ! trial target energy
real(8) :: xs_max ! max 0K xs over practical relative energies
real(8) :: xs_low ! 0K xs at lowest practical relative energy
real(8) :: xs_up ! 0K xs at highest practical relative energy
real(8) :: m ! slope for interpolation
real(8) :: R ! rejection criterion for DBRC / target speed
real(8) :: cdf_low ! xs cdf at lowest practical relative energy
real(8) :: cdf_up ! xs cdf at highest practical relative energy
real(8) :: cdf_rel ! trial xs cdf value
real(8) :: mu ! cosine between neutron and target velocities
integer :: i_E_low ! 0K index to lowest practical relative energy
integer :: i_E_up ! 0K index to highest practical relative energy
integer :: i_E_rel ! index to trial relative energy
integer :: n_grid ! number of energies on 0K grid
integer :: sampling_method ! method of target velocity sampling
awr = nuc % awr
! check if nuclide is a resonant scatterer
if (nuc % resonant) then
! sampling method to use
sampling_method = res_scat_method
! upper resonance scattering energy bound (target is at rest above this E)
if (E > res_scat_energy_max) then
v_target = ZERO
return
! lower resonance scattering energy bound (should be no resonances below)
else if (E < res_scat_energy_min) then
sampling_method = RES_SCAT_CXS
end if
! otherwise, use free gas model
else
if (E >= FREE_GAS_THRESHOLD * kT .and. awr > ONE) then
v_target = ZERO
return
else
sampling_method = RES_SCAT_CXS
end if
end if
! use appropriate target velocity sampling method
select case (sampling_method)
case (RES_SCAT_CXS)
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
case (RES_SCAT_WCM)
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
! adjust weight as prescribed by the weight correction method (wcm)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
xs_0K = elastic_xs_0K(E_rel, nuc)
wcf = xs_0K / xs_eff
wgt = wcf * wgt
case (RES_SCAT_DBRC, RES_SCAT_ARES)
E_red = sqrt(awr * E / kT)
E_low = max(ZERO, E_red - FOUR)**2 * kT / awr
E_up = (E_red + FOUR)**2 * kT / awr
! find lower and upper energy bound indices
! lower index
n_grid = size(nuc % energy_0K)
if (E_low < nuc % energy_0K(1)) then
i_E_low = 1
elseif (E_low > nuc % energy_0K(n_grid)) then
i_E_low = n_grid - 1
else
i_E_low = binary_search(nuc % energy_0K, n_grid, E_low)
end if
! upper index
if (E_up < nuc % energy_0K(1)) then
i_E_up = 1
elseif (E_up > nuc % energy_0K(n_grid)) then
i_E_up = n_grid - 1
else
i_E_up = binary_search(nuc % energy_0K, n_grid, E_up)
end if
if (i_E_up == i_E_low) then
! Handle degenerate case -- if the upper/lower bounds occur for the same
! index, then using cxs is probably a good approximation
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
else
if (sampling_method == RES_SCAT_DBRC) then
! interpolate xs since we're not exactly at the energy indices
xs_low = nuc % elastic_0K(i_E_low)
m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
xs_up = nuc % elastic_0K(i_E_up)
m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
! get max 0K xs value over range of practical relative energies
xs_max = max(xs_low, &
maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up)
DBRC_REJECT_LOOP: do
TARGET_ENERGY_LOOP: do
! sample target velocity with the constant cross section (cxs) approx.
call sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
if (E_rel < E_up) exit TARGET_ENERGY_LOOP
end do TARGET_ENERGY_LOOP
! perform Doppler broadening rejection correction (dbrc)
xs_0K = elastic_xs_0K(E_rel, nuc)
R = xs_0K / xs_max
if (prn() < R) exit DBRC_REJECT_LOOP
end do DBRC_REJECT_LOOP
elseif (sampling_method == RES_SCAT_ARES) then
! interpolate xs CDF since we're not exactly at the energy indices
! cdf value at lower bound attainable energy
m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
/ (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
cdf_low = nuc % xs_cdf(i_E_low - 1) &
+ m * (E_low - nuc % energy_0K(i_E_low))
if (E_low <= nuc % energy_0K(1)) cdf_low = ZERO
! cdf value at upper bound attainable energy
m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
/ (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
cdf_up = nuc % xs_cdf(i_E_up - 1) &
+ m * (E_up - nuc % energy_0K(i_E_up))
ARES_REJECT_LOOP: do
! directly sample Maxwellian
E_t = -kT * log(prn())
! sample a relative energy using the xs cdf
cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
i_E_up - i_E_low + 2, cdf_rel)
E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
- nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
/ (nuc % energy_0K(i_E_low + i_E_rel) &
- nuc % energy_0K(i_E_low + i_E_rel - 1))
E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
! perform rejection sampling on cosine between
! neutron and target velocities
mu = (E_t + awr * (E - E_rel)) / (TWO * sqrt(awr * E * E_t))
if (abs(mu) < ONE) then
! set and accept target velocity
E_t = E_t / awr
v_target = sqrt(E_t) * rotate_angle(uvw, mu)
exit ARES_REJECT_LOOP
end if
end do ARES_REJECT_LOOP
end if
end if
end select
end subroutine sample_target_velocity
!===============================================================================
! SAMPLE_CXS_TARGET_VELOCITY samples a target velocity based on the free gas
! scattering formulation, used by most Monte Carlo codes, in which cross section
! is assumed to be constant in energy. Excellent documentation for this method
! can be found in FRA-TM-123.
!===============================================================================
subroutine sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
type(Nuclide), intent(in) :: nuc ! target nuclide at temperature
real(8), intent(out) :: v_target(3)
real(8), intent(in) :: E
real(8), intent(in) :: uvw(3)
real(8), intent(in) :: kT ! equilibrium temperature of target in eV
real(8) :: awr ! target/neutron mass ratio
real(8) :: alpha ! probability of sampling f2 over f1
real(8) :: mu ! cosine of angle between neutron and target vel
real(8) :: r1, r2 ! pseudo-random numbers
real(8) :: c ! cosine used in maxwell sampling
real(8) :: accept_prob ! probability of accepting combination of vt and mu
real(8) :: beta_vn ! beta * speed of neutron
real(8) :: beta_vt ! beta * speed of target
real(8) :: beta_vt_sq ! (beta * speed of target)^2
real(8) :: vt ! speed of target
awr = nuc % awr
beta_vn = sqrt(awr * E / kT)
alpha = ONE/(ONE + sqrt(pi)*beta_vn/TWO)
do
! Sample two random numbers
r1 = prn()
r2 = prn()
if (prn() < alpha) then
! With probability alpha, we sample the distribution p(y) =
! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
! Carlo sampler
beta_vt_sq = -log(r1*r2)
else
! With probability 1-alpha, we sample the distribution p(y) = y^2 *
! e^(-y^2). This can be done with sampling scheme C61 from the Monte
! Carlo sampler
c = cos(PI/TWO * prn())
beta_vt_sq = -log(r1) - log(r2)*c*c
end if
! Determine beta * vt
beta_vt = sqrt(beta_vt_sq)
! Sample cosine of angle between neutron and target velocity
mu = TWO*prn() - ONE
! Determine rejection probability
accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
/(beta_vn + beta_vt)
! Perform rejection sampling on vt and mu
if (prn() < accept_prob) exit
end do
! Determine speed of target nucleus
vt = sqrt(beta_vt_sq*kT/awr)
! Determine velocity vector of target nucleus based on neutron's velocity
! and the sampled angle between them
v_target = vt * rotate_angle(uvw, mu)
end subroutine sample_cxs_target_velocity
!===============================================================================
! INELASTIC_SCATTER handles all reactions with a single secondary neutron (other
! than fission), i.e. level scattering, (n,np), (n,na), etc.
!===============================================================================
subroutine inelastic_scatter(nuc, rxn, p)
type(Nuclide), intent(in) :: nuc
type(Reaction), intent(in) :: rxn
type(Particle), intent(inout) :: p
integer :: i ! loop index
real(8) :: E ! energy in lab (incoming/outgoing)
real(8) :: mu ! cosine of scattering angle in lab
real(8) :: A ! atomic weight ratio of nuclide
real(8) :: E_in ! incoming energy
real(8) :: E_cm ! outgoing energy in center-of-mass
real(8) :: yield ! neutron yield
! copy energy of neutron
E_in = p % E
! sample outgoing energy and scattering cosine
call rxn % product_sample(1, E_in, E, mu)
! if scattering system is in center-of-mass, transfer cosine of scattering
! angle and outgoing energy from CM to LAB
if (rxn % scatter_in_cm) then
E_cm = E
! determine outgoing energy in lab
A = nuc%awr
E = E_cm + (E_in + TWO * mu * (A+ONE) * sqrt(E_in * E_cm)) &
/ ((A+ONE)*(A+ONE))
! determine outgoing angle in lab
mu = mu * sqrt(E_cm/E) + ONE/(A+ONE) * sqrt(E_in/E)
end if
! Because of floating-point roundoff, it may be possible for mu to be
! outside of the range [-1,1). In these cases, we just set mu to exactly -1
! or 1
if (abs(mu) > ONE) mu = sign(ONE,mu)
! Set outgoing energy and scattering angle
p % E = E
p % mu = mu
! change direction of particle
p % coord(1) % uvw = rotate_angle(p % coord(1) % uvw, mu)
! evaluate yield
yield = rxn % product_yield(1, E_in)
if (mod(yield, ONE) == ZERO) then
! If yield is integral, create exactly that many secondary particles
do i = 1, nint(yield) - 1
call particle_create_secondary(p, p % coord(1) % uvw, p % E, &
NEUTRON, run_CE=.true._C_BOOL)
end do
else
! Otherwise, change weight of particle based on yield
p % wgt = yield * p % wgt
end if
end subroutine inelastic_scatter
end module physics

721
src/physics.cpp
View file

@ -4,6 +4,7 @@
#include "openmc/constants.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/material.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/nuclide.h"
@ -11,12 +12,15 @@
#include "openmc/physics_common.h"
#include "openmc/random_lcg.h"
#include "openmc/reaction.h"
#include "openmc/secondary_uncorrelated.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/thermal.h"
#include "openmc/tallies/tally.h"
#include <algorithm> // for max, min
#include <cmath> // for sqrt, exp, log
#include <algorithm> // for max, min, max_element
#include <cmath> // for sqrt, exp, log, abs, copysign
#include <sstream>
namespace openmc {
@ -509,24 +513,24 @@ Reaction* sample_fission(int i_nuclide, double E)
// }
// Get grid index and interpolatoin factor and sample fission cdf
int i_temp = simulation::micro_xs[i_nuclide-1].index_temp;
int i_temp = simulation::micro_xs[i_nuclide-1].index_temp - 1;
int i_grid = simulation::micro_xs[i_nuclide-1].index_grid;
double f = simulation::micro_xs[i_nuclide-1].interp_factor;
double cutoff = prn() * simulation::micro_xs[i_nuclide-1].fission;
double prob = 0.0;
// Loop through each partial fission reaction type
for (auto& rx : nuc->reactions_) {
for (auto& rx : nuc->fission_rx_) {
// if energy is below threshold for this reaction, skip it
int threshold = rx->xs_[i_temp-1].threshold;
int threshold = rx->xs_[i_temp].threshold;
if (i_grid < threshold) continue;
// add to cumulative probability
prob += (1.0 - f) * rx->xs_[i_temp-1].value[i_grid - threshold]
prob += (1.0 - f) * rx->xs_[i_temp].value[i_grid - threshold]
+ f*rx->xs_[i_temp].value[i_grid - threshold + 1];
// Create fission bank sites if fission occurs
if (prob > cutoff) break;
if (prob > cutoff) return rx;
}
}
@ -596,173 +600,183 @@ void absorption(Particle* p, int i_nuclide)
}
}
// void scatter(Particle*, int i_nuclide, int i_nuc_mat)
// {
// // copy incoming direction
// uvw_old(:) = p->coord(1) % uvw
void scatter(Particle* p, int i_nuclide, int i_nuc_mat)
{
// copy incoming direction
Direction u_old {p->coord[0].uvw};
// // Get pointer to nuclide and grid index/interpolation factor
// nuc => nuclides(i_nuclide)
// i_temp = simulation::micro_xs[i_nuclide-1].index_temp
// i_grid = simulation::micro_xs[i_nuclide-1].index_grid
// f = simulation::micro_xs[i_nuclide-1].interp_factor
// Get pointer to nuclide and grid index/interpolation factor
const auto& nuc {data::nuclides[i_nuclide-1]};
const auto& micro {simulation::micro_xs[i_nuclide-1]};
int i_temp = micro.index_temp - 1;
int i_grid = micro.index_grid - 1;
double f = micro.interp_factor;
// // For tallying purposes, this routine might be called directly. In that
// // case, we need to sample a reaction via the cutoff variable
// cutoff = prn() * (micro_xs[i_nuclide-1].total - &
// simulation::micro_xs[i_nuclide-1].absorption)
// sampled = false
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
double cutoff = prn() * (micro.total - micro.absorption);
bool sampled = false;
// // Calculate elastic cross section if it wasn't precalculated
// if (micro_xs[i_nuclide-1].elastic == CACHE_INVALID) {
// nuc % calculate_elastic_xs(micro_xs(i_nuclide))
// }
// Calculate elastic cross section if it wasn't precalculated
if (micro.elastic == CACHE_INVALID) {
nuc->calculate_elastic_xs(i_nuclide);
}
// prob = simulation::micro_xs[i_nuclide-1].elastic - simulation::micro_xs[i_nuclide-1].thermal
// if (prob > cutoff) {
// // =======================================================================
// // NON-S(A,B) ELASTIC SCATTERING
double prob = micro.elastic - micro.thermal;
if (prob > cutoff) {
// =======================================================================
// NON-S(A,B) ELASTIC SCATTERING
// // Determine temperature
// if (nuc % mp_present) {
// kT = p->sqrtkT**2
// } else {
// kT = nuc % kTs(micro_xs[i_nuclide-1].index_temp)
// }
// Determine temperature
double kT;
// if (nuc % mp_present) {
// kT = p->sqrtkT**2
// } else {
kT = nuc->kTs_[i_temp];
// }
// // Perform collision physics for elastic scattering
// elastic_scatter(i_nuclide, nuc % reactions(1), kT, p->E, &
// p->coord(1) % uvw, p->mu, p->wgt)
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, nuc->reactions_[0].get(), kT,
&p->E, p->coord[0].uvw, &p->mu, &p->wgt);
// p->event_MT = ELASTIC
// sampled = true
// }
p->event_MT = ELASTIC;
sampled = true;
}
// prob = simulation::micro_xs[i_nuclide-1].elastic
// if (prob > cutoff && !sampled) {
// // =======================================================================
// // S(A,B) SCATTERING
prob = micro.elastic;
if (prob > cutoff && !sampled) {
// =======================================================================
// S(A,B) SCATTERING
// sab_scatter(i_nuclide, simulation::micro_xs[i_nuclide-1].index_sab, p->E, &
// p->coord(1) % uvw, p->mu)
sab_scatter(i_nuclide, micro.index_sab, &p->E, p->coord[0].uvw, &p->mu);
// p->event_MT = ELASTIC
// sampled = true
// }
p->event_MT = ELASTIC;
sampled = true;
}
// if (!sampled) {
// // =======================================================================
// // INELASTIC SCATTERING
if (!sampled) {
// =======================================================================
// INELASTIC SCATTERING
// j = 0
// do while (prob < cutoff)
// j = j + 1
// i = nuc % index_inelastic_scatter(j)
int j = 0;
int i;
while (prob < cutoff) {
i = nuc->index_inelastic_scatter_[j];
++j;
// // Check to make sure inelastic scattering reaction sampled
// if (i > size(nuc % reactions)) {
// particle_write_restart(p)
// fatal_error("Did not sample any reaction for nuclide " &
// &// trim(nuc % name))
// }
// Check to make sure inelastic scattering reaction sampled
if (i >= nuc->reactions_.size()) {
p->write_restart();
fatal_error("Did not sample any reaction for nuclide " + nuc->name_);
}
// associate (rx => nuc % reactions(i))
// // if energy is below threshold for this reaction, skip it
// threshold = rx % xs_threshold(i_temp)
// if (i_grid < threshold) cycle
// if energy is below threshold for this reaction, skip it
const auto& xs {nuc->reactions_[i]->xs_[i_temp]};
int threshold = xs.threshold - 1;
if (i_grid < threshold) continue;
// // add to cumulative probability
// prob = prob + ((1.0 - f)*rx % xs(i_temp, i_grid - threshold + 1) &
// + f*(rx % xs(i_temp, i_grid - threshold + 2)))
// end associate
// end do
// add to cumulative probability
prob += (1.0 - f)*xs.value[i_grid - threshold] +
f*xs.value[i_grid - threshold + 1];
}
// // Perform collision physics for inelastic scattering
// inelastic_scatter(nuc, nuc%reactions(i), p)
// p->event_MT = nuc % reactions(i) % MT
// Perform collision physics for inelastic scattering
const auto& rx {nuc->reactions_[i]};
inelastic_scatter(nuc.get(), rx.get(), p);
p->event_MT = rx->mt_;
}
// }
// Set event component
p->event = EVENT_SCATTER;
// // Set event component
// p->event = EVENT_SCATTER
// Sample new outgoing angle for isotropic-in-lab scattering
if (material_isotropic(p->material, i_nuc_mat)) {
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn() - 1.0;
double phi = 2.0*PI*prn();
Direction u_new;
u_new.x = mu;
u_new.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
u_new.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
// // Sample new outgoing angle for isotropic-in-lab scattering
// associate (mat => materials(p->material))
// if (mat % has_isotropic_nuclides) {
// if (materials(p->material) % p0(i_nuc_mat)) {
// // Sample isotropic-in-lab outgoing direction
// uvw_new(1) = 2.0 * prn() - 1.0
// phi = 2.0 * PI * prn()
// uvw_new(2) = std::cos(phi) * std::sqrt(1.0 - uvw_new(1)*uvw_new(1))
// uvw_new(3) = std::sin(phi) * std::sqrt(1.0 - uvw_new(1)*uvw_new(1))
// p->mu = dot_product(uvw_old, uvw_new)
p->mu = u_old.dot(u_new);
// // Change direction of particle
// p->coord(1) % uvw = uvw_new
// }
// }
// end associate
// }
// Change direction of particle
p->coord[0].uvw[0] = u_new.x;
p->coord[0].uvw[1] = u_new.y;
p->coord[0].uvw[2] = u_new.z;
}
}
// void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, double* E,
// Direction& u, double* mu_lab, double* wgt)
// {
// // get pointer to nuclide
// nuc => nuclides(i_nuclide)
void elastic_scatter(int i_nuclide, const Reaction* rx, double kT, double* E,
double* uvw, double* mu_lab, double* wgt)
{
// get pointer to nuclide
const auto& nuc {data::nuclides[i_nuclide-1]};
// vel = std::sqrt(E)
// awr = nuc % awr
double vel = std::sqrt(*E);
double awr = nuc->awr_;
// // Neutron velocity in LAB
// v_n = vel * uvw
// Neutron velocity in LAB
Direction u {uvw};
Direction v_n = vel*u;
// // Sample velocity of target nucleus
// if (!micro_xs[i_nuclide-1].use_ptable) {
// sample_target_velocity(nuc, v_t, E, uvw, v_n, wgt, &
// simulation::micro_xs[i_nuclide-1].elastic, kT)
// } else {
// v_t = 0.0
// }
// Sample velocity of target nucleus
Direction v_t {};
if (!simulation::micro_xs[i_nuclide-1].use_ptable) {
v_t = sample_target_velocity(nuc.get(), *E, u, v_n,
simulation::micro_xs[i_nuclide-1].elastic, kT, wgt);
}
// // Velocity of center-of-mass
// v_cm = (v_n + awr*v_t)/(awr + 1.0)
// Velocity of center-of-mass
Direction v_cm = (v_n + awr*v_t)/(awr + 1.0);
// // Transform to CM frame
// v_n = v_n - v_cm
// Transform to CM frame
v_n -= v_cm;
// // Find speed of neutron in CM
// vel = std::sqrt(dot_product(v_n, v_n))
// Find speed of neutron in CM
vel = v_n.norm();
// // Sample scattering angle
// mu_cm = rxn % sample_elastic_mu(E)
// Sample scattering angle, checking if it is an ncorrelated angle-energy
// distribution
double mu_cm;
auto& d = rx->products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) {
mu_cm = d_->angle().sample(*E);
} else {
mu_cm = 2.0*prn() - 1.0;
}
// // Determine direction cosines in CM
// uvw_cm = v_n/vel
// Determine direction cosines in CM
Direction u_cm = v_n/vel;
// // Rotate neutron velocity vector to new angle -- note that the speed of the
// // neutron in CM does not change in elastic scattering. However, the speed
// // will change when we convert back to LAB
// v_n = vel * rotate_angle(uvw_cm, mu_cm)
// Rotate neutron velocity vector to new angle -- note that the speed of the
// neutron in CM does not change in elastic scattering. However, the speed
// will change when we convert back to LAB
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr);
// // Transform back to LAB frame
// v_n = v_n + v_cm
// Transform back to LAB frame
v_n += v_cm;
// E = dot_product(v_n, v_n)
// vel = std::sqrt(E)
*E = v_n.dot(v_n);
vel = std::sqrt(*E);
// // compute cosine of scattering angle in LAB frame by taking dot product of
// // neutron's pre- and post-collision angle
// mu_lab = dot_product(uvw, v_n) / vel
// compute cosine of scattering angle in LAB frame by taking dot product of
// neutron's pre- and post-collision angle
*mu_lab = u.dot(v_n) / vel;
// // Set energy and direction of particle in LAB frame
// uvw = v_n / vel
// Set energy and direction of particle in LAB frame
u = v_n / vel;
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
// // Because of floating-point roundoff, it may be possible for mu_lab to be
// // outside of the range [-1,1). In these cases, we just set mu_lab to exactly
// // -1 or 1
// if (abs(mu_lab) > 1.0) mu_lab = sign(1.0,mu_lab)
// }
// Because of floating-point roundoff, it may be possible for mu_lab to be
// outside of the range [-1,1). In these cases, we just set mu_lab to exactly
// -1 or 1
if (std::abs(*mu_lab) > 1.0) *mu_lab = std::copysign(1.0, *mu_lab);
}
// void sab_scatter(int i_nuclide, int i_sab, double* E, Direction* u, double* mu)
// {
@ -775,212 +789,217 @@ void absorption(Particle* p, int i_nuclide)
// uvw = rotate_angle(uvw, mu)
// }
// void sample_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u,
// Direction v_neut, double* wgt, double xs_eff, double kT)
// {
// awr = nuc % awr
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_neut, double xs_eff, double kT, double* wgt)
{
// check if nuclide is a resonant scatterer
int sampling_method;
if (nuc->resonant_) {
// // check if nuclide is a resonant scatterer
// if (nuc % resonant) {
// sampling method to use
sampling_method = settings::res_scat_method;
// // sampling method to use
// sampling_method = res_scat_method
// upper resonance scattering energy bound (target is at rest above this E)
if (E > settings::res_scat_energy_max) {
return {};
// // upper resonance scattering energy bound (target is at rest above this E)
// if (E > res_scat_energy_max) {
// v_target = 0.0
// return
// lower resonance scattering energy bound (should be no resonances below)
} else if (E < settings::res_scat_energy_min) {
sampling_method = RES_SCAT_CXS;
}
// // lower resonance scattering energy bound (should be no resonances below)
// } else if (E < res_scat_energy_min) {
// sampling_method = RES_SCAT_CXS
// }
// otherwise, use free gas model
} else {
if (E >= FREE_GAS_THRESHOLD * kT && nuc->awr_ > 1.0) {
return {};
} else {
sampling_method = RES_SCAT_CXS;
}
}
// // otherwise, use free gas model
// } else {
// if (E >= FREE_GAS_THRESHOLD * kT && awr > 1.0) {
// v_target = 0.0
// return
// } else {
// sampling_method = RES_SCAT_CXS
// }
// }
// use appropriate target velocity sampling method
switch (sampling_method) {
case RES_SCAT_CXS:
// // use appropriate target velocity sampling method
// select case (sampling_method)
// case (RES_SCAT_CXS)
// sample target velocity with the constant cross section (cxs) approx.
return sample_cxs_target_velocity(nuc->awr_, E, u, kT);
// // sample target velocity with the constant cross section (cxs) approx.
// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
case RES_SCAT_WCM: {
// sample target velocity with the constant cross section (cxs) approx.
Direction v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT);
// case (RES_SCAT_WCM)
// adjust weight as prescribed by the weight correction method (wcm)
Direction v_rel = v_neut - v_target;
double E_rel = v_rel.dot(v_rel);
double xs_0K = nuc->elastic_xs_0K(E_rel);
*wgt *= xs_0K / xs_eff;
return v_target;
}
// // sample target velocity with the constant cross section (cxs) approx.
// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
case RES_SCAT_DBRC:
case RES_SCAT_ARES: {
double E_red = std::sqrt(nuc->awr_ * E / kT);
double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc->awr_;
double E_up = (E_red + 4.0)*(E_red + 4.0) * kT / nuc->awr_;
// // adjust weight as prescribed by the weight correction method (wcm)
// E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
// xs_0K = elastic_xs_0K(E_rel, nuc)
// wcf = xs_0K / xs_eff
// wgt = wcf * wgt
// find lower and upper energy bound indices
// lower index
int i_E_low;
if (E_low < nuc->energy_0K_.front()) {
i_E_low = 0;
} else if (E_low > nuc->energy_0K_.back()) {
i_E_low = nuc->energy_0K_.size() - 2;
} else {
i_E_low = lower_bound_index(nuc->energy_0K_.begin(),
nuc->energy_0K_.end(), E_low);
}
// case (RES_SCAT_DBRC, RES_SCAT_ARES)
// E_red = std::sqrt(awr * E / kT)
// E_low = std::max(0.0, E_red - FOUR)**2 * kT / awr
// E_up = (E_red + FOUR)**2 * kT / awr
// upper index
int i_E_up;
if (E_up < nuc->energy_0K_.front()) {
i_E_up = 0;
} else if (E_up > nuc->energy_0K_.back()) {
i_E_up = nuc->energy_0K_.size() - 2;
} else {
i_E_up = lower_bound_index(nuc->energy_0K_.begin(),
nuc->energy_0K_.end(), E_up);
}
// // find lower and upper energy bound indices
// // lower index
// n_grid = size(nuc % energy_0K)
// if (E_low < nuc % energy_0K(1)) {
// i_E_low = 1
// } else if (E_low > nuc % energy_0K(n_grid)) {
// i_E_low = n_grid - 1
// } else {
// i_E_low = binary_search(nuc % energy_0K, n_grid, E_low)
// }
if (i_E_up == i_E_low) {
// Handle degenerate case -- if the upper/lower bounds occur for the same
// index, then using cxs is probably a good approximation
return sample_cxs_target_velocity(nuc->awr_, E, u, kT);
}
// // upper index
// if (E_up < nuc % energy_0K(1)) {
// i_E_up = 1
// } else if (E_up > nuc % energy_0K(n_grid)) {
// i_E_up = n_grid - 1
// } else {
// i_E_up = binary_search(nuc % energy_0K, n_grid, E_up)
// }
if (sampling_method == RES_SCAT_DBRC) {
// interpolate xs since we're not exactly at the energy indices
double xs_low = nuc->elastic_0K_[i_E_low];
double m = (nuc->elastic_0K_[i_E_low + 1] - xs_low)
/ (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]);
xs_low += m * (E_low - nuc->energy_0K_[i_E_low]);
double xs_up = nuc->elastic_0K_[i_E_up];
m = (nuc->elastic_0K_[i_E_up + 1] - xs_up)
/ (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]);
xs_up += m * (E_up - nuc->energy_0K_[i_E_up]);
// if (i_E_up == i_E_low) {
// // Handle degenerate case -- if the upper/lower bounds occur for the same
// // index, then using cxs is probably a good approximation
// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
// get max 0K xs value over range of practical relative energies
double xs_max = *std::max_element(&nuc->elastic_0K_[i_E_low + 1],
&nuc->elastic_0K_[i_E_up + 1]);
xs_max = std::max({xs_low, xs_max, xs_up});
// } else {
// if (sampling_method == RES_SCAT_DBRC) {
// // interpolate xs since we're not exactly at the energy indices
// xs_low = nuc % elastic_0K(i_E_low)
// m = (nuc % elastic_0K(i_E_low + 1) - xs_low) &
// / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
// xs_low = xs_low + m * (E_low - nuc % energy_0K(i_E_low))
// xs_up = nuc % elastic_0K(i_E_up)
// m = (nuc % elastic_0K(i_E_up + 1) - xs_up) &
// / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
// xs_up = xs_up + m * (E_up - nuc % energy_0K(i_E_up))
while (true) {
double E_rel;
Direction v_target;
while (true) {
// sample target velocity with the constant cross section (cxs) approx.
v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT);
Direction v_rel = v_neut - v_target;
E_rel = v_rel.dot(v_rel);
if (E_rel < E_up) break;
}
// // get max 0K xs value over range of practical relative energies
// xs_max = std::max(xs_low, &
// maxval(nuc % elastic_0K(i_E_low + 1 : i_E_up)), xs_up)
// perform Doppler broadening rejection correction (dbrc)
double xs_0K = nuc->elastic_xs_0K(E_rel);
double R = xs_0K / xs_max;
if (prn() < R) return v_target;
}
// DBRC_REJECT_LOOP: do
// TARGET_ENERGY_LOOP: do
// // sample target velocity with the constant cross section (cxs) approx.
// sample_cxs_target_velocity(nuc, v_target, E, uvw, kT)
// E_rel = dot_product((v_neut - v_target), (v_neut - v_target))
// if (E_rel < E_up) exit TARGET_ENERGY_LOOP
// end do TARGET_ENERGY_LOOP
} else if (sampling_method == RES_SCAT_ARES) {
// interpolate xs CDF since we're not exactly at the energy indices
// cdf value at lower bound attainable energy
double m = (nuc->xs_cdf_[i_E_low] - nuc->xs_cdf_[i_E_low - 1])
/ (nuc->energy_0K_[i_E_low + 1] - nuc->energy_0K_[i_E_low]);
double cdf_low = nuc->xs_cdf_[i_E_low - 1]
+ m * (E_low - nuc->energy_0K_[i_E_low]);
if (E_low <= nuc->energy_0K_.front()) cdf_low = 0.0;
// // perform Doppler broadening rejection correction (dbrc)
// xs_0K = elastic_xs_0K(E_rel, nuc)
// R = xs_0K / xs_max
// if (prn() < R) exit DBRC_REJECT_LOOP
// end do DBRC_REJECT_LOOP
// cdf value at upper bound attainable energy
m = (nuc->xs_cdf_[i_E_up] - nuc->xs_cdf_[i_E_up - 1])
/ (nuc->energy_0K_[i_E_up + 1] - nuc->energy_0K_[i_E_up]);
double cdf_up = nuc->xs_cdf_[i_E_up - 1]
+ m*(E_up - nuc->energy_0K_[i_E_up]);
// } else if (sampling_method == RES_SCAT_ARES) {
// // interpolate xs CDF since we're not exactly at the energy indices
// // cdf value at lower bound attainable energy
// m = (nuc % xs_cdf(i_E_low) - nuc % xs_cdf(i_E_low - 1)) &
// / (nuc % energy_0K(i_E_low + 1) - nuc % energy_0K(i_E_low))
// cdf_low = nuc % xs_cdf(i_E_low - 1) &
// + m * (E_low - nuc % energy_0K(i_E_low))
// if (E_low <= nuc % energy_0K(1)) cdf_low = 0.0
while (true) {
// directly sample Maxwellian
double E_t = -kT * std::log(prn());
// // cdf value at upper bound attainable energy
// m = (nuc % xs_cdf(i_E_up) - nuc % xs_cdf(i_E_up - 1)) &
// / (nuc % energy_0K(i_E_up + 1) - nuc % energy_0K(i_E_up))
// cdf_up = nuc % xs_cdf(i_E_up - 1) &
// + m * (E_up - nuc % energy_0K(i_E_up))
// sample a relative energy using the xs cdf
double cdf_rel = cdf_low + prn()*(cdf_up - cdf_low);
int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1],
&nuc->xs_cdf_[i_E_up+1], cdf_rel);
double E_rel = nuc->energy_0K_[i_E_low + i_E_rel];
double m = (nuc->xs_cdf_[i_E_low + i_E_rel]
- nuc->xs_cdf_[i_E_low + i_E_rel - 1])
/ (nuc->energy_0K_[i_E_low + i_E_rel + 1]
- nuc->energy_0K_[i_E_low + i_E_rel]);
E_rel += (cdf_rel - nuc->xs_cdf_[i_E_low + i_E_rel - 1]) / m;
// ARES_REJECT_LOOP: do
// perform rejection sampling on cosine between
// neutron and target velocities
double mu = (E_t + nuc->awr_ * (E - E_rel)) /
(2.0 * std::sqrt(nuc->awr_ * E * E_t));
// // directly sample Maxwellian
// E_t = -kT * std::log(prn())
if (std::abs(mu) < 1.0) {
// set and accept target velocity
E_t /= nuc->awr_;
return std::sqrt(E_t) * rotate_angle(u, mu, nullptr);
}
}
}
}
}
}
// // sample a relative energy using the xs cdf
// cdf_rel = cdf_low + prn() * (cdf_up - cdf_low)
// i_E_rel = binary_search(nuc % xs_cdf(i_E_low-1:i_E_up), &
// i_E_up - i_E_low + 2, cdf_rel)
// E_rel = nuc % energy_0K(i_E_low + i_E_rel - 1)
// m = (nuc % xs_cdf(i_E_low + i_E_rel - 1) &
// - nuc % xs_cdf(i_E_low + i_E_rel - 2)) &
// / (nuc % energy_0K(i_E_low + i_E_rel) &
// - nuc % energy_0K(i_E_low + i_E_rel - 1))
// E_rel = E_rel + (cdf_rel - nuc % xs_cdf(i_E_low + i_E_rel - 2)) / m
Direction
sample_cxs_target_velocity(double awr, double E, Direction u, double kT)
{
double beta_vn = std::sqrt(awr * E / kT);
double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0);
// // perform rejection sampling on cosine between
// // neutron and target velocities
// mu = (E_t + awr * (E - E_rel)) / (2.0 * std::sqrt(awr * E * E_t))
double beta_vt_sq;
double mu;
while (true) {
// Sample two random numbers
double r1 = prn();
double r2 = prn();
// if (abs(mu) < 1.0) {
// // set and accept target velocity
// E_t = E_t / awr
// v_target = std::sqrt(E_t) * rotate_angle(uvw, mu)
// exit ARES_REJECT_LOOP
// }
// end do ARES_REJECT_LOOP
// }
// }
// end select
// }
if (prn() < alpha) {
// With probability alpha, we sample the distribution p(y) =
// y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
// Carlo sampler
// void sample_cxs_target_velocity(int i_nuclide, Direction* v_target, double E, Direction u,
// double kT)
// {
// awr = nuc % awr
beta_vt_sq = -std::log(r1*r2);
// beta_vn = std::sqrt(awr * E / kT)
// alpha = 1.0/(1.0 + std::sqrt(pi)*beta_vn/2.0)
} else {
// With probability 1-alpha, we sample the distribution p(y) = y^2 *
// e^(-y^2). This can be done with sampling scheme C61 from the Monte
// Carlo sampler
// do
// // Sample two random numbers
// r1 = prn()
// r2 = prn()
double c = std::cos(PI/2.0 * prn());
beta_vt_sq = -std::log(r1) - std::log(r2)*c*c;
}
// if (prn() < alpha) {
// // With probability alpha, we sample the distribution p(y) =
// // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
// // Carlo sampler
// Determine beta * vt
double beta_vt = std::sqrt(beta_vt_sq);
// beta_vt_sq = -std::log(r1*r2)
// Sample cosine of angle between neutron and target velocity
mu = 2.0*prn() - 1.0;
// } else {
// // With probability 1-alpha, we sample the distribution p(y) = y^2 *
// // e^(-y^2). This can be done with sampling scheme C61 from the Monte
// // Carlo sampler
// Determine rejection probability
double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq -
2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt);
// c = std::cos(PI/2.0 * prn())
// beta_vt_sq = -std::log(r1) - std::log(r2)*c*c
// }
// Perform rejection sampling on vt and mu
if (prn() < accept_prob) break;
}
// // Determine beta * vt
// beta_vt = std::sqrt(beta_vt_sq)
// Determine speed of target nucleus
double vt = std::sqrt(beta_vt_sq*kT/awr);
// // Sample cosine of angle between neutron and target velocity
// mu = 2.0*prn() - 1.0
// // Determine rejection probability
// accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
// /(beta_vn + beta_vt)
// // Perform rejection sampling on vt and mu
// if (prn() < accept_prob) exit
// end do
// // Determine speed of target nucleus
// vt = std::sqrt(beta_vt_sq*kT/awr)
// // Determine velocity vector of target nucleus based on neutron's velocity
// // and the sampled angle between them
// v_target = vt * rotate_angle(uvw, mu)
// }
// Determine velocity vector of target nucleus based on neutron's velocity
// and the sampled angle between them
return vt * rotate_angle(u, mu, nullptr);
}
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank* site)
{
@ -1073,53 +1092,55 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Bank
}
}
// void inelastic_scatter(int i_nuclide, const Reaction& rx, Particle* p)
// {
// // copy energy of neutron
// E_in = p->E
void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
{
// copy energy of neutron
double E_in = p->E;
// // sample outgoing energy and scattering cosine
// rxn % product_sample(1, E_in, E, mu)
// sample outgoing energy and scattering cosine
double E;
double mu;
rx->products_[0].sample(E_in, E, mu);
// // if scattering system is in center-of-mass, transfer cosine of scattering
// // angle and outgoing energy from CM to LAB
// if (rxn % scatter_in_cm) {
// E_cm = E
// if scattering system is in center-of-mass, transfer cosine of scattering
// angle and outgoing energy from CM to LAB
if (rx->scatter_in_cm_) {
double E_cm = E;
// // determine outgoing energy in lab
// A = nuc%awr
// E = E_cm + (E_in + 2.0 * mu * (A+1.0) * std::sqrt(E_in * E_cm)) &
// / ((A+1.0)*(A+1.0))
// determine outgoing energy in lab
double A = nuc->awr_;
E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm))
/ ((A + 1.0)*(A + 1.0));
// // determine outgoing angle in lab
// mu = mu * std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E)
// }
// determine outgoing angle in lab
mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E);
}
// // Because of floating-point roundoff, it may be possible for mu to be
// // outside of the range [-1,1). In these cases, we just set mu to exactly -1
// // or 1
// if (abs(mu) > 1.0) mu = sign(1.0,mu)
// Because of floating-point roundoff, it may be possible for mu to be
// outside of the range [-1,1). In these cases, we just set mu to exactly -1
// or 1
if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
// // Set outgoing energy and scattering angle
// p->E = E
// p->mu = mu
// Set outgoing energy and scattering angle
p->E = E;
p->mu = mu;
// // change direction of particle
// p->coord(1) % uvw = rotate_angle(p->coord(1) % uvw, mu)
// change direction of particle
rotate_angle_c(p->coord[0].uvw, mu, nullptr);
// // evaluate yield
// yield = rxn % product_yield(1, E_in)
// if (mod(yield, 1.0) == 0.0) {
// // If yield is integral, create exactly that many secondary particles
// do i = 1, nint(yield) - 1
// particle_create_secondary(p, p->coord(1) % uvw, p->E, &
// NEUTRON, run_CE=true)
// end do
// } else {
// // Otherwise, change weight of particle based on yield
// p->wgt = yield * p->wgt
// }
// }
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
if (std::floor(yield) == yield) {
// If yield is integral, create exactly that many secondary particles
for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
int neutron = static_cast<int>(ParticleType::neutron);
p->create_secondary(p->coord[0].uvw, p->E, neutron, true);
}
} else {
// Otherwise, change weight of particle based on yield
p->wgt *= yield;
}
}
void sample_secondary_photons(Particle* p, int i_nuclide)
{

18
src/position.cpp
View file

@ -60,4 +60,22 @@ Position::operator*=(double v)
return *this;
}
Position&
Position::operator/=(Position other)
{
x /= other.x;
y /= other.y;
z /= other.z;
return *this;
}
Position&
Position::operator/=(double v)
{
x /= v;
y /= v;
z /= v;
return *this;
}
} // namespace openmc

12
src/settings.cpp
View file

@ -806,18 +806,6 @@ void read_settings_xml()
//==============================================================================
extern "C" {
bool res_scat_nuclides_empty() {
return settings::res_scat_nuclides.empty();
}
int res_scat_nuclides_size() {
return settings::res_scat_nuclides.size();
}
bool res_scat_nuclides_cmp(int i, const char* name) {
return settings::res_scat_nuclides[i - 1] == name;
}
const char* path_cross_sections_c() {
return settings::path_cross_sections.c_str();
}