Add an Optical WLS (#1507)

* Add WavelengthShift

* Add a secondary photon store

* Fix a warning from the window build: conversion from celeritas::real_type to int

* Fix errors from unit and precision tests

* Rename ScalarPropertyWLS to WlsMaterialRecord

* Add a short description of WLS and changes for other review comments

* Add a note on the secondary memory allocation

* Add the unit for the expected average WLS time

* Use consistent float types to fix `float` build error

* Add todos

* Fix no-WLS construction

* Change the interaction action from_unchanged() to from_absorption() if there is no secondary photon generated from WLS

* Fixed the test result when there is no emitted WLS photon

---------

Co-authored-by: Soon Yung Jun <syjun@login29.chn.perlmutter.nersc.gov>
Co-authored-by: Soon Yung Jun <syjun@login20.chn.perlmutter.nersc.gov>
Co-authored-by: Soon Yung Jun <syjun@login31.chn.perlmutter.nersc.gov>
Co-authored-by: Seth R Johnson <johnsonsr@ornl.gov>
Co-authored-by: Soon Yung Jun <syjun@login07.chn.perlmutter.nersc.gov>

src/celeritas/CMakeLists.txt

@@ -85,6 +85,7 @@ list(APPEND SOURCES
   optical/MaterialParams.cc
   optical/TrackInitParams.cc
   optical/ScintillationParams.cc
+  optical/WavelengthShiftParams.cc
   optical/action/ActionGroups.cc
   optical/action/LocateVacanciesAction.cc
   optical/detail/OffloadParams.cc

src/celeritas/optical/CerenkovParams.cc

@@ -49,6 +49,7 @@ CerenkovParams::CerenkovParams(SPConstMaterial material)
         std::vector<real_type> integral(energy.size());
         for (size_type i = 1; i < energy.size(); ++i)
         {
+            // TODO: use trapezoidal integrator helper class
             integral[i] = integral[i - 1]
                           + real_type(0.5) * (energy[i] - energy[i - 1])
                                 * (1 / ipow<2>(refractive_index[i - 1])

src/celeritas/optical/Interaction.hh

@@ -8,7 +8,9 @@
 #pragma once
 
 #include "corecel/Macros.hh"
+#include "corecel/cont/Span.hh"
 #include "geocel/Types.hh"
+#include "celeritas/optical/TrackInitializer.hh"
 
 namespace celeritas
 {
@@ -29,10 +31,12 @@ struct Interaction
         scattered,  //!< Still alive, state has changed
         absorbed,  //!< Absorbed by the material
         unchanged,  //!< No state change, no secondaries
+        failed,  //!< Ran out of memory during sampling
     };
 
     Real3 direction;  //!< Post-interaction direction
     Real3 polarization;  //!< Post-interaction polarization
+    Span<TrackInitializer> secondaries;  //!< Emitted secondaries
     Action action{Action::scattered};  //!< Flags for interaction result
 
     //! Return an interaction respresenting an absorbed process
@@ -41,6 +45,9 @@ struct Interaction
     //! Return an interaction with no change in the track state
     static inline CELER_FUNCTION Interaction from_unchanged();
 
+    // Return an interaction representing a recoverable error
+    static inline CELER_FUNCTION Interaction from_failure();
+
     //! Whether the state changed but did not fail
     CELER_FUNCTION bool changed() const
     {
@@ -72,6 +79,17 @@ CELER_FUNCTION Interaction Interaction::from_unchanged()
     return result;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Indicate a failure to allocate memory for secondaries.
+ */
+CELER_FUNCTION Interaction Interaction::from_failure()
+{
+    Interaction result;
+    result.action = Action::failed;
+    return result;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace optical
 }  // namespace celeritas

src/celeritas/optical/WavelengthShiftData.hh

@@ -0,0 +1,80 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/optical/WavelengthShiftData.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Macros.hh"
+#include "corecel/Types.hh"
+#include "corecel/data/Collection.hh"
+#include "celeritas/Types.hh"
+#include "celeritas/grid/GenericGridData.hh"
+
+namespace celeritas
+{
+namespace optical
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Material dependent scalar property of wavelength shift (WLS).
+ */
+struct WlsMaterialRecord
+{
+    real_type mean_num_photons{};  //!< Mean number of reemitted photons
+    real_type time_constant{};  //!< Time delay of WLS [time]
+
+    //! Whether all data are assigned and valid
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return mean_num_photons > 0 && time_constant > 0;
+    }
+};
+
+//---------------------------------------------------------------------------//
+/*!
+ * Wavelength shift data
+ */
+template<Ownership W, MemSpace M>
+struct WavelengthShiftData
+{
+    template<class T>
+    using Items = Collection<T, W, M>;
+    template<class T>
+    using OpticalMaterialItems = Collection<T, W, M, OpticalMaterialId>;
+
+    //// MEMBER DATA ////
+
+    OpticalMaterialItems<WlsMaterialRecord> wls_record;
+
+    // Grid energy tabulated as a function of the cumulative probability.
+    OpticalMaterialItems<GenericGridRecord> energy_cdf;
+
+    // Backend data
+    Items<real_type> reals;
+
+    //// MEMBER FUNCTIONS ////
+
+    //! Whether all data are assigned and valid
+    explicit CELER_FUNCTION operator bool() const
+    {
+        return !wls_record.empty() && !energy_cdf.empty() && !reals.empty();
+    }
+
+    //! Assign from another set of data
+    template<Ownership W2, MemSpace M2>
+    WavelengthShiftData& operator=(WavelengthShiftData<W2, M2> const& other)
+    {
+        CELER_EXPECT(other);
+        wls_record = other.wls_record;
+        energy_cdf = other.energy_cdf;
+        reals = other.reals;
+        return *this;
+    }
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace optical
+}  // namespace celeritas

src/celeritas/optical/WavelengthShiftParams.cc

@@ -0,0 +1,108 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/optical/WavelengthShiftParams.cc
+//---------------------------------------------------------------------------//
+#include "WavelengthShiftParams.hh"
+
+#include <vector>
+
+#include "corecel/data/CollectionBuilder.hh"
+#include "celeritas/Types.hh"
+#include "celeritas/grid/GenericGridBuilder.hh"
+#include "celeritas/grid/GenericGridInserter.hh"
+#include "celeritas/io/ImportData.hh"
+
+namespace celeritas
+{
+namespace optical
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Construct wavelength shift (WLS) data with imported data.
+ */
+std::shared_ptr<WavelengthShiftParams>
+WavelengthShiftParams::from_import(ImportData const& data)
+{
+    CELER_EXPECT(!data.optical_materials.empty());
+
+    if (!std::any_of(
+            data.optical_materials.begin(),
+            data.optical_materials.end(),
+            [](auto const& iter) { return static_cast<bool>(iter.wls); }))
+    {
+        // No wavelength shift data present
+        return nullptr;
+    }
+
+    Input input;
+    for (auto const& mat : data.optical_materials)
+    {
+        input.data.push_back(mat.wls);
+    }
+    return std::make_shared<WavelengthShiftParams>(std::move(input));
+}
+
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with wavelength shift (WLS) input data.
+ */
+WavelengthShiftParams::WavelengthShiftParams(Input const& input)
+{
+    CELER_EXPECT(input.data.size() > 0);
+    HostVal<WavelengthShiftData> data;
+
+    CollectionBuilder wls_record{&data.wls_record};
+    GenericGridInserter insert_energy_cdf(&data.reals, &data.energy_cdf);
+    for (auto const& wls : input.data)
+    {
+        if (!wls)
+        {
+            // No WLS data for this material
+            wls_record.push_back({});
+            insert_energy_cdf();
+            continue;
+        }
+
+        // WLS material properties
+        WlsMaterialRecord record;
+        record.mean_num_photons = wls.mean_num_photons;
+        record.time_constant = wls.time_constant;
+        wls_record.push_back(record);
+
+        // Calculate the WLS cumulative probability of the emission spectrum
+        // Store WLS component tabulated as a function of photon energy
+        auto const& comp_vec = wls.component;
+        std::vector<double> cdf(comp_vec.x.size());
+
+        CELER_ASSERT(comp_vec.y[0] > 0);
+        // The value of cdf at the low edge is zero by default
+        cdf[0] = 0;
+        for (size_type i = 1; i < comp_vec.x.size(); ++i)
+        {
+            // TODO: use trapezoidal integrator helper class
+            cdf[i] = cdf[i - 1]
+                     + 0.5 * (comp_vec.x[i] - comp_vec.x[i - 1])
+                           * (comp_vec.y[i] + comp_vec.y[i - 1]);
+        }
+
+        // Normalize for the cdf probability
+        for (size_type i = 1; i < comp_vec.x.size(); ++i)
+        {
+            cdf[i] = cdf[i] / cdf.back();
+        }
+        // Note that energy and cdf are swapped for the inverse sampling
+        insert_energy_cdf(make_span(cdf), make_span(comp_vec.x));
+    }
+    CELER_ASSERT(data.energy_cdf.size() == input.data.size());
+    CELER_ASSERT(data.wls_record.size() == data.energy_cdf.size());
+
+    data_ = CollectionMirror<WavelengthShiftData>{std::move(data)};
+    CELER_ENSURE(data_);
+}
+
+//---------------------------------------------------------------------------//
+}  // namespace optical
+}  // namespace celeritas

src/celeritas/optical/WavelengthShiftParams.hh

@@ -0,0 +1,58 @@
+//----------------------------------*-C++-*----------------------------------//
+// Copyright 2024 UT-Battelle, LLC, and other Celeritas developers.
+// See the top-level COPYRIGHT file for details.
+// SPDX-License-Identifier: (Apache-2.0 OR MIT)
+//---------------------------------------------------------------------------//
+//! \file celeritas/optical/WavelengthShiftParams.hh
+//---------------------------------------------------------------------------//
+#pragma once
+
+#include "corecel/Types.hh"
+#include "corecel/data/CollectionMirror.hh"
+#include "corecel/data/ParamsDataInterface.hh"
+#include "celeritas/io/ImportOpticalMaterial.hh"
+
+#include "WavelengthShiftData.hh"
+
+namespace celeritas
+{
+struct ImportData;
+
+namespace optical
+{
+//---------------------------------------------------------------------------//
+/*!
+ * Build and manage wavelength shift (WLS) data.
+ */
+class WavelengthShiftParams final
+    : public ParamsDataInterface<WavelengthShiftData>
+{
+  public:
+    //! Material-dependent WLS data, indexed by \c OpticalMaterialId
+    struct Input
+    {
+        std::vector<ImportWavelengthShift> data;
+    };
+
+  public:
+    // Construct with imported data
+    static std::shared_ptr<WavelengthShiftParams>
+    from_import(ImportData const& data);
+
+    // Construct with WLS input data
+    explicit WavelengthShiftParams(Input const& input);
+
+    //! Access WLS data on the host
+    HostRef const& host_ref() const final { return data_.host_ref(); }
+
+    //! Access WLS data on the device
+    DeviceRef const& device_ref() const final { return data_.device_ref(); }
+
+  private:
+    // Host/device storage and reference
+    CollectionMirror<WavelengthShiftData> data_;
+};
+
+//---------------------------------------------------------------------------//
+}  // namespace optical
+}  // namespace celeritas