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 "MyTRIMRasterizer.h"
10 : #include "PKAGeneratorBase.h"
11 : #include "MooseRandom.h"
12 : #include "MooseMesh.h"
13 :
14 : // libmesh includes
15 : #include "libmesh/quadrature.h"
16 : #include "libmesh/parallel_algebra.h"
17 :
18 : #include <type_traits>
19 :
20 : // custom data load and data store methods for a struct with an std::vector member
21 : template <>
22 : inline void
23 57060 : dataStore(std::ostream & stream, MyTRIMRasterizer::AveragedData & ad, void * context)
24 : {
25 57060 : dataStore(stream, ad._elements, context);
26 57060 : dataStore(stream, ad._site_volume, context);
27 57060 : }
28 : template <>
29 : inline void
30 57060 : dataLoad(std::istream & stream, MyTRIMRasterizer::AveragedData & ad, void * context)
31 : {
32 57060 : dataLoad(stream, ad._elements, context);
33 57060 : dataLoad(stream, ad._site_volume, context);
34 57060 : }
35 :
36 : // custom data load and data store methods for a class with virtual members (vtable pointer must not
37 : // be (un)serialized)
38 : template <>
39 : inline void
40 872421 : dataStore(std::ostream & stream, MyTRIM_NS::IonBase & pl, void * context)
41 : {
42 872421 : dataStore(stream, pl._Z, context);
43 872421 : dataStore(stream, pl._m, context);
44 872421 : dataStore(stream, pl._E, context);
45 872421 : dataStore(stream, pl._dir, context);
46 872421 : dataStore(stream, pl._pos, context);
47 872421 : dataStore(stream, pl._seed, context);
48 872421 : dataStore(stream, pl._gen, context);
49 872421 : dataStore(stream, pl._id, context);
50 872421 : dataStore(stream, pl._tag, context);
51 872421 : dataStore(stream, pl._Ef, context);
52 872421 : dataStore(stream, pl._state, context);
53 872421 : }
54 : template <>
55 : inline void
56 872421 : dataLoad(std::istream & stream, MyTRIM_NS::IonBase & pl, void * context)
57 : {
58 872421 : dataLoad(stream, pl._Z, context);
59 872421 : dataLoad(stream, pl._m, context);
60 872421 : dataLoad(stream, pl._E, context);
61 872421 : dataLoad(stream, pl._dir, context);
62 872421 : dataLoad(stream, pl._pos, context);
63 872421 : dataLoad(stream, pl._seed, context);
64 872421 : dataLoad(stream, pl._gen, context);
65 872421 : dataLoad(stream, pl._id, context);
66 872421 : dataLoad(stream, pl._tag, context);
67 872421 : dataLoad(stream, pl._Ef, context);
68 872421 : dataLoad(stream, pl._state, context);
69 872421 : }
70 :
71 : registerMooseObject("MagpieApp", MyTRIMRasterizer);
72 :
73 : InputParameters
74 324 : MyTRIMRasterizer::validParams()
75 : {
76 324 : InputParameters params = ElementUserObject::validParams();
77 324 : params.addClassDescription("Gather the element distribution of the simulation domain for a TRIM "
78 : "binary collision Monte Carlo simulation");
79 :
80 648 : MooseEnum var_physical_meaning("STOICHIOMETRY NUMBER_DENSITY", "STOICHIOMETRY");
81 648 : params.addParam<MooseEnum>("var_physical_meaning",
82 : var_physical_meaning,
83 : "The physical meaning of the rasterizer variables.");
84 :
85 648 : params.addRequiredCoupledVar("var", "Variables to rasterize");
86 648 : params.addCoupledVar("periodic_var",
87 : "Optional variables that determines the periodicity. If not supplied the "
88 : "first argument of 'var' will be used.");
89 648 : params.addRequiredParam<std::vector<Real>>("M", "Element mass in amu");
90 648 : params.addRequiredParam<std::vector<Real>>("Z", "Nuclear charge in e");
91 648 : params.addParam<std::vector<Real>>("Mtol",
92 : "Tolerance on mass number for tagging PKAs with var id.");
93 648 : params.addParam<std::vector<Real>>("Ebind", {}, "Binding energy in eV");
94 648 : params.addParam<std::vector<Real>>("Edisp", {}, "Displacement threshold in eV");
95 648 : params.addParam<MaterialPropertyName>(
96 : "site_volume", "Lattice site volume in nm^3 (regardless of the chosen mesh units)");
97 648 : params.addRequiredParam<std::vector<UserObjectName>>("pka_generator",
98 : "List of PKA generating user objects");
99 324 : ExecFlagEnum setup_options(MooseUtils::getDefaultExecFlagEnum());
100 :
101 : // we run this object once a timestep
102 324 : setup_options = EXEC_TIMESTEP_BEGIN;
103 324 : params.set<ExecFlagEnum>("execute_on") = setup_options;
104 :
105 : // which TRIM Module to run for optional capabilities like energy deposition
106 648 : MooseEnum trim_module_options("CORE=0 ENERGY_DEPOSITION=1", "CORE");
107 648 : params.addParam<MooseEnum>("trim_module",
108 : trim_module_options,
109 : "TRIM Module to run for optional capabilities like energy deposition");
110 :
111 : // which units of length to use
112 648 : MooseEnum length_unit_options("ANGSTROM=0 NANOMETER MICROMETER", "ANGSTROM");
113 648 : params.addParam<MooseEnum>("length_unit",
114 : length_unit_options,
115 : "Length units of the MOOSE mesh. MyTRIM contains pretabulated "
116 : "crossection data with units so this option must be set correctly to "
117 : "obtain physical results.");
118 :
119 : // Advanced options
120 648 : params.addParam<unsigned int>("interval", 1, "The time step interval at which TRIM BCMC is run");
121 648 : params.addParam<Real>("analytical_energy_cutoff",
122 648 : 0.0,
123 : "Energy cutoff in eV below which recoils are not followed explicitly but "
124 : "effects are calculated analytically.");
125 648 : params.addParam<Real>("recoil_rate_scaling",
126 648 : 1.0,
127 : "A factor to scale computed reaction rates in the the PKAGenerator "
128 : "objects. This is useful to avoid extremely large PKA lists.");
129 648 : params.addParam<unsigned int>(
130 : "max_pka_count", "Desired number of PKAs to be run during each invocation of mytrim");
131 972 : params.addRangeCheckedParam<Real>("max_nrt_difference",
132 648 : 0.2,
133 : "max_nrt_difference > 0 & max_nrt_difference < 1",
134 : "The largest max-norm difference between number fractions for "
135 : "reusing existing polyatomic NRT.");
136 648 : params.addParamNamesToGroup("interval analytical_energy_cutoff max_pka_count", "Advanced");
137 :
138 648 : params.addParam<Real>("r_rec",
139 : "Recombination radius in Angstrom. Frenkel pairs with a separation "
140 : "distance lower than this will be removed from the cascade");
141 648 : params.addParamNamesToGroup("r_rec", "Recombination");
142 :
143 324 : return params;
144 648 : }
145 :
146 179 : MyTRIMRasterizer::MyTRIMRasterizer(const InputParameters & parameters)
147 : : ElementUserObject(parameters),
148 179 : _nvars(coupledComponents("var")),
149 179 : _dim(_mesh.dimension()),
150 179 : _var(_nvars),
151 179 : _site_volume_prop(nullptr),
152 358 : _pka_generator_names(getParam<std::vector<UserObjectName>>("pka_generator")),
153 179 : _pka_generators(),
154 358 : _periodic(isCoupled("periodic_var") ? coupled("periodic_var", 0) : coupled("var", 0)),
155 179 : _accumulated_time(0.0),
156 179 : _accumulated_time_old(0.0),
157 358 : _interval(getParam<unsigned int>("interval")),
158 537 : _perf_finalize(registerTimedSection("finalize", 2))
159 : {
160 179 : if (_nvars == 0)
161 0 : mooseError("Must couple variables to MyTRIMRasterizer.");
162 :
163 : // fill in trim parameters
164 179 : _trim_parameters.element_prototypes.resize(_nvars);
165 358 : _trim_parameters.analytical_cutoff = getParam<Real>("analytical_energy_cutoff");
166 358 : _trim_parameters.recoil_rate_scaling = getParam<Real>("recoil_rate_scaling");
167 358 : _trim_parameters.max_nrt_distance = getParam<Real>("max_nrt_difference");
168 179 : _trim_parameters.nrt_log_energy_spacing = 1.1;
169 358 : _trim_parameters.trim_module = getParam<MooseEnum>("trim_module").getEnum<TRIMModuleEnum>();
170 358 : if (isParamValid("max_pka_count"))
171 20 : _trim_parameters.desired_npka = getParam<unsigned int>("max_pka_count");
172 : else
173 169 : _trim_parameters.desired_npka = 0;
174 :
175 358 : auto trim_M = getParam<std::vector<Real>>("M");
176 537 : auto trim_Z = getParam<std::vector<Real>>("Z");
177 :
178 : std::vector<Real> mtol;
179 358 : if (isParamValid("Mtol"))
180 15 : mtol = getParam<std::vector<Real>>("Mtol");
181 : else
182 174 : mtol.assign(_nvars, 0.5);
183 :
184 179 : if (trim_M.size() != _nvars)
185 0 : mooseError("Parameter 'M' must have as many components as coupled variables.");
186 179 : if (trim_Z.size() != _nvars)
187 0 : mooseError("Parameter 'Z' must have as many components as coupled variables.");
188 179 : if (mtol.size() != _nvars)
189 0 : mooseError("Parameter 'mtol' must have as many components as coupled variables.");
190 :
191 : // error check masses and charges
192 436 : for (unsigned int i = 0; i < _nvars; ++i)
193 : {
194 257 : if (trim_Z[i] > trim_M[i])
195 0 : mooseError("Value of Z is larger than value of M for entry ", i);
196 257 : if (trim_Z[i] >= _pka_parameters._index_Z.size())
197 0 : mooseError("Value of Z is too large. Maximum Z supported is ",
198 0 : _pka_parameters._index_Z.size() - 1,
199 : " but one element has Z=",
200 : i);
201 : }
202 :
203 436 : for (unsigned int i = 0; i < _nvars; ++i)
204 : {
205 257 : _var[i] = &coupledValue("var", i);
206 257 : _trim_parameters.element_prototypes[i]._m = trim_M[i];
207 257 : _trim_parameters.element_prototypes[i]._Z = trim_Z[i];
208 : }
209 :
210 537 : auto trim_Ebind = getParam<std::vector<Real>>("Ebind");
211 179 : if (trim_Ebind.size() == _nvars)
212 39 : for (unsigned int i = 0; i < _nvars; ++i)
213 26 : _trim_parameters.element_prototypes[i]._Elbind = trim_Ebind[i];
214 166 : else if (trim_Ebind.empty())
215 397 : for (unsigned int i = 0; i < _nvars; ++i)
216 231 : _trim_parameters.element_prototypes[i]._Elbind = 3.0;
217 : else
218 0 : mooseError("Parameter 'Ebind' must have as many components as coupled variables (or left empty "
219 : "for a default of 3eV).");
220 :
221 537 : auto trim_Edisp = getParam<std::vector<Real>>("Edisp");
222 179 : if (trim_Edisp.size() == _nvars)
223 109 : for (unsigned int i = 0; i < _nvars; ++i)
224 61 : _trim_parameters.element_prototypes[i]._Edisp = trim_Edisp[i];
225 131 : else if (trim_Edisp.empty())
226 327 : for (unsigned int i = 0; i < _nvars; ++i)
227 196 : _trim_parameters.element_prototypes[i]._Edisp = 25.0;
228 : else
229 0 : mooseError("Parameter 'Edisp' must have as many components as coupled variables (or left empty "
230 : "for a default of 25eV).");
231 :
232 358 : if (isParamValid("r_rec"))
233 : {
234 10 : _trim_parameters.recombination = true;
235 20 : _trim_parameters.r_rec = getParam<Real>("r_rec");
236 : }
237 : else
238 169 : _trim_parameters.recombination = false;
239 :
240 : // fetch PKA Generators
241 358 : for (auto && name : _pka_generator_names)
242 179 : _pka_generators.push_back(&getUserObjectByName<PKAGeneratorBase>(name));
243 :
244 : // set up data for sample periodicity
245 623 : for (unsigned int i = 0; i < _dim; ++i)
246 : {
247 444 : _pbc[i] = _mesh.isRegularOrthogonal() && _mesh.isTranslatedPeriodic(_periodic, i);
248 :
249 444 : if (_pbc[i])
250 : {
251 229 : _min_dim(i) = _mesh.getMinInDimension(i);
252 229 : _max_dim(i) = _mesh.getMaxInDimension(i);
253 : }
254 : }
255 :
256 : // determine length scale factor for TRIM
257 358 : switch (getParam<MooseEnum>("length_unit").getEnum<Unit>())
258 : {
259 169 : case ANGSTROM:
260 169 : _trim_parameters.length_scale = 1.0;
261 169 : break;
262 :
263 5 : case NANOMETER:
264 5 : _trim_parameters.length_scale = 10.0;
265 5 : break;
266 :
267 5 : case MICROMETER:
268 5 : _trim_parameters.length_scale = 10000.0;
269 5 : break;
270 :
271 0 : default:
272 0 : mooseError("Unknown length unit.");
273 : }
274 :
275 537 : if (getParam<MooseEnum>("var_physical_meaning") == "STOICHIOMETRY")
276 : {
277 328 : if (!isParamValid("site_volume"))
278 0 : mooseError("Rasterizer variables are stoiciometric contents, site_volume must be provided.");
279 328 : _site_volume_prop = &getMaterialProperty<Real>("site_volume");
280 : }
281 : else
282 15 : _site_volume_conversion = Utility::pow<3>(_trim_parameters.length_scale) * 1e-3;
283 :
284 : // setup invariant PKA generation parameters
285 21659 : for (auto & nZ : _pka_parameters._index_Z)
286 : nZ = std::make_pair(0, 0);
287 179 : _pka_parameters._mass_charge_tuple.resize(_nvars);
288 179 : _pka_parameters._recoil_rate_scaling = _trim_parameters.recoil_rate_scaling;
289 436 : for (unsigned int i = 0; i < _nvars; ++i)
290 : {
291 257 : const auto Z = _trim_parameters.element_prototypes[i]._Z;
292 :
293 : // insert (mass, charge) pair
294 : _pka_parameters._mass_charge_tuple[i] =
295 : std::make_tuple(_trim_parameters.element_prototypes[i]._m, Z, mtol[i]);
296 :
297 : // increase the count of elements with the same Z
298 257 : auto & index_Z = _pka_parameters._index_Z[Z];
299 257 : index_Z.first++;
300 :
301 : // only set this to the first index (important for ionTag())
302 257 : if (index_Z.first == 1)
303 252 : index_Z.second = i;
304 : }
305 179 : }
306 :
307 : bool
308 541 : MyTRIMRasterizer::executeThisTimestep() const
309 : {
310 541 : return (_fe_problem.timeStep() - 1) % _interval == 0;
311 : }
312 :
313 : void
314 219 : MyTRIMRasterizer::initialize()
315 : {
316 219 : _execute_this_timestep = executeThisTimestep();
317 :
318 : // We roll back the accumulated time time if the preceeding timestep did
319 : // not converge
320 219 : if (!_fe_problem.nlConverged(/*nl_sys_num=*/0))
321 161 : _accumulated_time = _accumulated_time_old;
322 :
323 219 : if (_execute_this_timestep)
324 : {
325 219 : _trim_parameters.last_executed_dt = _fe_problem.dt();
326 : _material_map.clear();
327 219 : _pka_list.clear();
328 : }
329 :
330 : /// setup global PKA parameters for the current timestep
331 219 : _pka_parameters._dt = _accumulated_time + _fe_problem.dt();
332 219 : }
333 :
334 : void
335 185043 : MyTRIMRasterizer::execute()
336 : {
337 : // bail out early if not executing this timestep
338 185043 : if (!_execute_this_timestep)
339 0 : return;
340 :
341 : // average element concentrations
342 :
343 185043 : AveragedData average(_nvars);
344 : Real vol = 0.0;
345 :
346 : // average material data over elements
347 1160307 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
348 : {
349 975264 : const Real qpvol = _JxW[qp] * _coord[qp];
350 975264 : vol += qpvol;
351 :
352 : // average compositions on the element
353 2561064 : for (unsigned int i = 0; i < _nvars; ++i)
354 1585800 : average._elements[i] += qpvol * (*_var[i])[qp];
355 :
356 : // average site volume property
357 975264 : if (_site_volume_prop)
358 945264 : average._site_volume += qpvol * (*_site_volume_prop)[qp];
359 : else
360 30000 : average._site_volume += qpvol * _site_volume_conversion;
361 : }
362 :
363 : // divide by total element volume
364 185043 : if (vol > 0.0)
365 : {
366 450153 : for (unsigned int i = 0; i < _nvars; ++i)
367 265110 : average._elements[i] /= vol;
368 :
369 185043 : average._site_volume /= vol;
370 : }
371 :
372 : // store in map
373 185043 : _material_map[_current_elem->id()] = average;
374 :
375 : // update corrent element volume
376 185043 : _pka_parameters._volume = vol;
377 :
378 : // add PKAs for current element
379 370086 : for (auto && gen : _pka_generators)
380 185043 : gen->appendPKAs(_pka_list, _pka_parameters, average);
381 : }
382 :
383 : void
384 42 : MyTRIMRasterizer::threadJoin(const UserObject & y)
385 : {
386 : // if the map needs to be updated we merge the maps from all threads
387 42 : if (_execute_this_timestep)
388 : {
389 : const MyTRIMRasterizer & uo = static_cast<const MyTRIMRasterizer &>(y);
390 42 : _material_map.insert(uo._material_map.begin(), uo._material_map.end());
391 42 : _pka_list.insert(_pka_list.end(), uo._pka_list.begin(), uo._pka_list.end());
392 : }
393 42 : }
394 :
395 : void
396 177 : MyTRIMRasterizer::finalize()
397 : {
398 177 : TIME_SECTION(_perf_finalize);
399 :
400 : // save the accumulated time so that we can properly roll back if the step does not converge
401 177 : _accumulated_time_old = _accumulated_time;
402 :
403 : // bail out early if not executing this timestep
404 177 : if (!_execute_this_timestep)
405 : {
406 : // no BCMC done, so wee accumulate this step's time to be taken care of by a later BCMC run
407 0 : _accumulated_time += _fe_problem.dt();
408 : return;
409 : }
410 :
411 : // BCMC was run for the accumulated time - the debt is paid
412 177 : _accumulated_time = 0.0;
413 :
414 : // for single processor runs we do not need to do anything here
415 177 : if (_app.n_processors() > 1)
416 : {
417 : // create send buffer
418 : std::string send_buffer;
419 :
420 : // create byte buffers for the streams received from all processors
421 : std::vector<std::string> recv_buffers;
422 :
423 : // pack the complex datastructures into the string stream
424 84 : serialize(send_buffer);
425 :
426 : // broadcast serialized data to and receive from all processors
427 84 : _communicator.allgather(send_buffer, recv_buffers);
428 :
429 : // unpack the received data and merge it into the local data structures
430 84 : deserialize(recv_buffers);
431 84 : }
432 :
433 : // we will assign random seeds on proc 0 & reject on proc 0. To guarantee reproducibility
434 : // we need to sort the PKA list
435 177 : if (processor_id() == 0)
436 : {
437 135 : std::vector<unsigned int> pka_seeds(_pka_list.size());
438 :
439 : // only do this on proc 0 thread 0
440 1338621 : for (auto i = beginIndex(_pka_list); i < _pka_list.size(); ++i)
441 1338486 : pka_seeds[i] = MooseRandom::randl();
442 :
443 : // sort PKA list only on processor 0 & assign random number seeds
444 135 : std::sort(_pka_list.begin(),
445 : _pka_list.end(),
446 27322563 : [](MyTRIM_NS::IonBase a, MyTRIM_NS::IonBase b)
447 : {
448 28164169 : return (a._pos < b._pos) || (a._pos == b._pos && a._m < b._m) ||
449 27815824 : (a._pos == b._pos && a._m == b._m && a._E < b._E) ||
450 493261 : (a._pos == b._pos && a._m == b._m && a._E == b._E && a._Z < b._Z);
451 : });
452 :
453 : // store seeds in tag values
454 1338621 : for (auto i = beginIndex(_pka_list); i < _pka_list.size(); ++i)
455 1338486 : _pka_list[i]._seed = pka_seeds[i];
456 : }
457 :
458 : // rejection is performed in processor 0 only
459 177 : _trim_parameters.original_npka = _pka_list.size();
460 177 : if (_trim_parameters.desired_npka == 0 ||
461 : _trim_parameters.desired_npka > _trim_parameters.original_npka)
462 : {
463 169 : _trim_parameters.scaled_npka = _trim_parameters.original_npka;
464 169 : _trim_parameters.result_scaling_factor = 1.0 / _trim_parameters.recoil_rate_scaling;
465 : }
466 : else
467 : {
468 8 : if (processor_id() == 0)
469 : {
470 6 : Real acceptance_probability =
471 6 : Real(_trim_parameters.desired_npka) / Real(_trim_parameters.original_npka);
472 :
473 : // most straight-forward but probably inefficient implementation of rejection
474 6 : std::vector<MyTRIM_NS::IonBase> old_pka_list = _pka_list;
475 6 : _pka_list.resize(0);
476 54006 : for (auto & p : old_pka_list)
477 54000 : if (MooseRandom::rand() < acceptance_probability)
478 6309 : _pka_list.push_back(p);
479 :
480 : // save the size of the PKA list after rejection & the scaling factor
481 6 : _trim_parameters.scaled_npka = _pka_list.size();
482 6 : _trim_parameters.result_scaling_factor = Real(_trim_parameters.original_npka) /
483 6 : Real(_trim_parameters.scaled_npka) /
484 6 : _trim_parameters.recoil_rate_scaling;
485 6 : }
486 :
487 : // need to broadcast the size of the PKA list after rejection & result scaling factor
488 8 : _communicator.broadcast(_trim_parameters.scaled_npka);
489 8 : _communicator.broadcast(_trim_parameters.result_scaling_factor);
490 : }
491 :
492 : // communicate the PKA list if n_proc > 1
493 177 : if (_app.n_processors() > 1)
494 : {
495 : std::string pka_list_buffer;
496 84 : if (processor_id() == 0)
497 : {
498 : // pack the local _pka_list into a string buffer
499 42 : std::ostringstream oss;
500 42 : dataStore(oss, _pka_list, this);
501 42 : pka_list_buffer.assign(oss.str());
502 42 : }
503 :
504 : // communicate the pka list
505 84 : _communicator.broadcast(pka_list_buffer);
506 84 : if (processor_id() != 0)
507 : {
508 42 : _pka_list.resize(0);
509 42 : std::istringstream iss(pka_list_buffer);
510 42 : dataLoad(iss, _pka_list, this);
511 42 : }
512 : }
513 :
514 : // prune PKA list
515 177 : if (_app.n_processors() > 1)
516 : {
517 : // split PKAs into per-processor ranges
518 84 : std::vector<unsigned int> interval(_app.n_processors() + 1, 0);
519 252 : for (unsigned int i = 0; i < _app.n_processors(); ++i)
520 168 : interval[i + 1] = (_pka_list.size() - interval[i]) / (_app.n_processors() - i) + interval[i];
521 :
522 84 : auto begin = interval[processor_id()];
523 84 : auto end = interval[processor_id() + 1];
524 84 : std::vector<MyTRIM_NS::IonBase> own_pka_list(end - begin);
525 :
526 428346 : for (auto i = begin; i < end; ++i)
527 428262 : own_pka_list[i - begin] = _pka_list[i];
528 :
529 : _pka_list.swap(own_pka_list);
530 84 : }
531 : }
532 :
533 : const std::vector<Real> &
534 205490 : MyTRIMRasterizer::material(const Elem * elem) const
535 : {
536 205490 : auto i = _material_map.find(elem->id());
537 :
538 : // there should be data for every element in the mesh
539 205490 : if (i == _material_map.end())
540 0 : mooseError("Element not found in material map.");
541 :
542 205490 : return i->second._elements;
543 : }
544 :
545 : Real
546 201377 : MyTRIMRasterizer::siteVolume(const Elem * elem) const
547 : {
548 201377 : auto i = _material_map.find(elem->id());
549 :
550 : // there should be data for every element in the mesh
551 201377 : if (i == _material_map.end())
552 0 : mooseError("Element not found in material map.");
553 :
554 201377 : return i->second._site_volume;
555 : }
556 :
557 : Point
558 76249602 : MyTRIMRasterizer::periodicPoint(const Point & pos) const
559 : {
560 : // point to sample the material at
561 76249602 : Point p(pos(0), pos(1), _dim == 2 ? 0.0 : pos(2));
562 :
563 : // apply periodic boundary conditions
564 297736677 : for (unsigned int i = 0; i < _dim; ++i)
565 221487075 : if (_pbc[i])
566 : {
567 201485847 : const Real width = _max_dim(i) - _min_dim(i);
568 201485847 : p(i) -= std::floor((p(i) - _min_dim(i)) / width) * width;
569 : }
570 :
571 76249602 : return p;
572 : }
573 :
574 : bool
575 4113 : MyTRIMRasterizer::isTrackedSpecies(unsigned int atomic_number, Real mass_number) const
576 : {
577 : mooseAssert(_pka_generators[0], "PKA generator is not set.");
578 4113 : return _pka_generators[0]->ionTag(_pka_parameters, atomic_number, mass_number) != -1;
579 : }
580 :
581 : void
582 84 : MyTRIMRasterizer::serialize(std::string & serialized_buffer)
583 : {
584 : // stream for serializing the _material_map and _pka_list data structure to a byte stream
585 84 : std::ostringstream oss;
586 84 : dataStore(oss, _material_map, this);
587 84 : dataStore(oss, _pka_list, this);
588 :
589 : // Populate the passed in string pointer with the string stream's buffer contents
590 84 : serialized_buffer.assign(oss.str());
591 84 : }
592 :
593 : void
594 84 : MyTRIMRasterizer::deserialize(std::vector<std::string> & serialized_buffers)
595 : {
596 : mooseAssert(serialized_buffers.size() == _app.n_processors(),
597 : "Unexpected size of serialized_buffers: " << serialized_buffers.size());
598 :
599 : // The input string stream used for deserialization
600 84 : std::istringstream iss;
601 :
602 : // Loop over all datastructures for all procerssors to perfrom the gather operation
603 252 : for (unsigned int rank = 0; rank < serialized_buffers.size(); ++rank)
604 : {
605 : // skip the current processor (its data is already in the structures)
606 168 : if (rank == processor_id())
607 84 : continue;
608 :
609 : // populate the stream with a new buffer and reset stream state
610 : iss.str(serialized_buffers[rank]);
611 84 : iss.clear();
612 :
613 : // Load the communicated data into temporary structures
614 : MaterialMap other_material_map;
615 84 : dataLoad(iss, other_material_map, this);
616 : std::vector<MyTRIM_NS::IonBase> other_pka_list;
617 84 : dataLoad(iss, other_pka_list, this);
618 :
619 : // merge the data in with the current processor's data
620 84 : _material_map.insert(other_material_map.begin(), other_material_map.end());
621 :
622 : // merging the PKA lists
623 84 : _pka_list.insert(_pka_list.begin(), other_pka_list.begin(), other_pka_list.end());
624 84 : }
625 84 : }
|
