Add optical WLS secondary generation (#2339)

* Add WLS optical photon secondary generation

* Resample WLS polarization while not orthogonal and use more robust sampling method to avoid extremely high rejection rate

See also #2197 (comment)

* Fix missing update to core state counters num_pending in WLS generator action

* Save all generator accum stats in standalone optical runner

* Don't emit secondary photons if the incident photon energy is below the CDF sampling grid

* Add WLS generator test and geometry with WLS properties

* Clear distribution data buffer before invoking WLS interactor

* clang-tidy and copilot feedback

* Address initial feedback

* Address feedback

* Add WLS generator test

* Skip num_generated check with float build

---------

Co-authored-by: Seth R. Johnson <johnsonsr@ornl.gov>

Author: amandalund

src/celeritas/CMakeLists.txt

@@ -183,6 +183,7 @@ endif()
 
 if(CELERITAS_USE_CUDA OR CELERITAS_USE_HIP)
   list(APPEND SOURCES
+    optical/detail/Filler.cu
     track/detail/Filler.cu
     alongstep/detail/AlongStepKernels.cu
   )
@@ -407,6 +408,7 @@ celeritas_polysource(optical/gen/OffloadAction)
 celeritas_polysource(optical/gen/OffloadGatherAction)
 celeritas_polysource(optical/gen/PrimaryGeneratorAction)
 celeritas_polysource(optical/gen/DirectGeneratorAction)
+celeritas_polysource(optical/gen/WlsGeneratorAction)
 celeritas_polysource(optical/gen/detail/GeneratorAlgorithms)
 celeritas_polysource(optical/gen/detail/OffloadAlgorithms)
 celeritas_polysource(optical/surface/BoundaryAction)

src/celeritas/optical/CoreParams.cc

@@ -131,16 +131,13 @@ CoreParams::CoreParams(Input&& input) : input_(std::move(input))
     CP_VALIDATE_INPUT(surface);
     CP_VALIDATE_INPUT(detectors);
     CP_VALIDATE_INPUT(action_reg);
+    CP_VALIDATE_INPUT(aux_reg);
     CP_VALIDATE_INPUT(gen_reg);
     CP_VALIDATE_INPUT(max_streams);
 #undef CP_VALIDATE_INPUT
 
     CELER_EXPECT(input_);
 
-    if (!input_.aux_reg)
-    {
-        input_.aux_reg = std::make_shared<AuxParamsRegistry>();
-    }
     if (!input_.output_reg)
     {
         input_.output_reg = std::make_shared<OutputRegistry>();

src/celeritas/optical/CoreParams.hh

@@ -76,7 +76,7 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
         SPActionRegistry action_reg;
         SPOutputRegistry output_reg;
         SPGeneratorRegistry gen_reg;
-        SPAuxRegistry aux_reg;  //!< Optional, empty default
+        SPAuxRegistry aux_reg;
 
         // Problem definition and state
         SPConstCoreGeo geometry;
@@ -103,9 +103,9 @@ class CoreParams final : public ParamsDataInterface<CoreParamsData>
         explicit operator bool() const
         {
             return geometry && material && rng && sim && volume && surface
-                   && surface_physics && action_reg && gen_reg && max_streams
-                   && capacity.generators > 0 && capacity.tracks > 0
-                   && capacity.primaries > 0;
+                   && surface_physics && action_reg && aux_reg && gen_reg
+                   && max_streams && capacity.generators > 0
+                   && capacity.tracks > 0 && capacity.primaries > 0;
         }
     };
 

src/celeritas/optical/Interaction.hh

@@ -37,7 +37,6 @@ struct Interaction
     Real3 direction;  //!< Post-interaction direction
     Real3 polarization;  //!< Post-interaction polarization
     Action action{Action::scattered};  //!< Flags for interaction result
-    WlsDistributionData distribution;  //!< Data for generating WLS secondaries
 
     //! Return an interaction representing an absorbed process
     static inline CELER_FUNCTION Interaction from_absorption();

src/celeritas/optical/PhysicsParams.cc

@@ -6,13 +6,16 @@
 //---------------------------------------------------------------------------//
 #include "PhysicsParams.hh"
 
+#include "corecel/data/AuxParamsRegistry.hh"
 #include "corecel/io/Logger.hh"
 #include "corecel/sys/ActionRegistry.hh"
+#include "celeritas/phys/GeneratorRegistry.hh"
 
 #include "MaterialParams.hh"
 #include "MfpBuilder.hh"
 #include "Model.hh"
 #include "action/DiscreteSelectAction.hh"
+#include "gen/WlsGeneratorAction.hh"
 #include "model/AbsorptionModel.hh"
 #include "model/MieModel.hh"
 #include "model/RayleighModel.hh"
@@ -36,11 +39,16 @@ namespace optical
 PhysicsParams::PhysicsParams(inp::OpticalBulkPhysics const& input,
                              SPConstMaterials const& materials,
                              SPConstCoreMaterials const& core_materials,
-                             SPActionRegistry action_reg)
+                             SPActionRegistry action_reg,
+                             SPAuxRegistry aux_reg,
+                             SPGeneratorRegistry gen_reg,
+                             size_type gen_capacity)
 {
     CELER_EXPECT(materials);
     CELER_EXPECT(core_materials);
     CELER_EXPECT(action_reg);
+    CELER_EXPECT(aux_reg);
+    CELER_EXPECT(gen_reg);
 
     CELER_VALIDATE(input, << "no optical bulk models are present");
 
@@ -52,8 +60,13 @@ PhysicsParams::PhysicsParams(inp::OpticalBulkPhysics const& input,
         action_reg->insert(discrete_select_);
 
         // Build models
-        models_ = this->build_models(
-            input, materials, core_materials, *action_reg);
+        models_ = this->build_models(input,
+                                     materials,
+                                     core_materials,
+                                     *action_reg,
+                                     *aux_reg,
+                                     *gen_reg,
+                                     gen_capacity);
     }
 
     // Construct data
@@ -72,11 +85,17 @@ PhysicsParams::PhysicsParams(inp::OpticalBulkPhysics const& input,
 //---------------------------------------------------------------------------//
 /*!
  * Construct optical models and register them in the given registry.
+ *
+ * If WLS is enabled, a generator action is constructed for creating secondary
+ * photons.
  */
 auto PhysicsParams::build_models(inp::OpticalBulkPhysics const& input,
                                  SPConstMaterials materials,
                                  SPConstCoreMaterials core_materials,
-                                 ActionRegistry& action_reg) const -> VecModels
+                                 ActionRegistry& action_reg,
+                                 AuxParamsRegistry& aux_reg,
+                                 GeneratorRegistry& gen_reg,
+                                 size_type gen_capacity) const -> VecModels
 {
     VecModels models;
     if (input.absorption)
@@ -97,17 +116,44 @@ auto PhysicsParams::build_models(inp::OpticalBulkPhysics const& input,
             action_reg.next_id(), input.rayleigh, materials, core_materials));
         action_reg.insert(models.back());
     }
-    if (input.wls)
-    {
-        models.push_back(std::make_shared<WavelengthShiftModel>(
-            action_reg.next_id(), input.wls, materials, "wls"));
-        action_reg.insert(models.back());
-    }
-    if (input.wls2)
+    if (input.wls || input.wls2)
     {
-        models.push_back(std::make_shared<WavelengthShiftModel>(
-            action_reg.next_id(), input.wls2, materials, "wls2"));
-        action_reg.insert(models.back());
+        WlsGeneratorAction::Input gen_inp;
+        gen_inp.aux_id = aux_reg.next_id();
+        gen_inp.gen_id = gen_reg.next_id();
+        gen_inp.capacity = gen_capacity;
+
+        if (input.wls)
+        {
+            auto model
+                = std::make_shared<WavelengthShiftModel>(action_reg.next_id(),
+                                                         gen_inp.aux_id,
+                                                         input.wls,
+                                                         materials,
+                                                         GeneratorType::wls);
+            gen_inp.wls = model;
+            action_reg.insert(model);
+            models.push_back(std::move(model));
+        }
+        if (input.wls2)
+        {
+            auto model
+                = std::make_shared<WavelengthShiftModel>(action_reg.next_id(),
+                                                         gen_inp.aux_id,
+                                                         input.wls2,
+                                                         materials,
+                                                         GeneratorType::wls2);
+            gen_inp.wls2 = model;
+            action_reg.insert(model);
+            models.push_back(std::move(model));
+        }
+
+        // Action ID for the WLS generator must be greater than all model IDs
+        gen_inp.action_id = action_reg.next_id();
+        auto gen = std::make_shared<WlsGeneratorAction>(std::move(gen_inp));
+        action_reg.insert(gen);
+        aux_reg.insert(gen);
+        gen_reg.insert(gen);
     }
     return models;
 }

src/celeritas/optical/PhysicsParams.hh

@@ -17,6 +17,8 @@
 namespace celeritas
 {
 class ActionRegistry;
+class AuxParamsRegistry;
+class GeneratorRegistry;
 class MaterialParams;
 
 namespace optical
@@ -36,6 +38,8 @@ class PhysicsParams final : public ParamsDataInterface<PhysicsParamsData>
     //! \name Type aliases
     using ActionIdRange = Range<ActionId>;
     using SPActionRegistry = std::shared_ptr<ActionRegistry>;
+    using SPAuxRegistry = std::shared_ptr<AuxParamsRegistry>;
+    using SPGeneratorRegistry = std::shared_ptr<GeneratorRegistry>;
     using SPConstModel = std::shared_ptr<Model const>;
     using SPConstCoreMaterials
         = std::shared_ptr<celeritas::MaterialParams const>;
@@ -48,7 +52,10 @@ class PhysicsParams final : public ParamsDataInterface<PhysicsParamsData>
     explicit PhysicsParams(inp::OpticalBulkPhysics const&,
                            SPConstMaterials const&,
                            SPConstCoreMaterials const&,
-                           SPActionRegistry);
+                           SPActionRegistry,
+                           SPAuxRegistry,
+                           SPGeneratorRegistry,
+                           size_type);
 
     //! Number of optical models
     inline ModelId::size_type num_models() const { return models_.size(); }
@@ -81,7 +88,10 @@ class PhysicsParams final : public ParamsDataInterface<PhysicsParamsData>
     VecModels build_models(inp::OpticalBulkPhysics const&,
                            SPConstMaterials,
                            SPConstCoreMaterials,
-                           ActionRegistry&) const;
+                           ActionRegistry&,
+                           AuxParamsRegistry&,
+                           GeneratorRegistry&,
+                           size_type) const;
     void build_mfps(MaterialParams const& mats, HostValue& data) const;
     //!@}
 };

src/celeritas/optical/Runner.cc

@@ -11,6 +11,7 @@
 #include "corecel/io/OutputInterfaceAdapter.hh"
 #include "corecel/io/OutputRegistry.hh"
 #include "celeritas/inp/StandaloneInputIO.json.hh"
+#include "celeritas/phys/GeneratorRegistry.hh"
 #include "celeritas/setup/Problem.hh"
 
 #include "CoreParams.hh"
@@ -142,8 +143,13 @@ auto Runner::run() const -> Result
 
     Result result;
     result.counters = state_->accum();
-    result.counters.generators.push_back(
-        loaded_.problem.generator->counters(*state_->aux()).accum);
+    for (auto gen_id : range(GeneratorId(this->params()->gen_reg()->size())))
+    {
+        auto const gen = this->params()->gen_reg()->at(gen_id);
+        CELER_ASSERT(gen);
+        result.counters.generators.push_back(
+            gen->counters(*state_->aux()).accum);
+    }
     result.action_times
         = loaded_.problem.transporter->get_action_times(*state_->aux());
 

src/celeritas/optical/Types.cc

@@ -12,6 +12,21 @@ namespace celeritas
 {
 namespace optical
 {
+//---------------------------------------------------------------------------//
+/*!
+ * Get a string corresponding to the photon generating process.
+ */
+char const* to_cstring(GeneratorType value)
+{
+    static EnumStringMapper<GeneratorType> const to_cstring_impl{
+        "cherenkov",
+        "scintillation",
+        "wls",
+        "wls2",
+    };
+    return to_cstring_impl(value);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Get a string corresponding to an optical surface physics step.

src/celeritas/optical/Types.hh

@@ -33,6 +33,8 @@ enum class GeneratorType
 {
     cherenkov,
     scintillation,
+    wls,
+    wls2,
     size_
 };
 
@@ -95,6 +97,7 @@ enum class WlsDistribution
 // FREE FUNCTIONS
 //---------------------------------------------------------------------------//
 
+char const* to_cstring(GeneratorType);
 char const* to_cstring(SurfacePhysicsOrder);
 char const* to_cstring(ReflectionMode);
 char const* to_cstring(WlsDistribution);