Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #pragma once
11 :
12 : #include "GeneralUserObject.h"
13 :
14 : // Local Includes
15 : #include "ElemIndexHelper.h"
16 : #include "ParallelRayStudy.h"
17 : #include "RayTracingAttributes.h"
18 : #include "TraceData.h"
19 : #include "TraceRayBndElement.h"
20 :
21 : // MOOSE Includes
22 : #include "TheWarehouse.h"
23 :
24 : // libMesh includes
25 : #include "libmesh/mesh.h"
26 : #include "libmesh/elem_side_builder.h"
27 :
28 : // Forward Declarations
29 : class RayBoundaryConditionBase;
30 : class RayKernelBase;
31 : class TraceRay;
32 : class RayTracingObject;
33 :
34 : /**
35 : * Base class for Ray tracing studies that will generate Rays and then propagate
36 : * all of them to termination.
37 : *
38 : * Subclasses _must_ override generateRays()
39 : */
40 : class RayTracingStudy : public GeneralUserObject
41 : {
42 : public:
43 : RayTracingStudy(const InputParameters & parameters);
44 :
45 : /**
46 : * Key that is used for restricting access to acquireRayInternal().
47 : *
48 : * In general, we restrict access to the construction of Rays to only protected
49 : * methods within RayTracingStudy. However, there are a few classes/functions
50 : * that have a legimate reason to construct Rays. Here, the friends are said
51 : * classes/functions.
52 : */
53 : class AcquireRayInternalKey
54 : {
55 : friend class Parallel::Packing<std::shared_ptr<Ray>>;
56 : friend void dataLoad(std::istream & stream, std::shared_ptr<Ray> & ray, void * context);
57 : AcquireRayInternalKey() {}
58 : AcquireRayInternalKey(const AcquireRayInternalKey &) {}
59 : };
60 :
61 : /**
62 : * Key that is used for restricting access to moveRayToBufferDuringTrace() and
63 : * acquireRayDuringTrace().
64 : */
65 : class AcquireMoveDuringTraceKey
66 : {
67 : friend class RayBoundaryConditionBase;
68 : friend class RayKernelBase;
69 : AcquireMoveDuringTraceKey() {}
70 : AcquireMoveDuringTraceKey(const AcquireMoveDuringTraceKey &) {}
71 : };
72 :
73 : static InputParameters validParams();
74 :
75 : virtual void initialSetup() override;
76 : virtual void residualSetup() override;
77 : virtual void jacobianSetup() override;
78 : virtual void meshChanged() override;
79 : virtual void timestepSetup() override;
80 5463 : virtual void initialize() override {}
81 5363 : virtual void finalize() override {}
82 :
83 : /**
84 : * Executes the study (generates and propagates Rays)
85 : */
86 : virtual void execute() override;
87 :
88 : /**
89 : * Setup for on subdomain change or subdomain AND ray change during ray tracing
90 : * @param subdomain The subdomain changed to
91 : * @param tid Thread id
92 : * @param ray_id ID of the ray initiating the change
93 : */
94 : virtual void
95 : segmentSubdomainSetup(const SubdomainID subdomain, const THREAD_ID tid, const RayID ray_id);
96 :
97 : /**
98 : * Reinitialize objects for a Ray segment for ray tracing
99 : * @param elem The elem the segment is in
100 : * @param start Start point of the segment
101 : * @param end End point of the segment
102 : * @param length The length of the start -> end segment
103 : * @param tid Thread id
104 : */
105 : virtual void reinitSegment(
106 : const Elem * elem, const Point & start, const Point & end, const Real length, THREAD_ID tid);
107 :
108 : /**
109 : * Called at the end of a Ray segment
110 : * @param tid Thread id
111 : * @param ray The ray
112 : */
113 : virtual void postOnSegment(const THREAD_ID tid, const std::shared_ptr<Ray> & ray);
114 :
115 : /**
116 : * Called at the beginning of a trace for a ray
117 : * @param tid Thread id
118 : * @param ray The ray
119 : */
120 3267644 : virtual void preTrace(const THREAD_ID /* tid */, const std::shared_ptr<Ray> & /* ray */) {}
121 :
122 : /**
123 : * Method for executing the study so that it can be called out of the standard UO execute()
124 : */
125 : void executeStudy();
126 :
127 : /**
128 : * Total number of processor crossings for Rays that finished on this processor
129 : */
130 : unsigned long long int endingProcessorCrossings() const { return _ending_processor_crossings; }
131 : /**
132 : * Max number of total processor crossings for Rays that finished on this processor
133 : */
134 : unsigned int endingMaxProcessorCrossings() const { return _ending_max_processor_crossings; }
135 : /**
136 : * Total number of processor crossings
137 : */
138 8 : unsigned long long int totalProcessorCrossings() const { return _total_processor_crossings; }
139 : /**
140 : * Max number of processor crossings for all Rays
141 : */
142 8 : unsigned int maxProcessorCrossings() const { return _max_processor_crossings; }
143 :
144 : /**
145 : * Total number of Ray/element intersections for Rays that finished on this processor
146 : */
147 : unsigned long long int endingIntersections() const { return _ending_intersections; }
148 : /**
149 : * Max number of intersections for Rays that finished on this processor
150 : */
151 : unsigned int endingMaxIntersections() const { return _ending_max_intersections; }
152 : /**
153 : * Total number of Ray/element intersections
154 : */
155 : unsigned long long int totalIntersections() const { return _total_intersections; }
156 : /**
157 : * Max number of intersections for a Ray
158 : */
159 : unsigned int maxIntersections() const { return _max_intersections; }
160 : /**
161 : * Max number of trajectory changes for a Ray
162 : */
163 : unsigned int maxTrajectoryChanges() const { return _max_trajectory_changes; }
164 :
165 : /**
166 : * Total amount of distance traveled by the rays that end on this processor
167 : */
168 : Real endingDistance() const { return _ending_distance; }
169 : /**
170 : * Total distance traveled by all Rays
171 : */
172 1324 : Real totalDistance() const { return _total_distance; }
173 :
174 : unsigned long long int localTraceRayResult(const int result) const
175 : {
176 112 : return _local_trace_ray_results[result];
177 : }
178 :
179 : const ParallelRayStudy & parallelRayStudy() const { return *_parallel_ray_study; }
180 :
181 : /**
182 : * Max distance any Ray can travel
183 : */
184 23617495 : Real rayMaxDistance() const { return _ray_max_distance; }
185 :
186 : /**
187 : * Duration for execute() in seconds
188 : */
189 : Real executionTime() { return std::chrono::duration<Real>(_execution_time).count(); }
190 : /**
191 : * Duration for execute() in nanoseconds
192 : */
193 : Real executionTimeNano()
194 : {
195 : return std::chrono::duration<Real, std::nano>(_execution_time).count();
196 : }
197 : /**
198 : * Duration for creation of all Rays in seconds
199 : */
200 : Real generationTime() const { return std::chrono::duration<Real>(_generation_time).count(); }
201 : /**
202 : * Duration for creation of all Rays in seconds
203 : */
204 : Real propagationTime() const { return std::chrono::duration<Real>(_propagation_time).count(); }
205 :
206 : /**
207 : * Whether or not to tolerate failure
208 : */
209 8 : bool tolerateFailure() const { return _tolerate_failure; }
210 :
211 : /**
212 : * Whether or not to bank Rays on completion
213 : */
214 452 : bool bankRaysOnCompletion() const { return _bank_rays_on_completion; }
215 :
216 : /**
217 : * Whether or not to use Ray dependent subdomain setup
218 : */
219 3261117 : bool rayDependentSubdomainSetup() const { return _ray_dependent_subdomain_setup; }
220 :
221 : /**
222 : * Register a value to be filled in the data on a Ray with a given name.
223 : * @param name The name to register
224 : * \return The Ray data index for the registered value
225 : *
226 : * Note that this does not actually allocate the storage. It is simply a registration system to
227 : * keep track of the values stored in the Ray data.
228 : */
229 : RayDataIndex registerRayData(const std::string & name);
230 : /**
231 : * Register values to be filled in the data on a Ray with a given name.
232 : * @param names The names to register
233 : * \return The Ray data indices for the registered value
234 : *
235 : * Note that this does not actually allocate the storage. It is simply a registration system to
236 : * keep track of the values stored in the Ray data.
237 : */
238 : std::vector<RayDataIndex> registerRayData(const std::vector<std::string> & names);
239 : /**
240 : * Gets the index associated with a registered value in the Ray data
241 : * @param name The value name to get the index of
242 : * @param graceful Whether or not to exit gracefully if none is found (return
243 : * INVALID_RAY_DATA_INDEX)
244 : * \return The index for the value
245 : */
246 : RayDataIndex getRayDataIndex(const std::string & name, const bool graceful = false) const;
247 : /**
248 : * Gets the indices associated with registered values in the Ray data
249 : * @param names The value names to get the indices of
250 : * @param graceful Whether or not to exit gracefully if none is found (index is filled with
251 : * INVALID_RAY_DATA_INDEX)
252 : * \return The indices for the values
253 : */
254 : std::vector<RayDataIndex> getRayDataIndices(const std::vector<std::string> & names,
255 : const bool graceful = false) const;
256 : /**
257 : * Gets the name associated with a registered value in the Ray data
258 : * @param index The index to get the name of
259 : * \return The name associated with index
260 : */
261 : const std::string & getRayDataName(const RayDataIndex index) const;
262 : /**
263 : * Gets the names associated with registered values in the Ray data
264 : * @param indices The indices to get the names of
265 : * \return The associated names
266 : */
267 : std::vector<std::string> getRayDataNames(const std::vector<RayDataIndex> & indices) const;
268 :
269 : /**
270 : * The registered size of values in the Ray data.
271 : */
272 : std::size_t rayDataSize() const { return _ray_data_names.size(); }
273 : /**
274 : * Whether or not any Ray data are registered
275 : */
276 : bool hasRayData() const { return _ray_data_names.size(); }
277 : /**
278 : * The Ray data names
279 : */
280 118 : const std::vector<std::string> & rayDataNames() const { return _ray_data_names; }
281 :
282 : /**
283 : * Register a value to be filled in the aux data on a Ray with a given name.
284 : * @param name The name to register
285 : * \returns The Ray aux data index for the registered value
286 : *
287 : * Note that this does not actually allocate the storage. It is simply a registration system to
288 : * keep track of the values stored in the Ray aux data.
289 : */
290 : RayDataIndex registerRayAuxData(const std::string & name);
291 : /**
292 : * Register values to be filled in the aux data on a Ray with a given name.
293 : * @param names The names to register
294 : * \returns The Ray aux data indices for the registered values
295 : *
296 : * Note that this does not actually allocate the storage. It is simply a registration system to
297 : * keep track of the values stored in the Ray aux data.
298 : */
299 : std::vector<RayDataIndex> registerRayAuxData(const std::vector<std::string> & names);
300 : /**
301 : * Gets the index associated with a registered value in the Ray aux data
302 : * @param name The value name to get the index of
303 : * @param graceful Whether or not to exit gracefully if none is found (return
304 : * INVALID_RAY_DATA_INDEX)
305 : * \returns The index for the value
306 : */
307 : RayDataIndex getRayAuxDataIndex(const std::string & name, const bool graceful = false) const;
308 : /**
309 : * Gets the indices associated with registered values in the Ray aux data
310 : * @param names The value names to get the indices of
311 : * @param graceful Whether or not to exit gracefully if none is found (index is filled with
312 : * INVALID_RAY_DATA_INDEX)
313 : * \return The indices for the values
314 : */
315 : std::vector<RayDataIndex> getRayAuxDataIndices(const std::vector<std::string> & names,
316 : const bool graceful = false) const;
317 : /**
318 : * Gets the name associated with a registered value in the Ray aux data
319 : * @param index The index to get the name of
320 : * \return The name associated with index
321 : */
322 : const std::string & getRayAuxDataName(const RayDataIndex index) const;
323 : /**
324 : * Gets the names associated with registered values in the Ray aux data
325 : * @param indices The indices to get the names of
326 : * \return The associated names
327 : */
328 : std::vector<std::string> getRayAuxDataNames(const std::vector<RayDataIndex> & indices) const;
329 : /**
330 : * The registered size of values in the Ray aux data.
331 : */
332 : std::size_t rayAuxDataSize() const { return _ray_aux_data_names.size(); }
333 : /**
334 : * Whether or not any Ray aux data are registered
335 : */
336 : bool hasRayAuxData() const { return _ray_aux_data_names.size(); }
337 : /**
338 : * The Ray aux data names
339 : */
340 50 : const std::vector<std::string> & rayAuxDataNames() const { return _ray_aux_data_names; }
341 :
342 : /**
343 : * Whether or not there are currently any active RayKernel objects
344 : */
345 : bool hasRayKernels(const THREAD_ID tid);
346 : /**
347 : * Fills the active RayKernels associated with this study and a block into result
348 : */
349 : void getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid);
350 : /**
351 : * Fills the active RayKernels associated with this study into result
352 : */
353 : template <typename T>
354 16769 : void getRayKernels(std::vector<T *> & result, THREAD_ID tid)
355 : {
356 16769 : _fe_problem.theWarehouse()
357 : .query()
358 33538 : .condition<AttribRayTracingStudy>(this)
359 16769 : .condition<AttribSystem>("RayKernel")
360 16769 : .condition<AttribThread>(tid)
361 : .queryInto(result);
362 16769 : }
363 : /**
364 : * Fills the active RayKernels associeted with this study, block, and potentially Ray into result
365 : */
366 : void
367 : getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid, RayID ray_id);
368 : /**
369 : * Fills the active RayBCs associated with this study and a boundary into result
370 : */
371 : void getRayBCs(std::vector<RayBoundaryConditionBase *> & result, BoundaryID id, THREAD_ID tid);
372 : /**
373 : * Fills the active RayBCs associated with this study and boundaries result
374 : */
375 : template <typename T>
376 63645 : void getRayBCs(std::vector<T *> & result, const std::vector<BoundaryID> & ids, THREAD_ID tid)
377 : {
378 63645 : _fe_problem.theWarehouse()
379 : .query()
380 127290 : .condition<AttribRayTracingStudy>(this)
381 63645 : .condition<AttribSystem>("RayBoundaryCondition")
382 63645 : .condition<AttribBoundaries>(ids)
383 63645 : .condition<AttribThread>(tid)
384 : .queryInto(result);
385 63645 : }
386 : /**
387 : * Fills the active RayBCs associated with this study into result
388 : */
389 : template <typename T>
390 3810 : void getRayBCs(std::vector<T *> & result, THREAD_ID tid)
391 : {
392 3810 : _fe_problem.theWarehouse()
393 : .query()
394 7620 : .condition<AttribRayTracingStudy>(this)
395 3810 : .condition<AttribSystem>("RayBoundaryCondition")
396 3810 : .condition<AttribThread>(tid)
397 : .queryInto(result);
398 3810 : }
399 : /**
400 : * Fills the active RayBCs associated with thie study, boundary elements, and potentially Ray into
401 : * result.
402 : *
403 : * This is purposely virtual because it allows derived studies to optimize the retrieval of RayBCs
404 : * during the trace in TraceRay.
405 : */
406 : virtual void getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
407 : const std::vector<TraceRayBndElement> & bnd_elems,
408 : THREAD_ID tid,
409 : RayID ray_id);
410 :
411 : /**
412 : * Gets the current RayKernels for a thread, which are set in segmentSubdomainSetup()
413 : */
414 : const std::vector<RayKernelBase *> & currentRayKernels(THREAD_ID tid) const
415 : {
416 74414399 : return _threaded_current_ray_kernels[tid];
417 : }
418 :
419 : /**
420 : * Get the nodal bounding box for the domain
421 : */
422 2792189 : const BoundingBox & boundingBox() const { return _b_box; }
423 : /**
424 : * Get the loose nodal bounding box for the domain
425 : *
426 : * BoundingBox::contains_point does not have a tolerance, so we need to stretch the
427 : * box a little bit for some contains_point checks.
428 : */
429 9035212 : const BoundingBox & looseBoundingBox() const { return _loose_b_box; }
430 : /**
431 : * Get the inflated maximum length across the domain
432 : */
433 119019440 : Real domainMaxLength() const { return _domain_max_length; }
434 : /**
435 : * Get the current total volume of the domain
436 : */
437 7198 : Real totalVolume() const { return _total_volume; }
438 : /**
439 : * Whether or not the domain is rectangular (if it is prefectly encompassed by its bounding box)
440 : */
441 : bool isRectangularDomain() const;
442 :
443 : /**
444 : * Whether or not the local mesh has internal sidesets that have RayBCs on them
445 : *
446 : * NOTE: if useInternalSidesets() == false, this will be false even if the mesh does
447 : * have internal sidesets
448 : */
449 : bool hasInternalSidesets() const { return _internal_sidesets.size(); }
450 : /**
451 : * Get the internal sidesets (that have RayBC(s)) for each side for a given element
452 : *
453 : * This will be empty if the elem does not have any internal sidesets that have RayBC(s)
454 : */
455 : const std::vector<std::vector<BoundaryID>> & getInternalSidesets(const Elem * elem) const;
456 : /**
457 : * Gets the internal sidesets (that have RayBCs) within the local domain
458 : */
459 : const std::set<BoundaryID> & getInternalSidesets() const { return _internal_sidesets; }
460 : /**
461 : * Whether or not the side \s on elem \p elem is non-planar
462 : *
463 : * This is cached because checking whether or not a face is planar is costly
464 : */
465 9407924 : bool sideIsNonPlanar(const Elem * elem, const unsigned short s) const
466 : {
467 9407924 : return _has_non_planar_sides && _non_planar_sides[_elem_index_helper.getIndex(elem)][s] == 1;
468 : }
469 : /**
470 : * Whether or not the mesh has active elements of the same level
471 : *
472 : * Use this over sameLevelActiveElems(), which is for internally setting
473 : * _has_same_level_active_elems
474 : */
475 42628748 : bool hasSameLevelActiveElems() const { return _has_same_level_active_elems; }
476 :
477 : /**
478 : * INTERNAL methods for acquiring a Ray during a trace in RayKernels and RayBCs.
479 : *
480 : * You should not use these APIs directly. If you wish to acquire a Ray during generation
481 : * during generateRays()), use the protected RayTracingStudy::acquire{}Ray() methods.
482 : * If you wish to acquire a Ray during propagation in RayKernels and RayBC, use the
483 : * protected RayKernelBase::acquireRay() and RayBoundaryConditionBase::acquireRay(),
484 : * respectively.
485 : */
486 : ///@{
487 : std::shared_ptr<Ray> acquireRayDuringTrace(const THREAD_ID tid,
488 : const AcquireMoveDuringTraceKey &);
489 1407473 : std::shared_ptr<Ray> acquireRayInternal(const RayID id,
490 : const std::size_t data_size,
491 : const std::size_t aux_data_size,
492 : const bool reset,
493 : const AcquireRayInternalKey &)
494 : {
495 1407473 : return _parallel_ray_study->acquireParallelData(
496 2814946 : 0, this, id, data_size, aux_data_size, reset, Ray::ConstructRayKey());
497 : }
498 : ///@}
499 :
500 : /**
501 : * INTERNAL method for moving a Ray into the buffer during tracing.
502 : *
503 : * You should not and cannot use this API. It is protected by the
504 : * AcquireMoveDuringTraceKey, which is only constructable by
505 : * RayKernelBase and RayBoundaryConditionBase.
506 : */
507 : void moveRayToBufferDuringTrace(std::shared_ptr<Ray> & ray,
508 : const THREAD_ID tid,
509 : const AcquireMoveDuringTraceKey &);
510 : /**
511 : * Access to the libMesh MeshBase.
512 : *
513 : * This is needed for unpack routines for a Ray, which has a context of this study.
514 : */
515 1928236 : MeshBase & meshBase() const { return _mesh; }
516 :
517 : /**
518 : * @returns A reference to the MooseMesh associated with this study.
519 : */
520 1864 : MooseMesh & mesh() { return _mesh; }
521 :
522 : /**
523 : * Get the outward normal for a given element side.
524 : */
525 : virtual const Point &
526 : getSideNormal(const Elem * elem, const unsigned short side, const THREAD_ID tid);
527 : /**
528 : * Gets the outward normals for a given element.
529 : * Returns a pointer to the normal for the zeroth side.
530 : */
531 2 : virtual const Point * getElemNormals(const Elem * /* elem */, const THREAD_ID /* tid */)
532 : {
533 2 : mooseError("Unimplemented element normal caching in ", type(), "::getElemNormals()");
534 : }
535 :
536 : /**
537 : * Gets the data value for a banked ray with a given ID
538 : *
539 : * This will return the value replicated across all processors
540 : */
541 : RayData getBankedRayData(const RayID ray_id, const RayDataIndex index) const;
542 : /**
543 : * Gets the data value for a banked ray with a given ID
544 : *
545 : * This will return the value replicated across all processors
546 : */
547 : RayData getBankedRayAuxData(const RayID ray_id, const RayDataIndex index) const;
548 :
549 : /**
550 : * Gets the ID of a registered ray
551 : * @param name The name of said ray
552 : * @param graceful Whether or not to exit gracefully if none is found (with invalid_id)
553 : */
554 : RayID registeredRayID(const std::string & name, const bool graceful = false) const;
555 : /**
556 : * Gets the name of a registered ray
557 : * @param ray_id The ID of said ray
558 : */
559 : const std::string & registeredRayName(const RayID ray_id) const;
560 :
561 : /**
562 : * Whether or not ray registration is being used
563 : */
564 2958 : bool useRayRegistration() const { return _use_ray_registration; }
565 :
566 : /**
567 : * Whether or not to store the Ray data on the cached Ray traces
568 : */
569 47064 : bool dataOnCacheTraces() const { return _data_on_cache_traces; }
570 : /**
571 : * Whether or not to store the Ray aux data on the cached Ray traces
572 : */
573 47016 : bool auxDataOnCacheTraces() const { return _aux_data_on_cache_traces; }
574 : /**
575 : * Whether or not to cache individual element segments when _cache_traces = true
576 : */
577 40964 : bool segmentsOnCacheTraces() const { return _segments_on_cache_traces; }
578 :
579 : /**
580 : * Virtual that allows for selection in if a Ray should be cached or not (only used when
581 : * _cache_traces).
582 : */
583 3267644 : virtual bool shouldCacheTrace(const std::shared_ptr<Ray> & /* ray */) const
584 : {
585 3267644 : return _always_cache_traces;
586 : }
587 :
588 : /**
589 : * Initialize a Ray in the threaded cached trace map to be filled with segments
590 : */
591 : TraceData & initThreadedCachedTrace(const std::shared_ptr<Ray> & ray, THREAD_ID tid);
592 : /**
593 : * Get the cached trace data structure
594 : */
595 : const std::vector<TraceData> & getCachedTraces() const { return _cached_traces; }
596 :
597 : /**
598 : * Get the cached hmax for all elements in a subdomain
599 : *
600 : * Used for scaling tolerances in ray tracing.
601 : */
602 : Real subdomainHmax(const SubdomainID subdomain_id) const;
603 :
604 : /**
605 : * Entry point for acting on a ray when it is completed (shouldContinue() == false)
606 : */
607 : virtual void onCompleteRay(const std::shared_ptr<Ray> & ray);
608 :
609 : /**
610 : * Verifies that the Rays in the given range have unique Ray IDs.
611 : *
612 : * @param begin The beginning of the range of Rays to check
613 : * @param end The end of the range of Rays to check
614 : * @param global If true, this will complete the verification across all processors
615 : * @param error_suffix Entry point for additional information in the error message
616 : */
617 : void verifyUniqueRayIDs(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
618 : const std::vector<std::shared_ptr<Ray>>::const_iterator end,
619 : const bool global,
620 : const std::string & error_suffix) const;
621 :
622 : /**
623 : * Verifies that the Rays in the given range are unique. That is, that there are not multiple
624 : * shared_ptrs that point to the same Ray
625 : *
626 : * @param begin The beginning of the range of Rays to check
627 : * @param end The end of the range of Rays to check
628 : * @param error_suffix Entry point for additional information in the error message
629 : */
630 : void verifyUniqueRays(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
631 : const std::vector<std::shared_ptr<Ray>>::const_iterator end,
632 : const std::string & error_suffix);
633 :
634 : /**
635 : * Whether or not the study is propagating (tracing Rays)
636 : */
637 : bool currentlyPropagating() const { return _parallel_ray_study->currentlyExecuting(); }
638 : /**
639 : * Whether or not the study is generating
640 : */
641 : bool currentlyGenerating() const { return _parallel_ray_study->currentlyPreExecuting(); }
642 :
643 : /**
644 : * Gets the threaded TraceRay object for \p tid.
645 : */
646 : ///@{
647 20 : TraceRay & traceRay(const THREAD_ID tid) { return *_threaded_trace_ray[tid]; }
648 4119 : const TraceRay & traceRay(const THREAD_ID tid) const { return *_threaded_trace_ray[tid]; }
649 : ///@}
650 :
651 : /**
652 : * Whether or not to verify if Rays have valid information before being traced
653 : */
654 12314862 : bool verifyRays() const { return _verify_rays; }
655 : /**
656 : * Whether or not trace verification is enabled in devel/dbg modes
657 : */
658 : #ifndef NDEBUG
659 : bool verifyTraceIntersections() const { return _verify_trace_intersections; }
660 : #endif
661 :
662 : /**
663 : * Whether or not \p side is incoming on element \p elem in direction \p direction.
664 : */
665 : bool sideIsIncoming(const Elem * const elem,
666 : const unsigned short side,
667 : const Point & direction,
668 : const THREAD_ID tid);
669 :
670 : /**
671 : * Whether or not to produce a warning when interacting with a non-planar mesh.
672 : */
673 : bool warnNonPlanar() const { return _warn_non_planar; }
674 :
675 : /**
676 : * The underlying parallel study: used for the context for calling the packed range routines.
677 : */
678 : ParallelStudy<std::shared_ptr<Ray>, Ray> * parallelStudy() { return _parallel_ray_study.get(); }
679 :
680 : /**
681 : * Get an element's side pointer without excessive memory allocation
682 : *
683 : * @param elem The element to build a side for
684 : * @param s The side to build
685 : * @param tid The thread id
686 : * @return A pointer to the side element
687 : */
688 : const libMesh::Elem &
689 : elemSide(const libMesh::Elem & elem, const unsigned int s, const THREAD_ID tid = 0)
690 : {
691 1570032 : return _threaded_elem_side_builders[tid](elem, s);
692 : }
693 :
694 : protected:
695 : /**
696 : * Subclasses should override this to determine how to generate Rays.
697 : * This will be called within execute() and makes up the "generation phase"
698 : * of the algorithm.
699 : */
700 : virtual void generateRays() = 0;
701 :
702 : /**
703 : * Entry point before study execution
704 : */
705 5360 : virtual void preExecuteStudy() {}
706 : /**
707 : * Entry point after study execution
708 : */
709 5030 : virtual void postExecuteStudy() {}
710 :
711 : /**
712 : * Helper function for computing the total domain volume
713 : */
714 : Real computeTotalVolume();
715 :
716 : /**
717 : * Gets the writeable current RayKernels for a thread
718 : *
719 : * Allows for other ray studies to fill the current ray kernels in a custom manner
720 : */
721 : std::vector<RayKernelBase *> & currentRayKernelsWrite(THREAD_ID tid)
722 : {
723 : return _threaded_current_ray_kernels[tid];
724 : }
725 :
726 : /**
727 : * Reserve \p size entires in the Ray buffer.
728 : *
729 : * This can only be used within generateRays() and should be used when possible
730 : * to reserve entires in the buffer before adding Rays via moveRay(s)ToBuffer.
731 : */
732 : void reserveRayBuffer(const std::size_t size);
733 :
734 : /**
735 : * Determine whether or not the mesh currently has active elements that are
736 : * all the same level
737 : */
738 : bool sameLevelActiveElems() const;
739 :
740 : /**
741 : * Builds quadrature points for a given segment using the _segment_qrule
742 : * @param start Start point of the segment
743 : * @param end End point of the segment
744 : * @param length The lengh of the start -> end segment
745 : * @param points Points to fill into (should be sized ahead of time)
746 : * @param weights Weights to fill into (should be sized ahead of time)
747 : */
748 : virtual void buildSegmentQuadrature(const Point & start,
749 : const Point & end,
750 : const Real length,
751 : std::vector<Point> & points,
752 : std::vector<Real> & weights) const;
753 :
754 : /**
755 : * Get the Ray bank. This is the bank of Rays that have completed on this processor
756 : * after an execution of the study.
757 : *
758 : * This is only available when the private parameter _bank_rays_on_completion
759 : * is set to true.
760 : */
761 : const std::vector<std::shared_ptr<Ray>> & rayBank() const;
762 :
763 : /**
764 : * Gets the Ray with the ID \p ray_id from the Ray bank.
765 : *
766 : * If the Ray with \p ray_id is not found across all processors, this will error.
767 : *
768 : * This will ONLY return a valid Ray (not a null shared_ptr) on the processor that
769 : * has the Ray.
770 : */
771 : std::shared_ptr<Ray> getBankedRay(const RayID ray_id) const;
772 :
773 : /**
774 : * Resets the generation of unique RayIDs via generateUniqueRayID() to the beginning
775 : * of the range.
776 : */
777 : void resetUniqueRayIDs();
778 : /**
779 : * Resets the generation of unique replicated RayIDs accessed via generateReplicatedRayID().
780 : */
781 : void resetReplicatedRayIDs();
782 :
783 : /**
784 : * Gets all of the currently active RayTracingObjects
785 : */
786 : std::vector<RayTracingObject *> getRayTracingObjects();
787 :
788 : /**
789 : * Generates a unique RayID to be used for a Ray.
790 : *
791 : * This is used internally when acquiring new Rays.
792 : */
793 : virtual RayID generateUniqueRayID(const THREAD_ID tid);
794 : /**
795 : * Generates a Ray ID that is replicated across all processors.
796 : */
797 : RayID generateReplicatedRayID();
798 :
799 : /**
800 : * User APIs for constructing Rays within the RayTracingStudy.
801 : *
802 : * Rays can ONLY be constructed by users within the RayTracingStudy via the following methods.
803 : */
804 : ///@{
805 : /**
806 : * Acquire a Ray from the pool of Rays within generateRays().
807 : *
808 : * A unique ID is generated and assigned to the acquired Ray. The data and aux
809 : * data sizes are set according to the sizes required by the RayTracingStudy.
810 : */
811 : std::shared_ptr<Ray> acquireRay();
812 : /**
813 : * Acquire a Ray from the pool of Rays within generateRays(), without resizing
814 : * the data (sizes the data to zero).
815 : *
816 : * A unique ID is generated and assigned to the acquired Ray.
817 : */
818 : std::shared_ptr<Ray> acquireUnsizedRay();
819 : /**
820 : * Acquire a Ray from the pool of Rays within generateRays() in a replicated fashion.
821 : *
822 : * That is, this method must be called on all processors at the same time and the ID of the
823 : * resulting Ray is the same across all processors.
824 : *
825 : * The data and aux data sizes are set according to the sizes required by the RayTracingStudy.
826 : */
827 : std::shared_ptr<Ray> acquireReplicatedRay();
828 : /**
829 : * Acquires a Ray that that is copied from another Ray within generateRays().
830 : *
831 : * All of the information is copied except for the counters (intersections, processor crossings,
832 : * etc), which are reset.
833 : */
834 : std::shared_ptr<Ray> acquireCopiedRay(const Ray & ray);
835 : /**
836 : * Acquires a Ray with a given name within generateRays(). Used when ray registration is enabled,
837 : * that is, the private paramater '_use_ray_registration' == true.
838 : *
839 : * This method must be called on all processors at the same time with the same name.
840 : * This method can only be called on thread 0, which is why there is no thread argument.
841 : */
842 : std::shared_ptr<Ray> acquireRegisteredRay(const std::string & name);
843 : ///@}
844 :
845 : /**
846 : * Moves a ray to the buffer to be traced during generateRays().
847 : *
848 : * This moves the Ray into the buffer (with std::move), therefore \p ray will be nullptr after
849 : * this call.
850 : */
851 : void moveRayToBuffer(std::shared_ptr<Ray> & ray);
852 : /**
853 : * Moves rays to the buffer to be traced during generateRays().
854 : *
855 : * This moves the rays into the buffer (with std::move), therefore all valid rays
856 : * in \p rays will be nullptr after this call.
857 : */
858 : void moveRaysToBuffer(std::vector<std::shared_ptr<Ray>> & rays);
859 :
860 : /// The Mesh
861 : MooseMesh & _mesh;
862 : /// The Communicator
863 : const Parallel::Communicator & _comm;
864 : /// The rank of this processor (this actually takes time to lookup - so just do it once)
865 : const processor_id_type _pid;
866 :
867 : /// Whether or not to perform coverage checks on RayKernels
868 : const bool _ray_kernel_coverage_check;
869 : /// Whether not to warn if non-planar faces are found
870 : const bool _warn_non_planar;
871 : /// Whether or not to use Ray registration
872 : const bool _use_ray_registration;
873 : /// Whether or not to use the internal sidesets in ray tracing
874 : const bool _use_internal_sidesets;
875 : /// Whether or not to tolerate a Ray Tracing failure
876 : const bool _tolerate_failure;
877 : /// Whether or not to bank rays on completion
878 : const bool _bank_rays_on_completion;
879 : /// Whether or not subdomain setup is dependent on the Ray
880 : const bool _ray_dependent_subdomain_setup;
881 :
882 : /// Whether or not to cache traces on every trace execution
883 : bool _always_cache_traces;
884 : /// Whether or not to store the Ray data on the cache traces
885 : const bool _data_on_cache_traces;
886 : /// Whether or not to store the Ray aux data on the cache traces
887 : const bool _aux_data_on_cache_traces;
888 : /// Whether or not to cache individual element segments when caching
889 : const bool _segments_on_cache_traces;
890 : /// Max distance a Ray can travel before being killed (can change)
891 : const Real _ray_max_distance;
892 : /// Whether or not to verify if Rays have valid information before being traced
893 : const bool _verify_rays;
894 : /// Whether or not to verify the trace intersections in devel and dbg modes
895 : #ifndef NDEBUG
896 : const bool _verify_trace_intersections;
897 : #endif
898 :
899 : private:
900 : /**
901 : * Perform coverage checks (coverage of RayMaterials and RayKernels, if enabled)
902 : */
903 : void coverageChecks();
904 :
905 : /**
906 : * Perform checks to see if the listed dependencies in the RayTracingObjects exist
907 : */
908 : void dependencyChecks();
909 : /**
910 : * Verifies that the dependencies exist for a set of RayTracingObjects
911 : */
912 : void verifyDependenciesExist(const std::vector<RayTracingObject *> & rtos);
913 :
914 : /**
915 : * Check for if all of the element types in the mesh are supported by ray tracing
916 : */
917 : void traceableMeshChecks();
918 :
919 : /**
920 : * Check for overlapping PeriodicRayBC boundaries and check for cases in which
921 : * ghosting may not be sufficient with distributed mesh for overlapping
922 : * periodic boundaries.
923 : */
924 : void periodicBoundaryChecks();
925 :
926 : /**
927 : * Does the setup for internal sidesets. This includes:
928 : * - Setting if the local mesh has internal sidesets that have RayBCs and on which
929 : * boundary (stored in _internal_sidesets)
930 : * - Sets up the internal sideset mapping that maps element sides to the BoundaryIDs
931 : * on said side that have internal sidesets that have RayBCs
932 : * (stored in _internal_sidesets_map)
933 : *
934 : * If _use_internal_sidesets == false, this will still check to make sure that the
935 : * user does not have internal sidesets with RayBCs on them, and will report an error
936 : * if so.
937 : */
938 : void internalSidesetSetup();
939 :
940 : /**
941 : * Sets up the caching of whether or not each element side is non-planar, which is
942 : * stored in _non_planar_sides.
943 : *
944 : * This is done because this check is made within TraceRay, and determining whether
945 : * or not a side is non planar is quite expensive and we don't want to call it
946 : * mid-trace.
947 : */
948 : void nonPlanarSideSetup();
949 :
950 : /**
951 : * Sets up the _elem_index_helper, which is used for obtaining a contiguous index
952 : * for all elements that this processor knows about.
953 : *
954 : * Said index is used for quickly accessing things like internal sidesets and checking
955 : * if an element side is non-planar during tracing.
956 : */
957 : void localElemIndexSetup();
958 :
959 : /**
960 : * Sets up the maps from Ray to associated RayTracingObjects if _use_ray_registration.
961 : *
962 : * If ray registration is disabled, this makes sure no RayTracingObjects provided rays.
963 : */
964 : void registeredRaySetup();
965 :
966 : /**
967 : * Zero the AuxVariables that the registered AuxRayKernels contribute to.
968 : */
969 : void zeroAuxVariables();
970 :
971 : /**
972 : * Caches the hmax for all elements in each subdomain
973 : */
974 : void subdomainHMaxSetup();
975 :
976 : /**
977 : * Internal method for registering Ray data or Ray aux data with a name.
978 : */
979 : RayDataIndex registerRayDataInternal(const std::string & name, const bool aux);
980 : /**
981 : * Internal method for registering Ray data or Ray aux data with names.
982 : */
983 : std::vector<RayDataIndex> registerRayDataInternal(const std::vector<std::string> & names,
984 : const bool aux);
985 : /**
986 : * Internal method for getting the index of Ray data or Ray aux data.
987 : */
988 : RayDataIndex
989 : getRayDataIndexInternal(const std::string & name, const bool aux, const bool graceful) const;
990 : /**
991 : * Internal method for getting the indicies of Ray data or Ray aux data.
992 : */
993 : std::vector<RayDataIndex> getRayDataIndicesInternal(const std::vector<std::string> & names,
994 : const bool aux,
995 : const bool graceful) const;
996 :
997 : /**
998 : * Internal method for getting the name of Ray data or Ray aux data.
999 : */
1000 : const std::string & getRayDataNameInternal(const RayDataIndex index, const bool aux) const;
1001 :
1002 : /**
1003 : * Internal method for getting the value (replicated across all processors) in a Ray's data
1004 : * or aux data from the Ray banks.
1005 : */
1006 : RayData
1007 : getBankedRayDataInternal(const RayID ray_id, const RayDataIndex index, const bool aux) const;
1008 :
1009 : /**
1010 : * Registers a Ray with a given name.
1011 : *
1012 : * This is for internal use only and is an internal map for replicated Ray IDs -> names for
1013 : * easy access for studies that have a small number of rays to generate.
1014 : *
1015 : * @param name The name to register
1016 : * @return The allocated ID for the registered ray
1017 : */
1018 : RayID registerRay(const std::string & name);
1019 :
1020 : /// Threaded helpers for building element sides without extraneous allocation
1021 : std::vector<libMesh::ElemSideBuilder> _threaded_elem_side_builders;
1022 :
1023 : /// Timing
1024 : //@{
1025 : std::chrono::steady_clock::time_point _execution_start_time;
1026 : std::chrono::steady_clock::duration _execution_time;
1027 : std::chrono::steady_clock::duration _generation_time;
1028 : std::chrono::steady_clock::duration _propagation_time;
1029 : //@}
1030 :
1031 : /// The map from Ray data names to index
1032 : std::unordered_map<std::string, RayDataIndex> _ray_data_map;
1033 : /// The map from Ray aux data names to index
1034 : std::unordered_map<std::string, RayDataIndex> _ray_aux_data_map;
1035 :
1036 : /// The names for each Ray data entry
1037 : std::vector<std::string> _ray_data_names;
1038 : /// The names for each Ray aux data entry
1039 : std::vector<std::string> _ray_aux_data_names;
1040 :
1041 : /// Map from registered Ray name to ID
1042 : std::unordered_map<std::string, RayID> & _registered_ray_map;
1043 : /// Map from registered Ray ID to name
1044 : std::vector<std::string> & _reverse_registered_ray_map;
1045 :
1046 : /// Storage for the cached traces
1047 : std::vector<TraceData> _cached_traces;
1048 : /// The threaded storage for cached traces
1049 : std::vector<std::vector<TraceData>> _threaded_cached_traces;
1050 :
1051 : /// Number of currently cached objects for Jacobian/residual for each thread
1052 : std::vector<std::size_t> _num_cached;
1053 :
1054 : /// The BoundaryIDs on the local mesh that have internal RayBCs
1055 : std::set<BoundaryID> _internal_sidesets;
1056 : /// Internal sideset data, if internal sidesets exist (indexed with getLocalElemIndex())
1057 : std::vector<std::vector<std::vector<BoundaryID>>> _internal_sidesets_map;
1058 : /// Whether or not the local mesh has elements with non-planar sides
1059 : bool _has_non_planar_sides;
1060 : /// Non planar side data, which is for quick checking if an elem side is non-planar
1061 : /// We use unsigned short here to avoid a std::vector<bool>; 0 = false, otherwise true
1062 : std::vector<std::vector<unsigned short>> _non_planar_sides;
1063 : /// Whether or not the mesh has active elements of the same level
1064 : bool _has_same_level_active_elems;
1065 :
1066 : /// Nodal bounding box for the domain
1067 : libMesh::BoundingBox _b_box;
1068 : /// Loose nodal bounding box for the domain
1069 : libMesh::BoundingBox _loose_b_box;
1070 : /// An inflated max distance for the domain
1071 : Real _domain_max_length;
1072 : /// The total volume of the domain
1073 : Real _total_volume;
1074 :
1075 : /// Threaded cached subdomain query for RayKernelBase objects pertaining to this study
1076 : std::vector<TheWarehouse::QueryCache<AttribSubdomains>> _threaded_cache_ray_kernel;
1077 : /// Threaded cached boundary query for RayBC objects pertaining to this study
1078 : std::vector<TheWarehouse::QueryCache<AttribBoundaries>> _threaded_cache_ray_bc;
1079 : /// Threaded storage for all of the RayTracingObjects associated with a single Ray
1080 : std::vector<std::vector<std::set<const RayTracingObject *>>> _threaded_ray_object_registration;
1081 : /// The current RayKernel objects for each thread
1082 : std::vector<std::vector<RayKernelBase *>> _threaded_current_ray_kernels;
1083 : /// The TraceRay objects for each thread (they do the physical tracing)
1084 : std::vector<std::shared_ptr<TraceRay>> _threaded_trace_ray;
1085 : /// Face FE used for computing face normals for each thread
1086 : std::vector<std::unique_ptr<libMesh::FEBase>> _threaded_fe_face;
1087 : /// Face quadrature used for computing face normals for each thread
1088 : std::vector<std::unique_ptr<libMesh::QBase>> _threaded_q_face;
1089 : /// Threaded cache for side normals that have been computed already during tracing
1090 : std::vector<std::unordered_map<std::pair<const Elem *, unsigned short>, Point>>
1091 : _threaded_cached_normals;
1092 : /// Cumulative Ray bank - stored only when _bank_rays_on_completion
1093 : std::vector<std::shared_ptr<Ray>> _ray_bank;
1094 : /// Storage for the next available unique RayID, obtained via generateUniqueRayID()
1095 : std::vector<RayID> _threaded_next_ray_id;
1096 : /// Storage for the next available replicated RayID, obtained via generateReplicatedRayID()
1097 : RayID _replicated_next_ray_id;
1098 :
1099 : /// The study that used is to actually execute (trace) the Rays
1100 : const std::unique_ptr<ParallelRayStudy> _parallel_ray_study;
1101 :
1102 : /// Quadrature rule for laying points across a 1D ray segment
1103 : std::unique_ptr<libMesh::QBase> _segment_qrule;
1104 :
1105 : /// Total number of processor crossings for Rays that finished on this processor
1106 : unsigned long long int _ending_processor_crossings;
1107 : /// Max number of total processor crossings for Rays that finished on this processor
1108 : unsigned int _ending_max_processor_crossings;
1109 : /// Total number of processor crossings
1110 : unsigned long long int _total_processor_crossings;
1111 : /// Max number of processor crossings for all Rays
1112 : unsigned int _max_processor_crossings;
1113 :
1114 : /// Total number of Ray/element intersections for Rays that finished on this processor
1115 : unsigned long long int _ending_intersections;
1116 : /// Max number of intersections for Rays that finished on this processor
1117 : unsigned int _ending_max_intersections;
1118 : /// Max number of trajectory changes for Rays that finished on this processor
1119 : unsigned int _ending_max_trajectory_changes;
1120 : /// Total number of Ray/element intersections
1121 : unsigned long long int _total_intersections;
1122 : /// Max number of intersections for a single Ray
1123 : unsigned int _max_intersections;
1124 : /// Max number of trajectory changes for a single Ray
1125 : unsigned int _max_trajectory_changes;
1126 :
1127 : /// Total distance traveled by Rays that end on this processor
1128 : Real _ending_distance;
1129 : /// Total distance traveled by all Rays
1130 : Real _total_distance;
1131 :
1132 : /// Cumulative results on this processor from the threaded TraceRay objects
1133 : std::vector<unsigned long long int> _local_trace_ray_results;
1134 :
1135 : /// The cached hmax for all elements in a subdomain
1136 : std::unordered_map<SubdomainID, Real> _subdomain_hmax;
1137 :
1138 : /// Whether or not we've called initial setup - used to stop from late registration
1139 : bool _called_initial_setup;
1140 :
1141 : /// Helper for defining a local contiguous index for each element
1142 : ElemIndexHelper _elem_index_helper;
1143 :
1144 : /// Spin mutex object for locks
1145 : mutable Threads::spin_mutex _spin_mutex;
1146 : };
|
