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 : #include "RayTracingStudy.h"
11 :
12 : #include "AuxRayKernel.h"
13 : #include "RayBoundaryConditionBase.h"
14 : #include "RayKernel.h"
15 : #include "TraceRay.h"
16 : #include "TraceRayTools.h"
17 : #include "PeriodicRayBC.h"
18 :
19 : #include "AuxiliarySystem.h"
20 : #include "Assembly.h"
21 : #include "NonlinearSystemBase.h"
22 :
23 : #include "libmesh/enum_to_string.h"
24 : #include "libmesh/mesh_tools.h"
25 : #include "libmesh/parallel_sync.h"
26 : #include "libmesh/remote_elem.h"
27 : #include "libmesh/periodic_boundary.h"
28 : #include "libmesh/periodic_boundaries.h"
29 :
30 : using namespace libMesh;
31 :
32 : InputParameters
33 3980 : RayTracingStudy::validParams()
34 : {
35 3980 : auto params = GeneralUserObject::validParams();
36 :
37 : // Parameters for the execution of the rays
38 3980 : params += ParallelRayStudy::validParams();
39 :
40 11940 : params.addRangeCheckedParam<Real>("ray_distance",
41 7960 : std::numeric_limits<Real>::max(),
42 : "ray_distance > 0",
43 : "The maximum distance all Rays can travel");
44 :
45 7960 : params.addParam<bool>(
46 7960 : "tolerate_failure", false, "Whether or not to tolerate a ray tracing failure");
47 :
48 7960 : MooseEnum work_buffers("lifo circular", "circular");
49 7960 : params.addParam<MooseEnum>("work_buffer_type", work_buffers, "The work buffer type to use");
50 :
51 7960 : params.addParam<bool>(
52 7960 : "ray_kernel_coverage_check", true, "Whether or not to perform coverage checks on RayKernels");
53 7960 : params.addParam<bool>("warn_non_planar",
54 7960 : true,
55 : "Whether or not to produce a warning if any element faces are non-planar.");
56 :
57 7960 : params.addParam<bool>(
58 : "always_cache_traces",
59 7960 : false,
60 : "Whether or not to cache the Ray traces on every execution, primarily for use in output. "
61 : "Warning: this can get expensive very quick with a large number of rays!");
62 7960 : params.addParam<bool>("data_on_cache_traces",
63 7960 : false,
64 : "Whether or not to also cache the Ray's data when caching its traces");
65 7960 : params.addParam<bool>("aux_data_on_cache_traces",
66 7960 : false,
67 : "Whether or not to also cache the Ray's aux data when caching its traces");
68 7960 : params.addParam<bool>(
69 : "segments_on_cache_traces",
70 7960 : true,
71 : "Whether or not to cache individual segments when trace caching is enabled. If false, we "
72 : "will instead cache a segment for each part of the trace where the direction is the same. "
73 : "This minimizes the number of segments requied to represent the Ray's path, but removes the "
74 : "ability to show Ray field data on each segment through an element.");
75 :
76 7960 : params.addParam<bool>("use_internal_sidesets",
77 7960 : false,
78 : "Whether or not to use internal sidesets for RayBCs in ray tracing");
79 :
80 7960 : params.addParam<bool>("warn_subdomain_hmax",
81 7960 : true,
82 : "Whether or not to warn if the approximated hmax (constant on subdomain) "
83 : "varies significantly for an element");
84 :
85 7960 : params.addParam<bool>(
86 : "verify_rays",
87 7960 : true,
88 : "Whether or not to verify the generated Rays. This includes checking their "
89 : "starting information and the uniqueness of Rays before and after execution. This is also "
90 : "used by derived studies for more specific verification.");
91 7960 : params.addParam<bool>("verify_trace_intersections",
92 7960 : true,
93 : "Whether or not to verify the trace intersections in devel and dbg modes. "
94 : "Trace intersections are not verified regardless of this parameter in "
95 : "optimized modes (opt, oprof).");
96 :
97 7960 : params.addParam<bool>("allow_other_flags_with_prekernels",
98 7960 : false,
99 : "Whether or not to allow the list of execution flags to have PRE_KERNELS "
100 : "mixed with other flags. If this parameter is not set then if PRE_KERNELS "
101 : "is provided it must be the only execution flag.");
102 :
103 3980 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
104 3980 : exec_enum.addAvailableFlags(EXEC_PRE_KERNELS);
105 :
106 7960 : params.addParamNamesToGroup(
107 : "always_cache_traces data_on_cache_traces aux_data_on_cache_traces segments_on_cache_traces",
108 : "Trace cache");
109 7960 : params.addParamNamesToGroup("warn_non_planar warn_subdomain_hmax", "Tracing Warnings");
110 7960 : params.addParamNamesToGroup("ray_kernel_coverage_check verify_rays verify_trace_intersections",
111 : "Checks and verifications");
112 :
113 : // Whether or not each Ray must be registered using the registerRay() API
114 3980 : params.addPrivateParam<bool>("_use_ray_registration", true);
115 : // Whether or not to bank Rays on completion
116 3980 : params.addPrivateParam<bool>("_bank_rays_on_completion", true);
117 : /// Whether or not subdomain setup is dependent on the Ray
118 7960 : params.addPrivateParam<bool>("_ray_dependent_subdomain_setup", true);
119 :
120 : // Add a point neighbor relationship manager
121 7960 : params.addRelationshipManager("ElementPointNeighborLayers",
122 : Moose::RelationshipManagerType::GEOMETRIC |
123 : Moose::RelationshipManagerType::ALGEBRAIC,
124 5891 : [](const InputParameters &, InputParameters & rm_params)
125 5891 : { rm_params.set<unsigned short>("layers") = 1; });
126 :
127 3980 : return params;
128 3980 : }
129 :
130 2029 : RayTracingStudy::RayTracingStudy(const InputParameters & parameters)
131 : : GeneralUserObject(parameters),
132 2029 : _mesh(_fe_problem.mesh()),
133 2029 : _comm(_mesh.comm()),
134 2029 : _pid(_comm.rank()),
135 :
136 4058 : _ray_kernel_coverage_check(getParam<bool>("ray_kernel_coverage_check")),
137 4058 : _warn_non_planar(getParam<bool>("warn_non_planar")),
138 4058 : _use_ray_registration(getParam<bool>("_use_ray_registration")),
139 4058 : _use_internal_sidesets(getParam<bool>("use_internal_sidesets")),
140 4058 : _tolerate_failure(getParam<bool>("tolerate_failure")),
141 4058 : _bank_rays_on_completion(getParam<bool>("_bank_rays_on_completion")),
142 4058 : _ray_dependent_subdomain_setup(getParam<bool>("_ray_dependent_subdomain_setup")),
143 :
144 4058 : _always_cache_traces(getParam<bool>("always_cache_traces")),
145 4058 : _data_on_cache_traces(getParam<bool>("data_on_cache_traces")),
146 4058 : _aux_data_on_cache_traces(getParam<bool>("aux_data_on_cache_traces")),
147 4058 : _segments_on_cache_traces(getParam<bool>("segments_on_cache_traces")),
148 4058 : _ray_max_distance(getParam<Real>("ray_distance")),
149 4058 : _verify_rays(getParam<bool>("verify_rays")),
150 : #ifndef NDEBUG
151 : _verify_trace_intersections(getParam<bool>("verify_trace_intersections")),
152 : #endif
153 :
154 2029 : _threaded_elem_side_builders(libMesh::n_threads()),
155 :
156 2029 : _registered_ray_map(
157 2029 : declareRestartableData<std::unordered_map<std::string, RayID>>("registered_ray_map")),
158 2029 : _reverse_registered_ray_map(
159 4058 : declareRestartableData<std::vector<std::string>>("reverse_registered_ray_map")),
160 :
161 2029 : _threaded_cached_traces(libMesh::n_threads()),
162 :
163 2029 : _num_cached(libMesh::n_threads(), 0),
164 :
165 2029 : _has_non_planar_sides(true),
166 2029 : _has_same_level_active_elems(sameLevelActiveElems()),
167 :
168 2029 : _b_box(MeshTools::create_nodal_bounding_box(_mesh.getMesh())),
169 2029 : _domain_max_length(1.01 * (_b_box.max() - _b_box.min()).norm()),
170 2029 : _total_volume(computeTotalVolume()),
171 :
172 2029 : _threaded_cache_ray_kernel(libMesh::n_threads()),
173 2029 : _threaded_cache_ray_bc(libMesh::n_threads()),
174 2029 : _threaded_ray_object_registration(libMesh::n_threads()),
175 2029 : _threaded_current_ray_kernels(libMesh::n_threads()),
176 2029 : _threaded_trace_ray(libMesh::n_threads()),
177 2029 : _threaded_fe_face(libMesh::n_threads()),
178 2029 : _threaded_q_face(libMesh::n_threads()),
179 2029 : _threaded_cached_normals(libMesh::n_threads()),
180 2029 : _threaded_next_ray_id(libMesh::n_threads()),
181 :
182 2029 : _parallel_ray_study(std::make_unique<ParallelRayStudy>(*this, _threaded_trace_ray)),
183 :
184 2029 : _local_trace_ray_results(TraceRay::FAILED_TRACES + 1, 0),
185 :
186 2029 : _called_initial_setup(false),
187 :
188 8116 : _elem_index_helper(_mesh.getMesh(), name() + "_elem_index")
189 : {
190 : // Initialize a tracing object for each thread
191 4395 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
192 : {
193 : // Initialize a tracing object for each thread
194 4732 : _threaded_trace_ray[tid] = std::make_shared<TraceRay>(*this, tid);
195 :
196 : // Setup the face FEs for normal computation on the fly
197 4732 : _threaded_fe_face[tid] = FEBase::build(_mesh.dimension(), FEType(CONSTANT, MONOMIAL));
198 4732 : _threaded_q_face[tid] = QBase::build(QGAUSS, _mesh.dimension() - 1, CONSTANT);
199 : _threaded_fe_face[tid]->attach_quadrature_rule(_threaded_q_face[tid].get());
200 : _threaded_fe_face[tid]->get_normals();
201 : }
202 :
203 : // Evaluating on residual and Jacobian evaluation
204 2029 : if (_execute_enum.isValueSet(EXEC_PRE_KERNELS))
205 : {
206 282 : if (!getParam<bool>("allow_other_flags_with_prekernels") && _execute_enum.size() > 1)
207 2 : paramError("execute_on",
208 : "PRE_KERNELS cannot be mixed with any other execution flag.\nThat is, you cannot "
209 : "currently "
210 : "mix RayKernels that contribute to the Jacobian/residual with those that do not.");
211 :
212 92 : if (_app.useEigenvalue())
213 2 : mooseError("Execution on residual and Jacobian evaluation (execute_on = PRE_KERNELS)\n",
214 : "is not supported for an eigenvalue solve.");
215 : }
216 :
217 2025 : resetUniqueRayIDs();
218 2025 : resetReplicatedRayIDs();
219 :
220 : // Scale the bounding box for loose checking
221 : _loose_b_box = _b_box;
222 2025 : _loose_b_box.scale(TOLERANCE * TOLERANCE);
223 2025 : }
224 :
225 : void
226 1911 : RayTracingStudy::initialSetup()
227 : {
228 : // Keep track of initialSetup call to avoid registration of various things
229 1911 : _called_initial_setup = true;
230 :
231 : // Sets up a local index for each elem this proc knows about
232 1911 : localElemIndexSetup();
233 :
234 : // Check for RayKernel coverage
235 1911 : coverageChecks();
236 :
237 : // Make sure the dependencies exist, if any
238 1909 : dependencyChecks();
239 :
240 : // Check for traceable element types
241 1905 : traceableMeshChecks();
242 :
243 : // Check for sane periodic boundaries
244 1903 : periodicBoundaryChecks();
245 :
246 : // Setup for internal sidesets
247 1901 : internalSidesetSetup();
248 :
249 1897 : nonPlanarSideSetup();
250 :
251 : // Setup approximate hmax for each subdomain
252 1895 : subdomainHMaxSetup();
253 :
254 : // Call initial setup on all of the objects
255 5952 : for (auto & rto : getRayTracingObjects())
256 5952 : rto->initialSetup();
257 :
258 : // Check for proper exec flags with RayKernels
259 : std::vector<RayKernelBase *> ray_kernels;
260 1893 : getRayKernels(ray_kernels, 0);
261 3475 : for (const auto & rkb : ray_kernels)
262 1584 : if (dynamic_cast<RayKernel *>(rkb) && !_execute_enum.isValueSet(EXEC_PRE_KERNELS))
263 2 : mooseError("This study has RayKernel objects that contribute to residuals and Jacobians.",
264 : "\nIn this case, the study must use the execute_on = PRE_KERNELS");
265 :
266 : // Build 1D quadrature rule for along a segment
267 : _segment_qrule =
268 3782 : QBase::build(QGAUSS, 1, _fe_problem.getSystemBase(_sys.number()).getMinQuadratureOrder());
269 1891 : }
270 :
271 : void
272 3328 : RayTracingStudy::residualSetup()
273 : {
274 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
275 : mooseAssert(_num_cached[tid] == 0, "Cached residuals/Jacobians not empty");
276 :
277 8595 : for (auto & rto : getRayTracingObjects())
278 8595 : rto->residualSetup();
279 3328 : }
280 :
281 : void
282 180 : RayTracingStudy::jacobianSetup()
283 : {
284 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
285 : mooseAssert(_num_cached[tid] == 0, "Cached residuals/Jacobians not empty");
286 :
287 636 : for (auto & rto : getRayTracingObjects())
288 636 : rto->jacobianSetup();
289 180 : }
290 :
291 : void
292 1955 : RayTracingStudy::timestepSetup()
293 : {
294 5923 : for (auto & rto : getRayTracingObjects())
295 5923 : rto->timestepSetup();
296 1955 : }
297 :
298 : void
299 132 : RayTracingStudy::meshChanged()
300 : {
301 132 : localElemIndexSetup();
302 132 : internalSidesetSetup();
303 132 : nonPlanarSideSetup();
304 132 : subdomainHMaxSetup();
305 :
306 132 : _has_same_level_active_elems = sameLevelActiveElems();
307 :
308 282 : for (const auto & trace_ray : _threaded_trace_ray)
309 150 : trace_ray->meshChanged();
310 132 : }
311 :
312 : void
313 5360 : RayTracingStudy::execute()
314 : {
315 5360 : executeStudy();
316 5260 : }
317 :
318 : void
319 1911 : RayTracingStudy::coverageChecks()
320 : {
321 : // Check for coverage of RayKernels on domain
322 1911 : if (_ray_kernel_coverage_check)
323 : {
324 : std::vector<RayKernelBase *> ray_kernels;
325 1406 : getRayKernels(ray_kernels, 0);
326 :
327 : std::set<SubdomainID> ray_kernel_blocks;
328 2978 : for (const auto & rk : ray_kernels)
329 3144 : ray_kernel_blocks.insert(rk->blockIDs().begin(), rk->blockIDs().end());
330 :
331 : std::set<SubdomainID> missing;
332 2812 : std::set_difference(_mesh.meshSubdomains().begin(),
333 1406 : _mesh.meshSubdomains().end(),
334 : ray_kernel_blocks.begin(),
335 : ray_kernel_blocks.end(),
336 : std::inserter(missing, missing.begin()));
337 :
338 1406 : if (!missing.empty() && !ray_kernel_blocks.count(Moose::ANY_BLOCK_ID))
339 : {
340 2 : std::ostringstream error;
341 2 : error << "Subdomains { ";
342 2 : std::copy(missing.begin(), missing.end(), std::ostream_iterator<SubdomainID>(error, " "));
343 2 : error << "} do not have RayKernels defined!";
344 :
345 2 : mooseError(error.str());
346 0 : }
347 1404 : }
348 1909 : }
349 :
350 : void
351 1909 : RayTracingStudy::dependencyChecks()
352 : {
353 : std::vector<RayTracingObject *> ray_tracing_objects;
354 :
355 1909 : getRayKernels(ray_tracing_objects, 0);
356 1909 : verifyDependenciesExist(ray_tracing_objects);
357 :
358 1907 : getRayBCs(ray_tracing_objects, 0);
359 1907 : verifyDependenciesExist(ray_tracing_objects);
360 1905 : }
361 :
362 : void
363 3816 : RayTracingStudy::verifyDependenciesExist(const std::vector<RayTracingObject *> & rtos)
364 : {
365 7429 : for (const auto & rto : rtos)
366 3677 : for (const auto & dep_name : rto->getRequestedItems())
367 : {
368 : bool found = false;
369 136 : for (const auto & rto_search : rtos)
370 132 : if (rto_search->name() == dep_name)
371 : {
372 : found = true;
373 : break;
374 : }
375 :
376 64 : if (!found)
377 4 : rto->paramError("depends_on", "The ", rto->getBase(), " '", dep_name, "' does not exist");
378 : }
379 3812 : }
380 :
381 : void
382 1905 : RayTracingStudy::traceableMeshChecks()
383 : {
384 270851 : for (const auto & elem : *_mesh.getActiveLocalElementRange())
385 : {
386 268948 : if (_fe_problem.adaptivity().isOn())
387 : {
388 3310 : if (!TraceRayTools::isAdaptivityTraceableElem(elem))
389 2 : mooseError("Element type ",
390 2 : Utility::enum_to_string(elem->type()),
391 : " is not supported in ray tracing with adaptivity");
392 : }
393 265638 : else if (!TraceRayTools::isTraceableElem(elem))
394 0 : mooseError("Element type ",
395 0 : Utility::enum_to_string(elem->type()),
396 : " is not supported in ray tracing");
397 : }
398 1903 : }
399 :
400 : void
401 1903 : RayTracingStudy::periodicBoundaryChecks()
402 : {
403 : // Collect the PeriodicRayBCs
404 : std::vector<const RayBoundaryConditionBase *> rbc_ptrs;
405 1903 : getRayBCs(rbc_ptrs, 0);
406 : std::vector<const PeriodicRayBC *> prbc_ptrs;
407 3930 : for (const auto rbc_ptr : rbc_ptrs)
408 2027 : if (const auto prbc_ptr = dynamic_cast<const PeriodicRayBC *>(rbc_ptr))
409 27 : prbc_ptrs.push_back(prbc_ptr);
410 1903 : if (prbc_ptrs.empty())
411 : return;
412 :
413 : // Collect each of the periodic boundaries
414 : std::map<boundary_id_type,
415 : std::tuple<const PeriodicRayBC *,
416 : const libMesh::PeriodicBoundaryBase *,
417 : std::unordered_set<dof_id_type>>>
418 : boundary_map;
419 50 : for (const auto prbc_ptr : prbc_ptrs)
420 : {
421 123 : for (const auto & [bid, pb] : prbc_ptr->getPeriodicBoundaries())
422 : {
423 98 : const auto [it, inserted] = boundary_map.emplace(
424 : std::piecewise_construct,
425 98 : std::tuple{bid},
426 98 : std::forward_as_tuple(prbc_ptr, pb.get(), std::unordered_set<dof_id_type>()));
427 98 : if (!inserted)
428 2 : prbc_ptr->mooseError("The periodic boundary '",
429 2 : _mesh.getBoundaryString(bid),
430 : "' has been defined in both ",
431 2 : prbc_ptr->typeAndName(),
432 : " and ",
433 2 : std::get<0>(it->second)->typeAndName());
434 : }
435 : }
436 :
437 : // Because we don't have ghosting setup correctly yet, we need to check if any
438 : // of the periodic boundaries are neighbors with distributed mesh. If the
439 : // mesh is replicated, we don't need to check this. See #31280.
440 23 : if (comm().size() == 1 || !_mesh.isDistributedMesh())
441 : return;
442 :
443 : // Collect all of the nodes that are on each periodic boundary
444 2 : const auto & sideset_map = _mesh.getMesh().get_boundary_info().get_sideset_map();
445 6 : for (const auto & [elem, side_bid_pair] : sideset_map)
446 : {
447 4 : const auto [side, bid] = side_bid_pair;
448 4 : if (auto it = boundary_map.find(bid); it != boundary_map.end())
449 8 : for (const auto n : elem->nodes_on_side(side))
450 8 : std::get<2>(it->second).insert(elem->node_ref(n).id());
451 : }
452 :
453 : // Distributed meshes have distributed boundary information, so sync
454 6 : for (auto & bid_tuple_pair : boundary_map)
455 4 : comm().set_union(std::get<2>(bid_tuple_pair.second));
456 :
457 : // Check for periodic boundaries that share nodes
458 : std::map<std::pair<boundary_id_type, boundary_id_type>,
459 : std::pair<const PeriodicRayBC *, const PeriodicRayBC *>>
460 : warn_boundaries;
461 6 : for (auto it = boundary_map.begin(); it != boundary_map.end(); ++it)
462 : {
463 : const auto & [bid, tup] = *it;
464 : const auto [prbc_ptr, pb, node_ids] = tup;
465 :
466 6 : for (auto other_it = std::next(it); other_it != boundary_map.end(); ++other_it)
467 : {
468 : const auto & [other_bid, other_tup] = *other_it;
469 :
470 : // Don't check against boundaries that are paried together
471 2 : if (pb->pairedboundary == other_bid)
472 : continue;
473 :
474 0 : const auto other_prbc_ptr = std::get<0>(other_tup);
475 : const auto & other_node_ids = std::get<2>(other_tup);
476 0 : for (const auto node_id : node_ids)
477 : if (other_node_ids.count(node_id))
478 : {
479 0 : if (!warn_boundaries.count(std::make_pair(other_bid, bid)))
480 0 : warn_boundaries.emplace(std::make_pair(bid, other_bid),
481 0 : std::make_pair(prbc_ptr, other_prbc_ptr));
482 : break;
483 : }
484 : }
485 : }
486 :
487 2 : if (warn_boundaries.size())
488 : {
489 0 : std::ostringstream oss;
490 : oss << warn_boundaries.size()
491 0 : << " ray tracing periodic boundaries were found to be neighbors:\n\n";
492 0 : for (const auto & [bids_pair, prbc_ptrs_pair] : warn_boundaries)
493 : {
494 0 : const auto [bid, paired_bid] = bids_pair;
495 0 : const auto [prbc_ptr, paired_prbc_ptr] = prbc_ptrs_pair;
496 0 : oss << " '" << _mesh.getBoundaryString(bid) << "' (in " << prbc_ptr->typeAndName()
497 0 : << ") <-> '" << _mesh.getBoundaryString(paired_bid) << "' (in "
498 0 : << paired_prbc_ptr->typeAndName() << ")\n";
499 : }
500 : oss << "\nThe periodic propagation of rays at points where two or more periodic"
501 : << "\nboundaries meet is not fully supported with a distributed mesh."
502 0 : << "\n\nIf you encounter trace failures, you should use a replicated mesh.";
503 0 : mooseWarning(oss.str());
504 0 : }
505 1901 : }
506 :
507 : void
508 2043 : RayTracingStudy::localElemIndexSetup()
509 : {
510 : // TODO: We could probably minimize this to local active elements followed by
511 : // boundary point neighbors, but if using distibuted mesh it really shouldn't matter
512 2043 : _elem_index_helper.initialize(_mesh.getMesh().active_element_ptr_range());
513 2043 : }
514 :
515 : void
516 2033 : RayTracingStudy::internalSidesetSetup()
517 : {
518 : // Even if we have _use_internal_sidesets == false, we will make sure the user didn't add RayBCs
519 : // on internal boundaries
520 :
521 : // Clear the data structures and size the map based on the elements that we know about
522 : _internal_sidesets.clear();
523 2033 : _internal_sidesets_map.clear();
524 2033 : _internal_sidesets_map.resize(_elem_index_helper.maxIndex() + 1);
525 :
526 : // First, we are going to store all elements with internal sidesets (if any) that have active
527 : // RayBCs on them as elem -> vector of (side, vector of boundary ids)
528 148361 : for (const auto & bnd_elem : *_mesh.getBoundaryElementRange())
529 : {
530 146330 : Elem * elem = bnd_elem->_elem;
531 146330 : const unsigned int side = bnd_elem->_side;
532 146330 : const auto bnd_id = bnd_elem->_bnd_id;
533 :
534 : // Not internal
535 : const Elem * const neighbor = elem->neighbor_ptr(side);
536 146330 : if (!neighbor || neighbor == remote_elem)
537 142760 : continue;
538 :
539 : // No RayBCs on this sideset
540 : std::vector<RayBoundaryConditionBase *> result;
541 3570 : getRayBCs(result, bnd_id, 0);
542 3570 : if (result.empty())
543 : continue;
544 :
545 3570 : if (neighbor->subdomain_id() == elem->subdomain_id())
546 2 : mooseError("RayBCs exist on internal sidesets that are not bounded by a different",
547 : "\nsubdomain on each side.",
548 : "\n\nIn order to use RayBCs on internal sidesets, said sidesets must have",
549 : "\na different subdomain on each side.");
550 :
551 : // Mark that this boundary is an internal sideset with RayBC(s)
552 3568 : _internal_sidesets.insert(bnd_id);
553 :
554 : // Get elem's entry in the internal sidset data structure
555 : const auto index = _elem_index_helper.getIndex(elem);
556 : auto & entry = _internal_sidesets_map[index];
557 :
558 : // Initialize this elem's sides if they have not been already
559 3568 : if (entry.empty())
560 2772 : entry.resize(elem->n_sides(), std::vector<BoundaryID>());
561 :
562 : // Add the internal boundary to the side entry
563 3568 : entry[side].push_back(bnd_id);
564 3568 : }
565 :
566 2031 : if (!_use_internal_sidesets && _internal_sidesets.size())
567 2 : mooseError("RayBCs are defined on internal sidesets, but the study is not set to use ",
568 : "internal sidesets during tracing.",
569 : "\n\nSet the parameter use_internal_sidesets = true to enable this capability.");
570 2029 : }
571 :
572 : void
573 2029 : RayTracingStudy::nonPlanarSideSetup()
574 : {
575 2029 : _has_non_planar_sides = false;
576 2029 : bool warned = !_warn_non_planar;
577 :
578 : // Nothing to do here for 2D or 1D
579 2029 : if (_mesh.dimension() != 3)
580 : return;
581 :
582 : // Clear the data structure and size it based on the elements that we know about
583 765 : _non_planar_sides.clear();
584 765 : _non_planar_sides.resize(_elem_index_helper.maxIndex() + 1);
585 :
586 907772 : for (const Elem * elem : _mesh.getMesh().active_element_ptr_range())
587 : {
588 : const auto index = _elem_index_helper.getIndex(elem);
589 : auto & entry = _non_planar_sides[index];
590 302082 : entry.resize(elem->n_sides(), 0);
591 :
592 1713816 : for (const auto s : elem->side_index_range())
593 : {
594 1411736 : const auto & side = elemSide(*elem, s);
595 1411736 : if (side.n_vertices() < 4)
596 1018080 : continue;
597 :
598 393656 : if (!side.has_affine_map())
599 : {
600 6318 : entry[s] = 1;
601 6318 : _has_non_planar_sides = true;
602 :
603 6318 : if (!warned)
604 : {
605 2 : mooseWarning("The mesh contains non-planar faces.\n\n",
606 : "Ray tracing on non-planar faces is an approximation and may fail.\n\n",
607 : "Use at your own risk! You can disable this warning by setting the\n",
608 : "parameter 'warn_non_planar' to false.");
609 : warned = true;
610 : }
611 : }
612 : }
613 763 : }
614 : }
615 :
616 : void
617 2027 : RayTracingStudy::subdomainHMaxSetup()
618 : {
619 : // Setup map with subdomain keys
620 : _subdomain_hmax.clear();
621 4427 : for (const auto subdomain_id : _mesh.meshSubdomains())
622 2400 : _subdomain_hmax[subdomain_id] = std::numeric_limits<Real>::min();
623 :
624 : // Set local max for each subdomain
625 275965 : for (const auto & elem : *_mesh.getActiveLocalElementRange())
626 : {
627 273938 : auto & entry = _subdomain_hmax.at(elem->subdomain_id());
628 276655 : entry = std::max(entry, elem->hmax());
629 : }
630 :
631 : // Accumulate global max for each subdomain
632 2027 : _communicator.max(_subdomain_hmax);
633 :
634 4054 : if (getParam<bool>("warn_subdomain_hmax"))
635 : {
636 2027 : const auto warn_prefix = type() + " '" + name() + "': ";
637 2027 : const auto warn_suffix =
638 : "\n\nRay tracing uses an approximate element size for each subdomain to scale the\n"
639 : "tolerances used in computing ray intersections. This warning suggests that the\n"
640 : "approximate element size is not a good approximation. This is likely due to poor\n"
641 : "element aspect ratios.\n\n"
642 : "This warning is only output for the first element affected.\n"
643 : "To disable this warning, set warn_subdomain_hmax = false.\n";
644 :
645 275963 : for (const auto & elem : *_mesh.getActiveLocalElementRange())
646 : {
647 273938 : const auto hmin = elem->hmin();
648 273938 : const auto hmax = elem->hmax();
649 273938 : const auto max_hmax = subdomainHmax(elem->subdomain_id());
650 :
651 273938 : const auto hmax_rel = hmax / max_hmax;
652 273938 : if (hmax_rel < 1.e-2 || hmax_rel > 1.e2)
653 0 : mooseDoOnce(mooseWarning(warn_prefix,
654 : "Element hmax varies significantly from subdomain hmax.\n",
655 : warn_suffix,
656 : "First element affected:\n",
657 : Moose::stringify(*elem)););
658 :
659 273938 : const auto h_rel = max_hmax / hmin;
660 273938 : if (h_rel > 1.e2)
661 2 : mooseDoOnce(mooseWarning(warn_prefix,
662 : "Element hmin varies significantly from subdomain hmax.\n",
663 : warn_suffix,
664 : "First element affected:\n",
665 : Moose::stringify(*elem)););
666 : }
667 : }
668 2025 : }
669 :
670 : void
671 5294 : RayTracingStudy::registeredRaySetup()
672 : {
673 : // First, clear the objects associated with each Ray on each thread
674 5294 : const auto num_rays = _registered_ray_map.size();
675 11396 : for (auto & entry : _threaded_ray_object_registration)
676 : {
677 6102 : entry.clear();
678 6102 : entry.resize(num_rays);
679 : }
680 :
681 5294 : const auto rtos = getRayTracingObjects();
682 :
683 5294 : if (_use_ray_registration)
684 : {
685 : // All of the registered ray names - used when a RayTracingObject did not specify
686 : // any Rays so it should be associated with all Rays.
687 : std::vector<std::string> all_ray_names;
688 3563 : all_ray_names.reserve(_registered_ray_map.size());
689 11622 : for (const auto & pair : _registered_ray_map)
690 8059 : all_ray_names.push_back(pair.first);
691 :
692 8974 : for (auto & rto : rtos)
693 : {
694 : // The Ray names associated with this RayTracingObject
695 5413 : const auto & ray_names = rto->parameters().get<std::vector<std::string>>("rays");
696 : // The registration for RayTracingObjects for the thread rto is on
697 5413 : const auto tid = rto->parameters().get<THREAD_ID>("_tid");
698 5413 : auto & registration = _threaded_ray_object_registration[tid];
699 :
700 : // Register each Ray for this object in the registration
701 17197 : for (const auto & ray_name : (ray_names.empty() ? all_ray_names : ray_names))
702 : {
703 10401 : const auto id = registeredRayID(ray_name, /* graceful = */ true);
704 10401 : if (ray_names.size() && id == Ray::INVALID_RAY_ID)
705 2 : rto->paramError(
706 2 : "rays", "Supplied ray '", ray_name, "' is not a registered Ray in ", typeAndName());
707 10399 : registration[id].insert(rto);
708 : }
709 : }
710 3561 : }
711 : // Not using Ray registration
712 : else
713 : {
714 5205 : for (const auto & rto : rtos)
715 3476 : if (rto->parameters().get<std::vector<std::string>>("rays").size())
716 2 : rto->paramError(
717 : "rays",
718 : "Rays cannot be supplied when the study does not require Ray registration.\n\n",
719 : type(),
720 : " does not require Ray registration.");
721 : }
722 5290 : }
723 :
724 : void
725 5360 : RayTracingStudy::zeroAuxVariables()
726 : {
727 : std::set<std::string> vars_to_be_zeroed;
728 : std::vector<RayKernelBase *> ray_kernels;
729 5360 : getRayKernels(ray_kernels, 0);
730 10832 : for (auto & rk : ray_kernels)
731 : {
732 5472 : AuxRayKernel * aux_rk = dynamic_cast<AuxRayKernel *>(rk);
733 5472 : if (aux_rk)
734 40 : vars_to_be_zeroed.insert(aux_rk->variable().name());
735 : }
736 :
737 : std::vector<std::string> vars_to_be_zeroed_vec(vars_to_be_zeroed.begin(),
738 5360 : vars_to_be_zeroed.end());
739 5360 : _fe_problem.getAuxiliarySystem().zeroVariables(vars_to_be_zeroed_vec);
740 10720 : }
741 :
742 : void
743 2997087 : RayTracingStudy::segmentSubdomainSetup(const SubdomainID subdomain,
744 : const THREAD_ID tid,
745 : const RayID ray_id)
746 : {
747 : mooseAssert(currentlyPropagating(), "Should not call while not propagating");
748 :
749 : // Call subdomain setup on FE
750 2997087 : _fe_problem.subdomainSetup(subdomain, tid);
751 :
752 : std::set<MooseVariableFEBase *> needed_moose_vars;
753 : std::unordered_set<unsigned int> needed_mat_props;
754 :
755 : // Get RayKernels and their dependencies and call subdomain setup
756 2997087 : getRayKernels(_threaded_current_ray_kernels[tid], subdomain, tid, ray_id);
757 5998762 : for (auto & rkb : _threaded_current_ray_kernels[tid])
758 : {
759 3001675 : rkb->subdomainSetup();
760 :
761 3001675 : const auto & mv_deps = rkb->getMooseVariableDependencies();
762 3001675 : needed_moose_vars.insert(mv_deps.begin(), mv_deps.end());
763 :
764 3001675 : const auto & mp_deps = rkb->getMatPropDependencies();
765 3001675 : needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
766 : }
767 :
768 : // Prepare aux vars
769 3022169 : for (auto & var : needed_moose_vars)
770 25082 : if (var->kind() == Moose::VarKindType::VAR_AUXILIARY)
771 5156 : var->prepareAux();
772 :
773 2997087 : _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, tid);
774 2997087 : _fe_problem.prepareMaterials(needed_mat_props, subdomain, tid);
775 2997087 : }
776 :
777 : void
778 22560690 : RayTracingStudy::reinitSegment(
779 : const Elem * elem, const Point & start, const Point & end, const Real length, THREAD_ID tid)
780 : {
781 : mooseAssert(MooseUtils::absoluteFuzzyEqual((start - end).norm(), length), "Invalid length");
782 : mooseAssert(currentlyPropagating(), "Should not call while not propagating");
783 :
784 22560690 : _fe_problem.setCurrentSubdomainID(elem, tid);
785 :
786 : // If we have any variables or material properties that are active, we definitely need to reinit
787 45044792 : bool reinit = _fe_problem.hasActiveElementalMooseVariables(tid) ||
788 22484102 : _fe_problem.hasActiveMaterialProperties(tid);
789 : // If not, make sure that the RayKernels have not requested a reinit (this could happen when a
790 : // RayKernel doesn't have variables or materials but still does an integration and needs qps)
791 : if (!reinit)
792 44967508 : for (const RayKernelBase * rk : currentRayKernels(tid))
793 22483886 : if (rk->needSegmentReinit())
794 : {
795 : reinit = true;
796 : break;
797 : }
798 :
799 22560690 : if (reinit)
800 : {
801 77068 : _fe_problem.prepare(elem, tid);
802 :
803 : std::vector<Point> points;
804 : std::vector<Real> weights;
805 77068 : buildSegmentQuadrature(start, end, length, points, weights);
806 77068 : _fe_problem.reinitElemPhys(elem, points, tid);
807 77068 : _fe_problem.assembly(tid, _sys.number()).modifyArbitraryWeights(weights);
808 :
809 77068 : _fe_problem.reinitMaterials(elem->subdomain_id(), tid);
810 77068 : }
811 22560690 : }
812 :
813 : void
814 77068 : RayTracingStudy::buildSegmentQuadrature(const Point & start,
815 : const Point & end,
816 : const Real length,
817 : std::vector<Point> & points,
818 : std::vector<Real> & weights) const
819 : {
820 77068 : points.resize(_segment_qrule->n_points());
821 77068 : weights.resize(_segment_qrule->n_points());
822 :
823 : const Point diff = end - start;
824 : const Point sum = end + start;
825 : mooseAssert(MooseUtils::absoluteFuzzyEqual(length, diff.norm()), "Invalid length");
826 :
827 : // The standard quadrature rule should be on x = [-1, 1]
828 : // To scale the points, you...
829 : // - Scale to size of the segment in 3D
830 : // initial_scaled_qp = x_qp * 0.5 * (end - start) = 0.5 * x_qp * diff
831 : // - Shift quadrature midpoint to segment midpoint
832 : // final_qp = initial_scaled_qp + 0.5 * (end - start) = initial_scaled_qp + 0.5 * sum
833 : // = 0.5 * (x_qp * diff + sum)
834 216188 : for (unsigned int qp = 0; qp < _segment_qrule->n_points(); ++qp)
835 : {
836 139120 : points[qp] = 0.5 * (_segment_qrule->qp(qp)(0) * diff + sum);
837 139120 : weights[qp] = 0.5 * _segment_qrule->w(qp) * length;
838 : }
839 77068 : }
840 :
841 : void
842 23617483 : RayTracingStudy::postOnSegment(const THREAD_ID tid, const std::shared_ptr<Ray> & /* ray */)
843 : {
844 : mooseAssert(currentlyPropagating(), "Should not call while not propagating");
845 23617483 : if (!_fe_problem.currentlyComputingJacobian() && !_fe_problem.currentlyComputingResidual())
846 : mooseAssert(_num_cached[tid] == 0,
847 : "Values should only be cached when computing Jacobian/residual");
848 :
849 : // Fill into cached Jacobian/residuals if necessary
850 23617483 : if (_fe_problem.currentlyComputingJacobian())
851 : {
852 6096 : _fe_problem.cacheJacobian(tid);
853 :
854 6096 : if (++_num_cached[tid] == 20)
855 : {
856 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
857 250 : _fe_problem.addCachedJacobian(tid);
858 250 : _num_cached[tid] = 0;
859 : }
860 : }
861 23611387 : else if (_fe_problem.currentlyComputingResidual())
862 : {
863 51968 : _fe_problem.cacheResidual(tid);
864 :
865 51968 : if (++_num_cached[tid] == 20)
866 : {
867 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
868 1445 : _fe_problem.addCachedResidual(tid);
869 1445 : _num_cached[tid] = 0;
870 : }
871 : }
872 23617483 : }
873 :
874 : void
875 5360 : RayTracingStudy::executeStudy()
876 : {
877 10720 : TIME_SECTION("executeStudy", 2, "Executing Study");
878 :
879 : mooseAssert(_called_initial_setup, "Initial setup not called");
880 :
881 : // Reset ray start/complete timers
882 5360 : _total_intersections = 0;
883 5360 : _max_intersections = 0;
884 5360 : _max_trajectory_changes = 0;
885 :
886 : // Reset physical tracing stats
887 80400 : for (auto & val : _local_trace_ray_results)
888 75040 : val = 0;
889 :
890 : // Reset crossing and intersection
891 5360 : _ending_processor_crossings = 0;
892 5360 : _total_processor_crossings = 0;
893 5360 : _ending_max_processor_crossings = 0;
894 5360 : _ending_intersections = 0;
895 5360 : _ending_max_intersections = 0;
896 5360 : _ending_max_trajectory_changes = 0;
897 5360 : _ending_distance = 0;
898 5360 : _total_distance = 0;
899 :
900 : // Zero the AuxVariables that our AuxRayKernels contribute to before they accumulate
901 5360 : zeroAuxVariables();
902 :
903 5360 : preExecuteStudy();
904 14251 : for (auto & rto : getRayTracingObjects())
905 14251 : rto->preExecuteStudy();
906 :
907 5360 : _ray_bank.clear();
908 :
909 11561 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
910 : {
911 6201 : _threaded_trace_ray[tid]->preExecute();
912 : _threaded_cached_normals[tid].clear();
913 : }
914 :
915 5360 : _communicator.barrier();
916 5360 : _execution_start_time = std::chrono::steady_clock::now();
917 :
918 5360 : _parallel_ray_study->preExecute();
919 :
920 : {
921 : {
922 5360 : auto generation_start_time = std::chrono::steady_clock::now();
923 :
924 10720 : TIME_SECTION("generateRays", 2, "Generating Rays");
925 :
926 5360 : generateRays();
927 :
928 5298 : _generation_time = std::chrono::steady_clock::now() - generation_start_time;
929 5298 : }
930 :
931 : // At this point, nobody is working so this is good time to make sure
932 : // Rays are unique across all processors in the working buffer
933 5298 : if (verifyRays())
934 : {
935 5298 : verifyUniqueRays(_parallel_ray_study->workBuffer().begin(),
936 : _parallel_ray_study->workBuffer().end(),
937 : /* error_suffix = */ "after generateRays()");
938 :
939 10590 : verifyUniqueRayIDs(_parallel_ray_study->workBuffer().begin(),
940 : _parallel_ray_study->workBuffer().end(),
941 : /* global = */ true,
942 : /* error_suffix = */ "after generateRays()");
943 : }
944 :
945 5294 : registeredRaySetup();
946 :
947 : {
948 10580 : TIME_SECTION("propagateRays", 2, "Propagating Rays");
949 :
950 5290 : const auto propagation_start_time = std::chrono::steady_clock::now();
951 :
952 5290 : _parallel_ray_study->execute();
953 :
954 5264 : _propagation_time = std::chrono::steady_clock::now() - propagation_start_time;
955 5264 : }
956 : }
957 :
958 5264 : _execution_time = std::chrono::steady_clock::now() - _execution_start_time;
959 :
960 5264 : if (verifyRays())
961 : {
962 10528 : verifyUniqueRays(_parallel_ray_study->workBuffer().begin(),
963 : _parallel_ray_study->workBuffer().end(),
964 : /* error_suffix = */ "after tracing completed");
965 :
966 : #ifndef NDEBUG
967 : // Outside of debug, _ray_bank always holds all of the Rays that have ended on this processor
968 : // We can use this as a global point to check for unique IDs for every Ray that has traced
969 : verifyUniqueRayIDs(_ray_bank.begin(),
970 : _ray_bank.end(),
971 : /* global = */ true,
972 : /* error_suffix = */ "after tracing completed");
973 : #endif
974 : }
975 :
976 : // Update counters from the threaded trace objects
977 11321 : for (const auto & tr : _threaded_trace_ray)
978 90855 : for (std::size_t i = 0; i < _local_trace_ray_results.size(); ++i)
979 84798 : _local_trace_ray_results[i] += tr->results()[i];
980 :
981 : // Update local ending counters
982 5264 : _total_processor_crossings = _ending_processor_crossings;
983 5264 : _max_processor_crossings = _ending_max_processor_crossings;
984 5264 : _total_intersections = _ending_intersections;
985 5264 : _max_intersections = _ending_max_intersections;
986 5264 : _max_trajectory_changes = _ending_max_trajectory_changes;
987 5264 : _total_distance = _ending_distance;
988 : // ...and communicate the global values
989 5264 : _communicator.sum(_total_processor_crossings);
990 5264 : _communicator.max(_max_processor_crossings);
991 5264 : _communicator.sum(_total_intersections);
992 5264 : _communicator.max(_max_intersections);
993 5264 : _communicator.max(_max_trajectory_changes);
994 5264 : _communicator.sum(_total_distance);
995 :
996 : // Throw a warning with the number of failed (tolerated) traces
997 5264 : if (_tolerate_failure)
998 : {
999 6 : auto failures = _local_trace_ray_results[TraceRay::FAILED_TRACES];
1000 6 : _communicator.sum(failures);
1001 6 : if (failures)
1002 6 : mooseWarning(
1003 : type(), " '", name(), "': ", failures, " ray tracing failures were tolerated.\n");
1004 : }
1005 :
1006 : // Clear the current RayKernels
1007 11321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1008 6057 : _threaded_current_ray_kernels[tid].clear();
1009 :
1010 : // Move the threaded cache trace information into the full cached trace vector
1011 : // Here, we only clear the cached vectors so that we might not have to
1012 : // reallocate on future traces
1013 : std::size_t num_entries = 0;
1014 11321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1015 6057 : num_entries += _threaded_cached_traces[tid].size();
1016 5264 : _cached_traces.clear();
1017 5264 : _cached_traces.reserve(num_entries);
1018 11321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1019 : {
1020 14633 : for (const auto & entry : _threaded_cached_traces[tid])
1021 8576 : _cached_traces.emplace_back(std::move(entry));
1022 6057 : _threaded_cached_traces[tid].clear();
1023 : }
1024 :
1025 : // Add any stragglers that contribute to the Jacobian or residual
1026 11321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1027 6057 : if (_num_cached[tid] != 0)
1028 : {
1029 : mooseAssert(_fe_problem.currentlyComputingJacobian() ||
1030 : _fe_problem.currentlyComputingResidual(),
1031 : "Should not have cached values without Jacobian/residual computation");
1032 :
1033 3297 : if (_fe_problem.currentlyComputingJacobian())
1034 146 : _fe_problem.addCachedJacobian(tid);
1035 : else
1036 3151 : _fe_problem.addCachedResidual(tid);
1037 :
1038 3297 : _num_cached[tid] = 0;
1039 : }
1040 :
1041 : // AuxRayKernels may have modified AuxVariables
1042 5264 : _fe_problem.getAuxiliarySystem().solution().close();
1043 5264 : _fe_problem.getAuxiliarySystem().update();
1044 :
1045 : // Clear FE
1046 11321 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1047 : {
1048 6057 : _fe_problem.clearActiveElementalMooseVariables(tid);
1049 6057 : _fe_problem.clearActiveMaterialProperties(tid);
1050 : }
1051 :
1052 5264 : postExecuteStudy();
1053 14082 : for (auto & rto : getRayTracingObjects())
1054 14082 : rto->postExecuteStudy();
1055 5260 : }
1056 :
1057 : void
1058 2755265 : RayTracingStudy::onCompleteRay(const std::shared_ptr<Ray> & ray)
1059 : {
1060 : mooseAssert(currentlyPropagating(), "Should only be called during Ray propagation");
1061 :
1062 2755265 : _ending_processor_crossings += ray->processorCrossings();
1063 2755265 : _ending_max_processor_crossings =
1064 2758098 : std::max(_ending_max_processor_crossings, ray->processorCrossings());
1065 2755265 : _ending_intersections += ray->intersections();
1066 2766258 : _ending_max_intersections = std::max(_ending_max_intersections, ray->intersections());
1067 2755265 : _ending_max_trajectory_changes =
1068 2755562 : std::max(_ending_max_trajectory_changes, ray->trajectoryChanges());
1069 2755265 : _ending_distance += ray->distance();
1070 :
1071 : #ifdef NDEBUG
1072 : // In non-opt modes, we will always bank the Rays for debugging
1073 2755265 : if (_bank_rays_on_completion)
1074 : #endif
1075 2755263 : _ray_bank.emplace_back(ray);
1076 2755265 : }
1077 :
1078 : RayDataIndex
1079 1289 : RayTracingStudy::registerRayDataInternal(const std::string & name, const bool aux)
1080 : {
1081 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1082 :
1083 1289 : if (_called_initial_setup)
1084 4 : mooseError("Cannot register Ray ", (aux ? "aux " : ""), "data after initialSetup()");
1085 :
1086 1287 : auto & map = aux ? _ray_aux_data_map : _ray_data_map;
1087 : const auto find = map.find(name);
1088 1287 : if (find != map.end())
1089 27 : return find->second;
1090 :
1091 1260 : auto & other_map = aux ? _ray_data_map : _ray_aux_data_map;
1092 1260 : if (other_map.find(name) != other_map.end())
1093 4 : mooseError("Cannot register Ray aux data with name ",
1094 : name,
1095 : " because Ray ",
1096 2 : (aux ? "(non-aux)" : "aux"),
1097 : " data already exists with said name.");
1098 :
1099 : // Add into the name -> index map
1100 1258 : map.emplace(name, map.size());
1101 :
1102 : // Add into the index -> names vector
1103 1258 : auto & vector = aux ? _ray_aux_data_names : _ray_data_names;
1104 1258 : vector.push_back(name);
1105 :
1106 1258 : return map.size() - 1;
1107 : }
1108 :
1109 : std::vector<RayDataIndex>
1110 144 : RayTracingStudy::registerRayDataInternal(const std::vector<std::string> & names, const bool aux)
1111 : {
1112 144 : std::vector<RayDataIndex> indices(names.size());
1113 315 : for (std::size_t i = 0; i < names.size(); ++i)
1114 171 : indices[i] = registerRayDataInternal(names[i], aux);
1115 144 : return indices;
1116 0 : }
1117 :
1118 : RayDataIndex
1119 788 : RayTracingStudy::getRayDataIndexInternal(const std::string & name,
1120 : const bool aux,
1121 : const bool graceful) const
1122 : {
1123 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1124 :
1125 788 : const auto & map = aux ? _ray_aux_data_map : _ray_data_map;
1126 : const auto find = map.find(name);
1127 788 : if (find != map.end())
1128 774 : return find->second;
1129 :
1130 14 : if (graceful)
1131 : return Ray::INVALID_RAY_DATA_INDEX;
1132 :
1133 6 : const auto & other_map = aux ? _ray_data_map : _ray_aux_data_map;
1134 6 : if (other_map.find(name) != other_map.end())
1135 8 : mooseError("Ray data with name '",
1136 : name,
1137 : "' was not found.\n\n",
1138 : "However, Ray ",
1139 4 : (aux ? "non-aux" : "aux"),
1140 : " data with said name was found.\n",
1141 : "Did you mean to use ",
1142 6 : (aux ? "getRayDataIndex()/getRayDataIndices()?"
1143 : : "getRayAuxDataIndex()/getRayAuxDataIndices()"),
1144 : "?");
1145 :
1146 4 : mooseError("Unknown Ray ", (aux ? "aux " : ""), "data with name ", name);
1147 : }
1148 :
1149 : std::vector<RayDataIndex>
1150 326 : RayTracingStudy::getRayDataIndicesInternal(const std::vector<std::string> & names,
1151 : const bool aux,
1152 : const bool graceful) const
1153 : {
1154 326 : std::vector<RayDataIndex> indices(names.size());
1155 386 : for (std::size_t i = 0; i < names.size(); ++i)
1156 60 : indices[i] = getRayDataIndexInternal(names[i], aux, graceful);
1157 326 : return indices;
1158 0 : }
1159 :
1160 : const std::string &
1161 4898 : RayTracingStudy::getRayDataNameInternal(const RayDataIndex index, const bool aux) const
1162 : {
1163 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1164 :
1165 9796 : if ((aux ? rayAuxDataSize() : rayDataSize()) < index)
1166 4 : mooseError("Unknown Ray ", aux ? "aux " : "", "data with index ", index);
1167 9792 : return aux ? _ray_aux_data_names[index] : _ray_data_names[index];
1168 : }
1169 :
1170 : RayDataIndex
1171 494 : RayTracingStudy::registerRayData(const std::string & name)
1172 : {
1173 494 : return registerRayDataInternal(name, /* aux = */ false);
1174 : }
1175 :
1176 : std::vector<RayDataIndex>
1177 76 : RayTracingStudy::registerRayData(const std::vector<std::string> & names)
1178 : {
1179 76 : return registerRayDataInternal(names, /* aux = */ false);
1180 : }
1181 :
1182 : RayDataIndex
1183 632 : RayTracingStudy::getRayDataIndex(const std::string & name, const bool graceful /* = false */) const
1184 : {
1185 632 : return getRayDataIndexInternal(name, /* aux = */ false, graceful);
1186 : }
1187 :
1188 : std::vector<RayDataIndex>
1189 163 : RayTracingStudy::getRayDataIndices(const std::vector<std::string> & names,
1190 : const bool graceful /* = false */) const
1191 : {
1192 163 : return getRayDataIndicesInternal(names, /* aux = */ false, graceful);
1193 : }
1194 :
1195 : const std::string &
1196 2450 : RayTracingStudy::getRayDataName(const RayDataIndex index) const
1197 : {
1198 2450 : return getRayDataNameInternal(index, /* aux = */ false);
1199 : }
1200 :
1201 : RayDataIndex
1202 624 : RayTracingStudy::registerRayAuxData(const std::string & name)
1203 : {
1204 624 : return registerRayDataInternal(name, /* aux = */ true);
1205 : }
1206 :
1207 : std::vector<RayDataIndex>
1208 68 : RayTracingStudy::registerRayAuxData(const std::vector<std::string> & names)
1209 : {
1210 68 : return registerRayDataInternal(names, /* aux = */ true);
1211 : }
1212 :
1213 : RayDataIndex
1214 96 : RayTracingStudy::getRayAuxDataIndex(const std::string & name,
1215 : const bool graceful /* = false */) const
1216 : {
1217 96 : return getRayDataIndexInternal(name, /* aux = */ true, graceful);
1218 : }
1219 :
1220 : std::vector<RayDataIndex>
1221 163 : RayTracingStudy::getRayAuxDataIndices(const std::vector<std::string> & names,
1222 : const bool graceful /* = false */) const
1223 : {
1224 163 : return getRayDataIndicesInternal(names, /* aux = */ true, graceful);
1225 : }
1226 :
1227 : const std::string &
1228 2448 : RayTracingStudy::getRayAuxDataName(const RayDataIndex index) const
1229 : {
1230 2448 : return getRayDataNameInternal(index, /* aux = */ true);
1231 : }
1232 :
1233 : bool
1234 6201 : RayTracingStudy::hasRayKernels(const THREAD_ID tid)
1235 : {
1236 : std::vector<RayKernelBase *> result;
1237 6201 : getRayKernels(result, tid);
1238 6201 : return result.size();
1239 6201 : }
1240 :
1241 : void
1242 2997089 : RayTracingStudy::getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid)
1243 : {
1244 : // If the cache doesn't have any attributes yet, it means that we haven't set
1245 : // the conditions yet. We do this so that it can be generated on the fly on first use.
1246 2997089 : if (!_threaded_cache_ray_kernel[tid].numAttribs())
1247 : {
1248 1396 : if (!_called_initial_setup)
1249 2 : mooseError("Should not call getRayKernels() before initialSetup()");
1250 :
1251 1394 : auto query = _fe_problem.theWarehouse()
1252 1394 : .query()
1253 2788 : .condition<AttribRayTracingStudy>(this)
1254 1394 : .condition<AttribSystem>("RayKernel")
1255 1394 : .condition<AttribThread>(tid);
1256 1394 : _threaded_cache_ray_kernel[tid] = query.clone();
1257 1394 : }
1258 :
1259 : _threaded_cache_ray_kernel[tid].queryInto(result, id);
1260 2997087 : }
1261 :
1262 : void
1263 2997087 : RayTracingStudy::getRayKernels(std::vector<RayKernelBase *> & result,
1264 : SubdomainID id,
1265 : THREAD_ID tid,
1266 : RayID ray_id)
1267 : {
1268 : // No Ray registration: no need to sift through objects
1269 2997087 : if (!_use_ray_registration)
1270 : {
1271 2990007 : getRayKernels(result, id, tid);
1272 : }
1273 : // Has Ray registration: only pick the objects associated with ray_id
1274 : else
1275 : {
1276 : // Get all of the kernels on this block
1277 : std::vector<RayKernelBase *> rkbs;
1278 7080 : getRayKernels(rkbs, id, tid);
1279 :
1280 : // The RayTracingObjects associated with this ray
1281 7080 : const auto & ray_id_rtos = _threaded_ray_object_registration[tid][ray_id];
1282 :
1283 : // The result is the union of all of the kernels and the objects associated with this Ray
1284 7080 : result.clear();
1285 17088 : for (auto rkb : rkbs)
1286 10008 : if (ray_id_rtos.count(rkb))
1287 7168 : result.push_back(rkb);
1288 7080 : }
1289 2997087 : }
1290 :
1291 : void
1292 913198 : RayTracingStudy::getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
1293 : BoundaryID id,
1294 : THREAD_ID tid)
1295 : {
1296 : // If the cache doesn't have any attributes yet, it means that we haven't set
1297 : // the conditions yet. We do this so that it can be generated on the fly on first use.
1298 913198 : if (!_threaded_cache_ray_bc[tid].numAttribs())
1299 : {
1300 1327 : if (!_called_initial_setup)
1301 2 : mooseError("Should not call getRayBCs() before initialSetup()");
1302 :
1303 1325 : auto query = _fe_problem.theWarehouse()
1304 1325 : .query()
1305 2650 : .condition<AttribRayTracingStudy>(this)
1306 1325 : .condition<AttribSystem>("RayBoundaryCondition")
1307 1325 : .condition<AttribThread>(tid);
1308 1325 : _threaded_cache_ray_bc[tid] = query.clone();
1309 1325 : }
1310 :
1311 913196 : _threaded_cache_ray_bc[tid].queryInto(result, std::make_tuple(id, false));
1312 913196 : }
1313 :
1314 : void
1315 973271 : RayTracingStudy::getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
1316 : const std::vector<TraceRayBndElement> & bnd_elems,
1317 : THREAD_ID tid,
1318 : RayID ray_id)
1319 : {
1320 : // No Ray registration: no need to sift through objects
1321 973271 : if (!_use_ray_registration)
1322 : {
1323 971604 : if (bnd_elems.size() == 1)
1324 908200 : getRayBCs(result, bnd_elems[0].bnd_id, tid);
1325 : else
1326 : {
1327 63404 : std::vector<BoundaryID> bnd_ids(bnd_elems.size());
1328 195468 : for (MooseIndex(bnd_elems.size()) i = 0; i < bnd_elems.size(); ++i)
1329 132064 : bnd_ids[i] = bnd_elems[i].bnd_id;
1330 63404 : getRayBCs(result, bnd_ids, tid);
1331 63404 : }
1332 : }
1333 : // Has Ray registration: only pick the objects associated with ray_id
1334 : else
1335 : {
1336 : // Get all of the RayBCs on these boundaries
1337 : std::vector<RayBoundaryConditionBase *> rbcs;
1338 1667 : if (bnd_elems.size() == 1)
1339 1426 : getRayBCs(rbcs, bnd_elems[0].bnd_id, tid);
1340 : else
1341 : {
1342 241 : std::vector<BoundaryID> bnd_ids(bnd_elems.size());
1343 831 : for (MooseIndex(bnd_elems.size()) i = 0; i < bnd_elems.size(); ++i)
1344 590 : bnd_ids[i] = bnd_elems[i].bnd_id;
1345 241 : getRayBCs(rbcs, bnd_ids, tid);
1346 241 : }
1347 :
1348 : // The RayTracingObjects associated with this ray
1349 : mooseAssert(ray_id < _threaded_ray_object_registration[tid].size(), "Not in registration");
1350 1667 : const auto & ray_id_rtos = _threaded_ray_object_registration[tid][ray_id];
1351 :
1352 : // The result is the union of all of the kernels and the objects associated with this Ray
1353 1667 : result.clear();
1354 4646 : for (auto rbc : rbcs)
1355 2979 : if (ray_id_rtos.count(rbc))
1356 1693 : result.push_back(rbc);
1357 1667 : }
1358 973271 : }
1359 :
1360 : std::vector<RayTracingObject *>
1361 23270 : RayTracingStudy::getRayTracingObjects()
1362 : {
1363 : std::vector<RayTracingObject *> result;
1364 46540 : _fe_problem.theWarehouse().query().condition<AttribRayTracingStudy>(this).queryInto(result);
1365 23270 : return result;
1366 0 : }
1367 :
1368 : const std::vector<std::shared_ptr<Ray>> &
1369 664 : RayTracingStudy::rayBank() const
1370 : {
1371 664 : if (!_bank_rays_on_completion)
1372 2 : mooseError("The Ray bank is not available because the private parameter "
1373 : "'_bank_rays_on_completion' is set to false.");
1374 662 : if (currentlyGenerating() || currentlyPropagating())
1375 2 : mooseError("Cannot get the Ray bank during generation or propagation.");
1376 :
1377 660 : return _ray_bank;
1378 : }
1379 :
1380 : std::shared_ptr<Ray>
1381 426 : RayTracingStudy::getBankedRay(const RayID ray_id) const
1382 : {
1383 : // This is only a linear search - can be improved on with a map in the future
1384 : // if this is used on a larger scale
1385 426 : std::shared_ptr<Ray> ray;
1386 778 : for (const std::shared_ptr<Ray> & possible_ray : rayBank())
1387 670 : if (possible_ray->id() == ray_id)
1388 : {
1389 : ray = possible_ray;
1390 : break;
1391 : }
1392 :
1393 : // Make sure one and only one processor has the Ray
1394 426 : unsigned int have_ray = ray ? 1 : 0;
1395 426 : _communicator.sum(have_ray);
1396 426 : if (have_ray == 0)
1397 2 : mooseError("Could not find a Ray with the ID ", ray_id, " in the Ray banks.");
1398 :
1399 : // This should never happen... but let's make sure
1400 : mooseAssert(have_ray == 1, "Multiple rays with the same ID were found in the Ray banks");
1401 :
1402 424 : return ray;
1403 : }
1404 :
1405 : RayData
1406 426 : RayTracingStudy::getBankedRayDataInternal(const RayID ray_id,
1407 : const RayDataIndex index,
1408 : const bool aux) const
1409 : {
1410 : // Will be a nullptr shared_ptr if this processor doesn't own the Ray
1411 426 : const std::shared_ptr<Ray> ray = getBankedRay(ray_id);
1412 :
1413 424 : Real value = ray ? (aux ? ray->auxData(index) : ray->data(index)) : 0;
1414 424 : _communicator.sum(value);
1415 742 : return value;
1416 : }
1417 :
1418 : RayData
1419 394 : RayTracingStudy::getBankedRayData(const RayID ray_id, const RayDataIndex index) const
1420 : {
1421 394 : return getBankedRayDataInternal(ray_id, index, /* aux = */ false);
1422 : }
1423 :
1424 : RayData
1425 32 : RayTracingStudy::getBankedRayAuxData(const RayID ray_id, const RayDataIndex index) const
1426 : {
1427 32 : return getBankedRayDataInternal(ray_id, index, /* aux = */ true);
1428 : }
1429 :
1430 : RayID
1431 1602 : RayTracingStudy::registerRay(const std::string & name)
1432 : {
1433 : libmesh_parallel_only(comm());
1434 :
1435 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1436 :
1437 1602 : if (!_use_ray_registration)
1438 2 : mooseError("Cannot use registerRay() with Ray registration disabled");
1439 :
1440 : // This is parallel only for now. We could likely stagger the ID building like we do with
1441 : // the unique IDs, but it would require a sync point which isn't there right now
1442 : libmesh_parallel_only(comm());
1443 :
1444 1600 : const auto & it = _registered_ray_map.find(name);
1445 1600 : if (it != _registered_ray_map.end())
1446 0 : return it->second;
1447 :
1448 1600 : const auto id = _reverse_registered_ray_map.size();
1449 1600 : _registered_ray_map.emplace(name, id);
1450 1600 : _reverse_registered_ray_map.push_back(name);
1451 1600 : return id;
1452 : }
1453 :
1454 : RayID
1455 10801 : RayTracingStudy::registeredRayID(const std::string & name, const bool graceful /* = false */) const
1456 : {
1457 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1458 :
1459 10801 : if (!_use_ray_registration)
1460 2 : mooseError("Should not use registeredRayID() with Ray registration disabled");
1461 :
1462 10799 : const auto search = _registered_ray_map.find(name);
1463 10799 : if (search != _registered_ray_map.end())
1464 10791 : return search->second;
1465 :
1466 8 : if (graceful)
1467 : return Ray::INVALID_RAY_ID;
1468 :
1469 2 : mooseError("Attempted to obtain ID of registered Ray ",
1470 : name,
1471 : ", but a Ray with said name is not registered.");
1472 : }
1473 :
1474 : const std::string &
1475 40 : RayTracingStudy::registeredRayName(const RayID ray_id) const
1476 : {
1477 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1478 :
1479 40 : if (!_use_ray_registration)
1480 2 : mooseError("Should not use registeredRayName() with Ray registration disabled");
1481 :
1482 38 : if (_reverse_registered_ray_map.size() > ray_id)
1483 36 : return _reverse_registered_ray_map[ray_id];
1484 :
1485 2 : mooseError("Attempted to obtain name of registered Ray with ID ",
1486 : ray_id,
1487 : ", but a Ray with said ID is not registered.");
1488 : }
1489 :
1490 : Real
1491 2029 : RayTracingStudy::computeTotalVolume()
1492 : {
1493 2029 : Real volume = 0;
1494 271139 : for (const auto & elem : *_mesh.getActiveLocalElementRange())
1495 269110 : volume += elem->volume();
1496 2029 : _communicator.sum(volume);
1497 2029 : return volume;
1498 : }
1499 :
1500 : const std::vector<std::vector<BoundaryID>> &
1501 10274 : RayTracingStudy::getInternalSidesets(const Elem * elem) const
1502 : {
1503 : mooseAssert(_use_internal_sidesets, "Not using internal sidesets");
1504 : mooseAssert(hasInternalSidesets(), "Processor does not have internal sidesets");
1505 : mooseAssert(_internal_sidesets_map.size() > _elem_index_helper.getIndex(elem),
1506 : "Internal sideset map not initialized");
1507 :
1508 : const auto index = _elem_index_helper.getIndex(elem);
1509 10274 : return _internal_sidesets_map[index];
1510 : }
1511 :
1512 : TraceData &
1513 8576 : RayTracingStudy::initThreadedCachedTrace(const std::shared_ptr<Ray> & ray, THREAD_ID tid)
1514 : {
1515 : mooseAssert(shouldCacheTrace(ray), "Not caching trace");
1516 : mooseAssert(currentlyPropagating(), "Should only use while tracing");
1517 :
1518 8576 : _threaded_cached_traces[tid].emplace_back(ray);
1519 8576 : return _threaded_cached_traces[tid].back();
1520 : }
1521 :
1522 : void
1523 5732 : RayTracingStudy::verifyUniqueRayIDs(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
1524 : const std::vector<std::shared_ptr<Ray>>::const_iterator end,
1525 : const bool global,
1526 : const std::string & error_suffix) const
1527 : {
1528 : // Determine the unique set of Ray IDs on this processor,
1529 : // and if not locally unique throw an error. Once we build this set,
1530 : // we will send it to rank 0 to verify globally
1531 : std::set<RayID> local_rays;
1532 2708261 : for (const std::shared_ptr<Ray> & ray : as_range(begin, end))
1533 : {
1534 : mooseAssert(ray, "Null ray");
1535 :
1536 : // Try to insert into the set; the second entry in the pair
1537 : // will be false if it was not inserted
1538 2702533 : if (!local_rays.insert(ray->id()).second)
1539 : {
1540 4 : for (const std::shared_ptr<Ray> & other_ray : as_range(begin, end))
1541 4 : if (ray.get() != other_ray.get() && ray->id() == other_ray->id())
1542 8 : mooseError("Multiple Rays exist with ID ",
1543 4 : ray->id(),
1544 : " on processor ",
1545 4 : _pid,
1546 : " ",
1547 : error_suffix,
1548 : "\n\nOffending Ray information:\n\n",
1549 4 : ray->getInfo(),
1550 : "\n",
1551 4 : other_ray->getInfo());
1552 : }
1553 : }
1554 :
1555 : // Send IDs from all procs to rank 0 and verify on rank 0
1556 5728 : if (global)
1557 : {
1558 : // Package our local IDs and send to rank 0
1559 : std::map<processor_id_type, std::vector<RayID>> send_ids;
1560 5294 : if (local_rays.size())
1561 : send_ids.emplace(std::piecewise_construct,
1562 4149 : std::forward_as_tuple(0),
1563 4149 : std::forward_as_tuple(local_rays.begin(), local_rays.end()));
1564 : local_rays.clear();
1565 :
1566 : // Mapping on rank 0 from ID -> processor ID
1567 : std::map<RayID, processor_id_type> global_map;
1568 :
1569 : // Verify another processor's IDs against the global map on rank 0
1570 : const auto check_ids =
1571 4149 : [this, &global_map, &error_suffix](processor_id_type pid, const std::vector<RayID> & ids)
1572 : {
1573 2704482 : for (const RayID id : ids)
1574 : {
1575 2700333 : const auto emplace_pair = global_map.emplace(id, pid);
1576 :
1577 : // Means that this ID already exists in the map
1578 2700333 : if (!emplace_pair.second)
1579 0 : mooseError("Ray with ID ",
1580 : id,
1581 : " exists on ranks ",
1582 0 : emplace_pair.first->second,
1583 : " and ",
1584 : pid,
1585 : "\n",
1586 : error_suffix);
1587 : }
1588 9443 : };
1589 :
1590 5294 : Parallel::push_parallel_vector_data(_communicator, send_ids, check_ids);
1591 : }
1592 5728 : }
1593 :
1594 : void
1595 10562 : RayTracingStudy::verifyUniqueRays(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
1596 : const std::vector<std::shared_ptr<Ray>>::const_iterator end,
1597 : const std::string & error_suffix)
1598 : {
1599 : std::set<const Ray *> rays;
1600 2710901 : for (const std::shared_ptr<Ray> & ray : as_range(begin, end))
1601 2700341 : if (!rays.insert(ray.get()).second) // false if not inserted into rays
1602 2 : mooseError("Multiple shared_ptrs were found that point to the same Ray ",
1603 : error_suffix,
1604 : "\n\nOffending Ray:\n",
1605 2 : ray->getInfo());
1606 10560 : }
1607 :
1608 : void
1609 2699651 : RayTracingStudy::moveRayToBuffer(std::shared_ptr<Ray> & ray)
1610 : {
1611 : mooseAssert(currentlyGenerating(), "Can only use while generating");
1612 : mooseAssert(ray, "Null ray");
1613 : mooseAssert(ray->shouldContinue(), "Ray is not continuing");
1614 :
1615 2699651 : _parallel_ray_study->moveWorkToBuffer(ray, /* tid = */ 0);
1616 2699651 : }
1617 :
1618 : void
1619 180 : RayTracingStudy::moveRaysToBuffer(std::vector<std::shared_ptr<Ray>> & rays)
1620 : {
1621 : mooseAssert(currentlyGenerating(), "Can only use while generating");
1622 : #ifndef NDEBUG
1623 : for (const std::shared_ptr<Ray> & ray : rays)
1624 : {
1625 : mooseAssert(ray, "Null ray");
1626 : mooseAssert(ray->shouldContinue(), "Ray is not continuing");
1627 : }
1628 : #endif
1629 :
1630 180 : _parallel_ray_study->moveWorkToBuffer(rays, /* tid = */ 0);
1631 180 : }
1632 :
1633 : void
1634 54960 : RayTracingStudy::moveRayToBufferDuringTrace(std::shared_ptr<Ray> & ray,
1635 : const THREAD_ID tid,
1636 : const AcquireMoveDuringTraceKey &)
1637 : {
1638 : mooseAssert(ray, "Null ray");
1639 : mooseAssert(currentlyPropagating(), "Can only use while tracing");
1640 :
1641 54960 : _parallel_ray_study->moveWorkToBuffer(ray, tid);
1642 54960 : }
1643 :
1644 : void
1645 5042 : RayTracingStudy::reserveRayBuffer(const std::size_t size)
1646 : {
1647 5042 : if (!currentlyGenerating())
1648 2 : mooseError("Can only reserve in Ray buffer during generateRays()");
1649 :
1650 5040 : _parallel_ray_study->reserveBuffer(size);
1651 5040 : }
1652 :
1653 : const Point &
1654 4429805 : RayTracingStudy::getSideNormal(const Elem * elem, unsigned short side, const THREAD_ID tid)
1655 : {
1656 : std::unordered_map<std::pair<const Elem *, unsigned short>, Point> & cache =
1657 4429805 : _threaded_cached_normals[tid];
1658 :
1659 : // See if we've already cached this side normal
1660 4429805 : const auto elem_side_pair = std::make_pair(elem, side);
1661 : const auto search = cache.find(elem_side_pair);
1662 :
1663 : // Haven't cached this side normal: compute it and then cache it
1664 4429805 : if (search == cache.end())
1665 : {
1666 936906 : _threaded_fe_face[tid]->reinit(elem, side);
1667 : const auto & normal = _threaded_fe_face[tid]->get_normals()[0];
1668 : cache.emplace(elem_side_pair, normal);
1669 936906 : return normal;
1670 : }
1671 :
1672 : // Have cached this side normal: simply return it
1673 3492899 : return search->second;
1674 : }
1675 :
1676 : bool
1677 2161 : RayTracingStudy::sameLevelActiveElems() const
1678 : {
1679 2161 : unsigned int min_level = std::numeric_limits<unsigned int>::max();
1680 2161 : unsigned int max_level = std::numeric_limits<unsigned int>::min();
1681 :
1682 276647 : for (const auto & elem : *_mesh.getActiveLocalElementRange())
1683 : {
1684 274486 : const auto level = elem->level();
1685 274486 : min_level = std::min(level, min_level);
1686 274486 : max_level = std::max(level, max_level);
1687 : }
1688 :
1689 2161 : _communicator.min(min_level);
1690 2161 : _communicator.max(max_level);
1691 :
1692 2161 : return min_level == max_level;
1693 : }
1694 :
1695 : Real
1696 3544866 : RayTracingStudy::subdomainHmax(const SubdomainID subdomain_id) const
1697 : {
1698 : const auto find = _subdomain_hmax.find(subdomain_id);
1699 3544866 : if (find == _subdomain_hmax.end())
1700 2 : mooseError("Subdomain ", subdomain_id, " not found in subdomain hmax map");
1701 3544864 : return find->second;
1702 : }
1703 :
1704 : bool
1705 7198 : RayTracingStudy::isRectangularDomain() const
1706 : {
1707 : Real bbox_volume = 1;
1708 22771 : for (unsigned int d = 0; d < _mesh.dimension(); ++d)
1709 15573 : bbox_volume *= std::abs(_b_box.max()(d) - _b_box.min()(d));
1710 :
1711 7198 : return MooseUtils::absoluteFuzzyEqual(bbox_volume, totalVolume(), TOLERANCE);
1712 : }
1713 :
1714 : void
1715 2025 : RayTracingStudy::resetUniqueRayIDs()
1716 : {
1717 : libmesh_parallel_only(comm());
1718 :
1719 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1720 :
1721 : mooseAssert(!currentlyGenerating() && !currentlyPropagating(),
1722 : "Cannot be reset during generation or propagation");
1723 :
1724 4385 : for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1725 2360 : _threaded_next_ray_id[tid] = (RayID)_pid * (RayID)libMesh::n_threads() + (RayID)tid;
1726 2025 : }
1727 :
1728 : RayID
1729 2709316 : RayTracingStudy::generateUniqueRayID(const THREAD_ID tid)
1730 : {
1731 : // Get the current ID to return
1732 2709316 : const auto id = _threaded_next_ray_id[tid];
1733 :
1734 : // Advance so that the next call has the correct ID
1735 2709316 : _threaded_next_ray_id[tid] += (RayID)n_processors() * (RayID)libMesh::n_threads();
1736 :
1737 2709316 : return id;
1738 : }
1739 :
1740 : void
1741 2025 : RayTracingStudy::resetReplicatedRayIDs()
1742 : {
1743 : libmesh_parallel_only(comm());
1744 :
1745 : Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1746 :
1747 : mooseAssert(!currentlyGenerating() && !currentlyPropagating(),
1748 : "Cannot be reset during generation or propagation");
1749 :
1750 2025 : _replicated_next_ray_id = 0;
1751 2025 : }
1752 :
1753 : RayID
1754 632 : RayTracingStudy::generateReplicatedRayID()
1755 : {
1756 632 : return _replicated_next_ray_id++;
1757 : }
1758 :
1759 : bool
1760 2053515 : RayTracingStudy::sideIsIncoming(const Elem * const elem,
1761 : const unsigned short side,
1762 : const Point & direction,
1763 : const THREAD_ID tid)
1764 : {
1765 2053515 : const auto & normal = getSideNormal(elem, side, tid);
1766 : const auto dot = normal * direction;
1767 2053515 : return dot < TraceRayTools::TRACE_TOLERANCE;
1768 : }
1769 :
1770 : std::shared_ptr<Ray>
1771 2643796 : RayTracingStudy::acquireRay()
1772 : {
1773 : mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1774 :
1775 2643796 : return _parallel_ray_study->acquireParallelData(
1776 : /* tid = */ 0,
1777 5287592 : this,
1778 2643796 : generateUniqueRayID(/* tid = */ 0),
1779 5287592 : rayDataSize(),
1780 2643796 : rayAuxDataSize(),
1781 2643796 : /* reset = */ true,
1782 2643796 : Ray::ConstructRayKey());
1783 : }
1784 :
1785 : std::shared_ptr<Ray>
1786 10560 : RayTracingStudy::acquireUnsizedRay()
1787 : {
1788 : mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1789 :
1790 10560 : return _parallel_ray_study->acquireParallelData(/* tid = */ 0,
1791 21120 : this,
1792 10560 : generateUniqueRayID(/* tid = */ 0),
1793 21120 : /* data_size = */ 0,
1794 21120 : /* aux_data_size = */ 0,
1795 10560 : /* reset = */ true,
1796 10560 : Ray::ConstructRayKey());
1797 : }
1798 :
1799 : std::shared_ptr<Ray>
1800 632 : RayTracingStudy::acquireReplicatedRay()
1801 : {
1802 : mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1803 : libmesh_parallel_only(comm());
1804 :
1805 632 : return _parallel_ray_study->acquireParallelData(
1806 : /* tid = */ 0,
1807 1264 : this,
1808 632 : generateReplicatedRayID(),
1809 1264 : rayDataSize(),
1810 632 : rayAuxDataSize(),
1811 632 : /* reset = */ true,
1812 632 : Ray::ConstructRayKey());
1813 : }
1814 :
1815 : std::shared_ptr<Ray>
1816 1602 : RayTracingStudy::acquireRegisteredRay(const std::string & name)
1817 : {
1818 : mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1819 :
1820 : // Either register a Ray or get an already registered Ray id
1821 1602 : const RayID id = registerRay(name);
1822 :
1823 : // Acquire a Ray with the properly sized data initialized to zero
1824 1600 : return _parallel_ray_study->acquireParallelData(
1825 : /* tid = */ 0,
1826 3200 : this,
1827 : id,
1828 3200 : rayDataSize(),
1829 1600 : rayAuxDataSize(),
1830 1600 : /* reset = */ true,
1831 1600 : Ray::ConstructRayKey());
1832 : }
1833 :
1834 : std::shared_ptr<Ray>
1835 2698423 : RayTracingStudy::acquireCopiedRay(const Ray & ray)
1836 : {
1837 : mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1838 2698423 : return _parallel_ray_study->acquireParallelData(
1839 5396846 : /* tid = */ 0, &ray, Ray::ConstructRayKey());
1840 : }
1841 :
1842 : std::shared_ptr<Ray>
1843 54960 : RayTracingStudy::acquireRayDuringTrace(const THREAD_ID tid, const AcquireMoveDuringTraceKey &)
1844 : {
1845 : mooseAssert(currentlyPropagating(), "Can only use during propagation");
1846 54960 : return _parallel_ray_study->acquireParallelData(tid,
1847 109920 : this,
1848 54960 : generateUniqueRayID(tid),
1849 109920 : rayDataSize(),
1850 54960 : rayAuxDataSize(),
1851 54960 : /* reset = */ true,
1852 54960 : Ray::ConstructRayKey());
1853 : }
|
