mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Merge pull request #1092 from nelsonag/physics_mg_cpp
Convert physics_common.F90 and physics_mg.F90 to C++
This commit is contained in:
Paul Romano
GitHub
commit
fd12a55b73
26 changed files with 437 additions and 441 deletions
|
|
@ -345,7 +345,6 @@ add_library(libopenmc SHARED
|
|||
src/photon_physics.F90
|
||||
src/physics_common.F90
|
||||
src/physics.F90
|
||||
src/physics_mg.F90
|
||||
src/plot.F90
|
||||
src/plot_header.F90
|
||||
src/progress_header.F90
|
||||
|
|
@ -423,6 +422,8 @@ add_library(libopenmc SHARED
|
|||
src/nuclide.cpp
|
||||
src/output.cpp
|
||||
src/particle.cpp
|
||||
src/physics_common.cpp
|
||||
src/physics_mg.cpp
|
||||
src/plot.cpp
|
||||
src/position.cpp
|
||||
src/pugixml/pugixml_c.cpp
|
||||
|
|
|
|||
|
|
@ -149,6 +149,8 @@ extern "C" {
|
|||
extern int32_t n_tallies;
|
||||
extern int32_t n_universes;
|
||||
extern bool openmc_simulation_initialized;
|
||||
extern double global_tally_absorption;
|
||||
#pragma omp threadprivate(global_tally_absorption)
|
||||
|
||||
// Variables that are shared by necessity (can be removed from public header
|
||||
// later)
|
||||
|
|
|
|||
|
|
@ -26,6 +26,7 @@ class Material
|
|||
public:
|
||||
int32_t id; //!< Unique ID
|
||||
double volume_ {-1.0}; //!< Volume in [cm^3]
|
||||
bool fissionable {false}; //!< Does this material contain fissionable nuclides
|
||||
|
||||
//! \brief Default temperature for cells containing this material.
|
||||
//!
|
||||
|
|
|
|||
|
|
@ -22,7 +22,7 @@ namespace openmc {
|
|||
//! @return The requested percentile
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double normal_percentile_c(double p);
|
||||
extern "C" double normal_percentile(double p);
|
||||
|
||||
//==============================================================================
|
||||
//! Calculate the percentile of the Student's t distribution with a specified
|
||||
|
|
@ -58,7 +58,7 @@ extern "C" void calc_pn_c(int n, double x, double pnx[]);
|
|||
//! evaluated at x
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double evaluate_legendre_c(int n, const double data[], double x);
|
||||
extern "C" double evaluate_legendre(int n, const double data[], double x);
|
||||
|
||||
//==============================================================================
|
||||
//! Calculate the n-th order real spherical harmonics for a given angle (in
|
||||
|
|
@ -144,7 +144,7 @@ Direction rotate_angle(Direction u, double mu, double* phi);
|
|||
//! @result The sampled outgoing energy
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double maxwell_spectrum_c(double T);
|
||||
extern "C" double maxwell_spectrum(double T);
|
||||
|
||||
//==============================================================================
|
||||
//! Samples an energy from a Watt energy-dependent fission distribution.
|
||||
|
|
@ -159,7 +159,7 @@ extern "C" double maxwell_spectrum_c(double T);
|
|||
//! @result The sampled outgoing energy
|
||||
//==============================================================================
|
||||
|
||||
extern "C" double watt_spectrum_c(double a, double b);
|
||||
extern "C" double watt_spectrum(double a, double b);
|
||||
|
||||
//==============================================================================
|
||||
//! Doppler broadens the windowed multipole curvefit.
|
||||
|
|
|
|||
|
|
@ -18,6 +18,9 @@ extern std::vector<Mgxs> nuclides_MG;
|
|||
extern std::vector<Mgxs> macro_xs;
|
||||
extern "C" int num_energy_groups;
|
||||
|
||||
//TODO: When more of the Fortran is converted (input_xml, tallies, etc, also
|
||||
// bring over energy_bin_avg, energy_bins, etc, as vectors)
|
||||
|
||||
//==============================================================================
|
||||
// Mgxs data loading interface methods
|
||||
//==============================================================================
|
||||
|
|
@ -44,13 +47,6 @@ extern "C" void
|
|||
calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
|
||||
double& total_xs, double& abs_xs, double& nu_fiss_xs);
|
||||
|
||||
extern "C" void
|
||||
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
|
||||
double uvw[3]);
|
||||
|
||||
extern "C" void
|
||||
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout);
|
||||
|
||||
extern "C" double
|
||||
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg);
|
||||
|
||||
|
|
|
|||
16
include/openmc/physics_common.h
Normal file
16
include/openmc/physics_common.h
Normal file
|
|
@ -0,0 +1,16 @@
|
|||
//! \file physics_common.h
|
||||
//! A collection of physics methods common to MG, CE, photon, etc.
|
||||
|
||||
#ifndef OPENMC_PHYSICS_COMMON_H
|
||||
#define OPENMC_PHYSICS_COMMON_H
|
||||
|
||||
#include "openmc/particle.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//! \brief Performs the russian roulette operation for a particle
|
||||
extern "C" void
|
||||
russian_roulette(Particle* p);
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_PHYSICS_COMMON_H
|
||||
60
include/openmc/physics_mg.h
Normal file
60
include/openmc/physics_mg.h
Normal file
|
|
@ -0,0 +1,60 @@
|
|||
//! \file physics_mg.h
|
||||
//! Methods needed to perform the collision physics for multi-group mode
|
||||
|
||||
#ifndef OPENMC_PHYSICS_MG_H
|
||||
#define OPENMC_PHYSICS_MG_H
|
||||
|
||||
#include "openmc/capi.h"
|
||||
#include "openmc/particle.h"
|
||||
#include "openmc/nuclide.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//TODO: Remove energy_bin_avg and material_xs parameters when they reside on
|
||||
// the C-side this should happen after materials, physics, input, and tallies
|
||||
// are brought over
|
||||
|
||||
//! \brief samples particle behavior after a collision event.
|
||||
//! \param p Particle to operate on
|
||||
//! \param energy_bin_avg Average energy within each energy bin
|
||||
//! \param material_xs The cross section cache for the current material
|
||||
extern "C" void
|
||||
collision_mg(Particle* p, const double* energy_bin_avg,
|
||||
const MaterialMacroXS* material_xs);
|
||||
|
||||
//! \brief samples a reaction type.
|
||||
//!
|
||||
//! Note that there is special logic when suvival biasing is turned on since
|
||||
//! fission and disappearance are treated implicitly.
|
||||
//! \param p Particle to operate on
|
||||
//! \param energy_bin_avg Average energy within each energy bin
|
||||
//! \param material_xs The cross section cache for the current material
|
||||
void
|
||||
sample_reaction(Particle* p, const double* energy_bin_avg,
|
||||
const MaterialMacroXS* material_xs);
|
||||
|
||||
//! \brief Samples the scattering event
|
||||
//! \param p Particle to operate on
|
||||
//! \param energy_bin_avg Average energy within each energy bin
|
||||
void
|
||||
scatter(Particle* p, const double* energy_bin_avg);
|
||||
|
||||
//! \brief Determines the average total, prompt and delayed neutrons produced
|
||||
//! from fission and creates the appropriate bank sites.
|
||||
//! \param p Particle to operate on
|
||||
//! \param bank_array The particle bank to populate
|
||||
//! \param size_bank Number of particles currently in the bank
|
||||
//! \param bank_array_size Allocated size of the bank
|
||||
//! \param material_xs The cross section cache for the current material
|
||||
void
|
||||
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
|
||||
int64_t bank_array_size, const MaterialMacroXS* material_xs);
|
||||
|
||||
//! \brief Handles an absorption event
|
||||
//! \param p Particle to operate on
|
||||
//! \param material_xs The cross section cache for the current material
|
||||
void
|
||||
absorption(Particle* p, const MaterialMacroXS* material_xs);
|
||||
|
||||
} // namespace openmc
|
||||
#endif // OPENMC_PHYSICS_MG_H
|
||||
|
|
@ -12,8 +12,8 @@ _dll.t_percentile_c.argtypes = [c_double, c_int]
|
|||
_dll.calc_pn_c.restype = None
|
||||
_dll.calc_pn_c.argtypes = [c_int, c_double, ndpointer(c_double)]
|
||||
|
||||
_dll.evaluate_legendre_c.restype = c_double
|
||||
_dll.evaluate_legendre_c.argtypes = [c_int, POINTER(c_double), c_double]
|
||||
_dll.evaluate_legendre.restype = c_double
|
||||
_dll.evaluate_legendre.argtypes = [c_int, POINTER(c_double), c_double]
|
||||
|
||||
_dll.calc_rn_c.restype = None
|
||||
_dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
|
||||
|
|
@ -27,11 +27,11 @@ _dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
|
|||
_dll.rotate_angle_c.restype = None
|
||||
_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
|
||||
POINTER(c_double)]
|
||||
_dll.maxwell_spectrum_c.restype = c_double
|
||||
_dll.maxwell_spectrum_c.argtypes = [c_double]
|
||||
_dll.maxwell_spectrum.restype = c_double
|
||||
_dll.maxwell_spectrum.argtypes = [c_double]
|
||||
|
||||
_dll.watt_spectrum_c.restype = c_double
|
||||
_dll.watt_spectrum_c.argtypes = [c_double, c_double]
|
||||
_dll.watt_spectrum.restype = c_double
|
||||
_dll.watt_spectrum.argtypes = [c_double, c_double]
|
||||
|
||||
_dll.broaden_wmp_polynomials_c.restype = None
|
||||
_dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
|
||||
|
|
@ -100,9 +100,8 @@ def evaluate_legendre(data, x):
|
|||
"""
|
||||
|
||||
data_arr = np.array(data, dtype=np.float64)
|
||||
return _dll.evaluate_legendre_c(len(data),
|
||||
data_arr.ctypes.data_as(POINTER(c_double)),
|
||||
x)
|
||||
return _dll.evaluate_legendre(len(data),
|
||||
data_arr.ctypes.data_as(POINTER(c_double)), x)
|
||||
|
||||
|
||||
def calc_rn(n, uvw):
|
||||
|
|
@ -182,7 +181,7 @@ def calc_zn_rad(n, rho):
|
|||
zn_rad = np.zeros(num_bins, dtype=np.float64)
|
||||
_dll.calc_zn_rad_c(n, rho, zn_rad)
|
||||
return zn_rad
|
||||
|
||||
|
||||
|
||||
def rotate_angle(uvw0, mu, phi=None):
|
||||
""" Rotates direction cosines through a polar angle whose cosine is
|
||||
|
|
@ -231,7 +230,7 @@ def maxwell_spectrum(T):
|
|||
|
||||
"""
|
||||
|
||||
return _dll.maxwell_spectrum_c(T)
|
||||
return _dll.maxwell_spectrum(T)
|
||||
|
||||
|
||||
def watt_spectrum(a, b):
|
||||
|
|
@ -251,7 +250,7 @@ def watt_spectrum(a, b):
|
|||
|
||||
"""
|
||||
|
||||
return _dll.watt_spectrum_c(a, b)
|
||||
return _dll.watt_spectrum(a, b)
|
||||
|
||||
|
||||
def broaden_wmp_polynomials(E, dopp, n):
|
||||
|
|
|
|||
|
|
@ -90,7 +90,7 @@ Maxwell::Maxwell(pugi::xml_node node)
|
|||
|
||||
double Maxwell::sample() const
|
||||
{
|
||||
return maxwell_spectrum_c(theta_);
|
||||
return maxwell_spectrum(theta_);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
@ -110,7 +110,7 @@ Watt::Watt(pugi::xml_node node)
|
|||
|
||||
double Watt::sample() const
|
||||
{
|
||||
return watt_spectrum_c(a_, b_);
|
||||
return watt_spectrum(a_, b_);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
|
|
|||
|
|
@ -279,7 +279,7 @@ double MaxwellEnergy::sample(double E) const
|
|||
|
||||
while (true) {
|
||||
// Sample maxwell fission spectrum
|
||||
double E_out = maxwell_spectrum_c(theta);
|
||||
double E_out = maxwell_spectrum(theta);
|
||||
|
||||
// Accept energy based on restriction energy
|
||||
if (E_out <= E - u_) return E_out;
|
||||
|
|
@ -343,7 +343,7 @@ double WattEnergy::sample(double E) const
|
|||
|
||||
while (true) {
|
||||
// Sample energy-dependent Watt fission spectrum
|
||||
double E_out = watt_spectrum_c(a, b);
|
||||
double E_out = watt_spectrum(a, b);
|
||||
|
||||
// Accept energy based on restriction energy
|
||||
if (E_out <= E - u_) return E_out;
|
||||
|
|
|
|||
|
|
@ -2986,7 +2986,7 @@ contains
|
|||
|
||||
! Check if material is fissionable
|
||||
if (nuclides(materials(i) % nuclide(j)) % fissionable) then
|
||||
materials(i) % fissionable = .true.
|
||||
call materials(i) % set_fissionable(.true.)
|
||||
end if
|
||||
end do
|
||||
|
||||
|
|
|
|||
|
|
@ -122,6 +122,13 @@ extern "C" {
|
|||
material_map[id] = index - 1;
|
||||
}
|
||||
|
||||
bool material_fissionable(Material* mat) {return mat->fissionable;}
|
||||
|
||||
void material_set_fissionable(Material* mat, bool fissionable)
|
||||
{
|
||||
mat->fissionable = fissionable;
|
||||
}
|
||||
|
||||
void extend_materials_c(int32_t n)
|
||||
{
|
||||
materials.reserve(materials.size() + n);
|
||||
|
|
|
|||
|
|
@ -51,6 +51,20 @@ module material_header
|
|||
integer(C_INT32_T), intent(in), value :: index
|
||||
end subroutine material_set_id_c
|
||||
|
||||
function material_fissionable_c(mat_ptr) &
|
||||
bind(C, name='material_fissionable') result(fissionable)
|
||||
import C_PTR, C_BOOL
|
||||
type(C_PTR), intent(in), value :: mat_ptr
|
||||
logical(C_BOOL) :: fissionable
|
||||
end function material_fissionable_c
|
||||
|
||||
subroutine material_set_fissionable_c(mat_ptr, fissionable) &
|
||||
bind(C, name='material_set_fissionable')
|
||||
import C_PTR, C_BOOL
|
||||
type(C_PTR), intent(in), value :: mat_ptr
|
||||
logical(C_BOOL), intent(in), value :: fissionable
|
||||
end subroutine material_set_fissionable_c
|
||||
|
||||
subroutine extend_materials_c(n) bind(C)
|
||||
import C_INT32_T
|
||||
integer(C_INT32_T), intent(in), value :: n
|
||||
|
|
@ -95,7 +109,6 @@ module material_header
|
|||
character(20), allocatable :: sab_names(:) ! name of S(a,b) table
|
||||
|
||||
! Does this material contain fissionable nuclides? Is it depletable?
|
||||
logical :: fissionable = .false.
|
||||
logical :: depletable = .false.
|
||||
|
||||
! enforce isotropic scattering in lab for specific nuclides
|
||||
|
|
@ -105,6 +118,8 @@ module material_header
|
|||
contains
|
||||
procedure :: id => material_id
|
||||
procedure :: set_id => material_set_id
|
||||
procedure :: fissionable => material_fissionable
|
||||
procedure :: set_fissionable => material_set_fissionable
|
||||
procedure :: set_density => material_set_density
|
||||
procedure :: init_nuclide_index => material_init_nuclide_index
|
||||
procedure :: assign_sab_tables => material_assign_sab_tables
|
||||
|
|
@ -139,6 +154,18 @@ contains
|
|||
call material_set_id_c(this % ptr, id, index)
|
||||
end subroutine material_set_id
|
||||
|
||||
function material_fissionable(this) result(fissionable)
|
||||
class(Material), intent(in) :: this
|
||||
logical(C_BOOL) :: fissionable
|
||||
fissionable = material_fissionable_c(this % ptr)
|
||||
end function material_fissionable
|
||||
|
||||
subroutine material_set_fissionable(this, fissionable)
|
||||
class(Material),intent(in) :: this
|
||||
logical, intent(in) :: fissionable
|
||||
call material_set_fissionable_c(this % ptr, logical(fissionable, C_BOOL))
|
||||
end subroutine material_set_fissionable
|
||||
|
||||
function material_set_density(this, density) result(err)
|
||||
class(Material), intent(inout) :: this
|
||||
real(8), intent(in) :: density
|
||||
|
|
@ -685,7 +712,7 @@ contains
|
|||
integer(C_INT) :: err
|
||||
|
||||
if (index >= 1 .and. index <= size(materials)) then
|
||||
fissionable = materials(index) % fissionable
|
||||
fissionable = materials(index) % fissionable()
|
||||
err = 0
|
||||
else
|
||||
err = E_OUT_OF_BOUNDS
|
||||
|
|
|
|||
44
src/math.F90
44
src/math.F90
|
|
@ -12,10 +12,7 @@ module math
|
|||
public :: calc_rn
|
||||
public :: calc_zn
|
||||
public :: calc_zn_rad
|
||||
public :: evaluate_legendre
|
||||
public :: rotate_angle
|
||||
public :: maxwell_spectrum
|
||||
public :: watt_spectrum
|
||||
public :: faddeeva
|
||||
public :: w_derivative
|
||||
public :: broaden_wmp_polynomials
|
||||
|
|
@ -42,16 +39,6 @@ module math
|
|||
real(C_DOUBLE), intent(out) :: pnx(n + 1)
|
||||
end subroutine calc_pn
|
||||
|
||||
pure function evaluate_legendre_c_intfc(n, data, x) &
|
||||
bind(C, name='evaluate_legendre_c') result(val)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: n
|
||||
real(C_DOUBLE), intent(in) :: data(n)
|
||||
real(C_DOUBLE), value, intent(in) :: x
|
||||
real(C_DOUBLE) :: val
|
||||
end function evaluate_legendre_c_intfc
|
||||
|
||||
pure subroutine calc_rn(n, uvw, rn) bind(C, name='calc_rn_c')
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
|
|
@ -85,23 +72,6 @@ module math
|
|||
real(C_DOUBLE), optional, intent(in) :: phi
|
||||
end subroutine rotate_angle_c_intfc
|
||||
|
||||
function maxwell_spectrum(T) bind(C, name='maxwell_spectrum_c') &
|
||||
result(E_out)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
real(C_DOUBLE), value, intent(in) :: T
|
||||
real(C_DOUBLE) :: E_out
|
||||
end function maxwell_spectrum
|
||||
|
||||
function watt_spectrum(a, b) bind(C, name='watt_spectrum_c') &
|
||||
result(E_out)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
real(C_DOUBLE), value, intent(in) :: a
|
||||
real(C_DOUBLE), value, intent(in) :: b
|
||||
real(C_DOUBLE) :: E_out
|
||||
end function watt_spectrum
|
||||
|
||||
subroutine broaden_wmp_polynomials(E, dopp, n, factors) &
|
||||
bind(C, name='broaden_wmp_polynomials_c')
|
||||
use ISO_C_BINDING
|
||||
|
|
@ -157,20 +127,6 @@ module math
|
|||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! EVALUATE_LEGENDRE Find the value of f(x) given a set of Legendre coefficients
|
||||
! and the value of x
|
||||
!===============================================================================
|
||||
|
||||
pure function evaluate_legendre(data, x) result(val) bind(C)
|
||||
real(C_DOUBLE), intent(in) :: data(:)
|
||||
real(C_DOUBLE), intent(in) :: x
|
||||
real(C_DOUBLE) :: val
|
||||
|
||||
val = evaluate_legendre_c_intfc(size(data) - 1, data, x)
|
||||
|
||||
end function evaluate_legendre
|
||||
|
||||
!===============================================================================
|
||||
! ROTATE_ANGLE rotates direction cosines through a polar angle whose cosine is
|
||||
! mu and through an azimuthal angle sampled uniformly. Note that this is done
|
||||
|
|
|
|||
|
|
@ -6,7 +6,7 @@ namespace openmc {
|
|||
// Mathematical methods
|
||||
//==============================================================================
|
||||
|
||||
double normal_percentile_c(double p) {
|
||||
double normal_percentile(double p) {
|
||||
constexpr double p_low = 0.02425;
|
||||
constexpr double a[6] = {-3.969683028665376e1, 2.209460984245205e2,
|
||||
-2.759285104469687e2, 1.383577518672690e2,
|
||||
|
|
@ -79,7 +79,7 @@ double t_percentile_c(double p, int df){
|
|||
// 16 (4), pp. 1123-1132 (1987).
|
||||
double n = df;
|
||||
double k = 1. / (n - 2.);
|
||||
double z = normal_percentile_c(p);
|
||||
double z = normal_percentile(p);
|
||||
double z2 = z * z;
|
||||
t = std::sqrt(n * k) * (z + (z2 - 3.) * z * k / 4. + ((5. * z2 - 56.) * z2 +
|
||||
75.) * z * k * k / 96. + (((z2 - 27.) * 3. * z2 + 417.) * z2 - 315.) *
|
||||
|
|
@ -103,7 +103,7 @@ void calc_pn_c(int n, double x, double pnx[]) {
|
|||
}
|
||||
|
||||
|
||||
double evaluate_legendre_c(int n, const double data[], double x) {
|
||||
double evaluate_legendre(int n, const double data[], double x) {
|
||||
double pnx[n + 1];
|
||||
double val = 0.0;
|
||||
calc_pn_c(n, x, pnx);
|
||||
|
|
@ -658,7 +658,7 @@ Direction rotate_angle(Direction u, double mu, double* phi)
|
|||
}
|
||||
|
||||
|
||||
double maxwell_spectrum_c(double T) {
|
||||
double maxwell_spectrum(double T) {
|
||||
// Set the random numbers
|
||||
double r1 = prn();
|
||||
double r2 = prn();
|
||||
|
|
@ -674,8 +674,8 @@ double maxwell_spectrum_c(double T) {
|
|||
}
|
||||
|
||||
|
||||
double watt_spectrum_c(double a, double b) {
|
||||
double w = maxwell_spectrum_c(a);
|
||||
double watt_spectrum(double a, double b) {
|
||||
double w = maxwell_spectrum(a);
|
||||
double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);
|
||||
|
||||
return E_out;
|
||||
|
|
|
|||
|
|
@ -90,7 +90,8 @@ contains
|
|||
end if
|
||||
end do NUCLIDE_LOOP
|
||||
|
||||
mat % fissionable = query_fissionable_c(mat % n_nuclides, mat % nuclide)
|
||||
call mat % set_fissionable( &
|
||||
logical(query_fissionable_c(mat % n_nuclides, mat % nuclide)))
|
||||
|
||||
end do MATERIAL_LOOP
|
||||
|
||||
|
|
|
|||
|
|
@ -62,26 +62,6 @@ module mgxs_interface
|
|||
real(C_DOUBLE), intent(inout) :: nu_fiss_xs
|
||||
end subroutine calculate_xs_c
|
||||
|
||||
subroutine sample_scatter_c(i_mat, gin, gout, mu, wgt, uvw) bind(C)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: i_mat
|
||||
integer(C_INT), value, intent(in) :: gin
|
||||
integer(C_INT), intent(inout) :: gout
|
||||
real(C_DOUBLE), intent(inout) :: mu
|
||||
real(C_DOUBLE), intent(inout) :: wgt
|
||||
real(C_DOUBLE), intent(inout) :: uvw(1:3)
|
||||
end subroutine sample_scatter_c
|
||||
|
||||
subroutine sample_fission_energy_c(i_mat, gin, dg, gout) bind(C)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
integer(C_INT), value, intent(in) :: i_mat
|
||||
integer(C_INT), value, intent(in) :: gin
|
||||
integer(C_INT), intent(inout) :: dg
|
||||
integer(C_INT), intent(inout) :: gout
|
||||
end subroutine sample_fission_energy_c
|
||||
|
||||
subroutine get_name_c(index, name_len, name) bind(C)
|
||||
use ISO_C_BINDING
|
||||
implicit none
|
||||
|
|
@ -156,7 +136,7 @@ module mgxs_interface
|
|||
real(8), allocatable :: energy_bins(:)
|
||||
|
||||
! Midpoint of the energy group structure
|
||||
real(8), allocatable :: energy_bin_avg(:)
|
||||
real(C_DOUBLE), allocatable :: energy_bin_avg(:)
|
||||
|
||||
! Energy group structure with increasing energy
|
||||
real(C_DOUBLE), allocatable, target :: rev_energy_bins(:)
|
||||
|
|
|
|||
|
|
@ -97,36 +97,6 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
|
|||
|
||||
//==============================================================================
|
||||
|
||||
void
|
||||
sample_scatter_c(int i_mat, int gin, int& gout, double& mu, double& wgt,
|
||||
double uvw[3])
|
||||
{
|
||||
int gout_c = gout - 1;
|
||||
macro_xs[i_mat - 1].sample_scatter(gin - 1, gout_c, mu, wgt);
|
||||
|
||||
// adjust return value for fortran indexing
|
||||
gout = gout_c + 1;
|
||||
|
||||
// Rotate the angle
|
||||
rotate_angle_c(uvw, mu, nullptr);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
||||
void
|
||||
sample_fission_energy_c(int i_mat, int gin, int& dg, int& gout)
|
||||
{
|
||||
int dg_c = 0;
|
||||
int gout_c = 0;
|
||||
macro_xs[i_mat - 1].sample_fission_energy(gin - 1, dg_c, gout_c);
|
||||
|
||||
// adjust return values for fortran indexing
|
||||
dg = dg_c + 1;
|
||||
gout = gout_c + 1;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
|
||||
double
|
||||
get_nuclide_xs_c(int index, int xstype, int gin, int* gout, double* mu, int* dg)
|
||||
{
|
||||
|
|
|
|||
|
|
@ -1,33 +1,20 @@
|
|||
module physics_common
|
||||
|
||||
use constants
|
||||
use, intrinsic :: ISO_C_BINDING
|
||||
|
||||
use particle_header, only: Particle
|
||||
use random_lcg, only: prn
|
||||
use settings, only: weight_cutoff, weight_survive
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! RUSSIAN_ROULETTE
|
||||
! RUSSIAN_ROULETTE FROM C
|
||||
!===============================================================================
|
||||
|
||||
subroutine russian_roulette(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
if (p % wgt < weight_cutoff) then
|
||||
if (prn() < p % wgt / weight_survive) then
|
||||
p % wgt = weight_survive
|
||||
p % last_wgt = p % wgt
|
||||
else
|
||||
p % wgt = ZERO
|
||||
p % last_wgt = ZERO
|
||||
p % alive = .false.
|
||||
end if
|
||||
end if
|
||||
|
||||
end subroutine russian_roulette
|
||||
interface
|
||||
subroutine russian_roulette(p) bind(C)
|
||||
import Particle
|
||||
type(Particle), intent(inout) :: p
|
||||
end subroutine russian_roulette
|
||||
end interface
|
||||
|
||||
end module physics_common
|
||||
|
|
|
|||
26
src/physics_common.cpp
Normal file
26
src/physics_common.cpp
Normal file
|
|
@ -0,0 +1,26 @@
|
|||
#include "openmc/physics_common.h"
|
||||
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
//==============================================================================
|
||||
// RUSSIAN_ROULETTE
|
||||
//==============================================================================
|
||||
|
||||
void russian_roulette(Particle* p)
|
||||
{
|
||||
if (p->wgt < settings::weight_cutoff) {
|
||||
if (prn() < p->wgt / settings::weight_survive) {
|
||||
p->wgt = settings::weight_survive;
|
||||
p->last_wgt = p->wgt;
|
||||
} else {
|
||||
p->wgt = 0.;
|
||||
p->last_wgt = 0.;
|
||||
p->alive = false;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} //namespace openmc
|
||||
|
|
@ -1,275 +0,0 @@
|
|||
module physics_mg
|
||||
! This module contains the multi-group specific physics routines so as to not
|
||||
! hinder performance of the CE versions with multiple if-thens.
|
||||
|
||||
use bank_header
|
||||
use constants
|
||||
use error, only: fatal_error, warning, write_message
|
||||
use material_header, only: Material, materials
|
||||
use math, only: rotate_angle
|
||||
use mgxs_interface
|
||||
use message_passing
|
||||
use nuclide_header, only: material_xs
|
||||
use particle_header
|
||||
use physics_common
|
||||
use random_lcg, only: prn
|
||||
use settings
|
||||
use simulation_header
|
||||
use string, only: to_str
|
||||
use tally_header
|
||||
|
||||
implicit none
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
! COLLISION_MG samples a nuclide and reaction and then calls the appropriate
|
||||
! routine for that reaction
|
||||
!===============================================================================
|
||||
|
||||
subroutine collision_mg(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
! Add to collision counter for particle
|
||||
p % n_collision = p % n_collision + 1
|
||||
|
||||
! Sample nuclide/reaction for the material the particle is in
|
||||
call sample_reaction(p)
|
||||
|
||||
! Display information about collision
|
||||
if (verbosity >= 10 .or. trace) then
|
||||
call write_message(" " // "Energy Group = " // trim(to_str(p % g)))
|
||||
end if
|
||||
|
||||
end subroutine collision_mg
|
||||
|
||||
!===============================================================================
|
||||
! SAMPLE_REACTION samples a nuclide based on the macroscopic cross sections for
|
||||
! each nuclide within a material and then samples a reaction for that nuclide
|
||||
! and calls the appropriate routine to process the physics. Note that there is
|
||||
! special logic when suvival biasing is turned on since fission and
|
||||
! disappearance are treated implicitly.
|
||||
!===============================================================================
|
||||
|
||||
subroutine sample_reaction(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
type(Material), pointer :: mat
|
||||
|
||||
mat => materials(p % material)
|
||||
|
||||
! Create fission bank sites. Note that while a fission reaction is sampled,
|
||||
! it never actually "happens", i.e. the weight of the particle does not
|
||||
! change when sampling fission sites. The following block handles all
|
||||
! absorption (including fission)
|
||||
|
||||
if (mat % fissionable) then
|
||||
if (run_mode == MODE_EIGENVALUE) then
|
||||
call create_fission_sites(p, fission_bank, n_bank)
|
||||
elseif (run_mode == MODE_FIXEDSOURCE .and. create_fission_neutrons) then
|
||||
call create_fission_sites(p, p % secondary_bank, p % n_secondary)
|
||||
end if
|
||||
end if
|
||||
|
||||
! If survival biasing is being used, the following subroutine adjusts the
|
||||
! weight of the particle. Otherwise, it checks to see if absorption occurs
|
||||
|
||||
if (material_xs % absorption > ZERO) then
|
||||
call absorption(p)
|
||||
else
|
||||
p % absorb_wgt = ZERO
|
||||
end if
|
||||
if (.not. p % alive) return
|
||||
|
||||
! Sample a scattering reaction and determine the secondary energy of the
|
||||
! exiting neutron
|
||||
call scatter(p)
|
||||
|
||||
! Play russian roulette if survival biasing is turned on
|
||||
if (survival_biasing) then
|
||||
call russian_roulette(p)
|
||||
if (.not. p % alive) return
|
||||
end if
|
||||
|
||||
end subroutine sample_reaction
|
||||
|
||||
!===============================================================================
|
||||
! ABSORPTION
|
||||
!===============================================================================
|
||||
|
||||
subroutine absorption(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
if (survival_biasing) then
|
||||
! Determine weight absorbed in survival biasing
|
||||
p % absorb_wgt = (p % wgt * &
|
||||
material_xs % absorption / material_xs % total)
|
||||
|
||||
! Adjust weight of particle by probability of absorption
|
||||
p % wgt = p % wgt - p % absorb_wgt
|
||||
p % last_wgt = p % wgt
|
||||
|
||||
! Score implicit absorption estimate of keff
|
||||
!$omp atomic
|
||||
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
|
||||
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % absorb_wgt * &
|
||||
material_xs % nu_fission / material_xs % absorption
|
||||
else
|
||||
! See if disappearance reaction happens
|
||||
if (material_xs % absorption > prn() * material_xs % total) then
|
||||
! Score absorption estimate of keff
|
||||
!$omp atomic
|
||||
global_tallies(RESULT_VALUE, K_ABSORPTION) = &
|
||||
global_tallies(RESULT_VALUE, K_ABSORPTION) + p % wgt * &
|
||||
material_xs % nu_fission / material_xs % absorption
|
||||
|
||||
p % alive = .false.
|
||||
p % event = EVENT_ABSORB
|
||||
end if
|
||||
end if
|
||||
|
||||
end subroutine absorption
|
||||
|
||||
!===============================================================================
|
||||
! SCATTER
|
||||
!===============================================================================
|
||||
|
||||
subroutine scatter(p)
|
||||
|
||||
type(Particle), intent(inout) :: p
|
||||
|
||||
call sample_scatter_c(p % material, p % last_g, p % g, p % mu, &
|
||||
p % wgt, p % coord(1) % uvw)
|
||||
|
||||
! Update energy value for downstream compatability (in tallying)
|
||||
p % E = energy_bin_avg(p % g)
|
||||
|
||||
! Set event component
|
||||
p % event = EVENT_SCATTER
|
||||
|
||||
end subroutine scatter
|
||||
|
||||
!===============================================================================
|
||||
! CREATE_FISSION_SITES determines the average total, prompt, and delayed
|
||||
! neutrons produced from fission and creates appropriate bank sites.
|
||||
!===============================================================================
|
||||
|
||||
subroutine create_fission_sites(p, bank_array, size_bank)
|
||||
type(Particle), intent(inout) :: p
|
||||
type(Bank), intent(inout) :: bank_array(:)
|
||||
integer(8), intent(inout) :: size_bank
|
||||
|
||||
integer :: nu_d(MAX_DELAYED_GROUPS) ! number of delayed neutrons born
|
||||
integer :: i ! loop index
|
||||
integer :: dg ! delayed group
|
||||
integer :: gout ! group out
|
||||
integer :: nu ! actual number of neutrons produced
|
||||
real(8) :: nu_t ! total nu
|
||||
real(8) :: mu ! fission neutron angular cosine
|
||||
real(8) :: phi ! fission neutron azimuthal angle
|
||||
real(8) :: weight ! weight adjustment for ufs method
|
||||
|
||||
interface
|
||||
function ufs_get_weight(p) result(weight) bind(C)
|
||||
import Particle, C_DOUBLE
|
||||
type(Particle), intent(in) :: p
|
||||
real(C_DOUBLE) :: WEIGHT
|
||||
end function
|
||||
end interface
|
||||
|
||||
! TODO: Heat generation from fission
|
||||
|
||||
! If uniform fission source weighting is turned on, we increase of decrease
|
||||
! the expected number of fission sites produced
|
||||
|
||||
if (ufs) then
|
||||
weight = ufs_get_weight(p)
|
||||
else
|
||||
weight = ONE
|
||||
end if
|
||||
|
||||
! Determine expected number of neutrons produced
|
||||
nu_t = p % wgt / keff * weight * &
|
||||
material_xs % nu_fission / material_xs % total
|
||||
|
||||
! Sample number of neutrons produced
|
||||
if (prn() > nu_t - int(nu_t)) then
|
||||
nu = int(nu_t)
|
||||
else
|
||||
nu = int(nu_t) + 1
|
||||
end if
|
||||
|
||||
! Check for bank size getting hit. For fixed source calculations, this is a
|
||||
! fatal error. For eigenvalue calculations, it just means that k-effective
|
||||
! was too high for a single batch.
|
||||
if (size_bank + nu > size(bank_array)) then
|
||||
if (run_mode == MODE_FIXEDSOURCE) then
|
||||
call fatal_error("Secondary particle bank size limit reached. If you &
|
||||
&are running a subcritical multiplication problem, k-effective &
|
||||
&may be too close to one.")
|
||||
else
|
||||
if (master) call warning("Maximum number of sites in fission bank &
|
||||
&reached. This can result in irreproducible results using different &
|
||||
&numbers of processes/threads.")
|
||||
end if
|
||||
end if
|
||||
|
||||
! Bank source neutrons
|
||||
if (nu == 0 .or. size_bank == size(bank_array)) return
|
||||
|
||||
! Initialize counter of delayed neutrons encountered for each delayed group
|
||||
! to zero.
|
||||
nu_d(:) = 0
|
||||
|
||||
p % fission = .true. ! Fission neutrons will be banked
|
||||
do i = int(size_bank,4) + 1, int(min(size_bank + nu, int(size(bank_array),8)),4)
|
||||
! Bank source neutrons by copying particle data
|
||||
bank_array(i) % xyz = p % coord(1) % xyz
|
||||
|
||||
! Set particle as neutron
|
||||
bank_array(i) % particle = NEUTRON
|
||||
|
||||
! Set weight of fission bank site
|
||||
bank_array(i) % wgt = ONE/weight
|
||||
|
||||
! Sample cosine of angle -- fission neutrons are treated as being emitted
|
||||
! isotropically.
|
||||
mu = TWO * prn() - ONE
|
||||
|
||||
! Sample azimuthal angle uniformly in [0,2*pi)
|
||||
phi = TWO * PI * prn()
|
||||
bank_array(i) % uvw(1) = mu
|
||||
bank_array(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
|
||||
bank_array(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
|
||||
|
||||
! Sample secondary energy distribution for fission reaction and set energy
|
||||
! in fission bank
|
||||
call sample_fission_energy_c(p % material, p % g, dg, gout)
|
||||
|
||||
bank_array(i) % E = real(gout, 8)
|
||||
bank_array(i) % delayed_group = dg
|
||||
|
||||
! Set delayed group on particle too
|
||||
p % delayed_group = dg
|
||||
|
||||
! Increment the number of neutrons born delayed
|
||||
if (p % delayed_group > 0) then
|
||||
nu_d(p % delayed_group) = nu_d(p % delayed_group) + 1
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
! increment number of bank sites
|
||||
size_bank = min(size_bank + nu, int(size(bank_array),8))
|
||||
|
||||
! Store total weight banked for analog fission tallies
|
||||
p % n_bank = nu
|
||||
p % wgt_bank = nu/weight
|
||||
p % n_delayed_bank(:) = nu_d(:)
|
||||
|
||||
end subroutine create_fission_sites
|
||||
|
||||
end module physics_mg
|
||||
233
src/physics_mg.cpp
Normal file
233
src/physics_mg.cpp
Normal file
|
|
@ -0,0 +1,233 @@
|
|||
#include "openmc/physics_mg.h"
|
||||
|
||||
#include <stdexcept>
|
||||
#include <sstream>
|
||||
|
||||
#include "xtensor/xarray.hpp"
|
||||
|
||||
#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/mgxs_interface.h"
|
||||
#include "openmc/physics_common.h"
|
||||
#include "openmc/random_lcg.h"
|
||||
#include "openmc/settings.h"
|
||||
#include "openmc/simulation.h"
|
||||
|
||||
namespace openmc {
|
||||
|
||||
void
|
||||
collision_mg(Particle* p, const double* energy_bin_avg,
|
||||
const MaterialMacroXS* material_xs)
|
||||
{
|
||||
// Add to the collision counter for the particle
|
||||
p->n_collision++;
|
||||
|
||||
// Sample the reaction type
|
||||
sample_reaction(p, energy_bin_avg, material_xs);
|
||||
|
||||
// Display information about collision
|
||||
if ((settings::verbosity >= 10) || (openmc_trace)) {
|
||||
std::stringstream msg;
|
||||
msg << " Energy Group = " << p->g;
|
||||
write_message(msg, 1);
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
sample_reaction(Particle* p, const double* energy_bin_avg,
|
||||
const MaterialMacroXS* material_xs)
|
||||
{
|
||||
// Create fission bank sites. Note that while a fission reaction is sampled,
|
||||
// it never actually "happens", i.e. the weight of the particle does not
|
||||
// change when sampling fission sites. The following block handles all
|
||||
// absorption (including fission)
|
||||
|
||||
if (materials[p->material - 1]->fissionable) {
|
||||
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
|
||||
Bank* result_bank;
|
||||
int64_t result_bank_size;
|
||||
// Get pointer to fission bank from Fortran side
|
||||
openmc_fission_bank(&result_bank, &result_bank_size);
|
||||
create_fission_sites(p, result_bank, &n_bank, result_bank_size,
|
||||
material_xs);
|
||||
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
|
||||
(settings::create_fission_neutrons)) {
|
||||
create_fission_sites(p, p->secondary_bank, &(p->n_secondary),
|
||||
MAX_SECONDARY, material_xs);
|
||||
}
|
||||
}
|
||||
|
||||
// If survival biasing is being used, the following subroutine adjusts the
|
||||
// weight of the particle. Otherwise, it checks to see if absorption occurs.
|
||||
if (material_xs->absorption > 0.) {
|
||||
absorption(p, material_xs);
|
||||
} else {
|
||||
p->absorb_wgt = 0.;
|
||||
}
|
||||
if (!p->alive) return;
|
||||
|
||||
// Sample a scattering event to determine the energy of the exiting neutron
|
||||
scatter(p, energy_bin_avg);
|
||||
|
||||
// Play Russian roulette if survival biasing is turned on
|
||||
if (settings::survival_biasing) {
|
||||
russian_roulette(p);
|
||||
if (!p->alive) return;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
scatter(Particle* p, const double* energy_bin_avg)
|
||||
{
|
||||
// Adjust indices for Fortran to C++ indexing
|
||||
// TODO: Remove when no longer needed
|
||||
int gin = p->last_g - 1;
|
||||
int gout = p->g - 1;
|
||||
int i_mat = p->material - 1;
|
||||
macro_xs[i_mat].sample_scatter(gin, gout, p->mu, p->wgt);
|
||||
|
||||
// Adjust return value for fortran indexing
|
||||
// TODO: Remove when no longer needed
|
||||
p->g = gout + 1;
|
||||
|
||||
// Rotate the angle
|
||||
rotate_angle_c(p->coord[0].uvw, p->mu, nullptr);
|
||||
|
||||
// Update energy value for downstream compatability (in tallying)
|
||||
p->E = energy_bin_avg[gout];
|
||||
|
||||
// Set event component
|
||||
p->event = EVENT_SCATTER;
|
||||
}
|
||||
|
||||
void
|
||||
create_fission_sites(Particle* p, Bank* bank_array, int64_t* size_bank,
|
||||
int64_t bank_array_size, const MaterialMacroXS* material_xs)
|
||||
{
|
||||
// TODO: Heat generation from fission
|
||||
|
||||
// If uniform fission source weighting is turned on, we increase or decrease
|
||||
// the expected number of fission sites produced
|
||||
double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
|
||||
|
||||
// Determine the expected number of neutrons produced
|
||||
double nu_t = p->wgt / openmc_keff * weight * material_xs->nu_fission /
|
||||
material_xs->total;
|
||||
|
||||
// Sample the number of neutrons produced
|
||||
int nu = static_cast<int>(nu_t);
|
||||
if (prn() <= (nu_t - int(nu_t))) {
|
||||
nu++;
|
||||
}
|
||||
|
||||
// Check for the bank size getting hit. For fixed source calculations, this
|
||||
// is a fatal error; for eigenvalue calculations, it just means that k-eff
|
||||
// was too high for a single batch.
|
||||
if (*size_bank + nu > bank_array_size) {
|
||||
if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
|
||||
throw std::runtime_error{"Secondary particle bank size limit reached."
|
||||
" If you are running a subcritical multiplication problem,"
|
||||
" k-effective may be too close to one."};
|
||||
} else {
|
||||
if (mpi::master) {
|
||||
std::stringstream msg;
|
||||
msg << "Maximum number of sites in fission bank reached. This can"
|
||||
" result in irreproducible results using different numbers of"
|
||||
" processes/threads.";
|
||||
warning(msg);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Begin banking the source neutrons
|
||||
// First, if our bank is full then don't continue
|
||||
if ((nu == 0) || (*size_bank == bank_array_size)) return;
|
||||
|
||||
// Initialize the counter of delayed neutrons encountered for each delayed
|
||||
// group.
|
||||
double nu_d[MAX_DELAYED_GROUPS] = {0.};
|
||||
|
||||
p->fission = true;
|
||||
for (size_t i = static_cast<size_t>(*size_bank);
|
||||
i < static_cast<size_t>(std::min(*size_bank + nu, bank_array_size)); i++) {
|
||||
// Bank source neutrons by copying the particle data
|
||||
bank_array[i].xyz[0] = p->coord[0].xyz[0];
|
||||
bank_array[i].xyz[1] = p->coord[0].xyz[1];
|
||||
bank_array[i].xyz[2] = p->coord[0].xyz[2];
|
||||
|
||||
// Set that the bank particle is a neutron
|
||||
bank_array[i].particle = static_cast<int>(ParticleType::neutron);
|
||||
|
||||
// Set the weight of the fission bank site
|
||||
bank_array[i].wgt = 1. / weight;
|
||||
|
||||
// Sample the cosine of the angle, assuming fission neutrons are emitted
|
||||
// isotropically
|
||||
double mu = 2. * prn() - 1.;
|
||||
|
||||
// Sample the azimuthal angle uniformly in [0, 2.pi)
|
||||
double phi = 2. * PI * prn();
|
||||
bank_array[i].uvw[0] = mu;
|
||||
bank_array[i].uvw[1] = std::sqrt(1. - mu * mu) * std::cos(phi);
|
||||
bank_array[i].uvw[2] = std::sqrt(1. - mu * mu) * std::sin(phi);
|
||||
|
||||
// Sample secondary energy distribution for the fission reaction and set
|
||||
// the energy in the fission bank
|
||||
int dg;
|
||||
int gout;
|
||||
macro_xs[p->material - 1].sample_fission_energy(p->g - 1, dg, gout);
|
||||
bank_array[i].E = static_cast<double>(gout + 1);
|
||||
bank_array[i].delayed_group = dg + 1;
|
||||
|
||||
// Set the delayed group on the particle as well
|
||||
p->delayed_group = dg + 1;
|
||||
|
||||
// Increment the number of neutrons born delayed
|
||||
if (p->delayed_group > 0) {
|
||||
nu_d[dg]++;
|
||||
}
|
||||
}
|
||||
|
||||
// Increment number of bank sites
|
||||
*size_bank = std::min(*size_bank + nu, bank_array_size);
|
||||
|
||||
// Store the total weight banked for analog fission tallies
|
||||
p->n_bank = nu;
|
||||
p->wgt_bank = nu / weight;
|
||||
for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
|
||||
p->n_delayed_bank[d] = nu_d[d];
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
absorption(Particle* p, const MaterialMacroXS* material_xs)
|
||||
{
|
||||
if (settings::survival_biasing) {
|
||||
// Determine weight absorbed in survival biasing
|
||||
p->absorb_wgt = p->wgt * material_xs->absorption / material_xs->total;
|
||||
|
||||
// Adjust weight of particle by the probability of absorption
|
||||
p->wgt -= p->absorb_wgt;
|
||||
p->last_wgt = p->wgt;
|
||||
|
||||
// Score implicit absorpion estimate of keff
|
||||
#pragma omp atomic
|
||||
global_tally_absorption += p->absorb_wgt * material_xs->nu_fission /
|
||||
material_xs->absorption;
|
||||
} else {
|
||||
if (material_xs->absorption > prn() * material_xs->total) {
|
||||
#pragma omp atomic
|
||||
global_tally_absorption += p->wgt * material_xs->nu_fission /
|
||||
material_xs->absorption;
|
||||
p->alive = false;
|
||||
p->event = EVENT_ABSORB;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
} //namespace openmc
|
||||
|
|
@ -316,8 +316,8 @@ ScattDataLegendre::update_max_val()
|
|||
}
|
||||
|
||||
// Calculate probability
|
||||
double f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
|
||||
dist[gin][i_gout].data(), mu);
|
||||
double f = evaluate_legendre(dist[gin][i_gout].size() - 1,
|
||||
dist[gin][i_gout].data(), mu);
|
||||
|
||||
// if this is a new maximum, store it
|
||||
if (f > max_val[gin][i_gout]) max_val[gin][i_gout] = f;
|
||||
|
|
@ -339,8 +339,8 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
|
|||
f = 0.;
|
||||
} else {
|
||||
int i_gout = gout - gmin[gin];
|
||||
f = evaluate_legendre_c(dist[gin][i_gout].size() - 1,
|
||||
dist[gin][i_gout].data(), mu);
|
||||
f = evaluate_legendre(dist[gin][i_gout].size() - 1,
|
||||
dist[gin][i_gout].data(), mu);
|
||||
}
|
||||
return f;
|
||||
}
|
||||
|
|
@ -881,8 +881,8 @@ convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab,
|
|||
tab.fmu[gin][i_gout].resize(n_mu);
|
||||
for (int imu = 0; imu < n_mu; imu++) {
|
||||
tab.fmu[gin][i_gout][imu] =
|
||||
evaluate_legendre_c(leg.dist[gin][i_gout].size() - 1,
|
||||
leg.dist[gin][i_gout].data(), tab.mu[imu]);
|
||||
evaluate_legendre(leg.dist[gin][i_gout].size() - 1,
|
||||
leg.dist[gin][i_gout].data(), tab.mu[imu]);
|
||||
}
|
||||
|
||||
// Ensure positivity
|
||||
|
|
|
|||
|
|
@ -32,13 +32,13 @@ void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu) const
|
|||
double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
|
||||
|
||||
// x is essentially a Maxwellian distribution
|
||||
double x = maxwell_spectrum_c(1.0);
|
||||
double x = maxwell_spectrum(1.0);
|
||||
|
||||
double y;
|
||||
double r1, r2, r3, r4, r5, r6;
|
||||
switch (n_bodies_) {
|
||||
case 3:
|
||||
y = maxwell_spectrum_c(1.0);
|
||||
y = maxwell_spectrum(1.0);
|
||||
break;
|
||||
case 4:
|
||||
r1 = prn();
|
||||
|
|
|
|||
|
|
@ -137,7 +137,7 @@ module tally_header
|
|||
! invalidation. Thus, we use threadprivate variables to accumulate global
|
||||
! tallies and then reduce at the end of a generation.
|
||||
real(C_DOUBLE), public :: global_tally_collision = ZERO
|
||||
real(C_DOUBLE), public :: global_tally_absorption = ZERO
|
||||
real(C_DOUBLE), public, bind(C) :: global_tally_absorption = ZERO
|
||||
real(C_DOUBLE), public :: global_tally_tracklength = ZERO
|
||||
real(C_DOUBLE), public :: global_tally_leakage = ZERO
|
||||
!$omp threadprivate(global_tally_collision, global_tally_absorption, &
|
||||
|
|
|
|||
|
|
@ -17,7 +17,6 @@ module tracking
|
|||
use nuclide_header
|
||||
use particle_header
|
||||
use physics, only: collision
|
||||
use physics_mg, only: collision_mg
|
||||
use random_lcg, only: prn, prn_set_stream
|
||||
use settings
|
||||
use simulation_header
|
||||
|
|
@ -33,6 +32,16 @@ module tracking
|
|||
|
||||
implicit none
|
||||
|
||||
interface
|
||||
subroutine collision_mg(p, energy_bin_avg, material_xs) bind(C)
|
||||
import Particle, C_DOUBLE, MaterialMacroXS
|
||||
type(Particle), intent(inout) :: p
|
||||
real(C_DOUBLE), intent(in) :: energy_bin_avg(*)
|
||||
type(MaterialMacroXS), intent(in) :: material_xs
|
||||
end subroutine collision_mg
|
||||
|
||||
end interface
|
||||
|
||||
contains
|
||||
|
||||
!===============================================================================
|
||||
|
|
@ -226,7 +235,7 @@ contains
|
|||
if (run_CE) then
|
||||
call collision(p)
|
||||
else
|
||||
call collision_mg(p)
|
||||
call collision_mg(p, energy_bin_avg, material_xs)
|
||||
end if
|
||||
|
||||
! Score collision estimator tallies -- this is done after a collision
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue