RayTracingStudy.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "GeneralUserObject.h"

// Local Includes
#include "ElemIndexHelper.h"
#include "ParallelRayStudy.h"
#include "RayTracingAttributes.h"
#include "TraceData.h"
#include "TraceRayBndElement.h"

// MOOSE Includes
#include "TheWarehouse.h"

// libMesh includes
#include "libmesh/mesh.h"
#include "libmesh/elem_side_builder.h"

// Forward Declarations
class RayBoundaryConditionBase;
class RayKernelBase;
class TraceRay;
class RayTracingObject;

class RayTracingStudy : public GeneralUserObject
{
public:
  RayTracingStudy(const InputParameters & parameters);

  class AcquireRayInternalKey
  {
    friend class Parallel::Packing<std::shared_ptr<Ray>>;
    friend void dataLoad(std::istream & stream, std::shared_ptr<Ray> & ray, void * context);
    AcquireRayInternalKey() {}
    AcquireRayInternalKey(const AcquireRayInternalKey &) {}
  };

  class AcquireMoveDuringTraceKey
  {
    friend class RayBoundaryConditionBase;
    friend class RayKernelBase;
    AcquireMoveDuringTraceKey() {}
    AcquireMoveDuringTraceKey(const AcquireMoveDuringTraceKey &) {}
  };

  static InputParameters validParams();

  virtual void initialSetup() override;
  virtual void residualSetup() override;
  virtual void jacobianSetup() override;
  virtual void meshChanged() override;
  virtual void timestepSetup() override;
  virtual void initialize() override {}
  virtual void finalize() override {}

  virtual void execute() override;

  virtual void
  segmentSubdomainSetup(const SubdomainID subdomain, const THREAD_ID tid, const RayID ray_id);

  virtual void reinitSegment(
      const Elem * elem, const Point & start, const Point & end, const Real length, THREAD_ID tid);

  virtual void postOnSegment(const THREAD_ID tid, const std::shared_ptr<Ray> & ray);

  virtual void preTrace(const THREAD_ID /* tid */, const std::shared_ptr<Ray> & /* ray */) {}

  void executeStudy();

  unsigned long long int endingProcessorCrossings() const { return _ending_processor_crossings; }
  unsigned int endingMaxProcessorCrossings() const { return _ending_max_processor_crossings; }
  unsigned long long int totalProcessorCrossings() const { return _total_processor_crossings; }
  unsigned int maxProcessorCrossings() const { return _max_processor_crossings; }

  unsigned long long int endingIntersections() const { return _ending_intersections; }
  unsigned int endingMaxIntersections() const { return _ending_max_intersections; }
  unsigned long long int totalIntersections() const { return _total_intersections; }
  unsigned int maxIntersections() const { return _max_intersections; }
  unsigned int maxTrajectoryChanges() const { return _max_trajectory_changes; }

  Real endingDistance() const { return _ending_distance; }
  Real totalDistance() const { return _total_distance; }

  unsigned long long int localTraceRayResult(const int result) const
  {
    return _local_trace_ray_results[result];
  }

  const ParallelRayStudy & parallelRayStudy() const { return *_parallel_ray_study; }

  Real rayMaxDistance() const { return _ray_max_distance; }

  Real executionTime() { return std::chrono::duration<Real>(_execution_time).count(); }
  Real executionTimeNano()
  {
    return std::chrono::duration<Real, std::nano>(_execution_time).count();
  }
  Real generationTime() const { return std::chrono::duration<Real>(_generation_time).count(); }
  Real propagationTime() const { return std::chrono::duration<Real>(_propagation_time).count(); }

  bool tolerateFailure() const { return _tolerate_failure; }

  bool bankRaysOnCompletion() const { return _bank_rays_on_completion; }

  bool rayDependentSubdomainSetup() const { return _ray_dependent_subdomain_setup; }

  RayDataIndex registerRayData(const std::string & name);
  std::vector<RayDataIndex> registerRayData(const std::vector<std::string> & names);
  RayDataIndex getRayDataIndex(const std::string & name, const bool graceful = false) const;
  std::vector<RayDataIndex> getRayDataIndices(const std::vector<std::string> & names,
                                              const bool graceful = false) const;
  const std::string & getRayDataName(const RayDataIndex index) const;
  std::vector<std::string> getRayDataNames(const std::vector<RayDataIndex> & indices) const;

  std::size_t rayDataSize() const { return _ray_data_names.size(); }
  bool hasRayData() const { return _ray_data_names.size(); }
  const std::vector<std::string> & rayDataNames() const { return _ray_data_names; }

  RayDataIndex registerRayAuxData(const std::string & name);
  std::vector<RayDataIndex> registerRayAuxData(const std::vector<std::string> & names);
  RayDataIndex getRayAuxDataIndex(const std::string & name, const bool graceful = false) const;
  std::vector<RayDataIndex> getRayAuxDataIndices(const std::vector<std::string> & names,
                                                 const bool graceful = false) const;
  const std::string & getRayAuxDataName(const RayDataIndex index) const;
  std::vector<std::string> getRayAuxDataNames(const std::vector<RayDataIndex> & indices) const;
  std::size_t rayAuxDataSize() const { return _ray_aux_data_names.size(); }
  bool hasRayAuxData() const { return _ray_aux_data_names.size(); }
  const std::vector<std::string> & rayAuxDataNames() const { return _ray_aux_data_names; }

  bool hasRayKernels(const THREAD_ID tid);
  void getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid);
  template <typename T>
  void getRayKernels(std::vector<T *> & result, THREAD_ID tid)
  {
    _fe_problem.theWarehouse()
        .query()
        .condition<AttribRayTracingStudy>(this)
        .condition<AttribSystem>("RayKernel")
        .condition<AttribThread>(tid)
        .queryInto(result);
  }
  void
  getRayKernels(std::vector<RayKernelBase *> & result, SubdomainID id, THREAD_ID tid, RayID ray_id);
  void getRayBCs(std::vector<RayBoundaryConditionBase *> & result, BoundaryID id, THREAD_ID tid);
  template <typename T>
  void getRayBCs(std::vector<T *> & result, const std::vector<BoundaryID> & ids, THREAD_ID tid)
  {
    _fe_problem.theWarehouse()
        .query()
        .condition<AttribRayTracingStudy>(this)
        .condition<AttribSystem>("RayBoundaryCondition")
        .condition<AttribBoundaries>(ids)
        .condition<AttribThread>(tid)
        .queryInto(result);
  }
  template <typename T>
  void getRayBCs(std::vector<T *> & result, THREAD_ID tid)
  {
    _fe_problem.theWarehouse()
        .query()
        .condition<AttribRayTracingStudy>(this)
        .condition<AttribSystem>("RayBoundaryCondition")
        .condition<AttribThread>(tid)
        .queryInto(result);
  }
  virtual void getRayBCs(std::vector<RayBoundaryConditionBase *> & result,
                         const std::vector<TraceRayBndElement> & bnd_elems,
                         THREAD_ID tid,
                         RayID ray_id);

  const std::vector<RayKernelBase *> & currentRayKernels(THREAD_ID tid) const
  {
    return _threaded_current_ray_kernels[tid];
  }

  const BoundingBox & boundingBox() const { return _b_box; }
  const BoundingBox & looseBoundingBox() const { return _loose_b_box; }
  Real domainMaxLength() const { return _domain_max_length; }
  Real totalVolume() const { return _total_volume; }
  bool isRectangularDomain() const;

  bool hasInternalSidesets() const { return _internal_sidesets.size(); }
  const std::vector<std::vector<BoundaryID>> & getInternalSidesets(const Elem * elem) const;
  const std::set<BoundaryID> & getInternalSidesets() const { return _internal_sidesets; }
  bool sideIsNonPlanar(const Elem * elem, const unsigned short s) const
  {
    return _has_non_planar_sides && _non_planar_sides[_elem_index_helper.getIndex(elem)][s] == 1;
  }
  bool hasSameLevelActiveElems() const { return _has_same_level_active_elems; }

  std::shared_ptr<Ray> acquireRayDuringTrace(const THREAD_ID tid,
                                             const AcquireMoveDuringTraceKey &);
  std::shared_ptr<Ray> acquireRayInternal(const RayID id,
                                          const std::size_t data_size,
                                          const std::size_t aux_data_size,
                                          const bool reset,
                                          const AcquireRayInternalKey &)
  {
    return _parallel_ray_study->acquireParallelData(
        0, this, id, data_size, aux_data_size, reset, Ray::ConstructRayKey());
  }

  void moveRayToBufferDuringTrace(std::shared_ptr<Ray> & ray,
                                  const THREAD_ID tid,
                                  const AcquireMoveDuringTraceKey &);
  MeshBase & meshBase() const { return _mesh; }

  MooseMesh & mesh() { return _mesh; }

  virtual const Point &
  getSideNormal(const Elem * elem, const unsigned short side, const THREAD_ID tid);
  virtual const Point * getElemNormals(const Elem * /* elem */, const THREAD_ID /* tid */)
  {
    mooseError("Unimplemented element normal caching in ", type(), "::getElemNormals()");
  }

  RayData getBankedRayData(const RayID ray_id, const RayDataIndex index) const;
  RayData getBankedRayAuxData(const RayID ray_id, const RayDataIndex index) const;

  RayID registeredRayID(const std::string & name, const bool graceful = false) const;
  const std::string & registeredRayName(const RayID ray_id) const;

  bool useRayRegistration() const { return _use_ray_registration; }

  bool dataOnCacheTraces() const { return _data_on_cache_traces; }
  bool auxDataOnCacheTraces() const { return _aux_data_on_cache_traces; }
  bool segmentsOnCacheTraces() const { return _segments_on_cache_traces; }

  virtual bool shouldCacheTrace(const std::shared_ptr<Ray> & /* ray */) const
  {
    return _always_cache_traces;
  }

  TraceData & initThreadedCachedTrace(const std::shared_ptr<Ray> & ray, THREAD_ID tid);
  const std::vector<TraceData> & getCachedTraces() const { return _cached_traces; }

  Real subdomainHmax(const SubdomainID subdomain_id) const;

  virtual void onCompleteRay(const std::shared_ptr<Ray> & ray);

  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;

  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);

  bool currentlyPropagating() const { return _parallel_ray_study->currentlyExecuting(); }
  bool currentlyGenerating() const { return _parallel_ray_study->currentlyPreExecuting(); }

  TraceRay & traceRay(const THREAD_ID tid) { return *_threaded_trace_ray[tid]; }
  const TraceRay & traceRay(const THREAD_ID tid) const { return *_threaded_trace_ray[tid]; }

  bool verifyRays() const { return _verify_rays; }
#ifndef NDEBUG
  bool verifyTraceIntersections() const { return _verify_trace_intersections; }
#endif

  bool sideIsIncoming(const Elem * const elem,
                      const unsigned short side,
                      const Point & direction,
                      const THREAD_ID tid);

  bool warnNonPlanar() const { return _warn_non_planar; }

  ParallelStudy<std::shared_ptr<Ray>, Ray> * parallelStudy() { return _parallel_ray_study.get(); }

  const libMesh::Elem &
  elemSide(const libMesh::Elem & elem, const unsigned int s, const THREAD_ID tid = 0)
  {
    return _threaded_elem_side_builders[tid](elem, s);
  }

protected:
  virtual void generateRays() = 0;

  virtual void preExecuteStudy() {}
  virtual void postExecuteStudy() {}

  Real computeTotalVolume();

  std::vector<RayKernelBase *> & currentRayKernelsWrite(THREAD_ID tid)
  {
    return _threaded_current_ray_kernels[tid];
  }

  void reserveRayBuffer(const std::size_t size);

  bool sameLevelActiveElems() const;

  virtual void buildSegmentQuadrature(const Point & start,
                                      const Point & end,
                                      const Real length,
                                      std::vector<Point> & points,
                                      std::vector<Real> & weights) const;

  const std::vector<std::shared_ptr<Ray>> & rayBank() const;

  std::shared_ptr<Ray> getBankedRay(const RayID ray_id) const;

  void resetUniqueRayIDs();
  void resetReplicatedRayIDs();

  std::vector<RayTracingObject *> getRayTracingObjects();

  virtual RayID generateUniqueRayID(const THREAD_ID tid);
  RayID generateReplicatedRayID();


  std::shared_ptr<Ray> acquireRay();
  std::shared_ptr<Ray> acquireUnsizedRay();
  std::shared_ptr<Ray> acquireReplicatedRay();
  std::shared_ptr<Ray> acquireCopiedRay(const Ray & ray);
  std::shared_ptr<Ray> acquireRegisteredRay(const std::string & name);

  void moveRayToBuffer(std::shared_ptr<Ray> & ray);
  void moveRaysToBuffer(std::vector<std::shared_ptr<Ray>> & rays);

  MooseMesh & _mesh;
  const Parallel::Communicator & _comm;
  const processor_id_type _pid;

  const bool _ray_kernel_coverage_check;
  const bool _warn_non_planar;
  const bool _use_ray_registration;
  const bool _use_internal_sidesets;
  const bool _tolerate_failure;
  const bool _bank_rays_on_completion;
  const bool _ray_dependent_subdomain_setup;

  bool _always_cache_traces;
  const bool _data_on_cache_traces;
  const bool _aux_data_on_cache_traces;
  const bool _segments_on_cache_traces;
  const Real _ray_max_distance;
  const bool _verify_rays;
#ifndef NDEBUG
  const bool _verify_trace_intersections;
#endif

private:
  void coverageChecks();

  void dependencyChecks();
  void verifyDependenciesExist(const std::vector<RayTracingObject *> & rtos);

  void traceableMeshChecks();

  void internalSidesetSetup();

  void nonPlanarSideSetup();

  void localElemIndexSetup();

  void registeredRaySetup();

  void zeroAuxVariables();

  void subdomainHMaxSetup();

  RayDataIndex registerRayDataInternal(const std::string & name, const bool aux);
  std::vector<RayDataIndex> registerRayDataInternal(const std::vector<std::string> & names,
                                                    const bool aux);
  RayDataIndex
  getRayDataIndexInternal(const std::string & name, const bool aux, const bool graceful) const;
  std::vector<RayDataIndex> getRayDataIndicesInternal(const std::vector<std::string> & names,
                                                      const bool aux,
                                                      const bool graceful) const;

  const std::string & getRayDataNameInternal(const RayDataIndex index, const bool aux) const;

  RayData
  getBankedRayDataInternal(const RayID ray_id, const RayDataIndex index, const bool aux) const;

  RayID registerRay(const std::string & name);

  std::vector<libMesh::ElemSideBuilder> _threaded_elem_side_builders;


  std::chrono::steady_clock::time_point _execution_start_time;
  std::chrono::steady_clock::duration _execution_time;
  std::chrono::steady_clock::duration _generation_time;
  std::chrono::steady_clock::duration _propagation_time;

  std::unordered_map<std::string, RayDataIndex> _ray_data_map;
  std::unordered_map<std::string, RayDataIndex> _ray_aux_data_map;

  std::vector<std::string> _ray_data_names;
  std::vector<std::string> _ray_aux_data_names;

  std::unordered_map<std::string, RayID> & _registered_ray_map;
  std::vector<std::string> & _reverse_registered_ray_map;

  std::vector<TraceData> _cached_traces;
  std::vector<std::vector<TraceData>> _threaded_cached_traces;

  std::vector<std::size_t> _num_cached;

  std::set<BoundaryID> _internal_sidesets;
  std::vector<std::vector<std::vector<BoundaryID>>> _internal_sidesets_map;
  bool _has_non_planar_sides;
  std::vector<std::vector<unsigned short>> _non_planar_sides;
  bool _has_same_level_active_elems;

  libMesh::BoundingBox _b_box;
  libMesh::BoundingBox _loose_b_box;
  Real _domain_max_length;
  Real _total_volume;

  std::vector<TheWarehouse::QueryCache<AttribSubdomains>> _threaded_cache_ray_kernel;
  std::vector<TheWarehouse::QueryCache<AttribBoundaries>> _threaded_cache_ray_bc;
  std::vector<std::vector<std::set<const RayTracingObject *>>> _threaded_ray_object_registration;
  std::vector<std::vector<RayKernelBase *>> _threaded_current_ray_kernels;
  std::vector<std::shared_ptr<TraceRay>> _threaded_trace_ray;
  std::vector<std::unique_ptr<libMesh::FEBase>> _threaded_fe_face;
  std::vector<std::unique_ptr<libMesh::QBase>> _threaded_q_face;
  std::vector<std::unordered_map<std::pair<const Elem *, unsigned short>, Point>>
      _threaded_cached_normals;
  std::vector<std::shared_ptr<Ray>> _ray_bank;
  std::vector<RayID> _threaded_next_ray_id;
  RayID _replicated_next_ray_id;

  const std::unique_ptr<ParallelRayStudy> _parallel_ray_study;

  std::unique_ptr<libMesh::QBase> _segment_qrule;

  unsigned long long int _ending_processor_crossings;
  unsigned int _ending_max_processor_crossings;
  unsigned long long int _total_processor_crossings;
  unsigned int _max_processor_crossings;

  unsigned long long int _ending_intersections;
  unsigned int _ending_max_intersections;
  unsigned int _ending_max_trajectory_changes;
  unsigned long long int _total_intersections;
  unsigned int _max_intersections;
  unsigned int _max_trajectory_changes;

  Real _ending_distance;
  Real _total_distance;

  std::vector<unsigned long long int> _local_trace_ray_results;

  std::unordered_map<SubdomainID, Real> _subdomain_hmax;

  bool _called_initial_setup;

  ElemIndexHelper _elem_index_helper;

  mutable Threads::spin_mutex _spin_mutex;
};