Skip to content

Commit 607319d

Browse files
makecleanmakeclean
authored
Merge pull request #1 from Shimwell/first_draft
First draft
2 parents 11f0faf + 32db5fd commit 607319d

6 files changed

Lines changed: 553 additions & 0 deletions

File tree

‎README.md‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,5 @@
11
# parametric-plasma-source
22
Source and build files for parametric plasma source for use in fusion neutron transport calculations.
3+
4+
The plasma source is based on a paper by [C. Fausser et al](https://www.sciencedirect.com/science/article/pii/S0920379612000853)
5+

‎parametric_plasma_source/Makefile‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
2+
# Makefile to build dynamic sources for OpenMC
3+
# this assumes that your source sampling filename is
4+
# source_sampling.cpp
5+
#
6+
# you can add fortran, c and cpp dependencies to this source
7+
# adding them in FC_DEPS, C_DEPS, and CXX_DEPS accordingly
8+
9+
ifeq ($(FC), f77)
10+
FC = gfortran
11+
endif
12+
13+
default: all
14+
15+
ALL = source_sampling
16+
# add your fortran depencies here
17+
FC_DEPS =
18+
# add your c dependencies here
19+
C_DEPS =
20+
# add your cpp dependencies here
21+
CPP_DEPS = plasma_source.cpp
22+
DEPS = $(FC_DEPS) $(C_DEPS) $(CPP_DEPS)
23+
OPT_LEVEL = -O3
24+
FFLAGS = $(OPT_LEVEL) -fPIC
25+
C_FLAGS = -fPIC
26+
CXX_FLAGS = -fPIC
27+
28+
#this directory will need changing if your openmc path is different
29+
OPENMC_DIR = /opt/openmc
30+
OPENMC_INC_DIR = $(OPENMC_DIR)/include
31+
OPENMC_LIB_DIR = $(OPENMC_DIR)/lib
32+
# setting the so name is important
33+
LINK_FLAGS = $(OPT_LEVEL) -Wall -Wl,-soname,source_sampling.so
34+
# setting shared is important
35+
LINK_FLAGS += -L$(OPENMC_LIB_DIR) -lopenmc -lgfortran -fPIC -shared
36+
OPENMC_INCLUDES = -I$(OPENMC_INC_DIR) -I$(OPENMC_DIR)/vendor/pugixml
37+
38+
all: $(ALL)
39+
40+
source_sampling: $(DEPS)
41+
$(CXX) source_sampling.cpp $(DEPS) $(OPENMC_INCLUDES) $(LINK_FLAGS) -o $@.so
42+
# make any fortran objects
43+
%.o : %.F90
44+
$(FC) -c $(FFLAGS) $*.F90 -o $@
45+
# make any c objects
46+
%.o : %.c
47+
$(CC) -c $(FFLAGS) $*.c -o $@
48+
#make any cpp objects
49+
%.o : %.cpp
50+
$(CXX) -c $(FFLAGS) $*.cpp -o $@
51+
clean:
52+
rm -rf *.o *.mod

‎parametric_plasma_source/compile.sh‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
#!/bin/bash
2+
3+
make clean
4+
5+
make

‎parametric_plasma_source/plasma_source.cpp‎

Lines changed: 285 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,285 @@
1+
#include <iostream>
2+
#include <vector>
3+
#include <cmath>
4+
#include "plasma_source.hpp"
5+
#include <stdlib.h>
6+
#include "openmc/random_lcg.h"
7+
8+
#define RANDOM openmc::prn()
9+
10+
namespace plasma_source {
11+
12+
// default constructor
13+
PlasmaSource::PlasmaSource() {}
14+
15+
// large constructor
16+
PlasmaSource::PlasmaSource(const double ion_density_ped, const double ion_density_sep,
17+
const double ion_density_origin, const double ion_temp_ped,
18+
const double ion_temp_sep, const double ion_temp_origin,
19+
const double pedistal_rad, const double ion_density_peak,
20+
const double ion_temp_peak, const double minor_radius,
21+
const double major_radius, const double elongation,
22+
const double triangularity, const double shafranov,
23+
const std::string plasma_type, const int plasma_id,
24+
const int number_of_bins,
25+
const double min_toroidal_angle,
26+
const double max_toroidal_angle ) {
27+
28+
// set params
29+
ionDensityPedistal = ion_density_ped;
30+
ionDensitySeperatrix = ion_density_sep;
31+
ionDensityOrigin = ion_density_origin;
32+
ionTemperaturePedistal = ion_temp_ped;
33+
ionTemperatureSeperatrix = ion_temp_sep;
34+
ionTemperatureOrigin = ion_temp_origin;
35+
pedistalRadius = pedistal_rad;
36+
ionDensityPeaking = ion_density_peak;
37+
ionTemperaturePeaking = ion_temp_peak;
38+
minorRadius = minor_radius;
39+
majorRadius = major_radius;
40+
this->elongation = elongation;
41+
this->triangularity = triangularity;
42+
this->shafranov = shafranov;
43+
plasmaType = plasma_type;
44+
plasmaId = plasma_id;
45+
numberOfBins = number_of_bins;
46+
minToroidalAngle = min_toroidal_angle/180.*M_PI;
47+
maxToroidalAngle = max_toroidal_angle/180.*M_PI;
48+
49+
setup_plasma_source();
50+
}
51+
52+
// destructor
53+
PlasmaSource::~PlasmaSource(){}
54+
55+
// main master sample function
56+
void PlasmaSource::SampleSource(std::array<double,8> random_numbers,
57+
double &x,
58+
double &y,
59+
double &z,
60+
double &u,
61+
double &v,
62+
double &w,
63+
double &E) {
64+
double radius = 0.;
65+
int bin = 0;
66+
sample_source_radial(random_numbers[0],random_numbers[1],radius,bin);
67+
double r = 0.;
68+
convert_rad_to_rz(radius,random_numbers[2],r,z);
69+
convert_r_to_xy(r,random_numbers[3],x,y);
70+
sample_energy(bin,random_numbers[4],random_numbers[5],E);
71+
isotropic_direction(random_numbers[6],random_numbers[7],
72+
u,v,w);
73+
74+
}
75+
76+
/*
77+
* sample the pdf src_profile, to generate the sampled minor radius
78+
*/
79+
void PlasmaSource::sample_source_radial(double rn_store1, double rn_store2,
80+
double &sampled_radius, int &sampled_bin) {
81+
82+
for ( int i = 0 ; i < numberOfBins ; i++ ) {
83+
if ( rn_store1 <= source_profile[i] ) {
84+
if ( i > 0 ) {
85+
sampled_radius = (float(i-1)*(binWidth)) + (binWidth*(rn_store2));
86+
sampled_bin = i;
87+
return;
88+
} else {
89+
sampled_radius = binWidth*(rn_store2);
90+
sampled_bin = i;
91+
return;
92+
}
93+
}
94+
}
95+
96+
std::cerr << "error" << std::endl;
97+
std::cerr << "Sample position greater than plasma radius" << std::endl;
98+
exit(1);
99+
return;
100+
}
101+
102+
/*
103+
* sample the energy of the neutrons, updates energy neutron in mev
104+
*/
105+
void PlasmaSource::sample_energy(const int bin_number, double random_number1, double random_number2,
106+
double &energy_neutron) {
107+
// generate the normally distributed number
108+
const double twopi = 6.28318530718;
109+
double sample1 = std::sqrt(-2.0*std::log(random_number1));
110+
double sample2 = cos(twopi*(random_number2));
111+
energy_neutron = (5.59/2.35)*(ion_kt[bin_number])*sample1*sample2;
112+
energy_neutron += 14.08;
113+
// test energy limit
114+
// if (energy_neutron < 15.5){energy_neutron = 15.5} else {}
115+
return;
116+
}
117+
118+
/*
119+
* convert the sampled radius to an rz coordinate by using plasma parameters
120+
*/
121+
void PlasmaSource::convert_rad_to_rz( const double minor_sampled,
122+
const double rn_store,
123+
double &radius, double &height)
124+
{
125+
const double twopi = 6.28318530718;
126+
127+
double alpha = twopi*(rn_store);
128+
129+
double shift = shafranov*(1.0-std::pow(minor_sampled/(minorRadius),2));
130+
131+
radius = majorRadius + minor_sampled*cos(alpha+(triangularity*sin(alpha))) + shift;
132+
height = elongation*minor_sampled*sin(alpha);
133+
134+
135+
return;
136+
}
137+
138+
139+
/*
140+
* convert rz_to_xyz
141+
*/
142+
void PlasmaSource::convert_r_to_xy(const double r, const double rn_store,
143+
double &x, double &y)
144+
145+
{
146+
double toroidal_extent = maxToroidalAngle - minToroidalAngle;
147+
double toroidal_angle = toroidal_extent*rn_store + minToroidalAngle;
148+
x = r*sin(toroidal_angle);
149+
y = r*cos(toroidal_angle);
150+
return;
151+
}
152+
153+
/*
154+
* sets up the cumulatitive probability profile
155+
* on the basis of the ion temp and ion density
156+
* this portion is deterministic
157+
*/
158+
void PlasmaSource::setup_plasma_source()
159+
{
160+
double ion_d; // ion density
161+
double ion_t; // ion temp
162+
double sig_dt; // dt xs
163+
164+
std::vector<double> src_strength; // the source strength, n/m3
165+
double r;
166+
167+
binWidth = minorRadius/float(numberOfBins);
168+
double total = 0.0; // total source strength
169+
170+
for (int i = 0 ; i < numberOfBins ; i++) {
171+
r = binWidth * float(i);
172+
ion_d = ion_density(r);
173+
ion_t = ion_temperature(r);
174+
src_strength.push_back(std::pow(ion_d,2)*dt_xs(ion_t));
175+
ion_kt.push_back(sqrt(ion_t/1000.0)); // convert to sqrt(MeV)
176+
total += src_strength[i];
177+
}
178+
179+
// normalise the source profile
180+
double sum = 0 ;
181+
for ( int i = 0 ; i < numberOfBins ; i++) {
182+
sum += src_strength[i];
183+
source_profile.push_back(sum/total);
184+
}
185+
return;
186+
}
187+
188+
/*
189+
* function that returns the ion density given the
190+
* given the critical plasma parameters
191+
*/
192+
double PlasmaSource::ion_density(const double sample_radius)
193+
{
194+
double ion_dens = 0.0;
195+
196+
if( plasmaId == 0 ) {
197+
ion_dens = ionDensityOrigin*
198+
(1.0-std::pow(sample_radius/minorRadius,2));
199+
} else {
200+
if(sample_radius <= pedistalRadius) {
201+
ion_dens += ionDensityPedistal;
202+
double product;
203+
product = 1.0-std::pow(sample_radius/pedistalRadius,2);
204+
product = std::pow(product,ionDensityPeaking);
205+
ion_dens += (ionDensityOrigin-ionDensityPedistal)*
206+
(product);
207+
} else {
208+
ion_dens += ionDensitySeperatrix;
209+
double product;
210+
product = ionDensityPedistal - ionDensitySeperatrix;
211+
ion_dens += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius);
212+
}
213+
}
214+
215+
return ion_dens;
216+
}
217+
218+
/*
219+
* function that returns the ion density given the
220+
* given the critical plasma parameters
221+
*/
222+
double PlasmaSource::ion_temperature(const double sample_radius)
223+
{
224+
double ion_temp = 0.0;
225+
226+
if( plasmaId == 0 ) {
227+
ion_temp = ionTemperatureOrigin*
228+
(1.0-std::pow(sample_radius/minorRadius,
229+
ionTemperaturePeaking));
230+
} else {
231+
if(sample_radius <= pedistalRadius) {
232+
ion_temp += ionTemperaturePedistal;
233+
double product;
234+
product = 1.0-std::pow(sample_radius/pedistalRadius,2);
235+
product = std::pow(product,ionTemperaturePeaking);
236+
ion_temp += (ionTemperatureOrigin-
237+
ionTemperaturePedistal)*(product);
238+
} else {
239+
ion_temp += ionTemperatureSeperatrix;
240+
double product;
241+
product = ionTemperaturePedistal - ionTemperatureSeperatrix;
242+
ion_temp += product*(minorRadius-sample_radius)/(minorRadius-pedistalRadius);
243+
}
244+
}
245+
246+
return ion_temp;
247+
}
248+
249+
/*
250+
* returns the dt cross section for a given ion temp
251+
*/
252+
double PlasmaSource::dt_xs(double ion_temp)
253+
{
254+
double dt;
255+
double c[7]={2.5663271e-18,19.983026,2.5077133e-2,
256+
2.5773408e-3,6.1880463e-5,6.6024089e-2,
257+
8.1215505e-3};
258+
259+
double u = 1.0-ion_temp*(c[2]+ion_temp*(c[3]-c[4]*ion_temp))
260+
/(1.0+ion_temp*(c[5]+c[6]*ion_temp));
261+
262+
dt = c[0]/(std::pow(u,5./6.)*std::pow(ion_temp,2./3.0));
263+
dt *= exp(-1.*c[1]*std::pow(u/ion_temp,1./3.));
264+
265+
return dt;
266+
}
267+
268+
/*
269+
* returns the dt cross section for a given ion temp
270+
*/
271+
void PlasmaSource::isotropic_direction(const double random1,
272+
const double random2,
273+
double &u, double &v,
274+
double &w) {
275+
double t = 2*M_PI*random1;
276+
double p = acos(1. - 2.*random2);
277+
278+
u = sin(p)*cos(t);
279+
v = sin(p)*sin(t);
280+
w = cos(p);
281+
282+
return;
283+
}
284+
285+
} // end of namespace

0 commit comments

Comments
 (0)