I want to use the customized source to model more complicated source geometries in combination with giving a list of energies 。
You will be happy to hear that lots of energy distributions are available
Here is a link to the relevant part of the docs
https://docs.openmc.org/en/stable/pythonapi/stats.html#univariate-probability-distributions
Thank you very much！@ Shimwell
I used a custom source to compile the energy spectrum source。Custom sources are written in the c++ format, such as:
#include <memory> // for unique_ptr
#include <cmath> // for M_PI
#include "openmc/source.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/distribution_multi.h"
class CustomSource : public openmc::Source
{
openmc::SourceSite sample(uint64_t* seed) const
{
openmc::SourceSite particle;
// weight
particle.particle = openmc::ParticleType::neutron;
particle.wgt = 1.0;
// position
double angle = 2.0 * M_PI * openmc::prn(seed);
double r_min =5.0;
double r_max =7.5;
double radius_sq = r_min*r_min + openmc::prn(seed)*(r_max*r_max - r_min*r_min);
double radius = std::sqrt(radius_sq);
//double radius = 3.0;
particle.r.x = radius * std::cos(angle);
particle.r.y = radius * std::sin(angle);
particle.r.z = 0.0;
// angle
particle.u = {0.0, 0.0, 1.0};
particle.E = 14.08e6;
particle.delayed_group = 0;
return particle;
}
};
extern "C" std::unique_ptr<CustomSource> openmc_create_source(std::string parameters)
{
return std::make_unique<CustomSource>();
}
If I refer to the example , how can I use those energy distribution functions?
Or there are other ways to compile energy spectrum sources。
Ah sorry form the first question I had assumed you were using the python API.
Have you seen this c++ customized source by Andy Davis it has energy, radius and other sampling.
In your example you have 14MeV so this might do just what you are looking for.
Shimwell,Thank you very much!
I am mainly engaged in dose calculation of Boron neutron capture therapy （BNCT）。We wanted to use openmc as the dose engine，This is due to the relatively fast computational speed of openmc。
For BNCT, we need to get the energy, space, type and emission Angle of the particle， as the source term。
Since the openmc source project has not developed the function of rotation, the custom source method (c++) can only be used to write the source items after rotation 。
I want to write the rotated source term as a custom source, such as:
#include <memory> // for unique_ptr
#include <cmath> // for M_PI
#include "openmc/source.h"
#include "openmc/particle.h"
#include "openmc/random_lcg.h"
#include "openmc/distribution_multi.h"
class CustomSource : public openmc::Source
{
//public:
//CustomSource(double energy) : energy_{energy} { }
// Samples from an instance of this class.
openmc::SourceSite sample(uint64_t* seed) const
{
openmc::SourceSite particle;
// weight
particle.particle = openmc::ParticleType::neutron;
particle.wgt = 1.0;
openmc::SourceSite particle2;
particle2.particle = openmc::ParticleType::neutron;
particle2.wgt = 1.0;
// position
//沿x/y/z轴旋转
//定义
double Angle1,Angle2,Angle3,c3, s3,c2, s2,c1, s1;
Angle1=45;
Angle2=0;
//沿z轴旋转
Angle3=90;
//定义旋转矩阵A10
c3 = cos(Angle1/180*M_PI), s3=sin(Angle1/180*M_PI);
c2 = cos(Angle2/180*M_PI), s2=sin(Angle2/180*M_PI);
c1 = cos(Angle3/180*M_PI), s1=sin(Angle3/180*M_PI);
//定义圆面的中心点A0
double x0=0;
double y0=0;
double z0=-15;
double A0[3]={x0, y0, z0};
//r为1~7.5之间的任何值
double r_min =3;
double r_max =7.5;
double radius_sq = r_min*r_min + openmc::prn(seed)*(r_max*r_max - r_min*r_min);
double radius = std::sqrt(radius_sq);
//角度为0~360之间的任何值
double angle = 2* M_PI * openmc::prn(seed);
//初始源空间位置抽样
particle.r.x = radius * std::cos(angle);
particle.r.y = radius * std::sin(angle);
particle.r.z = 0.0;
//旋转后x/y/z的坐标
particle2.r.x= (particle.r.x+x0)*(c1*c2)+(particle.r.y+y0)*(-c2*s1)+(particle.r.z+z0)*(s2);
particle2.r.y= (particle.r.x+x0)*(c1*s2 + c3*s1)+(particle.r.y+y0)*(c1*c3 - s1*s2*s3)+(particle.r.z+z0)*(c2*s3);
particle2.r.z= (particle.r.x+x0)*(s1*s3 + c1*c3*s2)+(particle.r.y+y0)*(c3*s1*s2 + c1*s3)+(particle.r.z+z0)*(c2*c3);
//平移矩阵,x11、y11、z11平移参数
double x11=0;
double y11=20;
double z11=20;
//平移后x/y/z的坐标
particle.r.x=particle2.r.x+x11;
particle.r.y=particle2.r.y+y11;
particle.r.z=particle2.r.z+z11;
// angle
// //参考方向
// //两个点旋转平移后的矢量
double x1=0;
double y1=0;
double z1=0;
//矢量方向，源粒子入射参考方向
double x21=(s2)/sqrt((s2)*(s2)+(c2*s3)*(c2*s3)+(c2*c3)*(c2*c3));
double y21=(c2*s3)/sqrt((s2)*(s2)+(c2*s3)*(c2*s3)+(c2*c3)*(c2*c3));
double z21=(c2*c3)/sqrt((s2)*(s2)+(c2*s3)*(c2*s3)+(c2*c3)*(c2*c3));
particle.u = {x21, y21, z21};
particle.E = 14.08;
particle.delayed_group = 0;
return particle;
}
};
extern "C" std::unique_ptr<CustomSource> openmc_create_source(std::string parameter)
{
//double energy = std::stod(parameter);
return std::make_unique<CustomSource>();
}
However, only the spatial distribution, single energy and single direction emission of the rotated source term are solved.
I also need to solve the problem of when the energy of the particle is multiple energies and the emission direction is multiple directions.
At the same time, the particle emission direction should be written in the form of PolarAzimuthal() function.This is because before rotation, the emission direction of the particle is written in the following way：
# define 入射角度参数
mu = openmc.stats.Tabular([-1,
-8.3333E-01,
-6.6667E-01,
-5.0000E-01,
-3.3333E-01,
-1.6667E-01,
-1.6653E-01,
1.6667E-01,
3.3333E-01,
5.00E-01,
6.6667E-01,
8.3333E-01,
1.0000E+00],
[0,
0.003807914,
0.003064343,
0.002326964,
0.001658866,
0.001023727,
0.000365031,
0.013973863,
0.048526217,
0.094445023,
0.164759809,
0.265680708,
0.400367534
])
phi = openmc.stats.Uniform(0.0, 2*pi)
# define 柱坐标系参数，圆面的半径为5cm，中心点为（0，0，0）
radius=7.5
r11 = openmc.stats.PowerLaw(0.0, radius, 1.0)
phi1 = openmc.stats.Uniform(0.0, 2*pi)
z11 = openmc.stats.Discrete([0], [1.0])
# Define source
source = openmc.Source()
source.particle = 'neutron'
source.space = openmc.stats.CylindricalIndependent(r11, phi1, z11,origin=(0, 0, 0))
source.angle = openmc.stats.PolarAzimuthal(mu, phi, reference_uvw=(0, 0, 1))
source.energy = openmc.stats.Tabular([0,
1.000E-03,
1.239E-03,
1.535E-03,
1.901E-03,
2.355E-03,
2.918E-03
...],
[0.0000000E+00,
8.3152743E-04,
8.3131809E-04,
8.2065299E-04,
3.3259815E-03,
1.2242576E-03,
1.1891521E-03,
...0],
'histogram') #
Please, what should I do?
Or the same type of example.
thanks very mush!