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 :
9 : #ifdef GSL_ENABLED
10 :
11 : #include "DPAUserObjectBase.h"
12 : #include "PolyatomicDisplacementFunction.h"
13 :
14 : InputParameters
15 32 : DPAUserObjectBase::validParams()
16 : {
17 32 : InputParameters params = GeneralUserObject::validParams();
18 64 : params.addParam<Real>("irradiation_time", "Irradiation time used ");
19 32 : MultiMooseEnum damage_reaction_types("elastic inelastic");
20 64 : params.addParam<MultiMooseEnum>("damage_reaction_types",
21 : damage_reaction_types,
22 : "The neutron reaction causing radiation damage");
23 64 : params.addParam<std::vector<Real>>("Z", {}, "The atomic numbers of all present isotopes");
24 64 : params.addParam<std::vector<Real>>("A", {}, "The mass numbers of all present isotopes");
25 64 : params.addParam<std::vector<Real>>(
26 : "number_densities", {}, "The number densities of all present isotopes");
27 64 : params.addParam<std::vector<Real>>("energy_group_boundaries",
28 : {},
29 : "The neutron flux energy group boundaries in units of eV "
30 : "starting with the highest energy group");
31 64 : params.addParam<std::vector<Real>>("scalar_flux",
32 : {},
33 : "The values of the neutron scalar flux by energy group "
34 : "starting with the highest energy group");
35 64 : params.addParam<std::vector<std::vector<Real>>>(
36 : "cross_section",
37 : "The values of the cross sections. Each semicolon separated vector contains cross sections "
38 : "for a particular nuclide and reaction type. Each entry must be number of energy groups "
39 : "entries long. One vector must be provided for each "
40 : "nuclide/reaction type combination. The ordering is as follows: if there are reaction types"
41 : "a and b, and nuclides i, k, and l, the ordering will be xs_ai; xs_ak; xs_al; xs_bi; xs_bk, "
42 : "xs_bl");
43 64 : params.addParam<std::vector<std::vector<Real>>>(
44 : "Q", "The Q values by reaction type and then isotope. Assumed zero if not provided.");
45 32 : params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
46 32 : params.suppressParameter<ExecFlagEnum>("execute_on");
47 32 : return params;
48 32 : }
49 :
50 16 : DPAUserObjectBase::DPAUserObjectBase(const InputParameters & parameters)
51 : : GeneralUserObject(parameters),
52 32 : _is_transient_irradiation(!isParamValid("irradiation_time")),
53 32 : _irradiation_time(_is_transient_irradiation ? 0 : getParam<Real>("irradiation_time")),
54 32 : _neutron_reaction_types(getParam<MultiMooseEnum>("damage_reaction_types")),
55 16 : _nr(_neutron_reaction_types.size()),
56 32 : _atomic_numbers(getParam<std::vector<Real>>("Z")),
57 32 : _mass_numbers(getParam<std::vector<Real>>("A")),
58 32 : _number_densities(getParam<std::vector<Real>>("number_densities")),
59 32 : _energy_group_boundaries(getParam<std::vector<Real>>("energy_group_boundaries")),
60 64 : _scalar_flux(getParam<std::vector<Real>>("scalar_flux"))
61 : {
62 16 : if (_nr == 0)
63 0 : paramError("damage_reaction_types", "At least one damage mechanism must be provided");
64 :
65 32 : if (isParamValid("cross_section"))
66 : {
67 48 : std::vector<std::vector<Real>> xs = getParam<std::vector<std::vector<Real>>>("cross_section");
68 :
69 : // we know _nr so the length of cross sections should be divisible by _nr
70 16 : unsigned int n = xs.size();
71 16 : if (n % _nr != 0)
72 0 : paramError("cross_section",
73 : "The number of entries in the cross_section param must be divisible by number of "
74 : "reaction types ",
75 : _nr);
76 :
77 16 : unsigned int n_nuc = n / _nr;
78 16 : _cross_sections.resize(_nr);
79 : unsigned int p = 0;
80 40 : for (unsigned int r = 0; r < _nr; ++r)
81 : {
82 24 : _cross_sections[r].resize(n_nuc);
83 64 : for (unsigned int j = 0; j < n_nuc; ++j)
84 : {
85 40 : _cross_sections[r][j] = xs[p];
86 40 : ++p;
87 : }
88 : }
89 16 : }
90 16 : }
91 :
92 : void
93 16 : DPAUserObjectBase::prepare()
94 : {
95 : // This object support dynamic changes of number densities
96 : // appearance and disappearance of isotopes, and change in
97 : // cross sections. Therefore, this function needs to check
98 : // and redo a lot of things that would usually be done in
99 : // the constructor.
100 :
101 : // checks on scalar flux & energy groups
102 16 : _ng = _scalar_flux.size();
103 16 : if (_ng == 0)
104 0 : mooseError("The number of provided scalar fluxes is zero. Provide using parameter scalar_flux "
105 : "or set programmatically");
106 :
107 16 : if (_energy_group_boundaries.size() != _ng + 1)
108 0 : mooseError("The size of the energy group boundary is ",
109 0 : _energy_group_boundaries.size(),
110 : " but it must have exactly number of energy groups plus one entry ",
111 0 : _ng + 1);
112 :
113 120 : for (unsigned int j = 1; j < _energy_group_boundaries.size(); ++j)
114 104 : if (_energy_group_boundaries[j - 1] <= _energy_group_boundaries[j])
115 0 : mooseError("Entries of energy_group_boundaries must be strictly decreasing");
116 :
117 : // check on A, Z, number densities
118 16 : initAZNHelper();
119 :
120 : // checks on cross sections
121 16 : if (_cross_sections.size() != _nr)
122 0 : mooseError("Leading dimension of _cross_sections is ",
123 0 : _cross_sections.size(),
124 : " but should be equal to number of damage reaction types ",
125 0 : _nr);
126 40 : for (unsigned int i = 0; i < _nr; ++i)
127 : {
128 24 : if (_cross_sections[i].size() != _atomic_numbers.size())
129 0 : mooseError("Second dimension of _cross_sections for index ",
130 : i,
131 : " does not match number of isotope ",
132 0 : _atomic_numbers.size());
133 64 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
134 40 : if (_cross_sections[i][j].size() != _ng)
135 0 : mooseError("Third dimension of _cross_sections for indices ",
136 : i,
137 : " ",
138 : j,
139 : " has ",
140 0 : _cross_sections[i][j].size(),
141 : " entries when number of groups is ",
142 0 : _ng);
143 : }
144 :
145 : // get Q values, need to do this here to be sure that number of isotopes is set
146 32 : if (isParamValid("Q"))
147 : {
148 24 : _q_values = getParam<std::vector<std::vector<Real>>>("Q");
149 8 : if (_q_values.size() != _nr)
150 0 : paramError("Q", "Leading dimension must be of the same size as reaction types.");
151 24 : for (unsigned int j = 0; j < _nr; ++j)
152 16 : if (_q_values[j].size() != _atomic_numbers.size())
153 0 : paramError("Q", "Second dimension must be of same size as isotopes.");
154 : }
155 : else
156 : {
157 8 : _q_values.resize(_nr);
158 16 : for (unsigned int j = 0; j < _nr; ++j)
159 8 : _q_values[j].resize(_atomic_numbers.size());
160 : }
161 16 : }
162 :
163 : void
164 32 : DPAUserObjectBase::initAZNHelper()
165 : {
166 32 : unsigned int ns = getAtomicNumbers().size();
167 32 : if (ns == 0)
168 0 : mooseError("The size of the atomic number array is zero. Set Z parameter or set atomic number "
169 : "programmatically");
170 :
171 32 : if (_mass_numbers.size() != ns)
172 0 : mooseError("Size of mass number array does not match size of atomic number array. Size of mass "
173 : "numbers ",
174 0 : _mass_numbers.size(),
175 : ", size of atomic numbers ",
176 : ns);
177 :
178 32 : if (_number_densities.size() != ns)
179 0 : mooseError("Size of mass number array does not match size of atomic number array. Size of mass "
180 : "numbers ",
181 0 : _number_densities.size(),
182 : ", size of atomic numbers ",
183 : ns);
184 32 : }
185 :
186 : std::vector<unsigned int>
187 36 : DPAUserObjectBase::getAtomicNumbers() const
188 : {
189 : // in this kind we need to work a bit since
190 : // VPPs are Reals and we need to return unsigned int
191 36 : std::vector<unsigned int> Zs(_atomic_numbers.size());
192 92 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
193 : {
194 56 : unsigned int i = static_cast<unsigned int>(_atomic_numbers[j]);
195 56 : if (std::abs(_atomic_numbers[j] - i) > 1e-12)
196 0 : mooseError("Entry ", j, ":", _atomic_numbers[j], " is not a non-negative integer.");
197 56 : Zs[j] = i;
198 : }
199 36 : return Zs;
200 : }
201 :
202 : std::vector<Real>
203 4 : DPAUserObjectBase::getMassNumbers() const
204 : {
205 4 : return _mass_numbers;
206 : }
207 :
208 : std::vector<Real>
209 4 : DPAUserObjectBase::getNumberFractions() const
210 : {
211 : // need to normalize here
212 : Real sum = 0;
213 12 : for (unsigned int j = 0; j < _number_densities.size(); ++j)
214 8 : sum += _number_densities[j];
215 4 : std::vector<Real> nf(_number_densities.size());
216 12 : for (unsigned int j = 0; j < _number_densities.size(); ++j)
217 8 : nf[j] = _number_densities[j] / sum;
218 4 : return nf;
219 : }
220 :
221 : Real
222 28 : DPAUserObjectBase::getMaxEnergy() const
223 : {
224 : // between elastic & inelastic scattering, elastic scattering
225 : // allows the largest recoil energy = gamma * E_neutron
226 28 : Real max_neutron_energy = _energy_group_boundaries[0];
227 :
228 : // find maximum gamma
229 : Real min_mass = std::numeric_limits<Real>::max();
230 76 : for (unsigned int j = 0; j < _mass_numbers.size(); ++j)
231 48 : if (min_mass > _mass_numbers[j])
232 : min_mass = _mass_numbers[j];
233 28 : Real gamma = 4 * min_mass / (min_mass + 1) / (min_mass + 1);
234 28 : return gamma * max_neutron_energy;
235 : }
236 :
237 : bool
238 16 : DPAUserObjectBase::changed() const
239 : {
240 : // the size could have changed
241 0 : if (_atomic_numbers.size() != _atomic_numbers_old.size() ||
242 16 : _mass_numbers.size() != _mass_numbers_old.size() ||
243 : _number_densities.size() != _number_densities_old.size())
244 : return true;
245 :
246 : // the size is still the same
247 0 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
248 0 : if (!MooseUtils::absoluteFuzzyEqual(_atomic_numbers[j], _atomic_numbers_old[j], _tol) ||
249 0 : !MooseUtils::absoluteFuzzyEqual(_mass_numbers[j], _mass_numbers_old[j], _tol) ||
250 : !MooseUtils::absoluteFuzzyEqual(_number_densities[j], _number_densities_old[j], _tol))
251 : return true;
252 :
253 : return false;
254 : }
255 :
256 : void
257 16 : DPAUserObjectBase::accumulateDamage()
258 : {
259 16 : if (changed())
260 : {
261 16 : initAZNHelper();
262 16 : onCompositionChanged();
263 :
264 : // copy over A, Z, N to A, Z, N_old
265 16 : _atomic_numbers_old.resize(_atomic_numbers.size());
266 16 : _mass_numbers_old.resize(_atomic_numbers.size());
267 16 : _number_densities_old.resize(_atomic_numbers.size());
268 40 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
269 : {
270 24 : _atomic_numbers_old[j] = _atomic_numbers[j];
271 24 : _mass_numbers_old[j] = _mass_numbers[j];
272 24 : _number_densities_old[j] = _number_densities[j];
273 : }
274 : }
275 :
276 : // compute total number density
277 : Real total_number_density = 0;
278 40 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
279 24 : total_number_density += _number_densities[j];
280 :
281 : // damage function have been computed
282 : // so now we can start computing dpa
283 16 : if (!_is_transient_irradiation)
284 : {
285 16 : _dpa = 0;
286 : _partial_dpa.clear();
287 : }
288 :
289 16 : Real del_t = _is_transient_irradiation ? _dt : _irradiation_time;
290 :
291 40 : for (unsigned int x = 0; x < _nr; ++x)
292 : {
293 200 : for (unsigned int g = 0; g < _ng; ++g)
294 : {
295 176 : Real del_E = _energy_group_boundaries[g] - _energy_group_boundaries[g + 1];
296 :
297 496 : for (unsigned int i = 0; i < _atomic_numbers.size(); ++i)
298 928 : for (unsigned int j = 0; j < _atomic_numbers.size(); ++j)
299 : {
300 608 : Real v = del_t / del_E * _scalar_flux[g] * _number_densities[i] / total_number_density *
301 608 : _cross_sections[x][i][g] * neutronDamageEfficiency(i, j, g, x);
302 :
303 : // total dpa is incremented regardless
304 608 : _dpa += v;
305 :
306 : // partial dpa: if the entry exists, add to it, otherwise create it
307 1216 : const auto & p = _partial_dpa.find(zaidHelper(_atomic_numbers[i], _mass_numbers[i]));
308 608 : if (p == _partial_dpa.end())
309 24 : _partial_dpa[zaidHelper(_atomic_numbers[i], _mass_numbers[i])] = v;
310 : else
311 584 : p->second += v;
312 : }
313 : }
314 : }
315 16 : }
316 :
317 : Real
318 608 : DPAUserObjectBase::neutronDamageEfficiency(unsigned int i,
319 : unsigned int j,
320 : unsigned int g,
321 : unsigned int x) const
322 : {
323 : std::vector<Real> points;
324 : std::vector<Real> weights;
325 608 : PolyatomicDisplacementFunction::gslQuadRule(100, points, weights);
326 :
327 1216 : if (_neutron_reaction_types[x] == "elastic")
328 : {
329 320 : Real gamma = 4 * _mass_numbers[i] / (_mass_numbers[i] + 1) / (_mass_numbers[i] + 1);
330 320 : Real from_E = _energy_group_boundaries[g + 1];
331 320 : Real to_E = _energy_group_boundaries[g];
332 320 : Real delta_E = to_E - from_E;
333 :
334 : // integral over E
335 : Real integral = 0;
336 32320 : for (unsigned int qp = 0; qp < points.size(); ++qp)
337 : {
338 32000 : Real energy = 0.5 * (points[qp] + 1) * delta_E + from_E;
339 32000 : Real w = 0.5 * weights[qp] * delta_E;
340 32000 : integral += w * integralDamageFunction(gamma * energy, i, j) / energy;
341 : }
342 :
343 320 : return integral / gamma;
344 : }
345 576 : else if (_neutron_reaction_types[x] == "inelastic")
346 : {
347 : // warn if q value does not make sense
348 288 : if (_q_values[x][i] >= 0)
349 144 : mooseDoOnce(mooseWarning("Q value is greater or equal to zero for inelastic scattering."));
350 :
351 288 : Real d = std::abs(_q_values[x][i]) * (_mass_numbers[i] + 1) / _mass_numbers[i];
352 288 : Real gamma = 4 * _mass_numbers[i] / (_mass_numbers[i] + 1) / (_mass_numbers[i] + 1);
353 288 : Real from_E = _energy_group_boundaries[g + 1];
354 288 : Real to_E = _energy_group_boundaries[g];
355 288 : Real delta_E = to_E - from_E;
356 :
357 : // integral over E
358 : Real integral = 0;
359 29088 : for (unsigned int qp = 0; qp < points.size(); ++qp)
360 : {
361 28800 : Real energy = 0.5 * (points[qp] + 1) * delta_E + from_E;
362 28800 : Real Delta = d / energy;
363 :
364 : // make sure threshold energy is honored
365 28800 : if (Delta > 1)
366 13440 : continue;
367 :
368 15360 : Real sr = std::sqrt(1 - Delta);
369 15360 : Real Tmin = gamma / 2 * energy * (1 - sr - 0.5 * Delta);
370 15360 : Real Tmax = gamma / 2 * energy * (1 + sr - 0.5 * Delta);
371 15360 : Real w = 0.5 * weights[qp] * delta_E;
372 15360 : integral += w * (integralDamageFunction(Tmax, i, j) - integralDamageFunction(Tmin, i, j)) /
373 15360 : (energy * sr);
374 : }
375 :
376 288 : return integral / gamma;
377 : }
378 :
379 0 : mooseError("Neutron reaction type not recognized. Should never get here.");
380 : return 0;
381 : }
382 :
383 : std::string
384 632 : DPAUserObjectBase::zaidHelper(unsigned int Z, Real A) const
385 : {
386 632 : std::stringstream ss;
387 : ss << std::setprecision(4) << Z << A;
388 632 : return ss.str();
389 632 : }
390 :
391 : Real
392 0 : DPAUserObjectBase::getPartialDPA(unsigned int Z, Real A) const
393 : {
394 0 : const auto & p = _partial_dpa.find(zaidHelper(Z, A));
395 0 : if (p == _partial_dpa.end())
396 : return 0;
397 0 : return p->second;
398 : }
399 :
400 : #endif
|
