Line data Source code
1 : /**********************************************************************/
2 : /* DO NOT MODIFY THIS HEADER */
3 : /* MAGPIE - Mesoscale Atomistic Glue Program for Integrated Execution */
4 : /* */
5 : /* Copyright 2017 Battelle Energy Alliance, LLC */
6 : /* ALL RIGHTS RESERVED */
7 : /**********************************************************************/
8 : #ifdef GSL_ENABLED
9 :
10 : #include "ScatteringRecoilCrossSection.h"
11 : #include "MagpieUtils.h"
12 : #include "Function.h"
13 :
14 : // gsl includes
15 : #include "gsl/gsl_sf_legendre.h"
16 : #include "gsl/gsl_integration.h"
17 :
18 : InputParameters
19 96 : ScatteringRecoilCrossSection::validParams()
20 : {
21 96 : InputParameters params = GeneralUserObject::validParams();
22 96 : params.addClassDescription(
23 : "Base class for computing recoil scattering cross sections for a single isotope."
24 : "It outputs the coefficients for the Legendre expansion of the cross section up to the "
25 : "specified Legendre expansion order"
26 : "polynomials l, neutron energy groups g and recoil energy bins t. It also outputs the "
27 : "maximum and mininum"
28 : "cosines of the recoil atom in the laboratory frame.");
29 192 : params.addParam<unsigned int>("quadrature_order", 400, "Quadrature order for integration");
30 192 : params.addParam<Real>("atomic_mass", 1, "Atomic Mass of the isotope. Default to Hydrogen A = 1");
31 192 : params.addRequiredParam<std::vector<Real>>(
32 : "neutron_energy_limits",
33 : "Energy limits of the incident neutron in [eV] and descending order");
34 192 : params.addRequiredParam<std::vector<Real>>(
35 : "recoil_energy_limits", "Energy limits of the recoil atom in [eV] and descending order");
36 192 : params.addRequiredParam<FunctionName>("neutron_spectrum",
37 : "Function representing the reactor neutron spectrum");
38 192 : params.addRequiredParam<std::vector<FunctionName>>("scattering_xs",
39 : "Functions representing the neutron cross "
40 : "sections. NOTE: this is a vector to allow "
41 : "separate inputs for inelastic modes.");
42 192 : params.addParam<unsigned int>(
43 192 : "legendre_order", 10, "Order of Legendre polynomials where n = 0, ..., 10. Default to P10");
44 192 : params.addParam<std::string>(
45 : "cross_section_output_filename",
46 : "Name of the output file with the cross section coefficients (.csv)");
47 192 : params.addParam<std::string>("mu_L_output_filename",
48 : "Name of the output file with the mininum and maximum mu_L (.csv)");
49 96 : return params;
50 0 : }
51 :
52 48 : ScatteringRecoilCrossSection::ScatteringRecoilCrossSection(const InputParameters & parameters)
53 : : GeneralUserObject(parameters),
54 48 : _csv_tolerance(1.0e-10),
55 48 : _quad_order(getParam<unsigned int>("quadrature_order")),
56 48 : _neutron_spectrum(getFunction("neutron_spectrum")),
57 96 : _L(getParam<unsigned int>("legendre_order")),
58 96 : _atomic_mass(getParam<Real>("atomic_mass")),
59 96 : _neutron_energy_limits(getParam<std::vector<Real>>("neutron_energy_limits")),
60 192 : _recoil_energy_limits(getParam<std::vector<Real>>("recoil_energy_limits"))
61 : {
62 144 : if (isParamValid("cross_section_output_filename") ^ isParamValid("mu_L_output_filename"))
63 0 : mooseError("cross_section_output_filename and mu_L_output_filename must either both be present "
64 : "or absent");
65 :
66 : // get scattering cross sections
67 144 : std::vector<FunctionName> xs_names = getParam<std::vector<FunctionName>>("scattering_xs");
68 48 : _scattering_cross_section.resize(xs_names.size());
69 96 : for (unsigned int j = 0; j < xs_names.size(); ++j)
70 48 : _scattering_cross_section[j] = &getFunctionByName(xs_names[j]);
71 :
72 : // set up integration rule
73 48 : auto * qp_table = gsl_integration_glfixed_table_alloc(_quad_order);
74 48 : _quad_points.resize(_quad_order);
75 48 : _quad_weights.resize(_quad_order);
76 20848 : for (std::size_t j = 0; j < _quad_order; ++j)
77 : {
78 : double point, weight;
79 20800 : gsl_integration_glfixed_point(-1.0, 1.0, j, &point, &weight, qp_table);
80 20800 : _quad_points[j] = point;
81 20800 : _quad_weights[j] = weight;
82 : }
83 48 : gsl_integration_glfixed_table_free(qp_table);
84 :
85 48 : _alpha = std::pow(((_atomic_mass - 1) / (_atomic_mass + 1)), 2);
86 48 : _gamma = 4 * _atomic_mass / std::pow((_atomic_mass + 1), 2);
87 48 : _G = _neutron_energy_limits.size() - 1;
88 48 : _T = _recoil_energy_limits.size() - 1;
89 48 : }
90 :
91 : void
92 12 : ScatteringRecoilCrossSection::initialize()
93 : {
94 : // Calculate neutron spectrum over group g
95 12 : _xi_g.resize(_G);
96 24 : for (unsigned int g = 0; g < _G; ++g)
97 : {
98 12 : Real E_u = _neutron_energy_limits[g];
99 12 : Real E_l = _neutron_energy_limits[g + 1];
100 :
101 4812 : for (unsigned int p = 0; p < _quad_points.size(); ++p)
102 : {
103 4800 : Real E = 0.5 * (E_u - E_l) * _quad_points[p] + 0.5 * (E_u + E_l);
104 4800 : Real w_E = 0.5 * _quad_weights[p] * (E_u - E_l);
105 4800 : _xi_g[g] += w_E * _neutron_spectrum.value(E, Point());
106 : }
107 : }
108 :
109 : // Size the cross section array
110 12 : _recoil_coefficient.resize(_L + 1);
111 56 : for (unsigned int l = 0; l < _L + 1; ++l)
112 : {
113 44 : _recoil_coefficient[l].resize(_T);
114 228 : for (unsigned int t = 0; t < _T; ++t)
115 184 : _recoil_coefficient[l][t].assign(_G, 0.0);
116 : }
117 :
118 : // Size the lab frame cosine array and initialize to -1, both lab cosines -1 means
119 : // this g->t combination is impossible
120 12 : _mu_L.resize(_T);
121 68 : for (unsigned int t = 0; t < _T; ++t)
122 : {
123 56 : _mu_L[t].resize(_G);
124 112 : for (unsigned int g = 0; g < _G; ++g)
125 56 : _mu_L[t][g].assign(2, -1.0);
126 : }
127 12 : }
128 :
129 : void
130 12 : ScatteringRecoilCrossSection::finalize()
131 : {
132 48 : if (!(isParamValid("cross_section_output_filename") && isParamValid("mu_L_output_filename")))
133 0 : return;
134 24 : _recoil_xs_file_name = getParam<std::string>("cross_section_output_filename");
135 24 : _mu_L_file_name = getParam<std::string>("mu_L_output_filename");
136 :
137 : /*
138 : * Write the elastic recoil cross section (g -> t) output file
139 : *
140 : * t = 0 t = 1 t = ...
141 : * l = 0, 1, 2, ... l = 0, 1, 2, ... l = 0, 1, 2, ...
142 : * g 0
143 : * 1
144 : * 2
145 : * .
146 : * .
147 : * .
148 : */
149 12 : std::ofstream output_file;
150 12 : output_file.open(_recoil_xs_file_name);
151 :
152 : // print header
153 12 : output_file << "Neutron group,";
154 68 : for (unsigned int t = 0; t < _T; ++t)
155 240 : for (unsigned int l = 0; l < _L + 1; ++l)
156 : {
157 368 : output_file << "Recoil group t = " << t << " Moment order l = " << l;
158 184 : if (!(t == _T - 1 && l == _L))
159 172 : output_file << ",";
160 : }
161 : output_file << std::endl;
162 :
163 : // print data
164 24 : for (unsigned int g = 0; g < _G; ++g)
165 : {
166 12 : output_file << g << ",";
167 68 : for (unsigned int t = 0; t < _T; ++t)
168 240 : for (unsigned int l = 0; l < _L + 1; ++l)
169 : {
170 184 : output_file << csvPrint(_recoil_coefficient[l][t][g]);
171 184 : if (!(t == _T - 1 && l == _L))
172 172 : output_file << ',';
173 : }
174 : output_file << std::endl;
175 : }
176 12 : output_file.close();
177 :
178 : /*
179 : * Writes output file with maximum and mininum cosines in the Lab frame (mu_L)
180 : * It follows same structure as ERXS outfile file, but saves the mu_L_min and
181 : * mu_L_max per g -> t combination.
182 : */
183 12 : std::ofstream output_file2;
184 12 : output_file2.open(_mu_L_file_name);
185 : // print header
186 12 : output_file2 << "Neutron group,";
187 68 : for (unsigned int t = 0; t < _T; ++t)
188 : {
189 56 : output_file2 << "Recoil group t = " << t << " mu_min,"
190 112 : << "Recoil group t = " << t << " mu_max";
191 56 : if (t != _T - 1)
192 44 : output_file2 << ",";
193 : }
194 : output_file2 << std::endl;
195 :
196 : // print data
197 24 : for (unsigned int g = 0; g < _G; ++g)
198 : {
199 12 : output_file2 << g << ",";
200 68 : for (unsigned int t = 0; t < _T; ++t)
201 : {
202 112 : output_file2 << csvPrint(_mu_L[t][g][0]) << ',' << csvPrint(_mu_L[t][g][1]);
203 56 : if (t != _T - 1)
204 44 : output_file2 << ',';
205 : }
206 : output_file2 << std::endl;
207 : }
208 12 : output_file2.close();
209 12 : }
210 :
211 : // Find the neutron energy group given a neutron energy
212 : unsigned int
213 0 : ScatteringRecoilCrossSection::findNeutronEnergyGroup(Real energy)
214 : {
215 0 : for (unsigned int g = 0; g < _G; ++g)
216 : {
217 0 : if (energy < _neutron_energy_limits[g] && energy > _neutron_energy_limits[g + 1])
218 0 : return g;
219 : }
220 0 : mooseError("Should never get here");
221 : }
222 :
223 : Real
224 296 : ScatteringRecoilCrossSection::csvPrint(Real value) const
225 : {
226 296 : if (std::abs(value) < _csv_tolerance)
227 52 : return 0;
228 : return value;
229 : }
230 :
231 : Real
232 0 : ScatteringRecoilCrossSection::getSigmaRecoil(unsigned int g, unsigned int t, unsigned int l) const
233 : {
234 : mooseAssert(g < _G, "g is larger than _G [indexed from g = 0, ..., G - 1]");
235 : mooseAssert(t < _T, "t is larger than _T [indexed from t = 0, ..., T - 1]");
236 0 : if (l <= _L)
237 0 : return _recoil_coefficient[l][t][g];
238 : return 0.0;
239 : }
240 :
241 : Real
242 0 : ScatteringRecoilCrossSection::getMaxRecoilCosine(unsigned int g, unsigned t) const
243 : {
244 : mooseAssert(g < _G, "g is larger than _G [indexed from g = 0, ..., G - 1]");
245 : mooseAssert(t < _T, "t is larger than _T [indexed from t = 0, ..., T - 1]");
246 0 : return _mu_L[t][g][1];
247 : }
248 :
249 : Real
250 0 : ScatteringRecoilCrossSection::getMinRecoilCosine(unsigned int g, unsigned t) const
251 : {
252 : mooseAssert(g < _G, "g is larger than _G [indexed from g = 0, ..., G - 1]");
253 : mooseAssert(t < _T, "t is larger than _T [indexed from t = 0, ..., T - 1]");
254 0 : return _mu_L[t][g][0];
255 : }
256 :
257 : bool
258 0 : ScatteringRecoilCrossSection::isRecoilPossible(unsigned int g, unsigned int t) const
259 : {
260 0 : return !(_mu_L[t][g][0] == -1 && _mu_L[t][g][1] == -1);
261 : }
262 :
263 : #endif
|
