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 : #include "ThreadedRecoilLoopBase.h"
10 : #include "MooseMyTRIMCore.h"
11 : #include "MooseMyTRIMEnergyDeposition.h"
12 : #include "MooseMyTRIMSample.h"
13 : #include "MooseMesh.h"
14 :
15 : // libmesh includes
16 : #include "libmesh/quadrature.h"
17 :
18 : // Remove after idaholab/moose#26102
19 : #include "MagpieNanoflann.h"
20 :
21 : // Specialization for PointListAdaptor<MyTRIMDefectBufferItem>
22 : template <>
23 : inline const Point &
24 : PointListAdaptor<ThreadedRecoilLoopBase::MyTRIMDefectBufferItem>::getPoint(
25 : const ThreadedRecoilLoopBase::MyTRIMDefectBufferItem & item) const
26 : {
27 : return item.point;
28 : }
29 :
30 161 : ThreadedRecoilLoopBase::ThreadedRecoilLoopBase(const MyTRIMRasterizer & rasterizer,
31 161 : const MooseMesh & mesh)
32 161 : : _rasterizer(rasterizer),
33 161 : _trim_parameters(_rasterizer.getTrimParameters()),
34 161 : _nvars(_trim_parameters.nVars()),
35 161 : _mesh(mesh),
36 161 : _dim(_mesh.dimension())
37 : {
38 161 : _simconf.setLengthScale(_trim_parameters.length_scale);
39 161 : }
40 :
41 : // Splitting Constructor
42 37 : ThreadedRecoilLoopBase::ThreadedRecoilLoopBase(const ThreadedRecoilLoopBase & x,
43 37 : Threads::split /*split*/)
44 37 : : _rasterizer(x._rasterizer),
45 37 : _trim_parameters(x._trim_parameters),
46 37 : _nvars(x._nvars),
47 37 : _mesh(x._mesh),
48 37 : _dim(x._dim)
49 : {
50 37 : _simconf.setLengthScale(_trim_parameters.length_scale);
51 37 : }
52 :
53 : void
54 198 : ThreadedRecoilLoopBase::operator()(const PKARange & pka_list)
55 : {
56 : // fetch a point locator
57 396 : _pl = _mesh.getPointLocator();
58 :
59 : // permit querying points that are potentially outside the mesh
60 198 : _pl->enable_out_of_mesh_mode();
61 :
62 : // create a new sample class to bridge the MOOSE mesh and the MyTRIM domain
63 198 : MooseMyTRIMSample sample(_rasterizer, _mesh, &_simconf);
64 :
65 : // create a FIFO for recoils
66 198 : std::queue<MyTRIM_NS::IonBase *> recoils;
67 :
68 : // create a list potentially used for energy deposition
69 : std::list<std::pair<Point, Real>> edep_list;
70 :
71 : // build the requested TRIM module
72 198 : std::unique_ptr<MooseMyTRIMCore> TRIM;
73 198 : switch (_trim_parameters.trim_module)
74 : {
75 : // basic module with interstitial and vacancy generation
76 149 : case MyTRIMRasterizer::MYTRIM_CORE:
77 149 : TRIM.reset(new MooseMyTRIMCore(&_simconf, &sample));
78 : break;
79 :
80 : // record energy deposited to the lattice
81 49 : case MyTRIMRasterizer::MYTRIM_ENERGY_DEPOSITION:
82 49 : TRIM.reset(new MooseMyTRIMEnergyDeposition(&_simconf, &sample, edep_list));
83 : break;
84 :
85 0 : default:
86 0 : mooseError("Unknown TRIM module.");
87 : }
88 :
89 : // for instantaneous recombination
90 : const unsigned int requested_results = 10;
91 198 : std::vector<std::size_t> return_index(requested_results);
92 198 : std::vector<Real> return_dist_sqr(requested_results);
93 :
94 : // copy the pka list into the recoil queue
95 90924 : for (auto && pka : pka_list)
96 : {
97 : // seed the RNG with the seed assigned to the primary knock on atom (for parallel
98 : // reproducibility)
99 90726 : _simconf.seed(pka._seed);
100 :
101 : // push primary knock on atom onto the recoil queue
102 90726 : recoils.push(new MyTRIM_NS::IonBase(pka));
103 :
104 : // clear cascade recombination buffers
105 90726 : _interstitial_buffer.clear();
106 90726 : _vacancy_buffer.clear();
107 :
108 : MyTRIM_NS::IonBase * recoil;
109 5685663 : while (!recoils.empty())
110 : {
111 5594937 : recoil = recoils.front();
112 : recoils.pop();
113 5594937 : sample.averages(recoil);
114 :
115 : // project into xy plane
116 5594937 : if (_dim == 2)
117 589620 : recoil->_pos(2) = 0.0;
118 :
119 : // deal with what this recoil left behind
120 5594937 : if (recoil->_state == MyTRIM_NS::IonBase::VACANCY)
121 4270494 : _vacancy_buffer.push_back(MyTRIMDefectBufferItem(recoil->_pos, recoil->_tag));
122 1324443 : else if (recoil->_state == MyTRIM_NS::IonBase::REPLACEMENT)
123 1233717 : addDefectToResult(
124 2467434 : _rasterizer.periodicPoint(recoil->_pos), recoil->_tag, 1, REPLACEMENT_OUT);
125 : // ignore like-for-like substitutions
126 :
127 : // remove zero energy recoils
128 5594937 : if (recoil->_E > 0.0)
129 : {
130 : // full recoil or analytical approximation
131 5599050 : if (recoil->_E < _trim_parameters.analytical_cutoff &&
132 4113 : _rasterizer.isTrackedSpecies(recoil->_Z, recoil->_m))
133 : {
134 4113 : mooseDoOnce(mooseWarning("Skipping detailed cascade calculation below cutoff energy."));
135 4113 : const auto pp = _rasterizer.periodicPoint(recoil->_pos);
136 :
137 : #ifdef GSL_ENABLED
138 4113 : const auto elem = (*_pl)(pp);
139 :
140 : // the actual number density is immaterial for NRT, only number fractions are important
141 : // so we go ahead and first get it from the rasterizer and then normalize it to sum to 1
142 4113 : auto number_fractions = _rasterizer.material(elem);
143 :
144 : Real total_number_density = 0;
145 8916 : for (unsigned int j = 0; j < number_fractions.size(); ++j)
146 4803 : total_number_density += number_fractions[j];
147 :
148 8916 : for (unsigned int j = 0; j < number_fractions.size(); ++j)
149 4803 : number_fractions[j] /= total_number_density;
150 :
151 : unsigned int index;
152 : Real distance;
153 :
154 4113 : findBestNRTMatch(number_fractions, index, distance);
155 4113 : if (distance > _trim_parameters.max_nrt_distance)
156 21 : index = addNRTEntry(number_fractions);
157 :
158 : // look up entry for recoil Z & m
159 4113 : unsigned int species_index = _pa_nrt[index]->findSpeciesIndex(recoil->_Z, recoil->_m);
160 :
161 : // interpolate the replacement counts for all considered target species
162 : // vacancy energy is the energy of the cascade that goes into creating vacancies
163 : // vacancy_energy = sum_{targets: j} Ebind_j * n_{projectile, j}
164 : Real vacancy_energy = 0;
165 4113 : if (species_index != libMesh::invalid_uint)
166 8916 : for (unsigned int target_var = 0; target_var < _nvars; ++target_var)
167 : {
168 4803 : unsigned int target_species_index = _pa_nrt[index]->findSpeciesIndex(
169 4803 : _trim_parameters.element_prototypes[target_var]._Z,
170 4803 : _trim_parameters.element_prototypes[target_var]._m);
171 : // using linear perturbation theory estimate of g_ij(number_fractions)
172 4803 : Real w = _pa_nrt[index]->linearInterpolation(
173 : recoil->_E, species_index, target_species_index);
174 10986 : for (unsigned int l = 0; l < _pa_derivative_nrt[index]->nSpecies(); ++l)
175 6183 : w += _pa_derivative_nrt[index]->linearInterpolation(
176 6183 : recoil->_E, species_index, target_species_index, l) *
177 6183 : (number_fractions[l] - _pa_nrt[index]->numberFraction(l));
178 :
179 : // increment energy, vacancy and interstitial buffers
180 4803 : vacancy_energy += w * _trim_parameters.element_prototypes[target_var]._Elbind;
181 4803 : _vacancy_buffer.emplace_back(MyTRIMDefectBufferItem(recoil->_pos, target_var, w));
182 4803 : _interstitial_buffer.emplace_back(
183 4803 : MyTRIMDefectBufferItem(recoil->_pos, target_var, w));
184 : }
185 : else
186 0 : mooseError("NRT treatment did not find recoil with ", recoil->_Z, " ", recoil->_m);
187 :
188 : // add remaining recoil energy
189 4113 : if (_trim_parameters.trim_module == MyTRIMRasterizer::MYTRIM_ENERGY_DEPOSITION)
190 3423 : addEnergyToResult(pp, recoil->_E - vacancy_energy);
191 : #else
192 : mooseDoOnce(
193 : mooseWarning("GSL is not enabled and cutoff energy != 0. Polyatomic NRT requires "
194 : "GSL, current settings do not account for skipped ions."));
195 : // add remaining recoil energy
196 : if (_trim_parameters.trim_module == MyTRIMRasterizer::MYTRIM_ENERGY_DEPOSITION)
197 : addEnergyToResult(pp, recoil->_E);
198 : #endif
199 : }
200 : else
201 : {
202 : // follow this ion's trajectory and store recoils
203 5590824 : TRIM->trim(recoil, recoils);
204 :
205 : // are we tracking atoms of this type?
206 5590824 : if (recoil->_tag >= 0)
207 : {
208 5540343 : if (recoil->_state == MyTRIM_NS::IonBase::INTERSTITIAL)
209 4314117 : _interstitial_buffer.emplace_back(MyTRIMDefectBufferItem(recoil->_pos, recoil->_tag));
210 1226226 : else if (recoil->_state == MyTRIM_NS::IonBase::REPLACEMENT)
211 934383 : addDefectToResult(
212 1868766 : _rasterizer.periodicPoint(recoil->_pos), recoil->_tag, 1, REPLACEMENT_IN);
213 : // BUG: untracked PKAs in replacement collisions will cause a REPLACEMENT_OUT without a
214 : // REPLACEMENT_IN!
215 : }
216 :
217 : // store energy deposition
218 6349179 : for (auto & edep : edep_list)
219 758355 : addEnergyToResult(_rasterizer.periodicPoint(edep.first), edep.second);
220 : edep_list.clear();
221 : }
222 : }
223 :
224 : // done with this recoil
225 5594937 : delete recoil;
226 : }
227 :
228 : //
229 : // Process instantaneous recombination of this PKA's defects
230 : // Recombination only takes place within the cascade of an individual PKA.
231 : // Cascades are assumed to be non-overlapping in time and space.
232 : ///
233 90726 : if (_trim_parameters.recombination && !_vacancy_buffer.empty())
234 : {
235 : // 1. build kd-tree for the vacancies
236 : const unsigned int max_leaf_size = 50; // slightly affects runtime
237 : auto point_list =
238 : PointListAdaptor<MyTRIMDefectBufferItem>(_vacancy_buffer.begin(), _vacancy_buffer.end());
239 : auto kd_tree = std::make_unique<KDTreeType>(
240 1851 : LIBMESH_DIM, point_list, nanoflann::KDTreeSingleIndexAdaptorParams(max_leaf_size));
241 :
242 : mooseAssert(kd_tree != nullptr, "KDTree was not properly initialized.");
243 1851 : kd_tree->buildIndex();
244 :
245 1851 : const Real r_rec2 = _trim_parameters.r_rec * _trim_parameters.r_rec;
246 : nanoflann::SearchParameters params;
247 :
248 : // 2. iterate over interstitials and recombine them if they are with r_rec of a vacancy
249 : std::vector<nanoflann::ResultItem<std::size_t, Real>> ret_matches;
250 70002 : for (auto & i : _interstitial_buffer)
251 : {
252 : ret_matches.clear();
253 68151 : std::size_t n_result = kd_tree->radiusSearch(&(i.point(0)), r_rec2, ret_matches, params);
254 :
255 68517 : for (std::size_t j = 0; j < n_result; ++j)
256 : {
257 4974 : auto & v = _vacancy_buffer[ret_matches[j].first];
258 :
259 : // only allow interstitial to go into vacancy of the same type
260 4974 : if (v.variable_id == i.variable_id)
261 : {
262 : // mark vacancy-interstitial pair for deletion
263 4608 : i.variable_id = libMesh::invalid_uint;
264 4608 : v.variable_id = libMesh::invalid_uint;
265 4608 : break;
266 : }
267 : }
268 : }
269 1851 : }
270 :
271 : // add remaining defects to result
272 4409646 : for (auto & i : _interstitial_buffer)
273 4318920 : if (i.variable_id != libMesh::invalid_uint)
274 4314312 : addDefectToResult(
275 8628624 : _rasterizer.periodicPoint(i.point), i.variable_id, i.weight, INTERSTITIAL);
276 4366023 : for (auto & v : _vacancy_buffer)
277 4275297 : if (v.variable_id != libMesh::invalid_uint)
278 4270689 : addDefectToResult(_rasterizer.periodicPoint(v.point), v.variable_id, v.weight, VACANCY);
279 : }
280 396 : }
281 :
282 : #ifdef GSL_ENABLED
283 : unsigned int
284 21 : ThreadedRecoilLoopBase::addNRTEntry(const std::vector<Real> & number_fractions)
285 : {
286 21 : if (number_fractions.size() != _nvars)
287 0 : mooseError("number_fractions has wrong size");
288 :
289 : std::vector<MyTRIM_NS::Element> poly_mat;
290 48 : for (unsigned int j = 0; j < _nvars; ++j)
291 : {
292 27 : MyTRIM_NS::Element element;
293 27 : element._Z = _trim_parameters.element_prototypes[j]._Z;
294 27 : element._m = _trim_parameters.element_prototypes[j]._m;
295 27 : element._t = number_fractions[j];
296 27 : element._Edisp = _trim_parameters.element_prototypes[j]._Edisp;
297 27 : element._Elbind = _trim_parameters.element_prototypes[j]._Elbind;
298 27 : poly_mat.push_back(element);
299 : }
300 :
301 : /**
302 : * NOTE: using the net displacement rate here defined as
303 : * "[...] atoms displaced and not recaptured in subsequent replacement collisions"
304 : * TODO: this should be checked for accuracy, also net displacement rate is not as well
305 : * verified as total displacement rate
306 : */
307 42 : _pa_nrt.push_back(std::make_unique<PolyatomicDisplacementFunction>(poly_mat, NET));
308 42 : _pa_derivative_nrt.push_back(std::make_unique<PolyatomicDisplacementDerivativeFunction>(
309 42 : poly_mat, NET_DERIVATIVE, _pa_nrt.back().get()));
310 :
311 : // integrate P&K's equation to the analytical cutoff
312 : Real energy = _pa_nrt.back()->minEnergy();
313 :
314 : // some baby steps [10 eV] to get initial evolution right
315 21 : Real upper_linear_range = std::max(2 * energy, 100.0);
316 21 : unsigned int n_initial_steps = std::ceil((upper_linear_range - energy) / 10);
317 177 : for (unsigned int j = 0; j < n_initial_steps; ++j)
318 : {
319 156 : energy += 10;
320 156 : _pa_nrt.back()->advanceDisplacements(energy);
321 : }
322 :
323 : // constant in log(E) stepping
324 : for (;;)
325 : {
326 462 : energy *= _trim_parameters.nrt_log_energy_spacing;
327 462 : if (energy > _trim_parameters.analytical_cutoff)
328 : {
329 21 : _pa_nrt.back()->advanceDisplacements(_trim_parameters.analytical_cutoff);
330 : break;
331 : }
332 441 : _pa_nrt.back()->advanceDisplacements(energy);
333 : }
334 :
335 639 : for (unsigned int n = 1; n < _pa_nrt.back()->nEnergySteps(); ++n)
336 618 : _pa_derivative_nrt.back()->advanceDisplacements(_pa_nrt.back()->energyPoint(n));
337 :
338 21 : return _pa_nrt.size() - 1;
339 : }
340 :
341 : void
342 4113 : ThreadedRecoilLoopBase::findBestNRTMatch(const std::vector<Real> & number_fractions,
343 : unsigned int & index,
344 : Real & distance) const
345 : {
346 4113 : index = 0;
347 4113 : distance = std::numeric_limits<Real>::max();
348 :
349 4113 : if (_pa_nrt.size() == 0)
350 : return;
351 :
352 8853 : for (unsigned int i = 0; i < _pa_nrt.size(); ++i)
353 : {
354 : auto & nrt_candidate = _pa_nrt[i];
355 :
356 : mooseAssert(number_fractions.size() == nrt_candidate->nSpecies(),
357 : "Number densities have different sizes");
358 :
359 : // find the infinity norm difference
360 : Real inf_norm_difference = 0;
361 10866 : for (unsigned int j = 0; j < number_fractions.size(); ++j)
362 6108 : if (std::abs(number_fractions[j] - nrt_candidate->numberFraction(j)) > inf_norm_difference)
363 : inf_norm_difference = std::abs(number_fractions[j] - nrt_candidate->numberFraction(j));
364 :
365 4758 : if (distance > inf_norm_difference)
366 : {
367 4344 : distance = inf_norm_difference;
368 4344 : index = i;
369 : }
370 : }
371 : }
372 : #endif
|
