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

Merge remote-tracking branch 'upstream/develop' into cpp_geometry

Browse source
This commit is contained in:
Sterling Harper 2018-08-27 20:03:02 -04:00
commit 452de390ba
37 changed files with 792 additions and 1252 deletions
Download patch file Download diff file Expand all files Collapse all files

6
CMakeLists.txt
View file

@ -299,8 +299,6 @@ add_library(libopenmc SHARED
src/cmfd_solver.F90
src/constants.F90
src/dict_header.F90
src/distribution_multivariate.F90
src/distribution_univariate.F90
src/eigenvalue.F90
src/endf.F90
src/endf_header.F90
@ -341,8 +339,6 @@ add_library(libopenmc SHARED
src/settings.F90
src/simulation_header.F90
src/simulation.F90
src/source.F90
src/source_header.F90
src/state_point.F90
src/stl_vector.F90
src/string.F90
@ -401,6 +397,7 @@ add_library(libopenmc SHARED
src/message_passing.cpp
src/mgxs.cpp
src/mgxs_interface.cpp
src/nuclide.cpp
src/output.cpp
src/particle.cpp
src/plot.cpp
@ -416,6 +413,7 @@ add_library(libopenmc SHARED
src/scattdata.cpp
src/settings.cpp
src/simulation.cpp
src/source.cpp
src/state_point.cpp
src/string_functions.cpp
src/summary.cpp

3
include/openmc/capi.h
View file

@ -30,7 +30,6 @@ extern "C" {
int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_sources(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_filter_get_id(int32_t index, int32_t* id);
int openmc_filter_get_type(int32_t index, char* type);
@ -57,6 +56,7 @@ extern "C" {
int openmc_material_add_nuclide(int32_t index, const char name[], double density);
int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n);
int openmc_material_get_id(int32_t index, int32_t* id);
int openmc_material_get_fissionable(int32_t index, bool* fissionable);
int openmc_material_get_volume(int32_t index, double* volume);
int openmc_material_set_density(int32_t index, double density);
int openmc_material_set_densities(int32_t index, int n, const char** name, const double* density);
@ -84,7 +84,6 @@ extern "C" {
int openmc_simulation_finalize();
int openmc_simulation_init();
int openmc_source_bank(struct Bank** ptr, int64_t* n);
int openmc_source_set_strength(int32_t index, double strength);
int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max);
int openmc_spatial_legendre_filter_set_order(int32_t index, int order);

4
include/openmc/distribution.h
View file

@ -36,6 +36,10 @@ public:
//! Sample a value from the distribution
//! \return Sampled value
double sample() const;
// Properties
const std::vector<double>& x() const { return x_; }
const std::vector<double>& p() const { return p_; }
private:
std::vector<double> x_; //!< Possible outcomes
std::vector<double> p_; //!< Probability of each outcome

11
include/openmc/distribution_multi.h
View file

@ -3,6 +3,8 @@
#include <memory>
#include "pugixml.hpp"
#include "openmc/distribution.h"
#include "openmc/position.h"
@ -16,14 +18,15 @@ namespace openmc {
class UnitSphereDistribution {
public:
UnitSphereDistribution() { };
explicit UnitSphereDistribution(Direction u) : u_ref{u} { };
explicit UnitSphereDistribution(Direction u) : u_ref_{u} { };
explicit UnitSphereDistribution(pugi::xml_node node);
virtual ~UnitSphereDistribution() = default;
//! Sample a direction from the distribution
//! \return Direction sampled
virtual Direction sample() const = 0;
Direction u_ref {0.0, 0.0, 1.0}; //!< reference direction
Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
};
//==============================================================================
@ -33,6 +36,7 @@ public:
class PolarAzimuthal : public UnitSphereDistribution {
public:
PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi);
explicit PolarAzimuthal(pugi::xml_node node);
//! Sample a direction from the distribution
//! \return Direction sampled
@ -62,12 +66,15 @@ public:
class Monodirectional : public UnitSphereDistribution {
public:
Monodirectional(Direction u) : UnitSphereDistribution{u} { };
explicit Monodirectional(pugi::xml_node node) : UnitSphereDistribution{node} { };
//! Sample a direction from the distribution
//! \return Sampled direction
Direction sample() const;
};
using UPtrAngle = std::unique_ptr<UnitSphereDistribution>;
} // namespace openmc
#endif // DISTRIBUTION_MULTI_H

13
include/openmc/distribution_spatial.h
View file

@ -1,4 +1,4 @@
#ifndef OPENMC_DISTRIBTUION_SPATIAL_H
#ifndef OPENMC_DISTRIBUTION_SPATIAL_H
#define OPENMC_DISTRIBUTION_SPATIAL_H
#include "pugixml.hpp"
@ -43,15 +43,18 @@ private:
class SpatialBox : public SpatialDistribution {
public:
explicit SpatialBox(pugi::xml_node node);
explicit SpatialBox(pugi::xml_node node, bool fission=false);
//! Sample a position from the distribution
//! \return Sampled position
Position sample() const;
// Properties
bool only_fissionable() const { return only_fissionable_; }
private:
Position lower_left_; //!< Lower-left coordinates of box
Position upper_right_; //!< Upper-right coordinates of box
bool only_fissionable {false}; //!< Only accept sites in fissionable region?
bool only_fissionable_ {false}; //!< Only accept sites in fissionable region?
};
//==============================================================================
@ -60,6 +63,8 @@ private:
class SpatialPoint : public SpatialDistribution {
public:
SpatialPoint() : r_{} { };
SpatialPoint(Position r) : r_{r} { };
explicit SpatialPoint(pugi::xml_node node);
//! Sample a position from the distribution
@ -69,6 +74,8 @@ private:
Position r_; //!< Single position at which sites are generated
};
using UPtrSpace = std::unique_ptr<SpatialDistribution>;
} // namespace openmc
#endif // OPENMC_DISTRIBUTION_SPATIAL_H

20
include/openmc/file_utils.h Normal file
View file

@ -0,0 +1,20 @@
#ifndef OPENMC_FILE_UTILS_H
#define OPENMC_FILE_UTILS_H
#include <fstream> // for ifstream
#include <string>
namespace openmc {
//! Determine if a file exists
//! \param[in] filename Path to file
//! \return Whether file exists
inline bool file_exists(const std::string& filename)
{
std::ifstream s {filename};
return s.good();
}
} // namespace openmc
#endif // OPENMC_FILE_UTILS_H

1
include/openmc/mgxs_interface.h
View file

@ -16,6 +16,7 @@ namespace openmc {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
//==============================================================================
// Mgxs data loading interface methods

14
include/openmc/nuclide.h
View file

@ -1,10 +1,24 @@
//! \file nuclide.h
//! \brief Nuclide type and other associated types/data
#ifndef OPENMC_NUCLIDE_H
#define OPENMC_NUCLIDE_H
#include <array>
#include "openmc/constants.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
// Minimum/maximum transport energy for each particle type. Order corresponds to
// that of the ParticleType enum
extern std::array<double, 2> energy_min;
extern std::array<double, 2> energy_max;
//===============================================================================
//! Cached microscopic cross sections for a particular nuclide at the current
//! energy

2
include/openmc/particle.h
View file

@ -35,7 +35,7 @@ constexpr double REL_MAX_LOST_PARTICLES {1.0e-6};
//! Particle types
enum class ParticleType {
neutron, photon, electron, positron
neutron = 1, photon = 2, electron = 3, positron = 4
};
extern "C" {

4
include/openmc/settings.h
View file

@ -18,10 +18,12 @@ namespace openmc {
// Defined on Fortran side
extern "C" bool openmc_check_overlaps;
extern "C" bool openmc_particle_restart_run;
extern "C" bool openmc_photon_transport;
extern "C" bool openmc_restart_run;
extern "C" bool openmc_run_CE;
extern "C" int openmc_verbosity;
extern "C" bool openmc_write_all_tracks;
extern "C" double temperature_default;
extern "C" bool openmc_write_initial_source;
// Defined in .cpp
// TODO: Make strings instead of char* once Fortran is gone

25
include/openmc/simulation.h
View file

@ -1,14 +1,39 @@
//! \file simulation.h
//! \brief Variables/functions related to a running simulation
#ifndef OPENMC_SIMULATION_H
#define OPENMC_SIMULATION_H
#include <cstdint>
#include <vector>
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
extern "C" int openmc_current_batch;
extern "C" int openmc_current_gen;
extern "C" int64_t openmc_current_work;
extern "C" int openmc_n_lost_particles;
extern "C" int openmc_total_gen;
extern "C" bool openmc_trace;
extern std::vector<int64_t> work_index;
#pragma omp threadprivate(openmc_current_work, openmc_trace)
//==============================================================================
// Functions
//==============================================================================
//! Initialize simulation
extern "C" void openmc_simulation_init_c();
//! Determine number of particles to transport per process
void calculate_work();
} // namespace openmc
#endif // OPENMC_SIMULATION_H

63
include/openmc/source.h Normal file
View file

@ -0,0 +1,63 @@
//! \file source.h
//! \brief External source distributions
#ifndef OPENMC_SOURCE_H
#define OPENMC_SOURCE_H
#include <memory>
#include <vector>
#include "pugixml.hpp"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/capi.h"
#include "openmc/particle.h"
namespace openmc {
//==============================================================================
//! External source distribution
//==============================================================================
class SourceDistribution {
public:
// Constructors
SourceDistribution(UPtrSpace space, UPtrAngle angle, UPtrDist energy);
explicit SourceDistribution(pugi::xml_node node);
//! Sample from the external source distribution
//! \return Sampled site
Bank sample() const;
// Properties
double strength() const { return strength_; }
private:
ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
double strength_ {1.0}; //!< Source strength
UPtrSpace space_; //!< Spatial distribution
UPtrAngle angle_; //!< Angular distribution
UPtrDist energy_; //!< Energy distribution
};
//==============================================================================
// Global variables
//==============================================================================
extern std::vector<SourceDistribution> external_sources;
//==============================================================================
// Functions
//==============================================================================
//! Initialize source bank from file/distribution
extern "C" void initialize_source();
//! Sample a site from all external source distributions in proportion to their
//! source strength
//! \return Sampled source site
extern "C" Bank sample_external_source();
} // namespace openmc
#endif // OPENMC_SOURCE_H

13
openmc/filter_expansion.py
View file

@ -19,7 +19,7 @@ class ExpansionFilter(Filter):
if type(self) is not type(other):
return False
else:
return self.bins == other.bins
return hash(self) == hash(other)
@property
def order(self):
@ -325,8 +325,8 @@ class ZernikeFilter(ExpansionFilter):
This filter allows scores to be multiplied by Zernike polynomials of the
particle's position normalized to a given unit circle, up to a
user-specified order. The standard Zernike polynomials follow the definition by
Born and Wolf, *Principles of Optics* and are defined as
user-specified order. The standard Zernike polynomials follow the
definition by Born and Wolf, *Principles of Optics* and are defined as
.. math::
Z_n^m(\rho, \theta) = R_n^m(\rho) \cos (m\theta), \quad m > 0
@ -342,7 +342,7 @@ class ZernikeFilter(ExpansionFilter):
\frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where
is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where
:math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise.
Specifying a filter with order N tallies moments for all :math:`n` from 0
@ -390,6 +390,9 @@ class ZernikeFilter(ExpansionFilter):
def __hash__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tX', self.x)
string += '{: <16}=\t{}\n'.format('\tY', self.y)
string += '{: <16}=\t{}\n'.format('\tR', self.r)
return hash(string)
def __repr__(self):
@ -485,7 +488,7 @@ class ZernikeRadialFilter(ZernikeFilter):
is :math:`\frac{\pi}{n+1}`.
If there is only radial dependency, the polynomials are integrated over
the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta)
the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta)
= R_n^{0}(\rho)`. Note that :math:`n` could only be even orders.
Therefore, for a radial Zernike polynomials up to order of :math:`n`,
there are :math:`\frac{n}{2} + 1` terms in total. The indexing is from the

10
openmc/polynomial.py
View file

@ -1,9 +1,10 @@
import numpy as np
import openmc
import openmc.capi as capi
from collections.abc import Iterable
def legendre_from_expcoef(coef, domain= (-1,1)):
def legendre_from_expcoef(coef, domain=(-1, 1)):
"""Return a Legendre series object based on expansion coefficients.
Given a list of coefficients from FET tally and a array of down, return
@ -73,5 +74,8 @@ class ZernikeRadial(Polynomial):
return self._order
def __call__(self, r):
zn_rad = capi.calc_zn_rad(self.order, r)
return np.sum(self._norm_coef * zn_rad)
if isinstance(r, Iterable):
return [np.sum(self._norm_coef * capi.calc_zn_rad(self.order, r_i))
for r_i in r]
else:
return np.sum(self._norm_coef * capi.calc_zn_rad(self.order, r))

9
src/api.F90
View file

@ -20,7 +20,6 @@ module openmc_api
use random_lcg, only: openmc_get_seed, openmc_set_seed
use settings
use simulation_header
use source_header, only: openmc_extend_sources, openmc_source_set_strength
use state_point, only: openmc_statepoint_write
use tally_header
use tally_filter_header
@ -43,7 +42,6 @@ module openmc_api
public :: openmc_extend_filters
public :: openmc_extend_cells
public :: openmc_extend_materials
public :: openmc_extend_sources
public :: openmc_extend_tallies
public :: openmc_filter_get_id
public :: openmc_filter_get_type
@ -82,7 +80,6 @@ module openmc_api
public :: openmc_simulation_finalize
public :: openmc_simulation_init
public :: openmc_source_bank
public :: openmc_source_set_strength
public :: openmc_tally_allocate
public :: openmc_tally_get_estimator
public :: openmc_tally_get_id
@ -300,12 +297,16 @@ contains
use plot_header
use sab_header
use settings
use source_header
use surface_header
use tally_derivative_header
use trigger_header
use volume_header
interface
subroutine free_memory_source() bind(C)
end subroutine
end interface
call free_memory_geometry()
call free_memory_surfaces()
call free_memory_material()

43
src/distribution_multi.cpp
View file

@ -4,11 +4,30 @@
#include <cmath> // for sqrt, sin, cos, max
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/xml_interface.h"
namespace openmc {
//==============================================================================
// UnitSphereDistribution implementation
//==============================================================================
UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
{
// Read reference directional unit vector
if (check_for_node(node, "reference_uvw")) {
auto u_ref = get_node_array<double>(node, "reference_uvw");
if (u_ref.size() != 3)
fatal_error("Angular distribution reference direction must have "
"three parameters specified.");
u_ref_ = Direction(u_ref.data());
}
}
//==============================================================================
// PolarAzimuthal implementation
//==============================================================================
@ -16,15 +35,33 @@ namespace openmc {
PolarAzimuthal::PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi) :
UnitSphereDistribution{u}, mu_{std::move(mu)}, phi_{std::move(phi)} { }
PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
: UnitSphereDistribution{node}
{
if (check_for_node(node, "mu")) {
pugi::xml_node node_dist = node.child("mu");
mu_ = distribution_from_xml(node_dist);
} else {
mu_ = UPtrDist{new Uniform(-1., 1.)};
}
if (check_for_node(node, "phi")) {
pugi::xml_node node_dist = node.child("phi");
phi_ = distribution_from_xml(node_dist);
} else {
phi_ = UPtrDist{new Uniform(0.0, 2.0*PI)};
}
}
Direction PolarAzimuthal::sample() const
{
// Sample cosine of polar angle
double mu = mu_->sample();
if (mu == 1.0) return u_ref;
if (mu == 1.0) return u_ref_;
// Sample azimuthal angle
double phi = phi_->sample();
return rotate_angle(u_ref, mu, &phi);
return rotate_angle(u_ref_, mu, &phi);
}
//==============================================================================
@ -45,7 +82,7 @@ Direction Isotropic::sample() const
Direction Monodirectional::sample() const
{
return u_ref;
return u_ref_;
}
} // namespace openmc

257
src/distribution_multivariate.F90
View file

@ -1,257 +0,0 @@
module distribution_multivariate
use constants, only: ONE, TWO, PI
use distribution_univariate
use error, only: fatal_error
use random_lcg, only: prn
use math, only: rotate_angle
use xml_interface
implicit none
!===============================================================================
! UNITSPHEREDISTRIBUTION type defines a probability density function for points
! on the unit sphere. Extensions of this type are used to sample angular
! distributions for starting sources
!===============================================================================
type, abstract :: UnitSphereDistribution
real(8) :: reference_uvw(3)
contains
procedure(unitsphere_distribution_sample_), deferred :: sample
end type UnitSphereDistribution
abstract interface
function unitsphere_distribution_sample_(this) result(uvw)
import UnitSphereDistribution
class(UnitSphereDistribution), intent(in) :: this
real(8) :: uvw(3)
end function unitsphere_distribution_sample_
end interface
!===============================================================================
! Derived classes of UnitSphereDistribution
!===============================================================================
! Explicit distribution of polar and azimuthal angles
type, extends(UnitSphereDistribution) :: PolarAzimuthal
class(Distribution), allocatable :: mu
class(Distribution), allocatable :: phi
contains
procedure :: sample => polar_azimuthal_sample
end type PolarAzimuthal
! Uniform distribution on the unit sphere
type, extends(UnitSphereDistribution) :: Isotropic
contains
procedure :: sample => isotropic_sample
end type Isotropic
! Monodirectional distribution
type, extends(UnitSphereDistribution) :: Monodirectional
contains
procedure :: sample => monodirectional_sample
end type Monodirectional
!===============================================================================
! SPATIALDISTRIBUTION type defines a probability density function for arbitrary
! points in Euclidean space.
!===============================================================================
type, abstract :: SpatialDistribution
contains
procedure(spatial_distribution_from_xml_), deferred :: from_xml
procedure(spatial_distribution_sample_), deferred :: sample
end type SpatialDistribution
abstract interface
subroutine spatial_distribution_from_xml_(this, node)
import SpatialDistribution, XMLNode
class(SpatialDistribution), intent(inout) :: this
type(XMLNode), intent(in) :: node
end subroutine spatial_distribution_from_xml_
function spatial_distribution_sample_(this) result(xyz)
import SpatialDistribution
class(SpatialDistribution), intent(in) :: this
real(8) :: xyz(3)
end function spatial_distribution_sample_
end interface
type, extends(SpatialDistribution) :: CartesianIndependent
class(Distribution), allocatable :: x
class(Distribution), allocatable :: y
class(Distribution), allocatable :: z
contains
procedure :: from_xml => cartesian_independent_from_xml
procedure :: sample => cartesian_independent_sample
end type CartesianIndependent
type, extends(SpatialDistribution) :: SpatialBox
real(8) :: lower_left(3)
real(8) :: upper_right(3)
logical :: only_fissionable = .false.
contains
procedure :: from_xml => spatial_box_from_xml
procedure :: sample => spatial_box_sample
end type SpatialBox
type, extends(SpatialDistribution) :: SpatialPoint
real(8) :: xyz(3)
contains
procedure :: from_xml => spatial_point_from_xml
procedure :: sample => spatial_point_sample
end type SpatialPoint
contains
function polar_azimuthal_sample(this) result(uvw)
class(PolarAzimuthal), intent(in) :: this
real(8) :: uvw(3)
real(8) :: mu ! cosine of polar angle
real(8) :: phi ! azimuthal angle
! Sample cosine of polar angle
mu = this % mu % sample()
if (mu == ONE) then
uvw(:) = this % reference_uvw
else
! Sample azimuthal angle
phi = this % phi % sample()
uvw = rotate_angle(this % reference_uvw, mu, phi)
end if
end function polar_azimuthal_sample
function isotropic_sample(this) result(uvw)
class(Isotropic), intent(in) :: this
real(8) :: uvw(3)
real(8) :: phi
real(8) :: mu
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
uvw(1) = mu
uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
end function isotropic_sample
function monodirectional_sample(this) result(uvw)
class(Monodirectional), intent(in) :: this
real(8) :: uvw(3)
uvw(:) = this % reference_uvw
end function monodirectional_sample
subroutine cartesian_independent_from_xml(this, node)
class(CartesianIndependent), intent(inout) :: this
type(XMLNode), intent(in) :: node
type(XMLNode) :: node_dist
! Read distribution for x coordinate
if (check_for_node(node, "x")) then
node_dist = node % child("x")
call distribution_from_xml(this % x, node_dist)
else
allocate(Discrete :: this % x)
select type (dist => this % x)
type is (Discrete)
allocate(dist % x(1), dist % p(1))
dist % x(1) = ZERO
dist % p(1) = ONE
end select
end if
! Read distribution for y coordinate
if (check_for_node(node, "y")) then
node_dist = node % child("y")
call distribution_from_xml(this % y, node_dist)
else
allocate(Discrete :: this % y)
select type (dist => this % y)
type is (Discrete)
allocate(dist % x(1), dist % p(1))
dist % x(1) = ZERO
dist % p(1) = ONE
end select
end if
if (check_for_node(node, "z")) then
node_dist = node % child("z")
call distribution_from_xml(this % z, node_dist)
else
allocate(Discrete :: this % z)
select type (dist => this % z)
type is (Discrete)
allocate(dist % x(1), dist % p(1))
dist % x(1) = ZERO
dist % p(1) = ONE
end select
end if
end subroutine cartesian_independent_from_xml
function cartesian_independent_sample(this) result(xyz)
class(CartesianIndependent), intent(in) :: this
real(8) :: xyz(3)
xyz(1) = this % x % sample()
xyz(2) = this % y % sample()
xyz(3) = this % z % sample()
end function cartesian_independent_sample
subroutine spatial_box_from_xml(this, node)
class(SpatialBox), intent(inout) :: this
type(XMLNode), intent(in) :: node
real(8), allocatable :: temp_real(:)
! Make sure correct number of parameters are given
if (node_word_count(node, "parameters") /= 6) then
call fatal_error('Box/fission spatial source must have &
&six parameters specified.')
end if
! Read lower-right/upper-left coordinates
allocate(temp_real(6))
call get_node_array(node, "parameters", temp_real)
this % lower_left(:) = temp_real(1:3)
this % upper_right(:) = temp_real(4:6)
deallocate(temp_real)
end subroutine spatial_box_from_xml
function spatial_box_sample(this) result(xyz)
class(SpatialBox), intent(in) :: this
real(8) :: xyz(3)
integer :: i
real(8) :: r(3)
r = [ (prn(), i = 1,3) ]
xyz(:) = this % lower_left + r*(this % upper_right - this % lower_left)
end function spatial_box_sample
subroutine spatial_point_from_xml(this, node)
class(SpatialPoint), intent(inout) :: this
type(XMLNode), intent(in) :: node
! Make sure correct number of parameters are given
if (node_word_count(node, "parameters") /= 3) then
call fatal_error('Point spatial source must have &
&three parameters specified.')
end if
! Read location of point source
call get_node_array(node, "parameters", this % xyz)
end subroutine spatial_point_from_xml
function spatial_point_sample(this) result(xyz)
class(SpatialPoint), intent(in) :: this
real(8) :: xyz(3)
xyz(:) = this % xyz
end function spatial_point_sample
end module distribution_multivariate

3
src/distribution_spatial.cpp
View file

@ -55,7 +55,8 @@ Position CartesianIndependent::sample() const
// SpatialBox implementation
//==============================================================================
SpatialBox::SpatialBox(pugi::xml_node node)
SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
: only_fissionable_{fission}
{
// Read lower-right/upper-left coordinates
auto params = get_node_array<double>(node, "parameters");

373
src/distribution_univariate.F90
View file

@ -1,373 +0,0 @@
module distribution_univariate
use constants, only: ZERO, ONE, HALF, HISTOGRAM, LINEAR_LINEAR, &
MAX_LINE_LEN, MAX_WORD_LEN
use error, only: fatal_error
use random_lcg, only: prn
use math, only: maxwell_spectrum, watt_spectrum
use pugixml
use string, only: to_lower
use xml_interface
implicit none
!===============================================================================
! DISTRIBUTION type defines a probability density function
!===============================================================================
type, abstract :: Distribution
contains
procedure(distribution_sample_), deferred :: sample
end type Distribution
type DistributionContainer
class(Distribution), allocatable :: obj
end type DistributionContainer
abstract interface
function distribution_sample_(this) result(x)
import Distribution
class(Distribution), intent(in) :: this
real(8) :: x
end function distribution_sample_
end interface
!===============================================================================
! Derived classes of Distribution
!===============================================================================
! Discrete distribution
type, extends(Distribution) :: Discrete
real(8), allocatable :: x(:)
real(8), allocatable :: p(:)
contains
procedure :: sample => discrete_sample
procedure :: initialize => discrete_initialize
end type Discrete
! Uniform distribution over the interval [a,b]
type, extends(Distribution) :: Uniform
real(8) :: a
real(8) :: b
contains
procedure :: sample => uniform_sample
end type Uniform
! Maxwellian distribution of form c*E*exp(-E/a)
type, extends(Distribution) :: Maxwell
real(8) :: theta
contains
procedure :: sample => maxwell_sample
end type Maxwell
! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E))
type, extends(Distribution) :: Watt
real(8) :: a
real(8) :: b
contains
procedure :: sample => watt_sample
end type Watt
! Histogram or linear-linear interpolated tabular distribution
type, extends(Distribution) :: Tabular
integer :: interpolation
real(8), allocatable :: x(:) ! tabulated independent variable
real(8), allocatable :: p(:) ! tabulated probability density
real(8), allocatable :: c(:) ! cumulative distribution at tabulated values
contains
procedure :: sample => tabular_sample
procedure :: initialize => tabular_initialize
end type Tabular
type, extends(Distribution) :: Equiprobable
real(8), allocatable :: x(:)
contains
procedure :: sample => equiprobable_sample
end type Equiprobable
contains
function discrete_sample(this) result(x)
class(Discrete), intent(in) :: this
real(8) :: x
integer :: i ! loop counter
integer :: n ! size of distribution
real(8) :: c ! cumulative frequency
real(8) :: xi ! sampled CDF value
n = size(this%x)
if (n > 1) then
xi = prn()
c = ZERO
do i = 1, size(this%x)
c = c + this%p(i)
if (xi < c) exit
end do
x = this%x(i)
else
x = this%x(1)
end if
end function discrete_sample
subroutine discrete_initialize(this, x, p)
class(Discrete), intent(inout) :: this
real(8), intent(in) :: x(:)
real(8), intent(in) :: p(:)
integer :: n
! Check length of x, p arrays
if (size(x) /= size(p)) then
call fatal_error('Tabulated probabilities not of same length as &
&independent variable.')
end if
! Copy probability density function
n = size(x)
allocate(this%x(n), this%p(n))
this%x(:) = x(:)
this%p(:) = p(:)
! Normalize density function
this%p(:) = this%p(:)/sum(this%p)
end subroutine
function uniform_sample(this) result(x)
class(Uniform), intent(in) :: this
real(8) :: x
x = this%a + prn()*(this%b - this%a)
end function uniform_sample
function maxwell_sample(this) result(x)
class(Maxwell), intent(in) :: this
real(8) :: x
x = maxwell_spectrum(this%theta)
end function maxwell_sample
function watt_sample(this) result(x)
class(Watt), intent(in) :: this
real(8) :: x
x = watt_spectrum(this%a, this%b)
end function watt_sample
function tabular_sample(this) result(x)
class(Tabular), intent(in) :: this
real(8) :: x
integer :: i
real(8) :: c ! sampled cumulative frequency
real(8) :: m ! slope of PDF
real(8) :: x_i, x_i1 ! i-th and (i+1)th x values
real(8) :: c_i, c_i1 ! i-th and (i+1)th cumulative distribution values
real(8) :: p_i, p_i1 ! i-th and (i+1)th probability density values
! Sample value of CDF
c = prn()
! Find first CDF bin which is above the sampled value
c_i = this%c(1)
do i = 1, size(this%c) - 1
c_i1 = this%c(i + 1)
if (c <= c_i1) exit
c_i = c_i1
end do
! Determine bounding PDF values
x_i = this%x(i)
p_i = this%p(i)
if (this%interpolation == HISTOGRAM) then
! Histogram interpolation
if (p_i > ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i
end if
else
! Linear-linear interpolation
x_i1 = this%x(i + 1)
p_i1 = this%p(i + 1)
m = (p_i1 - p_i)/(x_i1 - x_i)
if (m == ZERO) then
x = x_i + (c - c_i)/p_i
else
x = x_i + (sqrt(max(ZERO, p_i*p_i + 2*m*(c - c_i))) - p_i)/m
end if
end if
end function tabular_sample
subroutine tabular_initialize(this, x, p, interp)
class(Tabular), intent(inout) :: this
real(8), intent(in) :: x(:)
real(8), intent(in) :: p(:)
integer, intent(in) :: interp
integer :: i
integer :: n
! Check interpolation parameter
if (interp /= HISTOGRAM .and. interp /= LINEAR_LINEAR) then
call fatal_error('Only histogram and linear-linear interpolation for tabular &
&distribution is supported.')
end if
! Check length of x, p arrays
if (size(x) /= size(p)) then
call fatal_error('Tabulated probabilities not of same length as &
&independent variable.')
end if
! Copy probability density function and interpolation parameter
n = size(x)
allocate(this%x(n), this%p(n), this%c(n))
this%interpolation = interp
this%x(:) = x(:)
this%p(:) = p(:)
! Calculate cumulative distribution function
this%c(1) = ZERO
do i = 2, n
if (this%interpolation == HISTOGRAM) then
this%c(i) = this%c(i-1) + this%p(i-1)*(this%x(i) - this%x(i-1))
elseif (this%interpolation == LINEAR_LINEAR) then
this%c(i) = this%c(i-1) + HALF*(this%p(i-1) + this%p(i)) * &
(this%x(i) - this%x(i-1))
end if
end do
! Normalize density and distribution functions
this%p(:) = this%p(:)/this%c(n)
this%c(:) = this%c(:)/this%c(n)
end subroutine tabular_initialize
function equiprobable_sample(this) result(x)
class(Equiprobable), intent(in) :: this
real(8) :: x
integer :: i
integer :: n
real(8) :: r
real(8) :: xl, xr
n = size(this%x)
r = prn()
i = 1 + int((n - 1)*r)
xl = this%x(i)
xr = this%x(i+1)
x = xl + ((n - 1)*r - i + ONE) * (xr - xl)
end function equiprobable_sample
subroutine distribution_from_xml(dist, node_dist)
class(Distribution), allocatable, intent(inout) :: dist
type(XMLNode), intent(in) :: node_dist
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: temp_str
integer :: n
integer :: temp_int
real(8), allocatable :: temp_real(:)
if (check_for_node(node_dist, "type")) then
! Determine type of distribution
call get_node_value(node_dist, "type", type)
! Determine number of parameters specified
if (check_for_node(node_dist, "parameters")) then
n = node_word_count(node_dist, "parameters")
else
n = 0
end if
! Allocate extension of Distribution
select case (to_lower(type))
case ('uniform')
allocate(Uniform :: dist)
if (n /= 2) then
call fatal_error('Uniform distribution must have two &
&parameters specified.')
end if
case ('maxwell')
allocate(Maxwell :: dist)
if (n /= 1) then
call fatal_error('Maxwell energy distribution must have one &
&parameter specified.')
end if
case ('watt')
allocate(Watt :: dist)
if (n /= 2) then
call fatal_error('Watt energy distribution must have two &
&parameters specified.')
end if
case ('discrete')
allocate(Discrete :: dist)
case ('tabular')
allocate(Tabular :: dist)
case default
call fatal_error('Invalid distribution type: ' // trim(type) // '.')
end select
! Read parameters and interpolation for distribution
select type (dist)
type is (Uniform)
allocate(temp_real(2))
call get_node_array(node_dist, "parameters", temp_real)
dist%a = temp_real(1)
dist%b = temp_real(2)
deallocate(temp_real)
type is (Maxwell)
call get_node_value(node_dist, "parameters", dist%theta)
type is (Watt)
allocate(temp_real(2))
call get_node_array(node_dist, "parameters", temp_real)
dist%a = temp_real(1)
dist%b = temp_real(2)
deallocate(temp_real)
type is (Discrete)
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n))
deallocate(temp_real)
type is (Tabular)
! Read interpolation
if (check_for_node(node_dist, "interpolation")) then
call get_node_value(node_dist, "interpolation", temp_str)
select case(to_lower(temp_str))
case ('histogram')
temp_int = HISTOGRAM
case ('linear-linear')
temp_int = LINEAR_LINEAR
case default
call fatal_error("Unknown interpolation type for distribution: " &
// trim(temp_str))
end select
else
temp_int = HISTOGRAM
end if
! Read and initialize tabular distribution
allocate(temp_real(n))
call get_node_array(node_dist, "parameters", temp_real)
call dist%initialize(temp_real(1:n/2), temp_real(n/2+1:n), temp_int)
deallocate(temp_real)
end select
end if
end subroutine distribution_from_xml
end module distribution_univariate

3
src/error.F90
View file

@ -270,8 +270,9 @@ contains
character(kind=C_CHAR), intent(in) :: message(message_len)
integer(C_INT), intent(in), value :: level
character(message_len+1) :: message_out
! Using * in the internal write adds an extra space at the beginning
write(message_out, *) message
call write_message(message_out, level)
call write_message(message_out(2:), level)
end subroutine write_message_from_c
end module error

1
src/initialize.cpp
View file

@ -2,7 +2,6 @@
#include <cstddef>
#include <cstring>
#include <iostream>
#include <sstream>
#include <string>

56
src/input_xml.F90
View file

@ -7,8 +7,6 @@ module input_xml
use cmfd_header, only: cmfd_mesh
use constants
use dict_header, only: DictIntInt, DictCharInt, DictEntryCI
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use error, only: fatal_error, warning, write_message, openmc_err_msg
use geometry, only: neighbor_lists
@ -28,7 +26,6 @@ module input_xml
use surface_header
use set_header, only: SetChar
use settings
use source_header
use stl_vector, only: VectorInt, VectorReal, VectorChar
use string, only: to_lower, to_str, str_to_int, str_to_real, &
starts_with, ends_with, split_string, &
@ -101,6 +98,13 @@ module input_xml
integer(C_INT32_T), intent(in), value :: univ
integer(C_INT) :: n
end function maximum_levels
subroutine set_particle_energy_bounds(particle, E_min, E_max) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: particle
real(C_DOUBLE), value :: E_min
real(C_DOUBLE), value :: E_max
end subroutine
end interface
contains
@ -211,7 +215,6 @@ contains
integer, allocatable :: temp_int_array(:)
integer :: n_tracks
logical :: file_exists
character(MAX_WORD_LEN) :: type
character(MAX_LINE_LEN) :: filename
type(XMLDocument) :: doc
type(XMLNode) :: root
@ -226,7 +229,6 @@ contains
type(XMLNode) :: node_vol
type(XMLNode) :: node_tab_leg
type(XMLNode), allocatable :: node_mesh_list(:)
type(XMLNode), allocatable :: node_source_list(:)
type(XMLNode), allocatable :: node_vol_list(:)
! Check if settings.xml exists
@ -455,46 +457,13 @@ contains
! ==========================================================================
! EXTERNAL SOURCE
! Get point to list of <source> elements and make sure there is at least one
call get_node_list(root, "source", node_source_list)
n = size(node_source_list)
if (n == 0) then
! Default source is isotropic point source at origin with Watt spectrum
allocate(external_source(1))
external_source % particle = NEUTRON
external_source % strength = ONE
allocate(SpatialPoint :: external_source(1) % space)
select type (space => external_source(1) % space)
type is (SpatialPoint)
space % xyz(:) = [ZERO, ZERO, ZERO]
end select
allocate(Isotropic :: external_source(1) % angle)
external_source(1) % angle % reference_uvw(:) = [ZERO, ZERO, ONE]
allocate(Watt :: external_source(1) % energy)
select type(energy => external_source(1) % energy)
type is (Watt)
energy % a = 0.988e6_8
energy % b = 2.249e-6_8
end select
else
! Allocate array for sources
allocate(external_source(n))
end if
! Handled on C++ side
! Check if we want to write out source
if (check_for_node(root, "write_initial_source")) then
call get_node_value(root, "write_initial_source", write_initial_source)
end if
! Read each source
do i = 1, n
call external_source(i) % from_xml(node_source_list(i), path_source)
end do
! Survival biasing
if (check_for_node(root, "survival_biasing")) then
call get_node_value(root, "survival_biasing", survival_biasing)
@ -1337,7 +1306,6 @@ contains
character(3) :: element ! name of element, e.g. Zr
character(MAX_WORD_LEN) :: units ! units on density
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
character(MAX_WORD_LEN), allocatable :: sarray(:)
real(8) :: val ! value entered for density
real(8) :: temp_dble ! temporary double prec. real
@ -3343,6 +3311,8 @@ contains
! Get the minimum and maximum energies
energy_min(NEUTRON) = energy_bins(num_energy_groups + 1)
energy_max(NEUTRON) = energy_bins(1)
call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), &
energy_max(NEUTRON))
! Get the datasets present in the library
call get_groups(file_id, names)
@ -3503,6 +3473,8 @@ contains
nuclides(i_nuclide) % grid(1) % energy(1))
energy_max(NEUTRON) = min(energy_max(NEUTRON), nuclides(i_nuclide) % &
grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy)))
call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), &
energy_max(NEUTRON))
end if
! Add name and alias to dictionary
@ -3536,6 +3508,8 @@ contains
energy_max(PHOTON) = min(energy_max(PHOTON), &
exp(elements(i_element) % energy(size(elements(i_element) &
% energy))))
call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), &
energy_max(PHOTON))
end if
! Add element to set
@ -3577,6 +3551,8 @@ contains
if (size(ttb_e_grid) >= 1) then
energy_min(PHOTON) = max(energy_min(PHOTON), ttb_e_grid(1))
energy_max(PHOTON) = min(energy_max(PHOTON), ttb_e_grid(size(ttb_e_grid)))
call set_particle_energy_bounds(PHOTON, energy_min(PHOTON), &
energy_max(PHOTON))
end if
! Take logarithm of energies since they are log-log interpolated

16
src/material_header.F90
View file

@ -678,6 +678,22 @@ contains
end function openmc_material_get_id
function openmc_material_get_fissionable(index, fissionable) result(err) bind(C)
! returns whether a material is fissionable
integer(C_INT32_T), value :: index
logical(C_BOOL), intent(out) :: fissionable
integer(C_INT) :: err
if (index >= 1 .and. index <= size(materials)) then
fissionable = materials(index) % fissionable
err = 0
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in materials array is out of bounds.")
end if
end function openmc_material_get_fissionable
function openmc_material_set_id(index, id) result(err) bind(C)
! Set the ID of a material
integer(C_INT32_T), value, intent(in) :: index

13
src/mgxs_interface.F90
View file

@ -147,7 +147,7 @@ module mgxs_interface
end interface
! Number of energy groups
integer(C_INT) :: num_energy_groups
integer(C_INT), bind(C) :: num_energy_groups
! Number of delayed groups
integer(C_INT) :: num_delayed_groups
@ -159,6 +159,13 @@ module mgxs_interface
real(8), allocatable :: energy_bin_avg(:)
! Energy group structure with increasing energy
real(8), allocatable :: rev_energy_bins(:)
real(C_DOUBLE), allocatable, target :: rev_energy_bins(:)
end module mgxs_interface
contains
function rev_energy_bins_ptr() result(ptr) bind(C)
type(C_PTR) :: ptr
ptr = C_LOC(rev_energy_bins(1))
end function
end module mgxs_interface

23
src/nuclide.cpp Normal file
View file

@ -0,0 +1,23 @@
#include "openmc/nuclide.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
std::array<double, 2> energy_min {0.0, 0.0};
std::array<double, 2> energy_max {INFTY, INFTY};
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void
set_particle_energy_bounds(int particle, double E_min, double E_max)
{
energy_min[particle - 1] = E_min;
energy_max[particle - 1] = E_max;
}
} // namespace openmc

8
src/settings.F90
View file

@ -9,7 +9,7 @@ module settings
! ============================================================================
! ENERGY TREATMENT RELATED VARIABLES
logical(C_BOOL) :: run_CE = .true. ! Run in CE mode?
logical(C_BOOL), bind(C, name='openmc_run_CE') :: run_CE = .true. ! Run in CE mode?
! ============================================================================
! CONTINUOUS-ENERGY CROSS SECTION RELATED VARIABLES
@ -26,7 +26,7 @@ module settings
integer :: n_log_bins ! number of bins for logarithmic grid
logical :: photon_transport = .false.
logical(C_BOOL), bind(C, name='openmc_photon_transport') :: photon_transport = .false.
integer :: electron_treatment = ELECTRON_TTB
! ============================================================================
@ -104,7 +104,8 @@ module settings
particle_restart_run = .false.
! Write out initial source
logical :: write_initial_source = .false.
logical(C_BOOL), bind(C, name='openmc_write_initial_source') :: &
write_initial_source = .false.
! Whether create fission neutrons or not. Only applied for MODE_FIXEDSOURCE
logical :: create_fission_neutrons = .true.
@ -120,7 +121,6 @@ module settings
character(MAX_FILE_LEN) :: path_input ! Path to input file
character(MAX_FILE_LEN) :: path_cross_sections = '' ! Path to cross_sections.xml
character(MAX_FILE_LEN) :: path_multipole ! Path to wmp library
character(MAX_FILE_LEN) :: path_source = '' ! Path to binary source
character(MAX_FILE_LEN) :: path_state_point ! Path to binary state point
character(MAX_FILE_LEN) :: path_source_point ! Path to binary source point
character(MAX_FILE_LEN) :: path_particle_restart ! Path to particle restart

22
src/settings.cpp
View file

@ -2,7 +2,11 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/distribution.h"
#include "openmc/distribution_multi.h"
#include "openmc/distribution_spatial.h"
#include "openmc/error.h"
#include "openmc/source.h"
#include "openmc/string_utils.h"
#include "openmc/xml_interface.h"
@ -98,6 +102,24 @@ void read_settings(pugi::xml_node* root)
temperature_range[0] = range[0];
temperature_range[1] = range[1];
}
// ==========================================================================
// EXTERNAL SOURCE
// Get point to list of <source> elements and make sure there is at least one
for (pugi::xml_node node : root->children("source")) {
external_sources.emplace_back(node);
}
// If no source specified, default to isotropic point source at origin with Watt spectrum
if (external_sources.empty()) {
SourceDistribution source {
UPtrSpace{new SpatialPoint({0.0, 0.0, 0.0})},
UPtrAngle{new Isotropic()},
UPtrDist{new Watt(0.988, 2.249e-6)}
};
external_sources.push_back(std::move(source));
}
}
} // namespace openmc

30
src/simulation.F90
View file

@ -30,7 +30,6 @@ module simulation
use random_lcg, only: set_particle_seed
use settings
use simulation_header
use source, only: initialize_source, sample_external_source
use state_point, only: openmc_statepoint_write, write_source_point, load_state_point
use string, only: to_str
use tally, only: accumulate_tallies, setup_active_tallies, &
@ -52,6 +51,19 @@ module simulation
integer(C_INT), parameter :: STATUS_EXIT_MAX_BATCH = 1
integer(C_INT), parameter :: STATUS_EXIT_ON_TRIGGER = 2
interface
subroutine openmc_simulation_init_c() bind(C)
end subroutine
subroutine initialize_source() bind(C)
end subroutine
function sample_external_source() result(site) bind(C)
import Bank
type(Bank) :: site
end function
end interface
contains
!===============================================================================
@ -241,7 +253,10 @@ contains
subroutine finalize_generation()
integer(8) :: i
interface
subroutine fill_source_bank_fixedsource() bind(C)
end subroutine
end interface
! Update global tallies with the omp private accumulation variables
!$omp parallel
@ -294,13 +309,7 @@ contains
elseif (run_mode == MODE_FIXEDSOURCE) then
! For fixed-source mode, we need to sample the external source
if (path_source == '') then
do i = 1, work
call set_particle_seed((total_gen + overall_generation()) * &
n_particles + work_index(rank) + i)
call sample_external_source(source_bank(i))
end do
end if
call fill_source_bank_fixedsource()
end if
end subroutine finalize_generation
@ -425,6 +434,9 @@ contains
! Skip if simulation has already been initialized
if (simulation_initialized) return
! Call initialization on C++ side
call openmc_simulation_init_c()
! Set up tally procedure pointers
call init_tally_routines()

47
src/simulation.cpp
View file

@ -1,4 +1,7 @@
#include "openmc/simulation.h"
#include "openmc/capi.h"
#include "openmc/message_passing.h"
// OPENMC_RUN encompasses all the main logic where iterations are performed
// over the batches, generations, and histories in a fixed source or k-eigenvalue
@ -16,3 +19,47 @@ int openmc_run() {
openmc_simulation_finalize();
return err;
}
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
std::vector<int64_t> work_index;
//==============================================================================
// Functions
//==============================================================================
void openmc_simulation_init_c()
{
// Determine how much work each process should do
calculate_work();
}
void calculate_work()
{
// Determine minimum amount of particles to simulate on each processor
int64_t min_work = n_particles/mpi::n_procs;
// Determine number of processors that have one extra particle
int64_t remainder = n_particles % mpi::n_procs;
int64_t i_bank = 0;
work_index.reserve(mpi::n_procs);
work_index.push_back(0);
for (int i = 0; i < mpi::n_procs; ++i) {
// Number of particles for rank i
int64_t work_i = i < remainder ? min_work + 1 : min_work;
// Set number of particles
if (mpi::rank == i) openmc_work = work_i;
// Set index into source bank for rank i
i_bank += work_i;
work_index.push_back(i_bank);
}
}
} // namespace openmc

6
src/simulation_header.F90
View file

@ -22,7 +22,7 @@ module simulation_header
integer(C_INT), bind(C, name='openmc_current_batch') :: current_batch ! current batch
integer(C_INT), bind(C, name='openmc_current_gen') :: current_gen ! current generation within a batch
integer :: total_gen = 0 ! total number of generations simulated
integer(C_INT), bind(C, name='openmc_total_gen') :: total_gen = 0 ! total number of generations simulated
logical(C_BOOL), bind(C, name='openmc_simulation_initialized') :: &
simulation_initialized = .false.
logical :: need_depletion_rx ! need to calculate depletion reaction rx?
@ -77,8 +77,8 @@ contains
! OVERALL_GENERATION determines the overall generation number
!===============================================================================
pure function overall_generation() result(gen)
integer :: gen
pure function overall_generation() result(gen) bind(C)
integer(C_INT) :: gen
gen = gen_per_batch*(current_batch - 1) + current_gen
end function overall_generation

141
src/source.F90
View file

@ -1,141 +0,0 @@
module source
#ifdef OPENMC_MPI
use message_passing
#endif
use algorithm, only: binary_search
use bank_header, only: Bank, source_bank
use constants
use distribution_univariate, only: Discrete
use distribution_multivariate, only: SpatialBox
use error, only: fatal_error
use geometry, only: find_cell
use hdf5_interface
use math
use message_passing, only: rank
use mgxs_interface, only: rev_energy_bins, num_energy_groups
use output, only: write_message
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_set_stream
use settings
use simulation_header
use source_header, only: external_source
use string, only: to_str
use state_point, only: read_source_bank, write_source_bank
implicit none
contains
!===============================================================================
! INITIALIZE_SOURCE initializes particles in the source bank
!===============================================================================
subroutine initialize_source()
integer(8) :: i ! loop index over bank sites
integer(8) :: id ! particle id
integer(HID_T) :: file_id
character(MAX_WORD_LEN) :: filetype
character(MAX_FILE_LEN) :: filename
type(Bank), pointer :: src ! source bank site
call write_message("Initializing source particles...", 5)
if (path_source /= '') then
! Read the source from a binary file instead of sampling from some
! assumed source distribution
call write_message('Reading source file from ' // trim(path_source) &
&// '...', 6)
! Open the binary file
file_id = file_open(path_source, 'r', parallel=.true.)
! Read the file type
call read_attribute(filetype, file_id, "filetype")
! Check to make sure this is a source file
if (filetype /= 'source' .and. filetype /= 'statepoint') then
call fatal_error("Specified starting source file not a source file &
&type.")
end if
! Read in the source bank
call read_source_bank(file_id, work_index, source_bank)
! Close file
call file_close(file_id)
else
! Generation source sites from specified distribution in user input
do i = 1, work
! Get pointer to source bank site
src => source_bank(i)
! initialize random number seed
id = total_gen*n_particles + work_index(rank) + i
call set_particle_seed(id)
! sample external source distribution
call sample_external_source(src)
end do
end if
! Write out initial source
if (write_initial_source) then
call write_message('Writing out initial source...', 5)
filename = trim(path_output) // 'initial_source.h5'
file_id = file_open(filename, 'w', parallel=.true.)
call write_source_bank(file_id, work_index, source_bank)
call file_close(file_id)
end if
end subroutine initialize_source
!===============================================================================
! SAMPLE_EXTERNAL_SOURCE samples the user-specified external source and stores
! the position, angle, and energy in a Bank type.
!===============================================================================
subroutine sample_external_source(site)
type(Bank), intent(inout) :: site ! source site
integer :: i ! dummy loop index
integer :: n_source ! number of source distributions
real(8) :: c ! cumulative frequency
real(8) :: xi
! Set the random number generator to the source stream.
call prn_set_stream(STREAM_SOURCE)
! Sample from among multiple source distributions
n_source = size(external_source)
if (n_source > 1) then
xi = prn()*sum(external_source(:) % strength)
c = ZERO
do i = 1, n_source
c = c + external_source(i) % strength
if (xi < c) exit
end do
else
i = 1
end if
! Sample source site from i-th source distribution
site = external_source(i) % sample()
! If running in MG, convert site % E to group
if (.not. run_CE) then
site % E = real(binary_search(rev_energy_bins, num_energy_groups + 1, &
site % E), 8)
site % E = num_energy_groups + 1 - site % E
end if
! Set the random number generator back to the tracking stream.
call prn_set_stream(STREAM_TRACKING)
end subroutine sample_external_source
end module source

378
src/source.cpp Normal file
View file

@ -0,0 +1,378 @@
#include "openmc/source.h"
#include <algorithm> // for move
#include <sstream> // for stringstream
#include "xtensor/xadapt.hpp"
#include "openmc/cell.h"
#include "openmc/error.h"
#include "openmc/file_utils.h"
#include "openmc/hdf5_interface.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/capi.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/state_point.h"
#include "openmc/xml_interface.h"
namespace openmc {
//==============================================================================
// Global variables
//==============================================================================
std::vector<SourceDistribution> external_sources;
//==============================================================================
// SourceDistribution implementation
//==============================================================================
SourceDistribution::SourceDistribution(UPtrSpace space, UPtrAngle angle, UPtrDist energy)
: space_{std::move(space)}, angle_{std::move(angle)}, energy_{std::move(energy)} { }
SourceDistribution::SourceDistribution(pugi::xml_node node)
{
// Check for particle type
if (check_for_node(node, "particle")) {
auto temp_str = get_node_value(node, "particle", true, true);
if (temp_str == "neutron") {
particle_ = ParticleType::neutron;
} else if (temp_str == "photon") {
particle_ = ParticleType::photon;
openmc_photon_transport = true;
} else {
fatal_error(std::string("Unknown source particle type: ") + temp_str);
}
}
// Check for source strength
if (check_for_node(node, "strength")) {
strength_ = std::stod(get_node_value(node, "strength"));
}
// Check for external source file
if (check_for_node(node, "file")) {
// Copy path of source file
path_source = get_node_value(node, "file", false, true);
// Check if source file exists
if (!file_exists(path_source)) {
std::stringstream msg;
msg << "Source file '" << path_source << "' does not exist.";
fatal_error(msg);
}
} else {
// Spatial distribution for external source
if (check_for_node(node, "space")) {
// Get pointer to spatial distribution
pugi::xml_node node_space = node.child("space");
// Check for type of spatial distribution and read
std::string type;
if (check_for_node(node_space, "type"))
type = get_node_value(node_space, "type", true, true);
if (type == "cartesian") {
space_ = UPtrSpace{new CartesianIndependent(node_space)};
} else if (type == "box") {
space_ = UPtrSpace{new SpatialBox(node_space)};
} else if (type == "fission") {
space_ = UPtrSpace{new SpatialBox(node_space, true)};
} else if (type == "point") {
space_ = UPtrSpace{new SpatialPoint(node_space)};
} else {
std::stringstream msg;
msg << "Invalid spatial distribution for external source: " << type;
fatal_error(msg);
}
} else {
// If no spatial distribution specified, make it a point source
space_ = UPtrSpace{new SpatialPoint()};
}
// Determine external source angular distribution
if (check_for_node(node, "angle")) {
// Get pointer to angular distribution
pugi::xml_node node_angle = node.child("angle");
// Check for type of angular distribution
std::string type;
if (check_for_node(node_angle, "type"))
type = get_node_value(node_angle, "type", true, true);
if (type == "isotropic") {
angle_ = UPtrAngle{new Isotropic()};
} else if (type == "monodirectional") {
angle_ = UPtrAngle{new Monodirectional(node_angle)};
} else if (type == "mu-phi") {
angle_ = UPtrAngle{new PolarAzimuthal(node_angle)};
} else {
std::stringstream msg;
msg << "Invalid angular distribution for external source: " << type;
fatal_error(msg);
}
} else {
angle_ = UPtrAngle{new Isotropic()};
}
// Determine external source energy distribution
if (check_for_node(node, "energy")) {
pugi::xml_node node_dist = node.child("energy");
energy_ = distribution_from_xml(node_dist);
} else {
// Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
energy_ = UPtrDist{new Watt(0.988e6, 2.249e-6)};
}
}
}
Bank SourceDistribution::sample() const
{
Bank site;
// Set weight to one by default
site.wgt = 1.0;
// Repeat sampling source location until a good site has been found
bool found = false;
int n_reject = 0;
static int n_accept = 0;
while (!found) {
// Set particle type
site.particle = static_cast<int>(particle_);
// Sample spatial distribution
Position r = space_->sample();
site.xyz[0] = r.x;
site.xyz[1] = r.y;
site.xyz[2] = r.z;
// Now search to see if location exists in geometry
int32_t cell_index, instance;
int err = openmc_find_cell(site.xyz, &cell_index, &instance);
found = (err != OPENMC_E_GEOMETRY);
// Check if spatial site is in fissionable material
if (found) {
auto space_box = dynamic_cast<SpatialBox*>(space_.get());
if (space_box) {
if (space_box->only_fissionable()) {
// Determine material
auto c = cells[cell_index - 1];
int32_t mat_index = c->material_[instance];
auto m = materials[mat_index];
if (mat_index == MATERIAL_VOID) {
found = false;
} else {
bool fissionable;
openmc_material_get_fissionable(mat_index + 1, &fissionable);
if (!fissionable) found = false;
}
}
}
}
// Check for rejection
if (!found) {
++n_reject;
if (n_reject >= EXTSRC_REJECT_THRESHOLD &&
static_cast<double>(n_accept)/n_reject <= EXTSRC_REJECT_FRACTION) {
fatal_error("More than 95% of external source sites sampled were "
"rejected. Please check your external source definition.");
}
}
}
// Increment number of accepted samples
++n_accept;
// Sample angle
Direction u = angle_->sample();
site.uvw[0] = u.x;
site.uvw[1] = u.y;
site.uvw[2] = u.z;
// Check for monoenergetic source above maximum particle energy
auto p = static_cast<int>(particle_);
auto energy_ptr = dynamic_cast<Discrete*>(energy_.get());
if (energy_ptr) {
auto energies = xt::adapt(energy_ptr->x());
if (xt::any(energies > energy_max[p-1])) {
fatal_error("Source energy above range of energies of at least "
"one cross section table");
} else if (xt::any(energies < energy_min[p-1])) {
fatal_error("Source energy below range of energies of at least "
"one cross section table");
}
}
while (true) {
// Sample energy spectrum
site.E = energy_->sample();
// Resample if energy falls outside minimum or maximum particle energy
if (site.E < energy_max[p-1] && site.E > energy_min[p-1]) break;
}
// Set delayed group
site.delayed_group = 0;
return site;
}
//==============================================================================
// Non-member functions
//==============================================================================
void initialize_source()
{
write_message("Initializing source particles...", 5);
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
if (path_source != "") {
// Read the source from a binary file instead of sampling from some
// assumed source distribution
std::stringstream msg;
msg << "Reading source file from " << path_source << "...";
write_message(msg, 6);
// Open the binary file
hid_t file_id = file_open(path_source, 'r', true);
// Read the file type
std::string filetype;
read_attribute(file_id, "filetype", filetype);
// Check to make sure this is a source file
if (filetype != "source" && filetype != "statepoint") {
fatal_error("Specified starting source file not a source file type.");
}
// Read in the source bank
read_source_bank(file_id, work_index.data(), source_bank);
// Close file
file_close(file_id);
} else {
// Generation source sites from specified distribution in user input
for (int64_t i = 0; i < openmc_work; ++i) {
// initialize random number seed
int64_t id = openmc_total_gen*n_particles + work_index[openmc::mpi::rank] + i + 1;
set_particle_seed(id);
// sample external source distribution
source_bank[i] = sample_external_source();
}
}
// Write out initial source
if (openmc_write_initial_source) {
write_message("Writing out initial source...", 5);
std::string filename = path_output + "initial_source.h5";
hid_t file_id = file_open(filename, 'w', true);
write_source_bank(file_id, work_index.data(), source_bank);
file_close(file_id);
}
}
extern "C" double* rev_energy_bins_ptr();
Bank sample_external_source()
{
// Set the random number generator to the source stream.
prn_set_stream(STREAM_SOURCE);
// Determine total source strength
double total_strength = 0.0;
for (auto& s : external_sources)
total_strength += s.strength();
// Sample from among multiple source distributions
int i = 0;
if (external_sources.size() > 1) {
double xi = prn()*total_strength;
double c = 0.0;
for (; i < external_sources.size(); ++i) {
c += external_sources[i].strength();
if (xi < c) break;
}
}
// Sample source site from i-th source distribution
Bank site {external_sources[i].sample()};
// If running in MG, convert site % E to group
if (!openmc_run_CE) {
// Get pointer to rev_energy_bins array on Fortran side
double* rev_energy_bins = rev_energy_bins_ptr();
int n = num_energy_groups + 1;
site.E = lower_bound_index(rev_energy_bins, rev_energy_bins + n, site.E);
site.E = num_energy_groups - site.E;
}
// Set the random number generator back to the tracking stream.
prn_set_stream(STREAM_TRACKING);
return site;
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
extern "C" void free_memory_source()
{
external_sources.clear();
}
extern "C" double total_source_strength()
{
double strength = 0.0;
for (const auto& s : external_sources) {
strength += s.strength();
}
return strength;
}
// Needed in fill_source_bank_fixedsource
extern "C" int overall_generation();
//! Fill source bank at end of generation for fixed source simulations
extern "C" void fill_source_bank_fixedsource()
{
if (path_source.empty()) {
// Get pointer to source bank
Bank* source_bank;
int64_t n;
openmc_source_bank(&source_bank, &n);
for (int64_t i = 0; i < openmc_work; ++i) {
// initialize random number seed
int64_t id = (openmc_total_gen + overall_generation())*n_particles +
work_index[openmc::mpi::rank] + i + 1;
set_particle_seed(id);
// sample external source distribution
source_bank[i] = sample_external_source();
}
}
}
} // namespace openmc

388
src/source_header.F90
View file

@ -1,388 +0,0 @@
module source_header
use, intrinsic :: ISO_C_BINDING
use bank_header, only: Bank
use constants
use distribution_univariate
use distribution_multivariate
use error
use geometry, only: find_cell
use material_header, only: materials
use nuclide_header, only: energy_min, energy_max
use particle_header
use settings, only: photon_transport
use string, only: to_lower
use xml_interface
implicit none
private
public :: free_memory_source
public :: openmc_extend_sources
public :: openmc_source_set_strength
integer :: n_accept = 0 ! Number of samples accepted
integer :: n_reject = 0 ! Number of samples rejected
!===============================================================================
! SOURCEDISTRIBUTION describes an external source of particles for a
! fixed-source problem or for the starting source in a k eigenvalue problem
!===============================================================================
type, public :: SourceDistribution
integer :: particle ! particle type
real(8) :: strength = ONE ! source strength
class(SpatialDistribution), allocatable :: space ! spatial distribution
class(UnitSphereDistribution), allocatable :: angle ! angle distribution
class(Distribution), allocatable :: energy ! energy distribution
contains
procedure :: from_xml
procedure :: sample
end type SourceDistribution
! Number of external source distributions
integer(C_INT32_T), public, bind(C) :: n_sources = 0
! External source distributions
type(SourceDistribution), public, allocatable :: external_source(:)
contains
subroutine from_xml(this, node, path_source)
class(SourceDistribution), intent(inout) :: this
type(XMLNode), intent(in) :: node
character(MAX_FILE_LEN), intent(out) :: path_source
integer :: n
logical :: file_exists
character(MAX_WORD_LEN) :: type, temp_str
type(XMLNode) :: node_space
type(XMLNode) :: node_angle
type(XMLNode) :: node_dist
! Check for particle type
if (check_for_node(node, "particle")) then
call get_node_value(node, "particle", temp_str)
select case (to_lower(temp_str))
case ('neutron')
this % particle = NEUTRON
case ('photon')
this % particle = PHOTON
photon_transport = .true.
case default
call fatal_error('Unknown source particle type: ' // trim(temp_str))
end select
else
this % particle = NEUTRON
end if
! Check for source strength
if (check_for_node(node, "strength")) then
call get_node_value(node, "strength", this % strength)
end if
! Check for external source file
if (check_for_node(node, "file")) then
! Copy path of source file
call get_node_value(node, "file", path_source)
! Check if source file exists
inquire(FILE=path_source, EXIST=file_exists)
if (.not. file_exists) then
call fatal_error("Source file '" // trim(path_source) &
// "' does not exist!")
end if
else
! Spatial distribution for external source
if (check_for_node(node, "space")) then
! Get pointer to spatial distribution
node_space = node % child("space")
! Check for type of spatial distribution
type = ''
if (check_for_node(node_space, "type")) &
call get_node_value(node_space, "type", type)
select case (to_lower(type))
case ('cartesian')
allocate(CartesianIndependent :: this % space)
case ('box')
allocate(SpatialBox :: this % space)
case ('fission')
allocate(SpatialBox :: this % space)
select type(space => this % space)
type is (SpatialBox)
space % only_fissionable = .true.
end select
case ('point')
allocate(SpatialPoint :: this % space)
case default
call fatal_error("Invalid spatial distribution for external source: "&
// trim(type))
end select
! Read spatial distribution from XML
call this % space % from_xml(node_space)
else
! If no spatial distribution specified, make it a point source
allocate(SpatialPoint :: this % space)
select type (space => this % space)
type is (SpatialPoint)
space % xyz(:) = [ZERO, ZERO, ZERO]
end select
end if
! Determine external source angular distribution
if (check_for_node(node, "angle")) then
! Get pointer to angular distribution
node_angle = node % child("angle")
! Check for type of angular distribution
type = ''
if (check_for_node(node_angle, "type")) &
call get_node_value(node_angle, "type", type)
select case (to_lower(type))
case ('isotropic')
allocate(Isotropic :: this % angle)
case ('monodirectional')
allocate(Monodirectional :: this % angle)
case ('mu-phi')
allocate(PolarAzimuthal :: this % angle)
case default
call fatal_error("Invalid angular distribution for external source: "&
// trim(type))
end select
! Read reference directional unit vector
if (check_for_node(node_angle, "reference_uvw")) then
n = node_word_count(node_angle, "reference_uvw")
if (n /= 3) then
call fatal_error('Angular distribution reference direction must have &
&three parameters specified.')
end if
call get_node_array(node_angle, "reference_uvw", &
this % angle % reference_uvw)
else
! By default, set reference unit vector to be positive z-direction
this % angle % reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Read parameters for angle distribution
select type (angle => this % angle)
type is (Monodirectional)
call get_node_array(node_angle, "reference_uvw", &
this % angle % reference_uvw)
type is (PolarAzimuthal)
if (check_for_node(node_angle, "mu")) then
node_dist = node_angle % child("mu")
call distribution_from_xml(angle % mu, node_dist)
else
allocate(Uniform :: angle%mu)
select type (mu => angle%mu)
type is (Uniform)
mu % a = -ONE
mu % b = ONE
end select
end if
if (check_for_node(node_angle, "phi")) then
node_dist = node_angle % child("phi")
call distribution_from_xml(angle % phi, node_dist)
else
allocate(Uniform :: angle%phi)
select type (phi => angle%phi)
type is (Uniform)
phi % a = ZERO
phi % b = TWO*PI
end select
end if
end select
else
! Set default angular distribution isotropic
allocate(Isotropic :: this % angle)
this % angle % reference_uvw(:) = [ZERO, ZERO, ONE]
end if
! Determine external source energy distribution
if (check_for_node(node, "energy")) then
node_dist = node % child("energy")
call distribution_from_xml(this % energy, node_dist)
else
! Default to a Watt spectrum with parameters 0.988 MeV and 2.249 MeV^-1
allocate(Watt :: this % energy)
select type(energy => this % energy)
type is (Watt)
energy % a = 0.988e6_8
energy % b = 2.249e-6_8
end select
end if
end if
end subroutine from_xml
function sample(this) result(site)
class(SourceDistribution), intent(in) :: this
type(Bank) :: site
logical :: found ! Does the source particle exist within geometry?
type(Particle) :: p ! Temporary particle for using find_cell
! Set weight to one by default
site % wgt = ONE
! Repeat sampling source location until a good site has been found
found = .false.
do while (.not. found)
! Set particle defaults
call particle_initialize(p)
! Set particle type
site % particle = this % particle
! Sample spatial distribution
site % xyz(:) = this % space % sample()
! Fill p with needed data
p % coord(1) % xyz(:) = site % xyz
p % coord(1) % uvw(:) = [ ONE, ZERO, ZERO ]
! Now search to see if location exists in geometry
call find_cell(p, found)
! Check if spatial site is in fissionable material
if (found) then
select type (space => this % space)
type is (SpatialBox)
if (space % only_fissionable) then
if (p % material == MATERIAL_VOID) then
found = .false.
elseif (.not. materials(p % material) % fissionable) then
found = .false.
end if
end if
end select
end if
! Check for rejection
if (.not. found) then
n_reject = n_reject + 1
if (n_reject >= EXTSRC_REJECT_THRESHOLD .and. &
real(n_accept, 8)/n_reject <= EXTSRC_REJECT_FRACTION) then
call fatal_error("More than 95% of external source sites sampled &
&were rejected. Please check your external source definition.")
end if
end if
end do
! Increment number of accepted samples
n_accept = n_accept + 1
call particle_clear(p)
! Sample angle
site % uvw(:) = this % angle % sample()
! Check for monoenergetic source above maximum particle energy
select type (energy => this % energy)
type is (Discrete)
if (any(energy % x > energy_max(this % particle))) then
call fatal_error("Source energy above range of energies of at least &
&one cross section table")
else if (any(energy % x < energy_min(this % particle))) then
call fatal_error("Source energy below range of energies of at least &
&one cross section table")
end if
end select
do
! Sample energy spectrum
site % E = this % energy % sample()
! Resample if energy falls outside minimum or maximum particle energy
if (site % E < energy_max(this % particle) .and. &
site % E > energy_min(this % particle)) exit
end do
! Set delayed group
site % delayed_group = 0
end function sample
!===============================================================================
! FREE_MEMORY_SOURCE deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_source()
n_sources = 0
if (allocated(external_source)) deallocate(external_source)
end subroutine free_memory_source
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_extend_sources(n, index_start, index_end) result(err) bind(C)
! Extend the external_source array by n elements
integer(C_INT32_T), value, intent(in) :: n
integer(C_INT32_T), optional, intent(out) :: index_start
integer(C_INT32_T), optional, intent(out) :: index_end
integer(C_INT) :: err
type(SourceDistribution), allocatable :: temp(:) ! temporary array
if (n_sources == 0) then
! Allocate external_source array
allocate(external_source(n))
else
! Allocate external_source array with increased size
allocate(temp(n_sources + n))
! Copy original source array to temporary array
temp(1:n_sources) = external_source
! Move allocation from temporary array
call move_alloc(FROM=temp, TO=external_source)
end if
! Return indices in external_source array
if (present(index_start)) index_start = n_sources + 1
if (present(index_end)) index_end = n_sources + n
n_sources = n_sources + n
err = 0
end function openmc_extend_sources
function openmc_source_set_strength(index, strength) result(err) bind(C)
integer(C_INT32_T), value, intent(in) :: index
real(C_DOUBLE), value, intent(in) :: strength
integer(C_INT) :: err
if (index >= 1 .and. index <= n_sources) then
if (strength > ZERO) then
external_source(index) % strength = strength
err = 0
else
err = E_INVALID_ARGUMENT
call set_errmsg("Source strength must be positive.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in external source array is out of bounds.")
end if
end function openmc_source_set_strength
end module source_header

10
src/tallies/tally_header.F90
View file

@ -9,7 +9,6 @@ module tally_header
use message_passing, only: n_procs
use nuclide_header, only: nuclide_dict
use settings, only: reduce_tallies, run_mode
use source_header, only: external_source
use stl_vector, only: VectorInt
use string, only: to_lower, to_f_string, str_to_int, to_str
use tally_filter_header, only: TallyFilterContainer, filters, n_filters
@ -169,6 +168,13 @@ contains
real(C_DOUBLE) :: val
real(C_DOUBLE) :: total_source
interface
function total_source_strength() result(strength) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: strength
end function
end interface
! Increment number of realizations
if (reduce_tallies) then
this % n_realizations = this % n_realizations + 1
@ -178,7 +184,7 @@ contains
! Calculate total source strength for normalization
if (run_mode == MODE_FIXEDSOURCE) then
total_source = sum(external_source(:) % strength)
total_source = total_source_strength()
else
total_source = ONE
end if

11
src/xml_interface.F90
View file

@ -18,6 +18,7 @@ module xml_interface
interface get_node_value
module procedure get_node_value_bool
module procedure get_node_value_cbool
module procedure get_node_value_integer
module procedure get_node_value_long
module procedure get_node_value_double
@ -118,6 +119,16 @@ contains
end if
end subroutine get_node_value_bool
subroutine get_node_value_cbool(node, name, val)
type(XMLNode), intent(in) :: node
character(*), intent(in) :: name
logical(kind=C_BOOL), intent(out) :: val
logical :: val_
call get_node_value_bool(node, name, val_)
val = val_
end subroutine get_node_value_cbool
!===============================================================================
! GET_NODE_VALUE_INTEGER returns a integer value from an attribute or node
!===============================================================================

2
tests/unit_tests/test_capi.py
View file

@ -49,7 +49,7 @@ def capi_init(pincell_model):
@pytest.fixture(scope='module')
def capi_simulation_init(pincell_model):
def capi_simulation_init(capi_init):
openmc.capi.simulation_init()
yield

15
tests/unit_tests/test_polynomials.py
View file

@ -27,3 +27,18 @@ def test_zernike_radial():
test_vals = zn_rad(rho)
assert ref_vals == test_vals
rho = [0.2, 0.5]
# Reference solution from running the Fortran implementation
raw_zn1 = np.array([
1.00000000e+00, -9.20000000e-01, 7.69600000e-01,
-5.66720000e-01, 3.35219200e-01, -1.01747000e-01])
raw_zn2 = np.array([
1.00000000e+00, -5.00000000e-01, -1.25000000e-01,
4.37500000e-01, -2.89062500e-01, -8.98437500e-02])
ref_vals = [np.sum(norm_coeff * raw_zn1), np.sum(norm_coeff * raw_zn2)]
test_vals = zn_rad(rho)
assert np.allclose(ref_vals, test_vals)