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 "MDRunBase.h"
10 : #include "MooseMesh.h"
11 :
12 : // Remove after idaholab/moose#26102
13 : #include "MagpieNanoflann.h"
14 :
15 : #include "libmesh/mesh_tools.h"
16 :
17 : // custom data load and data store methods for MDParticles class
18 : template <>
19 : inline void
20 : dataStore(std::ostream & stream, MDRunBase::MDParticles & pl, void * context)
21 : {
22 : dataStore(stream, pl.pos, context);
23 : dataStore(stream, pl.id, context);
24 : dataStore(stream, pl.elem_id, context);
25 : dataStore(stream, pl.properties, context);
26 : }
27 :
28 : template <>
29 : inline void
30 : dataLoad(std::istream & stream, MDRunBase::MDParticles & pl, void * context)
31 : {
32 : dataLoad(stream, pl.pos, context);
33 : dataLoad(stream, pl.id, context);
34 : dataLoad(stream, pl.elem_id, context);
35 : dataLoad(stream, pl.properties, context);
36 : }
37 :
38 : InputParameters
39 48 : MDRunBase::validParams()
40 : {
41 48 : InputParameters params = GeneralUserObject::validParams();
42 48 : params.set<ExecFlagEnum>("execute_on") = EXEC_TIMESTEP_BEGIN;
43 48 : params.suppressParameter<ExecFlagEnum>("execute_on");
44 96 : params.addParam<MultiMooseEnum>("md_particle_properties",
45 96 : MDRunBase::mdParticleProperties(),
46 : "Properties of MD particles to be obtained and stored.");
47 96 : params.addRangeCheckedParam<Real>(
48 : "max_granular_radius",
49 : 0,
50 : "max_granular_radius>=0",
51 : "Maximum radius of granular particles. This can be important for partitioning.");
52 48 : params.addClassDescription(
53 : "Base class for execution of coupled molecular dynamics MOOSE calculations.");
54 48 : return params;
55 0 : }
56 :
57 24 : MDRunBase::MDRunBase(const InputParameters & parameters)
58 : : GeneralUserObject(parameters),
59 24 : _properties(getParam<MultiMooseEnum>("md_particle_properties")),
60 48 : _granular(_properties.contains("radius")),
61 24 : _mesh(_subproblem.mesh()),
62 24 : _dim(_mesh.dimension()),
63 48 : _nproc(_app.n_processors())
64 : {
65 48 : _max_granular_radius = getParam<Real>("max_granular_radius");
66 : // if the calculation is granular, max_particle_radius must be set
67 48 : if (_granular && !parameters.isParamSetByUser("max_granular_radius"))
68 0 : paramError("max_granular_radius", "max_granular_radius must be set for granular calculations");
69 :
70 : // set up the map from property ID to index
71 24 : _md_particles._prop_size = _properties.size();
72 44 : for (unsigned int j = 0; j < _properties.size(); ++j)
73 20 : _md_particles._map_props[_properties.get(j)] = j;
74 :
75 : // _md_particles._r_index is a short-cut for the "radius" index
76 24 : if (_granular)
77 12 : _md_particles._r_index = propIndex("radius");
78 24 : }
79 :
80 : void
81 24 : MDRunBase::initialSetup()
82 : {
83 24 : _bbox = MeshTools::create_processor_bounding_box(_mesh, processor_id());
84 :
85 : // inflate bounding box
86 96 : for (unsigned int d = 0; d < _dim; ++d)
87 : {
88 72 : _bbox.first(d) -= _max_granular_radius;
89 72 : _bbox.second(d) += _max_granular_radius;
90 : }
91 24 : }
92 :
93 : void
94 260 : MDRunBase::timestepSetup()
95 : {
96 : // update/init the processor bounding box
97 260 : _bbox = MeshTools::create_processor_bounding_box(_mesh, processor_id());
98 :
99 : // inflate bounding box
100 1040 : for (unsigned int d = 0; d < _dim; ++d)
101 : {
102 780 : _bbox.first(d) -= _max_granular_radius;
103 780 : _bbox.second(d) += _max_granular_radius;
104 : }
105 :
106 : // callback for updating md particle list
107 260 : updateParticleList();
108 260 : }
109 :
110 : void
111 260 : MDRunBase::updateKDTree()
112 : {
113 520 : _kd_tree = std::make_unique<KDTree>(_md_particles.pos, 50);
114 260 : }
115 :
116 : void
117 35856 : MDRunBase::elemParticles(unique_id_type elem_id, std::vector<unsigned int> & elem_particles) const
118 : {
119 35856 : if (_elem_particles.find(elem_id) != _elem_particles.end())
120 4752 : elem_particles = _elem_particles.find(elem_id)->second;
121 : else
122 31104 : elem_particles = {};
123 35856 : }
124 :
125 : void
126 504 : MDRunBase::granularElementVolumes(unique_id_type elem_id,
127 : std::vector<std::pair<unsigned int, Real>> & gran_vol) const
128 : {
129 : mooseAssert(_granular,
130 : "Radius must be provided as MD property to allow granular volume computation.");
131 504 : if (_elem_granular_volumes.find(elem_id) != _elem_granular_volumes.end())
132 504 : gran_vol = _elem_granular_volumes.find(elem_id)->second;
133 : else
134 0 : gran_vol = {};
135 504 : }
136 :
137 : Real
138 432 : MDRunBase::particleProperty(unsigned int j, unsigned int prop_id) const
139 : {
140 : // ensure that entry exists
141 : mooseAssert(j < _md_particles.properties.size(),
142 : "Particle index " << j << " not found in _md_particles. properties vector has length "
143 : << _md_particles.properties.size());
144 432 : return _md_particles.properties[j][propIndex(prop_id)];
145 : }
146 :
147 : void
148 260 : MDRunBase::mapMDParticles()
149 : {
150 : // clear data structures
151 : _elem_particles.clear();
152 260 : _md_particles.elem_id.assign(_md_particles.pos.size(), libMesh::DofObject::invalid_unique_id);
153 :
154 : // loop over semi-local elements to ensure consistent handling of
155 : // points on processor boundaries
156 260 : const libMesh::MeshBase & mesh_base = _mesh.getMesh();
157 520 : for (const auto & elem : as_range(mesh_base.active_semilocal_elements_begin(),
158 5408 : mesh_base.active_semilocal_elements_end()))
159 : {
160 : // find all points within an inflated bounding box
161 : std::vector<nanoflann::ResultItem<std::size_t, Real>> indices_dist;
162 2184 : BoundingBox bbox = elem->loose_bounding_box();
163 : Point center = 0.5 * (bbox.min() + bbox.max());
164 2184 : Real radius = (bbox.max() - center).norm();
165 2184 : _kd_tree->radiusSearch(center, radius, indices_dist);
166 :
167 2932 : for (auto & p : indices_dist)
168 : {
169 748 : Point candidate = _md_particles.pos[p.first];
170 :
171 : // avoid double counting of elements on element boundaries, smallest
172 : // elem id gets to own the point, this is consistent across processors
173 : // since we loop over semi-local elements
174 748 : if (_md_particles.elem_id[p.first] != libMesh::DofObject::invalid_unique_id &&
175 168 : elem->unique_id() > _md_particles.elem_id[p.first])
176 168 : continue;
177 :
178 : // contains_point performs cheap bounding box test, hence no need to do it before,
179 : // serious candidates need to go through the expensive test
180 580 : if (elem->contains_point(candidate))
181 : {
182 : // convenience variable for the MD particle we deal with
183 424 : unsigned int pp = p.first;
184 :
185 : // remove this particle from the entry for the old element to avoid double
186 : // counting
187 424 : if (_md_particles.elem_id[pp] != libMesh::DofObject::invalid_unique_id)
188 : {
189 : // get the stored unique id of the _md_particle
190 0 : unique_id_type old_elem_id = _md_particles.elem_id[pp];
191 :
192 : // entry in _elem_particles should exist but guard w/ assert
193 : mooseAssert(_elem_particles.find(old_elem_id) != _elem_particles.end(),
194 : "Entry " << old_elem_id << " in _elem_particles should exist.");
195 0 : std::vector<unsigned int> id_list = _elem_particles.find(old_elem_id)->second;
196 :
197 : // go through this list and find the entry with the right MD particle entry
198 : // and remove it
199 0 : for (unsigned int i = 0; i < id_list.size(); ++i)
200 0 : if (id_list[i] == pp)
201 : id_list.erase(id_list.begin() + i);
202 : }
203 :
204 : // insert entry for new element
205 424 : if (_elem_particles.find(elem->unique_id()) != _elem_particles.end())
206 112 : _elem_particles[elem->unique_id()].push_back(pp);
207 : else
208 312 : _elem_particles[elem->unique_id()] = {pp};
209 :
210 : // re-assigning the element id
211 424 : _md_particles.elem_id[pp] = elem->unique_id();
212 : }
213 : }
214 260 : }
215 260 : }
216 :
217 : void
218 16 : MDRunBase::updateElementGranularVolumes()
219 : {
220 : // clear the granular volume map
221 : _elem_granular_volumes.clear();
222 :
223 : /*
224 : // Ideally we want to compute _max_granular_radius automatically
225 : // but the value is needed for setting bounding boxes to partition
226 : // the particles on the processors. For now, it's an input parameter
227 : // until a path forward is determined.
228 : // The commented lines compute the largest granular radius
229 : _max_granular_radius = 0;
230 : for (auto & p : _md_particles.properties)
231 : if (p[7] > _max_granular_radius)
232 : _max_granular_radius = p[7];
233 : */
234 : /// do a sanity check _max_granular_radius
235 16 : Real mgr = 0;
236 80 : for (auto & p : _md_particles.properties)
237 64 : if (p[_md_particles._r_index] > mgr)
238 16 : mgr = p[_md_particles._r_index];
239 16 : if (mgr > _max_granular_radius)
240 0 : mooseError("Granular particle with radius: ",
241 : mgr,
242 : " exceeds max_granular_radius: ",
243 0 : _max_granular_radius);
244 :
245 : /// loop over all local elements
246 16 : ConstElemRange * active_local_elems = _mesh.getActiveLocalElementRange();
247 160 : for (const auto & elem : *active_local_elems)
248 : {
249 : // find all points within an inflated bounding box
250 : std::vector<nanoflann::ResultItem<std::size_t, Real>> indices_dist;
251 144 : BoundingBox bbox = elem->loose_bounding_box();
252 : Point center = 0.5 * (bbox.min() + bbox.max());
253 :
254 : // inflate the search sphere by the maximum granular radius
255 144 : Real radius = (bbox.max() - center).norm() + _max_granular_radius;
256 144 : _kd_tree->radiusSearch(center, radius, indices_dist);
257 :
258 : // prepare _elem_granular_candidates entry
259 144 : _elem_granular_volumes[elem->unique_id()] = {};
260 :
261 : // construct this element's overlap object
262 144 : ElemType t = elem->type();
263 144 : OVERLAP::Hexahedron hex = overlapUnitHex();
264 144 : OVERLAP::Tetrahedron tet = overlapUnitTet();
265 144 : if (t == HEX8)
266 72 : hex = overlapHex(elem);
267 72 : else if (t == TET4)
268 72 : tet = overlapTet(elem);
269 : else
270 0 : mooseError("Element type ", t, "not implemented");
271 :
272 : // loop through all MD particles that the search turned up and test overlap
273 360 : for (unsigned int j = 0; j < indices_dist.size(); ++j)
274 : {
275 : // construct OVERLAP::sphere object from MD granular particle
276 216 : unsigned int k = indices_dist[j].first;
277 : OVERLAP::Sphere sph(OVERLAP::vector_t{_md_particles.pos[k](0),
278 : _md_particles.pos[k](1),
279 : _md_particles.pos[k](2)},
280 216 : _md_particles.properties[k][_md_particles._r_index]);
281 :
282 : // compute the overlap
283 : Real ovlp = 0.0;
284 216 : if (t == HEX8)
285 144 : ovlp = OVERLAP::overlap(sph, hex);
286 72 : else if (t == TET4)
287 72 : ovlp = OVERLAP::overlap(sph, tet);
288 :
289 : // if the overlap is larger than 0, make entry in _elem_granular_volumes
290 216 : if (ovlp > 0.0)
291 432 : _elem_granular_volumes[elem->unique_id()].push_back(std::pair<unsigned int, Real>(k, ovlp));
292 : }
293 : }
294 16 : }
295 :
296 : MultiMooseEnum
297 78 : MDRunBase::mdParticleProperties()
298 : {
299 156 : return MultiMooseEnum("vel_x=0 vel_y=1 vel_z=2 force_x=3 force_y=4 force_z=5 charge=6 radius=7");
300 : }
301 :
302 : unsigned int
303 648 : MDRunBase::propIndex(unsigned int prop_id) const
304 : {
305 : auto it = _md_particles._map_props.find(prop_id);
306 648 : if (it == _md_particles._map_props.end())
307 0 : mooseError("Property id ", prop_id, " is not present in _map_props map.");
308 648 : return it->second;
309 : }
310 :
311 : unsigned int
312 12 : MDRunBase::propIndex(const std::string & prop_name) const
313 : {
314 12 : unsigned int prop_id = mdParticleProperties().find(prop_name)->id();
315 12 : return propIndex(prop_id);
316 : }
317 :
318 : OVERLAP::Hexahedron
319 72 : MDRunBase::overlapHex(const Elem * elem) const
320 : {
321 : Point p;
322 72 : p = elem->point(0);
323 : OVERLAP::vector_t v0{p(0), p(1), p(2)};
324 72 : p = elem->point(1);
325 : OVERLAP::vector_t v1{p(0), p(1), p(2)};
326 72 : p = elem->point(2);
327 : OVERLAP::vector_t v2{p(0), p(1), p(2)};
328 72 : p = elem->point(3);
329 : OVERLAP::vector_t v3{p(0), p(1), p(2)};
330 72 : p = elem->point(4);
331 : OVERLAP::vector_t v4{p(0), p(1), p(2)};
332 72 : p = elem->point(5);
333 : OVERLAP::vector_t v5{p(0), p(1), p(2)};
334 72 : p = elem->point(6);
335 : OVERLAP::vector_t v6{p(0), p(1), p(2)};
336 72 : p = elem->point(7);
337 : OVERLAP::vector_t v7{p(0), p(1), p(2)};
338 72 : return OVERLAP::Hexahedron{v0, v1, v2, v3, v4, v5, v6, v7};
339 : }
340 :
341 : OVERLAP::Hexahedron
342 144 : MDRunBase::overlapUnitHex() const
343 : {
344 : OVERLAP::vector_t v0{-1, -1, -1};
345 : OVERLAP::vector_t v1{1, -1, -1};
346 : OVERLAP::vector_t v2{1, 1, -1};
347 : OVERLAP::vector_t v3{-1, 1, -1};
348 : OVERLAP::vector_t v4{-1, -1, 1};
349 : OVERLAP::vector_t v5{1, -1, 1};
350 : OVERLAP::vector_t v6{1, 1, 1};
351 : OVERLAP::vector_t v7{-1, 1, 1};
352 144 : return OVERLAP::Hexahedron{v0, v1, v2, v3, v4, v5, v6, v7};
353 : }
354 :
355 : OVERLAP::Tetrahedron
356 72 : MDRunBase::overlapTet(const Elem * elem) const
357 : {
358 : Point p;
359 72 : p = elem->point(0);
360 : OVERLAP::vector_t v0{p(0), p(1), p(2)};
361 72 : p = elem->point(1);
362 : OVERLAP::vector_t v1{p(0), p(1), p(2)};
363 72 : p = elem->point(2);
364 : OVERLAP::vector_t v2{p(0), p(1), p(2)};
365 72 : p = elem->point(3);
366 : OVERLAP::vector_t v3{p(0), p(1), p(2)};
367 72 : return OVERLAP::Tetrahedron{v0, v1, v2, v3};
368 : }
369 :
370 : OVERLAP::Tetrahedron
371 144 : MDRunBase::overlapUnitTet() const
372 : {
373 : OVERLAP::vector_t v0{0, 0, 0};
374 : OVERLAP::vector_t v1{1, 0, 0};
375 : OVERLAP::vector_t v2{0, 1, 0};
376 : OVERLAP::vector_t v3{0, 0, 1};
377 144 : return OVERLAP::Tetrahedron{v0, v1, v2, v3};
378 : }
|
