https://mooseframework.inl.gov
RayTracingStudy.C
Go to the documentation of this file.
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 // Local Includes
13 #include "AuxRayKernel.h"
15 #include "RayKernel.h"
16 #include "TraceRay.h"
17 #include "TraceRayTools.h"
18 
19 // MOOSE Includes
20 #include "AuxiliarySystem.h"
21 #include "Assembly.h"
22 #include "NonlinearSystemBase.h"
23 
24 // libMesh Includes
25 #include "libmesh/enum_to_string.h"
26 #include "libmesh/mesh_tools.h"
27 #include "libmesh/parallel_sync.h"
28 #include "libmesh/remote_elem.h"
29 
30 using namespace libMesh;
31 
34 {
35  auto params = GeneralUserObject::validParams();
36 
37  // Parameters for the execution of the rays
39 
40  params.addRangeCheckedParam<Real>("ray_distance",
41  std::numeric_limits<Real>::max(),
42  "ray_distance > 0",
43  "The maximum distance all Rays can travel");
44 
45  params.addParam<bool>(
46  "tolerate_failure", false, "Whether or not to tolerate a ray tracing failure");
47 
48  MooseEnum work_buffers("lifo circular", "circular");
49  params.addParam<MooseEnum>("work_buffer_type", work_buffers, "The work buffer type to use");
50 
51  params.addParam<bool>(
52  "ray_kernel_coverage_check", true, "Whether or not to perform coverage checks on RayKernels");
53  params.addParam<bool>("warn_non_planar",
54  true,
55  "Whether or not to produce a warning if any element faces are non-planar.");
56 
57  params.addParam<bool>(
58  "always_cache_traces",
59  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  params.addParam<bool>("data_on_cache_traces",
63  false,
64  "Whether or not to also cache the Ray's data when caching its traces");
65  params.addParam<bool>("aux_data_on_cache_traces",
66  false,
67  "Whether or not to also cache the Ray's aux data when caching its traces");
68  params.addParam<bool>(
69  "segments_on_cache_traces",
70  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  params.addParam<bool>("use_internal_sidesets",
77  false,
78  "Whether or not to use internal sidesets for RayBCs in ray tracing");
79 
80  params.addParam<bool>("warn_subdomain_hmax",
81  true,
82  "Whether or not to warn if the approximated hmax (constant on subdomain) "
83  "varies significantly for an element");
84 
85  params.addParam<bool>(
86  "verify_rays",
87  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  params.addParam<bool>("verify_trace_intersections",
92  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  params.addParam<bool>("allow_other_flags_with_prekernels",
98  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  ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
105 
106  params.addParamNamesToGroup(
107  "always_cache_traces data_on_cache_traces aux_data_on_cache_traces segments_on_cache_traces",
108  "Trace cache");
109  params.addParamNamesToGroup("warn_non_planar warn_subdomain_hmax", "Tracing Warnings");
110  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  params.addPrivateParam<bool>("_use_ray_registration", true);
115  // Whether or not to bank Rays on completion
116  params.addPrivateParam<bool>("_bank_rays_on_completion", true);
118  params.addPrivateParam<bool>("_ray_dependent_subdomain_setup", true);
119 
120  // Add a point neighbor relationship manager
121  params.addRelationshipManager("ElementPointNeighborLayers",
124  [](const InputParameters &, InputParameters & rm_params)
125  { rm_params.set<unsigned short>("layers") = 1; });
126 
127  return params;
128 }
129 
131  : GeneralUserObject(parameters),
132  _mesh(_fe_problem.mesh()),
133  _comm(_mesh.comm()),
134  _pid(_comm.rank()),
135 
136  _ray_kernel_coverage_check(getParam<bool>("ray_kernel_coverage_check")),
137  _warn_non_planar(getParam<bool>("warn_non_planar")),
138  _use_ray_registration(getParam<bool>("_use_ray_registration")),
139  _use_internal_sidesets(getParam<bool>("use_internal_sidesets")),
140  _tolerate_failure(getParam<bool>("tolerate_failure")),
141  _bank_rays_on_completion(getParam<bool>("_bank_rays_on_completion")),
142  _ray_dependent_subdomain_setup(getParam<bool>("_ray_dependent_subdomain_setup")),
143 
144  _always_cache_traces(getParam<bool>("always_cache_traces")),
145  _data_on_cache_traces(getParam<bool>("data_on_cache_traces")),
146  _aux_data_on_cache_traces(getParam<bool>("aux_data_on_cache_traces")),
147  _segments_on_cache_traces(getParam<bool>("segments_on_cache_traces")),
148  _ray_max_distance(getParam<Real>("ray_distance")),
149  _verify_rays(getParam<bool>("verify_rays")),
150 #ifndef NDEBUG
151  _verify_trace_intersections(getParam<bool>("verify_trace_intersections")),
152 #endif
153 
154  _threaded_elem_side_builders(libMesh::n_threads()),
155 
156  _registered_ray_map(
157  declareRestartableData<std::unordered_map<std::string, RayID>>("registered_ray_map")),
158  _reverse_registered_ray_map(
159  declareRestartableData<std::vector<std::string>>("reverse_registered_ray_map")),
160 
161  _threaded_cached_traces(libMesh::n_threads()),
162 
163  _num_cached(libMesh::n_threads(), 0),
164 
165  _has_non_planar_sides(true),
166  _has_same_level_active_elems(sameLevelActiveElems()),
167 
168  _b_box(MeshTools::create_nodal_bounding_box(_mesh.getMesh())),
169  _domain_max_length(1.01 * (_b_box.max() - _b_box.min()).norm()),
170  _total_volume(computeTotalVolume()),
171 
172  _threaded_cache_ray_kernel(libMesh::n_threads()),
173  _threaded_cache_ray_bc(libMesh::n_threads()),
174  _threaded_ray_object_registration(libMesh::n_threads()),
175  _threaded_current_ray_kernels(libMesh::n_threads()),
176  _threaded_trace_ray(libMesh::n_threads()),
177  _threaded_fe_face(libMesh::n_threads()),
178  _threaded_q_face(libMesh::n_threads()),
179  _threaded_cached_normals(libMesh::n_threads()),
180  _threaded_next_ray_id(libMesh::n_threads()),
181 
182  _parallel_ray_study(std::make_unique<ParallelRayStudy>(*this, _threaded_trace_ray)),
183 
184  _local_trace_ray_results(TraceRay::FAILED_TRACES + 1, 0),
185 
186  _called_initial_setup(false),
187 
188  _elem_index_helper(_mesh.getMesh(), name() + "_elem_index")
189 {
190  // Initialize a tracing object for each thread
191  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
192  {
193  // Initialize a tracing object for each thread
194  _threaded_trace_ray[tid] = std::make_shared<TraceRay>(*this, tid);
195 
196  // Setup the face FEs for normal computation on the fly
197  _threaded_fe_face[tid] = FEBase::build(_mesh.dimension(), FEType(CONSTANT, MONOMIAL));
198  _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
205  {
206  if (!getParam<bool>("allow_other_flags_with_prekernels") && _execute_enum.size() > 1)
207  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  if (_app.useEigenvalue())
213  mooseError("Execution on residual and Jacobian evaluation (execute_on = PRE_KERNELS)\n",
214  "is not supported for an eigenvalue solve.");
215  }
216 
219 
220  // Scale the bounding box for loose checking
222  _loose_b_box.scale(TOLERANCE * TOLERANCE);
223 }
224 
225 void
227 {
228  // Keep track of initialSetup call to avoid registration of various things
229  _called_initial_setup = true;
230 
231  // Sets up a local index for each elem this proc knows about
233 
234  // Check for RayKernel coverage
235  coverageChecks();
236 
237  // Make sure the dependencies exist, if any
239 
240  // Check for traceable element types
242 
243  // Setup for internal sidesets
245 
247 
248  // Setup approximate hmax for each subdomain
250 
251  // Call initial setup on all of the objects
252  for (auto & rto : getRayTracingObjects())
253  rto->initialSetup();
254 
255  // Check for proper exec flags with RayKernels
256  std::vector<RayKernelBase *> ray_kernels;
257  getRayKernels(ray_kernels, 0);
258  for (const auto & rkb : ray_kernels)
259  if (dynamic_cast<RayKernel *>(rkb) && !_execute_enum.isValueSet(EXEC_PRE_KERNELS))
260  mooseError("This study has RayKernel objects that contribute to residuals and Jacobians.",
261  "\nIn this case, the study must use the execute_on = PRE_KERNELS");
262 
263  // Build 1D quadrature rule for along a segment
265  QBase::build(QGAUSS, 1, _fe_problem.getSystemBase(_sys.number()).getMinQuadratureOrder());
266 }
267 
268 void
270 {
271  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
272  mooseAssert(_num_cached[tid] == 0, "Cached residuals/Jacobians not empty");
273 
274  for (auto & rto : getRayTracingObjects())
275  rto->residualSetup();
276 }
277 
278 void
280 {
281  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
282  mooseAssert(_num_cached[tid] == 0, "Cached residuals/Jacobians not empty");
283 
284  for (auto & rto : getRayTracingObjects())
285  rto->jacobianSetup();
286 }
287 
288 void
290 {
291  for (auto & rto : getRayTracingObjects())
292  rto->timestepSetup();
293 }
294 
295 void
297 {
302 
304 
305  for (const auto & trace_ray : _threaded_trace_ray)
306  trace_ray->meshChanged();
307 }
308 
309 void
311 {
312  executeStudy();
313 }
314 
315 void
317 {
318  // Check for coverage of RayKernels on domain
320  {
321  std::vector<RayKernelBase *> ray_kernels;
322  getRayKernels(ray_kernels, 0);
323 
324  std::set<SubdomainID> ray_kernel_blocks;
325  for (const auto & rk : ray_kernels)
326  ray_kernel_blocks.insert(rk->blockIDs().begin(), rk->blockIDs().end());
327 
328  std::set<SubdomainID> missing;
329  std::set_difference(_mesh.meshSubdomains().begin(),
330  _mesh.meshSubdomains().end(),
331  ray_kernel_blocks.begin(),
332  ray_kernel_blocks.end(),
333  std::inserter(missing, missing.begin()));
334 
335  if (!missing.empty() && !ray_kernel_blocks.count(Moose::ANY_BLOCK_ID))
336  {
337  std::ostringstream error;
338  error << "Subdomains { ";
339  std::copy(missing.begin(), missing.end(), std::ostream_iterator<SubdomainID>(error, " "));
340  error << "} do not have RayKernels defined!";
341 
342  mooseError(error.str());
343  }
344  }
345 }
346 
347 void
349 {
350  std::vector<RayTracingObject *> ray_tracing_objects;
351 
352  getRayKernels(ray_tracing_objects, 0);
353  verifyDependenciesExist(ray_tracing_objects);
354 
355  getRayBCs(ray_tracing_objects, 0);
356  verifyDependenciesExist(ray_tracing_objects);
357 }
358 
359 void
360 RayTracingStudy::verifyDependenciesExist(const std::vector<RayTracingObject *> & rtos)
361 {
362  for (const auto & rto : rtos)
363  for (const auto & dep_name : rto->getRequestedItems())
364  {
365  bool found = false;
366  for (const auto & rto_search : rtos)
367  if (rto_search->name() == dep_name)
368  {
369  found = true;
370  break;
371  }
372 
373  if (!found)
374  rto->paramError("depends_on", "The ", rto->getBase(), " '", dep_name, "' does not exist");
375  }
376 }
377 
378 void
380 {
381 
382  for (const auto & elem : *_mesh.getActiveLocalElementRange())
383  {
384  if (_fe_problem.adaptivity().isOn())
385  {
387  mooseError("Element type ",
388  Utility::enum_to_string(elem->type()),
389  " is not supported in ray tracing with adaptivity");
390  }
391  else if (!TraceRayTools::isTraceableElem(elem))
392  mooseError("Element type ",
393  Utility::enum_to_string(elem->type()),
394  " is not supported in ray tracing");
395  }
396 }
397 
398 void
400 {
401  // TODO: We could probably minimize this to local active elements followed by
402  // boundary point neighbors, but if using distibuted mesh it really shouldn't matter
403  _elem_index_helper.initialize(_mesh.getMesh().active_element_ptr_range());
404 }
405 
406 void
408 {
409  // Even if we have _use_internal_sidesets == false, we will make sure the user didn't add RayBCs
410  // on internal boundaries
411 
412  // Clear the data structures and size the map based on the elements that we know about
413  _internal_sidesets.clear();
414  _internal_sidesets_map.clear();
416 
417  // First, we are going to store all elements with internal sidesets (if any) that have active
418  // RayBCs on them as elem -> vector of (side, vector of boundary ids)
419  for (const auto & bnd_elem : *_mesh.getBoundaryElementRange())
420  {
421  Elem * elem = bnd_elem->_elem;
422  const unsigned int side = bnd_elem->_side;
423  const auto bnd_id = bnd_elem->_bnd_id;
424 
425  // Not internal
426  const Elem * const neighbor = elem->neighbor_ptr(side);
427  if (!neighbor || neighbor == remote_elem)
428  continue;
429 
430  // No RayBCs on this sideset
431  std::vector<RayBoundaryConditionBase *> result;
432  getRayBCs(result, bnd_id, 0);
433  if (result.empty())
434  continue;
435 
436  if (neighbor->subdomain_id() == elem->subdomain_id())
437  mooseError("RayBCs exist on internal sidesets that are not bounded by a different",
438  "\nsubdomain on each side.",
439  "\n\nIn order to use RayBCs on internal sidesets, said sidesets must have",
440  "\na different subdomain on each side.");
441 
442  // Mark that this boundary is an internal sideset with RayBC(s)
443  _internal_sidesets.insert(bnd_id);
444 
445  // Get elem's entry in the internal sidset data structure
446  const auto index = _elem_index_helper.getIndex(elem);
447  auto & entry = _internal_sidesets_map[index];
448 
449  // Initialize this elem's sides if they have not been already
450  if (entry.empty())
451  entry.resize(elem->n_sides(), std::vector<BoundaryID>());
452 
453  // Add the internal boundary to the side entry
454  entry[side].push_back(bnd_id);
455  }
456 
458  mooseError("RayBCs are defined on internal sidesets, but the study is not set to use ",
459  "internal sidesets during tracing.",
460  "\n\nSet the parameter use_internal_sidesets = true to enable this capability.");
461 }
462 
463 void
465 {
466  _has_non_planar_sides = false;
467  bool warned = !_warn_non_planar;
468 
469  // Nothing to do here for 2D or 1D
470  if (_mesh.dimension() != 3)
471  return;
472 
473  // Clear the data structure and size it based on the elements that we know about
474  _non_planar_sides.clear();
476 
477  for (const Elem * elem : _mesh.getMesh().active_element_ptr_range())
478  {
479  const auto index = _elem_index_helper.getIndex(elem);
480  auto & entry = _non_planar_sides[index];
481  entry.resize(elem->n_sides(), 0);
482 
483  for (const auto s : elem->side_index_range())
484  {
485  const auto & side = elemSide(*elem, s);
486  if (side.n_vertices() < 4)
487  continue;
488 
489  if (!side.has_affine_map())
490  {
491  entry[s] = 1;
492  _has_non_planar_sides = true;
493 
494  if (!warned)
495  {
496  mooseWarning("The mesh contains non-planar faces.\n\n",
497  "Ray tracing on non-planar faces is an approximation and may fail.\n\n",
498  "Use at your own risk! You can disable this warning by setting the\n",
499  "parameter 'warn_non_planar' to false.");
500  warned = true;
501  }
502  }
503  }
504  }
505 }
506 
507 void
509 {
510  // Setup map with subdomain keys
511  _subdomain_hmax.clear();
512  for (const auto subdomain_id : _mesh.meshSubdomains())
513  _subdomain_hmax[subdomain_id] = std::numeric_limits<Real>::min();
514 
515  // Set local max for each subdomain
516  for (const auto & elem : *_mesh.getActiveLocalElementRange())
517  {
518  auto & entry = _subdomain_hmax.at(elem->subdomain_id());
519  entry = std::max(entry, elem->hmax());
520  }
521 
522  // Accumulate global max for each subdomain
524 
525  if (getParam<bool>("warn_subdomain_hmax"))
526  {
527  const auto warn_prefix = type() + " '" + name() + "': ";
528  const auto warn_suffix =
529  "\n\nRay tracing uses an approximate element size for each subdomain to scale the\n"
530  "tolerances used in computing ray intersections. This warning suggests that the\n"
531  "approximate element size is not a good approximation. This is likely due to poor\n"
532  "element aspect ratios.\n\n"
533  "This warning is only output for the first element affected.\n"
534  "To disable this warning, set warn_subdomain_hmax = false.\n";
535 
536  for (const auto & elem : *_mesh.getActiveLocalElementRange())
537  {
538  const auto hmin = elem->hmin();
539  const auto hmax = elem->hmax();
540  const auto max_hmax = subdomainHmax(elem->subdomain_id());
541 
542  const auto hmax_rel = hmax / max_hmax;
543  if (hmax_rel < 1.e-2 || hmax_rel > 1.e2)
544  mooseDoOnce(mooseWarning(warn_prefix,
545  "Element hmax varies significantly from subdomain hmax.\n",
546  warn_suffix,
547  "First element affected:\n",
548  Moose::stringify(*elem)););
549 
550  const auto h_rel = max_hmax / hmin;
551  if (h_rel > 1.e2)
552  mooseDoOnce(mooseWarning(warn_prefix,
553  "Element hmin varies significantly from subdomain hmax.\n",
554  warn_suffix,
555  "First element affected:\n",
556  Moose::stringify(*elem)););
557  }
558  }
559 }
560 
561 void
563 {
564  // First, clear the objects associated with each Ray on each thread
565  const auto num_rays = _registered_ray_map.size();
566  for (auto & entry : _threaded_ray_object_registration)
567  {
568  entry.clear();
569  entry.resize(num_rays);
570  }
571 
572  const auto rtos = getRayTracingObjects();
573 
575  {
576  // All of the registered ray names - used when a RayTracingObject did not specify
577  // any Rays so it should be associated with all Rays.
578  std::vector<std::string> all_ray_names;
579  all_ray_names.reserve(_registered_ray_map.size());
580  for (const auto & pair : _registered_ray_map)
581  all_ray_names.push_back(pair.first);
582 
583  for (auto & rto : rtos)
584  {
585  // The Ray names associated with this RayTracingObject
586  const auto & ray_names = rto->parameters().get<std::vector<std::string>>("rays");
587  // The registration for RayTracingObjects for the thread rto is on
588  const auto tid = rto->parameters().get<THREAD_ID>("_tid");
589  auto & registration = _threaded_ray_object_registration[tid];
590 
591  // Register each Ray for this object in the registration
592  for (const auto & ray_name : (ray_names.empty() ? all_ray_names : ray_names))
593  {
594  const auto id = registeredRayID(ray_name, /* graceful = */ true);
595  if (ray_names.size() && id == Ray::INVALID_RAY_ID)
596  rto->paramError(
597  "rays", "Supplied ray '", ray_name, "' is not a registered Ray in ", typeAndName());
598  registration[id].insert(rto);
599  }
600  }
601  }
602  // Not using Ray registration
603  else
604  {
605  for (const auto & rto : rtos)
606  if (rto->parameters().get<std::vector<std::string>>("rays").size())
607  rto->paramError(
608  "rays",
609  "Rays cannot be supplied when the study does not require Ray registration.\n\n",
610  type(),
611  " does not require Ray registration.");
612  }
613 }
614 
615 void
617 {
618  std::set<std::string> vars_to_be_zeroed;
619  std::vector<RayKernelBase *> ray_kernels;
620  getRayKernels(ray_kernels, 0);
621  for (auto & rk : ray_kernels)
622  {
623  AuxRayKernel * aux_rk = dynamic_cast<AuxRayKernel *>(rk);
624  if (aux_rk)
625  vars_to_be_zeroed.insert(aux_rk->variable().name());
626  }
627 
628  std::vector<std::string> vars_to_be_zeroed_vec(vars_to_be_zeroed.begin(),
629  vars_to_be_zeroed.end());
630  _fe_problem.getAuxiliarySystem().zeroVariables(vars_to_be_zeroed_vec);
631 }
632 
633 void
635  const THREAD_ID tid,
636  const RayID ray_id)
637 {
638  mooseAssert(currentlyPropagating(), "Should not call while not propagating");
639 
640  // Call subdomain setup on FE
641  _fe_problem.subdomainSetup(subdomain, tid);
642 
643  std::set<MooseVariableFEBase *> needed_moose_vars;
644  std::unordered_set<unsigned int> needed_mat_props;
645 
646  // Get RayKernels and their dependencies and call subdomain setup
647  getRayKernels(_threaded_current_ray_kernels[tid], subdomain, tid, ray_id);
648  for (auto & rkb : _threaded_current_ray_kernels[tid])
649  {
650  rkb->subdomainSetup();
651 
652  const auto & mv_deps = rkb->getMooseVariableDependencies();
653  needed_moose_vars.insert(mv_deps.begin(), mv_deps.end());
654 
655  const auto & mp_deps = rkb->getMatPropDependencies();
656  needed_mat_props.insert(mp_deps.begin(), mp_deps.end());
657  }
658 
659  // Prepare aux vars
660  for (auto & var : needed_moose_vars)
661  if (var->kind() == Moose::VarKindType::VAR_AUXILIARY)
662  var->prepareAux();
663 
664  _fe_problem.setActiveElementalMooseVariables(needed_moose_vars, tid);
665  _fe_problem.prepareMaterials(needed_mat_props, subdomain, tid);
666 }
667 
668 void
670  const Elem * elem, const Point & start, const Point & end, const Real length, THREAD_ID tid)
671 {
672  mooseAssert(MooseUtils::absoluteFuzzyEqual((start - end).norm(), length), "Invalid length");
673  mooseAssert(currentlyPropagating(), "Should not call while not propagating");
674 
676 
677  // If we have any variables or material properties that are active, we definitely need to reinit
678  bool reinit = _fe_problem.hasActiveElementalMooseVariables(tid) ||
680  // If not, make sure that the RayKernels have not requested a reinit (this could happen when a
681  // RayKernel doesn't have variables or materials but still does an integration and needs qps)
682  if (!reinit)
683  for (const RayKernelBase * rk : currentRayKernels(tid))
684  if (rk->needSegmentReinit())
685  {
686  reinit = true;
687  break;
688  }
689 
690  if (reinit)
691  {
692  _fe_problem.prepare(elem, tid);
693 
694  std::vector<Point> points;
695  std::vector<Real> weights;
696  buildSegmentQuadrature(start, end, length, points, weights);
697  _fe_problem.reinitElemPhys(elem, points, tid);
698  _fe_problem.assembly(tid, _sys.number()).modifyArbitraryWeights(weights);
699 
701  }
702 }
703 
704 void
706  const Point & end,
707  const Real length,
708  std::vector<Point> & points,
709  std::vector<Real> & weights) const
710 {
711  points.resize(_segment_qrule->n_points());
712  weights.resize(_segment_qrule->n_points());
713 
714  const Point diff = end - start;
715  const Point sum = end + start;
716  mooseAssert(MooseUtils::absoluteFuzzyEqual(length, diff.norm()), "Invalid length");
717 
718  // The standard quadrature rule should be on x = [-1, 1]
719  // To scale the points, you...
720  // - Scale to size of the segment in 3D
721  // initial_scaled_qp = x_qp * 0.5 * (end - start) = 0.5 * x_qp * diff
722  // - Shift quadrature midpoint to segment midpoint
723  // final_qp = initial_scaled_qp + 0.5 * (end - start) = initial_scaled_qp + 0.5 * sum
724  // = 0.5 * (x_qp * diff + sum)
725  for (unsigned int qp = 0; qp < _segment_qrule->n_points(); ++qp)
726  {
727  points[qp] = 0.5 * (_segment_qrule->qp(qp)(0) * diff + sum);
728  weights[qp] = 0.5 * _segment_qrule->w(qp) * length;
729  }
730 }
731 
732 void
733 RayTracingStudy::postOnSegment(const THREAD_ID tid, const std::shared_ptr<Ray> & /* ray */)
734 {
735  mooseAssert(currentlyPropagating(), "Should not call while not propagating");
737  mooseAssert(_num_cached[tid] == 0,
738  "Values should only be cached when computing Jacobian/residual");
739 
740  // Fill into cached Jacobian/residuals if necessary
742  {
744 
745  if (++_num_cached[tid] == 20)
746  {
747  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
749  _num_cached[tid] = 0;
750  }
751  }
753  {
755 
756  if (++_num_cached[tid] == 20)
757  {
758  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
760  _num_cached[tid] = 0;
761  }
762  }
763 }
764 
765 void
767 {
768  TIME_SECTION("executeStudy", 2, "Executing Study");
769 
770  mooseAssert(_called_initial_setup, "Initial setup not called");
771 
772  // Reset ray start/complete timers
774  _max_intersections = 0;
776 
777  // Reset physical tracing stats
778  for (auto & val : _local_trace_ray_results)
779  val = 0;
780 
781  // Reset crossing and intersection
788  _ending_distance = 0;
789  _total_distance = 0;
790 
791  // Zero the AuxVariables that our AuxRayKernels contribute to before they accumulate
793 
794  preExecuteStudy();
795  for (auto & rto : getRayTracingObjects())
796  rto->preExecuteStudy();
797 
798  _ray_bank.clear();
799 
800  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
801  {
802  _threaded_trace_ray[tid]->preExecute();
803  _threaded_cached_normals[tid].clear();
804  }
805 
807  _execution_start_time = std::chrono::steady_clock::now();
808 
809  _parallel_ray_study->preExecute();
810 
811  {
812  {
813  auto generation_start_time = std::chrono::steady_clock::now();
814 
815  TIME_SECTION("generateRays", 2, "Generating Rays");
816 
817  generateRays();
818 
819  _generation_time = std::chrono::steady_clock::now() - generation_start_time;
820  }
821 
822  // At this point, nobody is working so this is good time to make sure
823  // Rays are unique across all processors in the working buffer
824  if (verifyRays())
825  {
826  verifyUniqueRays(_parallel_ray_study->workBuffer().begin(),
827  _parallel_ray_study->workBuffer().end(),
828  /* error_suffix = */ "after generateRays()");
829 
830  verifyUniqueRayIDs(_parallel_ray_study->workBuffer().begin(),
831  _parallel_ray_study->workBuffer().end(),
832  /* global = */ true,
833  /* error_suffix = */ "after generateRays()");
834  }
835 
837 
838  {
839  TIME_SECTION("propagateRays", 2, "Propagating Rays");
840 
841  const auto propagation_start_time = std::chrono::steady_clock::now();
842 
843  _parallel_ray_study->execute();
844 
845  _propagation_time = std::chrono::steady_clock::now() - propagation_start_time;
846  }
847  }
848 
849  _execution_time = std::chrono::steady_clock::now() - _execution_start_time;
850 
851  if (verifyRays())
852  {
853  verifyUniqueRays(_parallel_ray_study->workBuffer().begin(),
854  _parallel_ray_study->workBuffer().end(),
855  /* error_suffix = */ "after tracing completed");
856 
857 #ifndef NDEBUG
858  // Outside of debug, _ray_bank always holds all of the Rays that have ended on this processor
859  // We can use this as a global point to check for unique IDs for every Ray that has traced
861  _ray_bank.end(),
862  /* global = */ true,
863  /* error_suffix = */ "after tracing completed");
864 #endif
865  }
866 
867  // Update counters from the threaded trace objects
868  for (const auto & tr : _threaded_trace_ray)
869  for (std::size_t i = 0; i < _local_trace_ray_results.size(); ++i)
870  _local_trace_ray_results[i] += tr->results()[i];
871 
872  // Update local ending counters
879  // ...and communicate the global values
886 
887  // Throw a warning with the number of failed (tolerated) traces
888  if (_tolerate_failure)
889  {
891  _communicator.sum(failures);
892  if (failures)
893  mooseWarning(
894  type(), " '", name(), "': ", failures, " ray tracing failures were tolerated.\n");
895  }
896 
897  // Clear the current RayKernels
898  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
900 
901  // Move the threaded cache trace information into the full cached trace vector
902  // Here, we only clear the cached vectors so that we might not have to
903  // reallocate on future traces
904  std::size_t num_entries = 0;
905  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
906  num_entries += _threaded_cached_traces[tid].size();
907  _cached_traces.clear();
908  _cached_traces.reserve(num_entries);
909  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
910  {
911  for (const auto & entry : _threaded_cached_traces[tid])
912  _cached_traces.emplace_back(std::move(entry));
913  _threaded_cached_traces[tid].clear();
914  }
915 
916  // Add any stragglers that contribute to the Jacobian or residual
917  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
918  if (_num_cached[tid] != 0)
919  {
920  mooseAssert(_fe_problem.currentlyComputingJacobian() ||
922  "Should not have cached values without Jacobian/residual computation");
923 
926  else
928 
929  _num_cached[tid] = 0;
930  }
931 
932  // AuxRayKernels may have modified AuxVariables
935 
936  // Clear FE
937  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
938  {
941  }
942 
944  for (auto & rto : getRayTracingObjects())
945  rto->postExecuteStudy();
946 }
947 
948 void
949 RayTracingStudy::onCompleteRay(const std::shared_ptr<Ray> & ray)
950 {
951  mooseAssert(currentlyPropagating(), "Should only be called during Ray propagation");
952 
953  _ending_processor_crossings += ray->processorCrossings();
955  std::max(_ending_max_processor_crossings, ray->processorCrossings());
956  _ending_intersections += ray->intersections();
957  _ending_max_intersections = std::max(_ending_max_intersections, ray->intersections());
959  std::max(_ending_max_trajectory_changes, ray->trajectoryChanges());
960  _ending_distance += ray->distance();
961 
962 #ifdef NDEBUG
963  // In non-opt modes, we will always bank the Rays for debugging
965 #endif
966  _ray_bank.emplace_back(ray);
967 }
968 
970 RayTracingStudy::registerRayDataInternal(const std::string & name, const bool aux)
971 {
972  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
973 
975  mooseError("Cannot register Ray ", (aux ? "aux " : ""), "data after initialSetup()");
976 
977  auto & map = aux ? _ray_aux_data_map : _ray_data_map;
978  const auto find = map.find(name);
979  if (find != map.end())
980  return find->second;
981 
982  auto & other_map = aux ? _ray_data_map : _ray_aux_data_map;
983  if (other_map.find(name) != other_map.end())
984  mooseError("Cannot register Ray aux data with name ",
985  name,
986  " because Ray ",
987  (aux ? "(non-aux)" : "aux"),
988  " data already exists with said name.");
989 
990  // Add into the name -> index map
991  map.emplace(name, map.size());
992 
993  // Add into the index -> names vector
994  auto & vector = aux ? _ray_aux_data_names : _ray_data_names;
995  vector.push_back(name);
996 
997  return map.size() - 1;
998 }
999 
1000 std::vector<RayDataIndex>
1001 RayTracingStudy::registerRayDataInternal(const std::vector<std::string> & names, const bool aux)
1002 {
1003  std::vector<RayDataIndex> indices(names.size());
1004  for (std::size_t i = 0; i < names.size(); ++i)
1005  indices[i] = registerRayDataInternal(names[i], aux);
1006  return indices;
1007 }
1008 
1011  const bool aux,
1012  const bool graceful) const
1013 {
1014  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1015 
1016  const auto & map = aux ? _ray_aux_data_map : _ray_data_map;
1017  const auto find = map.find(name);
1018  if (find != map.end())
1019  return find->second;
1020 
1021  if (graceful)
1023 
1024  const auto & other_map = aux ? _ray_data_map : _ray_aux_data_map;
1025  if (other_map.find(name) != other_map.end())
1026  mooseError("Ray data with name '",
1027  name,
1028  "' was not found.\n\n",
1029  "However, Ray ",
1030  (aux ? "non-aux" : "aux"),
1031  " data with said name was found.\n",
1032  "Did you mean to use ",
1033  (aux ? "getRayDataIndex()/getRayDataIndices()?"
1034  : "getRayAuxDataIndex()/getRayAuxDataIndices()"),
1035  "?");
1036 
1037  mooseError("Unknown Ray ", (aux ? "aux " : ""), "data with name ", name);
1038 }
1039 
1040 std::vector<RayDataIndex>
1041 RayTracingStudy::getRayDataIndicesInternal(const std::vector<std::string> & names,
1042  const bool aux,
1043  const bool graceful) const
1044 {
1045  std::vector<RayDataIndex> indices(names.size());
1046  for (std::size_t i = 0; i < names.size(); ++i)
1047  indices[i] = getRayDataIndexInternal(names[i], aux, graceful);
1048  return indices;
1049 }
1050 
1051 const std::string &
1052 RayTracingStudy::getRayDataNameInternal(const RayDataIndex index, const bool aux) const
1053 {
1054  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1055 
1056  if ((aux ? rayAuxDataSize() : rayDataSize()) < index)
1057  mooseError("Unknown Ray ", aux ? "aux " : "", "data with index ", index);
1058  return aux ? _ray_aux_data_names[index] : _ray_data_names[index];
1059 }
1060 
1063 {
1064  return registerRayDataInternal(name, /* aux = */ false);
1065 }
1066 
1067 std::vector<RayDataIndex>
1068 RayTracingStudy::registerRayData(const std::vector<std::string> & names)
1069 {
1070  return registerRayDataInternal(names, /* aux = */ false);
1071 }
1072 
1074 RayTracingStudy::getRayDataIndex(const std::string & name, const bool graceful /* = false */) const
1075 {
1076  return getRayDataIndexInternal(name, /* aux = */ false, graceful);
1077 }
1078 
1079 std::vector<RayDataIndex>
1080 RayTracingStudy::getRayDataIndices(const std::vector<std::string> & names,
1081  const bool graceful /* = false */) const
1082 {
1083  return getRayDataIndicesInternal(names, /* aux = */ false, graceful);
1084 }
1085 
1086 const std::string &
1088 {
1089  return getRayDataNameInternal(index, /* aux = */ false);
1090 }
1091 
1094 {
1095  return registerRayDataInternal(name, /* aux = */ true);
1096 }
1097 
1098 std::vector<RayDataIndex>
1099 RayTracingStudy::registerRayAuxData(const std::vector<std::string> & names)
1100 {
1101  return registerRayDataInternal(names, /* aux = */ true);
1102 }
1103 
1106  const bool graceful /* = false */) const
1107 {
1108  return getRayDataIndexInternal(name, /* aux = */ true, graceful);
1109 }
1110 
1111 std::vector<RayDataIndex>
1112 RayTracingStudy::getRayAuxDataIndices(const std::vector<std::string> & names,
1113  const bool graceful /* = false */) const
1114 {
1115  return getRayDataIndicesInternal(names, /* aux = */ true, graceful);
1116 }
1117 
1118 const std::string &
1120 {
1121  return getRayDataNameInternal(index, /* aux = */ true);
1122 }
1123 
1124 bool
1126 {
1127  std::vector<RayKernelBase *> result;
1128  getRayKernels(result, tid);
1129  return result.size();
1130 }
1131 
1132 void
1133 RayTracingStudy::getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid)
1134 {
1135  // If the cache doesn't have any attributes yet, it means that we haven't set
1136  // the conditions yet. We do this so that it can be generated on the fly on first use.
1137  if (!_threaded_cache_ray_kernel[tid].numAttribs())
1138  {
1139  if (!_called_initial_setup)
1140  mooseError("Should not call getRayKernels() before initialSetup()");
1141 
1142  auto query = _fe_problem.theWarehouse()
1143  .query()
1144  .condition<AttribRayTracingStudy>(this)
1145  .condition<AttribSystem>("RayKernel")
1146  .condition<AttribThread>(tid);
1147  _threaded_cache_ray_kernel[tid] = query.clone();
1148  }
1149 
1150  _threaded_cache_ray_kernel[tid].queryInto(result, id);
1151 }
1152 
1153 void
1154 RayTracingStudy::getRayKernels(std::vector<RayKernelBase *> & result,
1155  SubdomainID id,
1156  THREAD_ID tid,
1157  RayID ray_id)
1158 {
1159  // No Ray registration: no need to sift through objects
1160  if (!_use_ray_registration)
1161  {
1162  getRayKernels(result, id, tid);
1163  }
1164  // Has Ray registration: only pick the objects associated with ray_id
1165  else
1166  {
1167  // Get all of the kernels on this block
1168  std::vector<RayKernelBase *> rkbs;
1169  getRayKernels(rkbs, id, tid);
1170 
1171  // The RayTracingObjects associated with this ray
1172  const auto & ray_id_rtos = _threaded_ray_object_registration[tid][ray_id];
1173 
1174  // The result is the union of all of the kernels and the objects associated with this Ray
1175  result.clear();
1176  for (auto rkb : rkbs)
1177  if (ray_id_rtos.count(rkb))
1178  result.push_back(rkb);
1179  }
1180 }
1181 
1182 void
1183 RayTracingStudy::getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
1184  BoundaryID id,
1185  THREAD_ID tid)
1186 {
1187  // If the cache doesn't have any attributes yet, it means that we haven't set
1188  // the conditions yet. We do this so that it can be generated on the fly on first use.
1189  if (!_threaded_cache_ray_bc[tid].numAttribs())
1190  {
1191  if (!_called_initial_setup)
1192  mooseError("Should not call getRayBCs() before initialSetup()");
1193 
1194  auto query = _fe_problem.theWarehouse()
1195  .query()
1196  .condition<AttribRayTracingStudy>(this)
1197  .condition<AttribSystem>("RayBoundaryCondition")
1198  .condition<AttribThread>(tid);
1199  _threaded_cache_ray_bc[tid] = query.clone();
1200  }
1201 
1202  _threaded_cache_ray_bc[tid].queryInto(result, std::make_tuple(id, false));
1203 }
1204 
1205 void
1206 RayTracingStudy::getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
1207  const std::vector<TraceRayBndElement> & bnd_elems,
1208  THREAD_ID tid,
1209  RayID ray_id)
1210 {
1211  // No Ray registration: no need to sift through objects
1212  if (!_use_ray_registration)
1213  {
1214  if (bnd_elems.size() == 1)
1215  getRayBCs(result, bnd_elems[0].bnd_id, tid);
1216  else
1217  {
1218  std::vector<BoundaryID> bnd_ids(bnd_elems.size());
1219  for (MooseIndex(bnd_elems.size()) i = 0; i < bnd_elems.size(); ++i)
1220  bnd_ids[i] = bnd_elems[i].bnd_id;
1221  getRayBCs(result, bnd_ids, tid);
1222  }
1223  }
1224  // Has Ray registration: only pick the objects associated with ray_id
1225  else
1226  {
1227  // Get all of the RayBCs on these boundaries
1228  std::vector<RayBoundaryConditionBase *> rbcs;
1229  if (bnd_elems.size() == 1)
1230  getRayBCs(rbcs, bnd_elems[0].bnd_id, tid);
1231  else
1232  {
1233  std::vector<BoundaryID> bnd_ids(bnd_elems.size());
1234  for (MooseIndex(bnd_elems.size()) i = 0; i < bnd_elems.size(); ++i)
1235  bnd_ids[i] = bnd_elems[i].bnd_id;
1236  getRayBCs(rbcs, bnd_ids, tid);
1237  }
1238 
1239  // The RayTracingObjects associated with this ray
1240  mooseAssert(ray_id < _threaded_ray_object_registration[tid].size(), "Not in registration");
1241  const auto & ray_id_rtos = _threaded_ray_object_registration[tid][ray_id];
1242 
1243  // The result is the union of all of the kernels and the objects associated with this Ray
1244  result.clear();
1245  for (auto rbc : rbcs)
1246  if (ray_id_rtos.count(rbc))
1247  result.push_back(rbc);
1248  }
1249 }
1250 
1251 std::vector<RayTracingObject *>
1253 {
1254  std::vector<RayTracingObject *> result;
1255  _fe_problem.theWarehouse().query().condition<AttribRayTracingStudy>(this).queryInto(result);
1256  return result;
1257 }
1258 
1259 const std::vector<std::shared_ptr<Ray>> &
1261 {
1263  mooseError("The Ray bank is not available because the private parameter "
1264  "'_bank_rays_on_completion' is set to false.");
1266  mooseError("Cannot get the Ray bank during generation or propagation.");
1267 
1268  return _ray_bank;
1269 }
1270 
1271 std::shared_ptr<Ray>
1273 {
1274  // This is only a linear search - can be improved on with a map in the future
1275  // if this is used on a larger scale
1276  std::shared_ptr<Ray> ray;
1277  for (const std::shared_ptr<Ray> & possible_ray : rayBank())
1278  if (possible_ray->id() == ray_id)
1279  {
1280  ray = possible_ray;
1281  break;
1282  }
1283 
1284  // Make sure one and only one processor has the Ray
1285  unsigned int have_ray = ray ? 1 : 0;
1286  _communicator.sum(have_ray);
1287  if (have_ray == 0)
1288  mooseError("Could not find a Ray with the ID ", ray_id, " in the Ray banks.");
1289 
1290  // This should never happen... but let's make sure
1291  mooseAssert(have_ray == 1, "Multiple rays with the same ID were found in the Ray banks");
1292 
1293  return ray;
1294 }
1295 
1296 RayData
1298  const RayDataIndex index,
1299  const bool aux) const
1300 {
1301  // Will be a nullptr shared_ptr if this processor doesn't own the Ray
1302  const std::shared_ptr<Ray> ray = getBankedRay(ray_id);
1303 
1304  Real value = ray ? (aux ? ray->auxData(index) : ray->data(index)) : 0;
1306  return value;
1307 }
1308 
1309 RayData
1310 RayTracingStudy::getBankedRayData(const RayID ray_id, const RayDataIndex index) const
1311 {
1312  return getBankedRayDataInternal(ray_id, index, /* aux = */ false);
1313 }
1314 
1315 RayData
1317 {
1318  return getBankedRayDataInternal(ray_id, index, /* aux = */ true);
1319 }
1320 
1321 RayID
1323 {
1324  libmesh_parallel_only(comm());
1325 
1326  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1327 
1328  if (!_use_ray_registration)
1329  mooseError("Cannot use registerRay() with Ray registration disabled");
1330 
1331  // This is parallel only for now. We could likely stagger the ID building like we do with
1332  // the unique IDs, but it would require a sync point which isn't there right now
1333  libmesh_parallel_only(comm());
1334 
1335  const auto & it = _registered_ray_map.find(name);
1336  if (it != _registered_ray_map.end())
1337  return it->second;
1338 
1339  const auto id = _reverse_registered_ray_map.size();
1340  _registered_ray_map.emplace(name, id);
1341  _reverse_registered_ray_map.push_back(name);
1342  return id;
1343 }
1344 
1345 RayID
1346 RayTracingStudy::registeredRayID(const std::string & name, const bool graceful /* = false */) const
1347 {
1348  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1349 
1350  if (!_use_ray_registration)
1351  mooseError("Should not use registeredRayID() with Ray registration disabled");
1352 
1353  const auto search = _registered_ray_map.find(name);
1354  if (search != _registered_ray_map.end())
1355  return search->second;
1356 
1357  if (graceful)
1358  return Ray::INVALID_RAY_ID;
1359 
1360  mooseError("Attempted to obtain ID of registered Ray ",
1361  name,
1362  ", but a Ray with said name is not registered.");
1363 }
1364 
1365 const std::string &
1367 {
1368  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1369 
1370  if (!_use_ray_registration)
1371  mooseError("Should not use registeredRayName() with Ray registration disabled");
1372 
1373  if (_reverse_registered_ray_map.size() > ray_id)
1374  return _reverse_registered_ray_map[ray_id];
1375 
1376  mooseError("Attempted to obtain name of registered Ray with ID ",
1377  ray_id,
1378  ", but a Ray with said ID is not registered.");
1379 }
1380 
1381 Real
1383 {
1384  Real volume = 0;
1385  for (const auto & elem : *_mesh.getActiveLocalElementRange())
1386  volume += elem->volume();
1388  return volume;
1389 }
1390 
1391 const std::vector<std::vector<BoundaryID>> &
1392 RayTracingStudy::getInternalSidesets(const Elem * elem) const
1393 {
1394  mooseAssert(_use_internal_sidesets, "Not using internal sidesets");
1395  mooseAssert(hasInternalSidesets(), "Processor does not have internal sidesets");
1396  mooseAssert(_internal_sidesets_map.size() > _elem_index_helper.getIndex(elem),
1397  "Internal sideset map not initialized");
1398 
1399  const auto index = _elem_index_helper.getIndex(elem);
1400  return _internal_sidesets_map[index];
1401 }
1402 
1403 TraceData &
1404 RayTracingStudy::initThreadedCachedTrace(const std::shared_ptr<Ray> & ray, THREAD_ID tid)
1405 {
1406  mooseAssert(shouldCacheTrace(ray), "Not caching trace");
1407  mooseAssert(currentlyPropagating(), "Should only use while tracing");
1408 
1409  _threaded_cached_traces[tid].emplace_back(ray);
1410  return _threaded_cached_traces[tid].back();
1411 }
1412 
1413 void
1414 RayTracingStudy::verifyUniqueRayIDs(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
1415  const std::vector<std::shared_ptr<Ray>>::const_iterator end,
1416  const bool global,
1417  const std::string & error_suffix) const
1418 {
1419  // Determine the unique set of Ray IDs on this processor,
1420  // and if not locally unique throw an error. Once we build this set,
1421  // we will send it to rank 0 to verify globally
1422  std::set<RayID> local_rays;
1423  for (const std::shared_ptr<Ray> & ray : as_range(begin, end))
1424  {
1425  mooseAssert(ray, "Null ray");
1426 
1427  // Try to insert into the set; the second entry in the pair
1428  // will be false if it was not inserted
1429  if (!local_rays.insert(ray->id()).second)
1430  {
1431  for (const std::shared_ptr<Ray> & other_ray : as_range(begin, end))
1432  if (ray.get() != other_ray.get() && ray->id() == other_ray->id())
1433  mooseError("Multiple Rays exist with ID ",
1434  ray->id(),
1435  " on processor ",
1436  _pid,
1437  " ",
1438  error_suffix,
1439  "\n\nOffending Ray information:\n\n",
1440  ray->getInfo(),
1441  "\n",
1442  other_ray->getInfo());
1443  }
1444  }
1445 
1446  // Send IDs from all procs to rank 0 and verify on rank 0
1447  if (global)
1448  {
1449  // Package our local IDs and send to rank 0
1450  std::map<processor_id_type, std::vector<RayID>> send_ids;
1451  if (local_rays.size())
1452  send_ids.emplace(std::piecewise_construct,
1453  std::forward_as_tuple(0),
1454  std::forward_as_tuple(local_rays.begin(), local_rays.end()));
1455  local_rays.clear();
1456 
1457  // Mapping on rank 0 from ID -> processor ID
1458  std::map<RayID, processor_id_type> global_map;
1459 
1460  // Verify another processor's IDs against the global map on rank 0
1461  const auto check_ids =
1462  [this, &global_map, &error_suffix](processor_id_type pid, const std::vector<RayID> & ids)
1463  {
1464  for (const RayID id : ids)
1465  {
1466  const auto emplace_pair = global_map.emplace(id, pid);
1467 
1468  // Means that this ID already exists in the map
1469  if (!emplace_pair.second)
1470  mooseError("Ray with ID ",
1471  id,
1472  " exists on ranks ",
1473  emplace_pair.first->second,
1474  " and ",
1475  pid,
1476  "\n",
1477  error_suffix);
1478  }
1479  };
1480 
1481  Parallel::push_parallel_vector_data(_communicator, send_ids, check_ids);
1482  }
1483 }
1484 
1485 void
1486 RayTracingStudy::verifyUniqueRays(const std::vector<std::shared_ptr<Ray>>::const_iterator begin,
1487  const std::vector<std::shared_ptr<Ray>>::const_iterator end,
1488  const std::string & error_suffix)
1489 {
1490  std::set<const Ray *> rays;
1491  for (const std::shared_ptr<Ray> & ray : as_range(begin, end))
1492  if (!rays.insert(ray.get()).second) // false if not inserted into rays
1493  mooseError("Multiple shared_ptrs were found that point to the same Ray ",
1494  error_suffix,
1495  "\n\nOffending Ray:\n",
1496  ray->getInfo());
1497 }
1498 
1499 void
1500 RayTracingStudy::moveRayToBuffer(std::shared_ptr<Ray> & ray)
1501 {
1502  mooseAssert(currentlyGenerating(), "Can only use while generating");
1503  mooseAssert(ray, "Null ray");
1504  mooseAssert(ray->shouldContinue(), "Ray is not continuing");
1505 
1506  _parallel_ray_study->moveWorkToBuffer(ray, /* tid = */ 0);
1507 }
1508 
1509 void
1510 RayTracingStudy::moveRaysToBuffer(std::vector<std::shared_ptr<Ray>> & rays)
1511 {
1512  mooseAssert(currentlyGenerating(), "Can only use while generating");
1513 #ifndef NDEBUG
1514  for (const std::shared_ptr<Ray> & ray : rays)
1515  {
1516  mooseAssert(ray, "Null ray");
1517  mooseAssert(ray->shouldContinue(), "Ray is not continuing");
1518  }
1519 #endif
1520 
1521  _parallel_ray_study->moveWorkToBuffer(rays, /* tid = */ 0);
1522 }
1523 
1524 void
1526  const THREAD_ID tid,
1527  const AcquireMoveDuringTraceKey &)
1528 {
1529  mooseAssert(ray, "Null ray");
1530  mooseAssert(currentlyPropagating(), "Can only use while tracing");
1531 
1532  _parallel_ray_study->moveWorkToBuffer(ray, tid);
1533 }
1534 
1535 void
1536 RayTracingStudy::reserveRayBuffer(const std::size_t size)
1537 {
1538  if (!currentlyGenerating())
1539  mooseError("Can only reserve in Ray buffer during generateRays()");
1540 
1541  _parallel_ray_study->reserveBuffer(size);
1542 }
1543 
1544 const Point &
1545 RayTracingStudy::getSideNormal(const Elem * elem, unsigned short side, const THREAD_ID tid)
1546 {
1547  std::unordered_map<std::pair<const Elem *, unsigned short>, Point> & cache =
1549 
1550  // See if we've already cached this side normal
1551  const auto elem_side_pair = std::make_pair(elem, side);
1552  const auto search = cache.find(elem_side_pair);
1553 
1554  // Haven't cached this side normal: compute it and then cache it
1555  if (search == cache.end())
1556  {
1557  _threaded_fe_face[tid]->reinit(elem, side);
1558  const auto & normal = _threaded_fe_face[tid]->get_normals()[0];
1559  cache.emplace(elem_side_pair, normal);
1560  return normal;
1561  }
1562 
1563  // Have cached this side normal: simply return it
1564  return search->second;
1565 }
1566 
1567 bool
1569 {
1570  unsigned int min_level = std::numeric_limits<unsigned int>::max();
1571  unsigned int max_level = std::numeric_limits<unsigned int>::min();
1572 
1573  for (const auto & elem : *_mesh.getActiveLocalElementRange())
1574  {
1575  const auto level = elem->level();
1576  min_level = std::min(level, min_level);
1577  max_level = std::max(level, max_level);
1578  }
1579 
1580  _communicator.min(min_level);
1581  _communicator.max(max_level);
1582 
1583  return min_level == max_level;
1584 }
1585 
1586 Real
1588 {
1589  const auto find = _subdomain_hmax.find(subdomain_id);
1590  if (find == _subdomain_hmax.end())
1591  mooseError("Subdomain ", subdomain_id, " not found in subdomain hmax map");
1592  return find->second;
1593 }
1594 
1595 bool
1597 {
1598  Real bbox_volume = 1;
1599  for (unsigned int d = 0; d < _mesh.dimension(); ++d)
1600  bbox_volume *= std::abs(_b_box.max()(d) - _b_box.min()(d));
1601 
1602  return MooseUtils::absoluteFuzzyEqual(bbox_volume, totalVolume(), TOLERANCE);
1603 }
1604 
1605 void
1607 {
1608  libmesh_parallel_only(comm());
1609 
1610  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1611 
1612  mooseAssert(!currentlyGenerating() && !currentlyPropagating(),
1613  "Cannot be reset during generation or propagation");
1614 
1615  for (THREAD_ID tid = 0; tid < libMesh::n_threads(); ++tid)
1617 }
1618 
1619 RayID
1621 {
1622  // Get the current ID to return
1623  const auto id = _threaded_next_ray_id[tid];
1624 
1625  // Advance so that the next call has the correct ID
1627 
1628  return id;
1629 }
1630 
1631 void
1633 {
1634  libmesh_parallel_only(comm());
1635 
1636  Threads::spin_mutex::scoped_lock lock(_spin_mutex);
1637 
1638  mooseAssert(!currentlyGenerating() && !currentlyPropagating(),
1639  "Cannot be reset during generation or propagation");
1640 
1642 }
1643 
1644 RayID
1646 {
1647  return _replicated_next_ray_id++;
1648 }
1649 
1650 bool
1652  const unsigned short side,
1653  const Point & direction,
1654  const THREAD_ID tid)
1655 {
1656  const auto & normal = getSideNormal(elem, side, tid);
1657  const auto dot = normal * direction;
1658  return dot < TraceRayTools::TRACE_TOLERANCE;
1659 }
1660 
1661 std::shared_ptr<Ray>
1663 {
1664  mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1665 
1666  return _parallel_ray_study->acquireParallelData(
1667  /* tid = */ 0,
1668  this,
1669  generateUniqueRayID(/* tid = */ 0),
1670  rayDataSize(),
1671  rayAuxDataSize(),
1672  /* reset = */ true,
1674 }
1675 
1676 std::shared_ptr<Ray>
1678 {
1679  mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1680 
1681  return _parallel_ray_study->acquireParallelData(/* tid = */ 0,
1682  this,
1683  generateUniqueRayID(/* tid = */ 0),
1684  /* data_size = */ 0,
1685  /* aux_data_size = */ 0,
1686  /* reset = */ true,
1688 }
1689 
1690 std::shared_ptr<Ray>
1692 {
1693  mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1694  libmesh_parallel_only(comm());
1695 
1696  return _parallel_ray_study->acquireParallelData(
1697  /* tid = */ 0,
1698  this,
1700  rayDataSize(),
1701  rayAuxDataSize(),
1702  /* reset = */ true,
1704 }
1705 
1706 std::shared_ptr<Ray>
1708 {
1709  mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1710 
1711  // Either register a Ray or get an already registered Ray id
1712  const RayID id = registerRay(name);
1713 
1714  // Acquire a Ray with the properly sized data initialized to zero
1715  return _parallel_ray_study->acquireParallelData(
1716  /* tid = */ 0,
1717  this,
1718  id,
1719  rayDataSize(),
1720  rayAuxDataSize(),
1721  /* reset = */ true,
1723 }
1724 
1725 std::shared_ptr<Ray>
1727 {
1728  mooseAssert(currentlyGenerating(), "Can only use during generateRays()");
1729  return _parallel_ray_study->acquireParallelData(
1730  /* tid = */ 0, &ray, Ray::ConstructRayKey());
1731 }
1732 
1733 std::shared_ptr<Ray>
1735 {
1736  mooseAssert(currentlyPropagating(), "Can only use during propagation");
1737  return _parallel_ray_study->acquireParallelData(tid,
1738  this,
1739  generateUniqueRayID(tid),
1740  rayDataSize(),
1741  rayAuxDataSize(),
1742  /* reset = */ true,
1744 }
bool hasRayKernels(const THREAD_ID tid)
Whether or not there are currently any active RayKernel objects.
bool currentlyGenerating() const
Whether or not the study is generating.
void getRayKernels(std::vector< RayKernelBase *> &result, SubdomainID id, THREAD_ID tid)
Fills the active RayKernels associated with this study and a block into result.
void internalSidesetSetup()
Does the setup for internal sidesets.
RayDataIndex registerRayData(const std::string &name)
Register a value to be filled in the data on a Ray with a given name.
RayData getBankedRayDataInternal(const RayID ray_id, const RayDataIndex index, const bool aux) const
Internal method for getting the value (replicated across all processors) in a Ray&#39;s data or aux data ...
virtual void postOnSegment(const THREAD_ID tid, const std::shared_ptr< Ray > &ray)
Called at the end of a Ray segment.
const bool _tolerate_failure
Whether or not to tolerate a Ray Tracing failure.
Real totalVolume() const
Get the current total volume of the domain.
void traceableMeshChecks()
Check for if all of the element types in the mesh are supported by ray tracing.
virtual void onCompleteRay(const std::shared_ptr< Ray > &ray)
Entry point for acting on a ray when it is completed (shouldContinue() == false)
libMesh::ConstElemRange * getActiveLocalElementRange()
const std::string & getRayDataNameInternal(const RayDataIndex index, const bool aux) const
Internal method for getting the name of Ray data or Ray aux data.
const std::vector< std::shared_ptr< Ray > > & rayBank() const
Get the Ray bank.
static const RayDataIndex INVALID_RAY_DATA_INDEX
Invalid index into a Ray&#39;s data.
Definition: Ray.h:200
std::vector< std::vector< std::set< const RayTracingObject * > > > _threaded_ray_object_registration
Threaded storage for all of the RayTracingObjects associated with a single Ray.
std::vector< std::string > & _reverse_registered_ray_map
Map from registered Ray ID to name.
unsigned long long int _total_intersections
Total number of Ray/element intersections.
RayID registerRay(const std::string &name)
Registers a Ray with a given name.
unsigned long long int _ending_intersections
Total number of Ray/element intersections for Rays that finished on this processor.
unsigned int _max_trajectory_changes
Max number of trajectory changes for a single Ray.
unsigned long int RayID
Type for a Ray&#39;s ID.
Definition: Ray.h:43
void verifyUniqueRayIDs(const std::vector< std::shared_ptr< Ray >>::const_iterator begin, const std::vector< std::shared_ptr< Ray >>::const_iterator end, const bool global, const std::string &error_suffix) const
Verifies that the Rays in the given range have unique Ray IDs.
unsigned int n_threads()
bool _has_same_level_active_elems
Whether or not the mesh has active elements of the same level.
auto norm() const -> decltype(std::norm(Real()))
bool sameLevelActiveElems() const
Determine whether or not the mesh currently has active elements that are all the same level...
bool absoluteFuzzyEqual(const T &var1, const T2 &var2, const T3 &tol=libMesh::TOLERANCE *libMesh::TOLERANCE)
virtual void initialSetup() override
void moveRayToBufferDuringTrace(std::shared_ptr< Ray > &ray, const THREAD_ID tid, const AcquireMoveDuringTraceKey &)
INTERNAL method for moving a Ray into the buffer during tracing.
RayDataIndex registerRayAuxData(const std::string &name)
Register a value to be filled in the aux data on a Ray with a given name.
Class that is used as a parameter to the public constructors/reset methods.
Definition: Ray.h:87
void paramError(const std::string &param, Args... args) const
MooseMesh & _mesh
The Mesh.
T & getMesh(MooseMesh &mesh)
function to cast mesh
Definition: SCM.h:35
Attribute for the RayTracingStudy a RayTracingObject is associated with.
static InputParameters validParams()
virtual void prepare(const Elem *elem, const THREAD_ID tid) override
NumericVector< Number > & solution()
std::vector< std::vector< TraceData > > _threaded_cached_traces
The threaded storage for cached traces.
const bool _use_internal_sidesets
Whether or not to use the internal sidesets in ray tracing.
Data structure that stores information for output of a partial trace of a Ray on a processor...
Definition: TraceData.h:42
QGAUSS
bool isTraceableElem(const Elem *elem)
bool isAdaptivityTraceableElem(const Elem *elem)
std::shared_ptr< Ray > acquireRay()
User APIs for constructing Rays within the RayTracingStudy.
std::size_t rayAuxDataSize() const
The registered size of values in the Ray aux data.
std::unordered_map< SubdomainID, Real > _subdomain_hmax
The cached hmax for all elements in a subdomain.
libMesh::dof_id_type maxIndex() const
Gets the maximum index generated using this object.
libMesh::BoundingBox _loose_b_box
Loose nodal bounding box for the domain.
T & set(const std::string &name, bool quiet_mode=false)
const std::set< BoundaryID > & getInternalSidesets() const
Gets the internal sidesets (that have RayBCs) within the local domain.
std::vector< std::vector< unsigned short > > _non_planar_sides
Non planar side data, which is for quick checking if an elem side is non-planar We use unsigned short...
RayData getBankedRayData(const RayID ray_id, const RayDataIndex index) const
Gets the data value for a banked ray with a given ID.
Base object for the RayKernel syntax.
Definition: RayKernelBase.h:27
if(subdm)
MeshBase & mesh
void initialize(const libMesh::SimpleRange< libMesh::MeshBase::element_iterator > elems)
Initializes the indices in a contiguous manner for the given element range.
Real _ending_distance
Total distance traveled by Rays that end on this processor.
virtual void reinitSegment(const Elem *elem, const Point &start, const Point &end, const Real length, THREAD_ID tid)
Reinitialize objects for a Ray segment for ray tracing.
unsigned int size() const
RayDataIndex registerRayDataInternal(const std::string &name, const bool aux)
Internal method for registering Ray data or Ray aux data with a name.
const std::vector< RayKernelBase * > & currentRayKernels(THREAD_ID tid) const
Gets the current RayKernels for a thread, which are set in segmentSubdomainSetup() ...
const Parallel::Communicator & comm() const
Key that is used for restricting access to moveRayToBufferDuringTrace() and acquireRayDuringTrace().
void addAvailableFlags(const ExecFlagType &flag, Args... flags)
std::shared_ptr< Ray > acquireReplicatedRay()
Acquire a Ray from the pool of Rays within generateRays() in a replicated fashion.
std::shared_ptr< Ray > acquireRayDuringTrace(const THREAD_ID tid, const AcquireMoveDuringTraceKey &)
INTERNAL methods for acquiring a Ray during a trace in RayKernels and RayBCs.
const Parallel::Communicator & _communicator
std::shared_ptr< Ray > acquireRegisteredRay(const std::string &name)
Acquires a Ray with a given name within generateRays().
void clearActiveMaterialProperties(const THREAD_ID tid)
The following methods are specializations for using the Parallel::packed_range_* routines for a vecto...
RayID _replicated_next_ray_id
Storage for the next available replicated RayID, obtained via generateReplicatedRayID() ...
std::set< BoundaryID > _internal_sidesets
The BoundaryIDs on the local mesh that have internal RayBCs.
virtual RayID generateUniqueRayID(const THREAD_ID tid)
Generates a unique RayID to be used for a Ray.
bool isOn()
std::vector< std::vector< RayKernelBase * > > _threaded_current_ray_kernels
The current RayKernel objects for each thread.
auto max(const L &left, const R &right)
std::vector< TraceData > _cached_traces
Storage for the cached traces.
RayDataIndex getRayAuxDataIndex(const std::string &name, const bool graceful=false) const
Gets the index associated with a registered value in the Ray aux data.
std::vector< std::shared_ptr< TraceRay > > _threaded_trace_ray
The TraceRay objects for each thread (they do the physical tracing)
virtual Assembly & assembly(const THREAD_ID tid, const unsigned int sys_num) override
Real subdomainHmax(const SubdomainID subdomain_id) const
Get the cached hmax for all elements in a subdomain.
const Real TRACE_TOLERANCE
The standard tolerance to use in tracing.
Definition: TraceRayTools.h:48
std::vector< RayID > _threaded_next_ray_id
Storage for the next available unique RayID, obtained via generateUniqueRayID()
bool isRectangularDomain() const
Whether or not the domain is rectangular (if it is prefectly encompassed by its bounding box) ...
std::vector< std::string > _ray_aux_data_names
The names for each Ray aux data entry.
virtual void setActiveElementalMooseVariables(const std::set< MooseVariableFEBase * > &moose_vars, const THREAD_ID tid) override
Threads::spin_mutex _spin_mutex
Spin mutex object for locks.
uint8_t processor_id_type
std::shared_ptr< Ray > getBankedRay(const RayID ray_id) const
Gets the Ray with the ID ray_id from the Ray bank.
unsigned long long int _total_processor_crossings
Total number of processor crossings.
processor_id_type n_processors() const
CONSTANT
void nonPlanarSideSetup()
Sets up the caching of whether or not each element side is non-planar, which is stored in _non_planar...
bool verifyRays() const
Whether or not to verify if Rays have valid information before being traced.
const std::string & name() const
unsigned int RayDataIndex
Type for the index into the data and aux data on a Ray.
Definition: Ray.h:51
void coverageChecks()
Perform coverage checks (coverage of RayMaterials and RayKernels, if enabled)
bool & useEigenvalue()
void resetReplicatedRayIDs()
Resets the generation of unique replicated RayIDs accessed via generateReplicatedRayID().
Real value(unsigned n, unsigned alpha, unsigned beta, Real x)
const bool _bank_rays_on_completion
Whether or not to bank rays on completion.
MeshBase & getMesh()
void min(const T &r, T &o, Request &req) const
unsigned int _ending_max_intersections
Max number of intersections for Rays that finished on this processor.
sideset clear()
TheWarehouse & theWarehouse() const
virtual void postExecuteStudy()
Entry point after study execution.
MooseVariableFE< Real > & variable()
Gets the variable this AuxRayKernel contributes to.
Definition: AuxRayKernel.h:38
unsigned int _max_processor_crossings
Max number of processor crossings for all Rays.
const std::unique_ptr< ParallelRayStudy > _parallel_ray_study
The study that used is to actually execute (trace) the Rays.
virtual unsigned int dimension() const
std::vector< RayTracingObject * > getRayTracingObjects()
Gets all of the currently active RayTracingObjects.
virtual void residualSetup() override
const Point & min() const
const bool & currentlyComputingResidual() const
boundary_id_type BoundaryID
const std::string name
Definition: Setup.h:20
std::vector< TheWarehouse::QueryCache< AttribSubdomains > > _threaded_cache_ray_kernel
Threaded cached subdomain query for RayKernelBase objects pertaining to this study.
virtual void execute() override
Executes the study (generates and propagates Rays)
SimpleRange< IndexType > as_range(const std::pair< IndexType, IndexType > &p)
unsigned int _max_intersections
Max number of intersections for a single Ray.
unsigned int _ending_max_trajectory_changes
Max number of trajectory changes for Rays that finished on this processor.
unsigned int _ending_max_processor_crossings
Max number of total processor crossings for Rays that finished on this processor. ...
const std::string & getRayDataName(const RayDataIndex index) const
Gets the name associated with a registered value in the Ray data.
void reserveRayBuffer(const std::size_t size)
Reserve size entires in the Ray buffer.
void zeroAuxVariables()
Zero the AuxVariables that the registered AuxRayKernels contribute to.
std::vector< std::unique_ptr< libMesh::FEBase > > _threaded_fe_face
Face FE used for computing face normals for each thread.
const std::string & type() const
void dependencyChecks()
Perform checks to see if the listed dependencies in the RayTracingObjects exist.
static InputParameters validParams()
const std::string & getRayAuxDataName(const RayDataIndex index) const
Gets the name associated with a registered value in the Ray aux data.
virtual void cacheResidual(const THREAD_ID tid) override
Basic datastructure for a ray that will traverse the mesh.
Definition: Ray.h:56
std::chrono::steady_clock::time_point _execution_start_time
Timing.
std::vector< TheWarehouse::QueryCache< AttribBoundaries > > _threaded_cache_ray_bc
Threaded cached boundary query for RayBC objects pertaining to this study.
SystemBase & _sys
virtual void jacobianSetup() override
std::string typeAndName() const
bool isValueSet(const std::string &value) const
unsigned int number() const
auto norm(const T &a) -> decltype(std::abs(a))
std::string stringify(const T &t)
virtual void preExecuteStudy()
Entry point before study execution.
std::unordered_map< std::string, RayDataIndex > _ray_aux_data_map
The map from Ray aux data names to index.
AuxiliarySystem & getAuxiliarySystem()
Real volume(const MeshBase &mesh, unsigned int dim=libMesh::invalid_uint)
const bool _warn_non_planar
Whether not to warn if non-planar faces are found.
libMesh::BoundingBox _b_box
Nodal bounding box for the domain.
virtual void meshChanged() override
libMesh::dof_id_type getIndex(const libMesh::Elem *elem) const
Get the index associated with the element elem.
virtual void setCurrentSubdomainID(const Elem *elem, const THREAD_ID tid) override
virtual void generateRays()=0
Subclasses should override this to determine how to generate Rays.
RayID generateReplicatedRayID()
Generates a Ray ID that is replicated across all processors.
virtual void close()=0
MONOMIAL
const ExecFlagEnum & _execute_enum
RayDataIndex getRayDataIndexInternal(const std::string &name, const bool aux, const bool graceful) const
Internal method for getting the index of Ray data or Ray aux data.
void verifyDependenciesExist(const std::vector< RayTracingObject *> &rtos)
Verifies that the dependencies exist for a set of RayTracingObjects.
unsigned long long int _ending_processor_crossings
Total number of processor crossings for Rays that finished on this processor.
query_obj query
ElemIndexHelper _elem_index_helper
Helper for defining a local contiguous index for each element.
RayTracingStudy(const InputParameters &parameters)
virtual unsigned int n_sides() const=0
std::vector< std::unordered_map< std::pair< const Elem *, unsigned short >, Point > > _threaded_cached_normals
Threaded cache for side normals that have been computed already during tracing.
TraceData & initThreadedCachedTrace(const std::shared_ptr< Ray > &ray, THREAD_ID tid)
Initialize a Ray in the threaded cached trace map to be filled with segments.
virtual const SystemBase & getSystemBase(const unsigned int sys_num) const
const Elem * neighbor_ptr(unsigned int i) const
virtual const Point & getSideNormal(const Elem *elem, const unsigned short side, const THREAD_ID tid)
Get the outward normal for a given element side.
const SubdomainID ANY_BLOCK_ID
unsigned int level() const
std::vector< RayDataIndex > getRayAuxDataIndices(const std::vector< std::string > &names, const bool graceful=false) const
Gets the indices associated with registered values in the Ray aux data.
std::unordered_map< std::string, RayDataIndex > _ray_data_map
The map from Ray data names to index.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::vector< RayDataIndex > getRayDataIndices(const std::vector< std::string > &names, const bool graceful=false) const
Gets the indices associated with registered values in the Ray data.
bool currentlyPropagating() const
Whether or not the study is propagating (tracing Rays)
subdomain_id_type subdomain_id() const
bool _has_non_planar_sides
Whether or not the local mesh has elements with non-planar sides.
virtual void subdomainSetup(SubdomainID subdomain, const THREAD_ID tid)
const libMesh::Elem & elemSide(const libMesh::Elem &elem, const unsigned int s, const THREAD_ID tid=0)
Get an element&#39;s side pointer without excessive memory allocation.
void max(const T &r, T &o, Request &req) const
void reinitMaterials(SubdomainID blk_id, const THREAD_ID tid, bool swap_stateful=true)
FEProblemBase & _fe_problem
void getRayBCs(std::vector< RayBoundaryConditionBase *> &result, BoundaryID id, THREAD_ID tid)
Fills the active RayBCs associated with this study and a boundary into result.
std::vector< std::size_t > _num_cached
Number of currently cached objects for Jacobian/residual for each thread.
std::chrono::steady_clock::duration _propagation_time
Query query()
virtual void segmentSubdomainSetup(const SubdomainID subdomain, const THREAD_ID tid, const RayID ray_id)
Setup for on subdomain change or subdomain AND ray change during ray tracing.
void mooseWarning(Args &&... args) const
RayData getBankedRayAuxData(const RayID ray_id, const RayDataIndex index) const
Gets the data value for a banked ray with a given ID.
bool _called_initial_setup
Whether or not we&#39;ve called initial setup - used to stop from late registration.
void mooseError(Args &&... args) const
void verifyUniqueRays(const std::vector< std::shared_ptr< Ray >>::const_iterator begin, const std::vector< std::shared_ptr< Ray >>::const_iterator end, const std::string &error_suffix)
Verifies that the Rays in the given range are unique.
std::vector< RayDataIndex > getRayDataIndicesInternal(const std::vector< std::string > &names, const bool aux, const bool graceful) const
Internal method for getting the indicies of Ray data or Ray aux data.
Real _total_distance
Total distance traveled by all Rays.
const processor_id_type _pid
The rank of this processor (this actually takes time to lookup - so just do it once) ...
std::unordered_map< std::string, RayID > & _registered_ray_map
Map from registered Ray name to ID.
std::vector< std::string > _ray_data_names
The names for each Ray data entry.
float RayData
Type for a Ray&#39;s data.
Definition: Ray.h:46
const std::string & registeredRayName(const RayID ray_id) const
Gets the name of a registered ray.
void moveRayToBuffer(std::shared_ptr< Ray > &ray)
Moves a ray to the buffer to be traced during generateRays().
const ExecFlagType EXEC_PRE_KERNELS
RayID registeredRayID(const std::string &name, const bool graceful=false) const
Gets the ID of a registered ray.
virtual void buildSegmentQuadrature(const Point &start, const Point &end, const Real length, std::vector< Point > &points, std::vector< Real > &weights) const
Builds quadrature points for a given segment using the _segment_qrule.
libMesh::StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement *> * getBoundaryElementRange()
const bool _use_ray_registration
Whether or not to use Ray registration.
virtual bool hasActiveElementalMooseVariables(const THREAD_ID tid) const
const bool _ray_kernel_coverage_check
Whether or not to perform coverage checks on RayKernels.
void localElemIndexSetup()
Sets up the _elem_index_helper, which is used for obtaining a contiguous index for all elements that ...
const bool & currentlyComputingJacobian() const
std::vector< std::shared_ptr< Ray > > _ray_bank
Cumulative Ray bank - stored only when _bank_rays_on_completion.
libMesh::BoundingBox create_nodal_bounding_box(const MeshBase &mesh)
const Point & max() const
std::vector< unsigned long long int > _local_trace_ray_results
Cumulative results on this processor from the threaded TraceRay objects.
std::chrono::steady_clock::duration _generation_time
void registeredRaySetup()
Sets up the maps from Ray to associated RayTracingObjects if _use_ray_registration.
void scale(const Real factor)
void moveRaysToBuffer(std::vector< std::shared_ptr< Ray >> &rays)
Moves rays to the buffer to be traced during generateRays().
void executeStudy()
Method for executing the study so that it can be called out of the standard UO execute() ...
void prepareMaterials(const std::unordered_set< unsigned int > &consumer_needed_mat_props, const SubdomainID blk_id, const THREAD_ID tid)
Adaptivity & adaptivity()
bool hasActiveMaterialProperties(const THREAD_ID tid) const
virtual void cacheJacobian(const THREAD_ID tid) override
auto min(const L &left, const R &right)
virtual void reinitElemPhys(const Elem *elem, const std::vector< Point > &phys_points_in_elem, const THREAD_ID tid) override
std::shared_ptr< Ray > acquireUnsizedRay()
Acquire a Ray from the pool of Rays within generateRays(), without resizing the data (sizes the data ...
RayDataIndex getRayDataIndex(const std::string &name, const bool graceful=false) const
Gets the index associated with a registered value in the Ray data.
void resetUniqueRayIDs()
Resets the generation of unique RayIDs via generateUniqueRayID() to the beginning of the range...
virtual void addCachedResidual(const THREAD_ID tid) override
std::chrono::steady_clock::duration _execution_time
virtual void zeroVariables(std::vector< std::string > &vars_to_be_zeroed)
virtual void clearActiveElementalMooseVariables(const THREAD_ID tid) override
std::vector< std::unique_ptr< libMesh::QBase > > _threaded_q_face
Face quadrature used for computing face normals for each thread.
std::shared_ptr< Ray > acquireCopiedRay(const Ray &ray)
Acquires a Ray that that is copied from another Ray within generateRays().
Traces Rays through the mesh on a single processor.
Definition: TraceRay.h:46
std::vector< std::vector< std::vector< BoundaryID > > > _internal_sidesets_map
Internal sideset data, if internal sidesets exist (indexed with getLocalElemIndex()) ...
void subdomainHMaxSetup()
Caches the hmax for all elements in each subdomain.
bool sideIsIncoming(const Elem *const elem, const unsigned short side, const Point &direction, const THREAD_ID tid)
Whether or not side is incoming on element elem in direction direction.
static const RayID INVALID_RAY_ID
Invalid Ray ID.
Definition: Ray.h:202
std::unique_ptr< libMesh::QBase > _segment_qrule
Quadrature rule for laying points across a 1D ray segment.
bool hasInternalSidesets() const
Whether or not the local mesh has internal sidesets that have RayBCs on them.
std::size_t rayDataSize() const
The registered size of values in the Ray data.
unsigned int THREAD_ID
virtual void timestepSetup() override
const std::set< SubdomainID > & meshSubdomains() const
virtual void addCachedJacobian(const THREAD_ID tid) override
virtual bool shouldCacheTrace(const std::shared_ptr< Ray > &) const
Virtual that allows for selection in if a Ray should be cached or not (only used when _cache_traces)...
const RemoteElem * remote_elem
Real computeTotalVolume()
Helper function for computing the total domain volume.