Multiple WLS models and shared actions

[hhollenb]

In Geant there are two WLS models which are identical copies of each other. They reflect possibly different WLS properties of a given molecule, or multiple types of WLS molecules in a given material. Similar issues arise with multiple scintillation models, so discussing how to support multiple WLS models is a good proving ground for other cases.

The Problem

The imported WLS data consists of:

• The MFP table for all optical materials in ImportOpticalModel
• The ImportWavelengthShift for each optical material
  All WLS models use the same executor / interactor. The only difference is that different models pass different WavelengthShiftData to the executor.

Ideas

0. Don't support multiple models (right now)

1. Create copies of each model that are completely independent

This is the most straightforward, but requires some extra logic for how to map and build multiple identical WLS models. The main concern (which might not be important) is that these models will all be separate actions and therefore won't launch with the same kernels. Because the executor and interactor code is the same (just different references to WavelengthShiftData) they could be readily combined into the same kernel launch.

Whether this is an actual concern depends on how many WLS models there are, how much they fill each warp / block, and how much they overlap. This is probably something worth testing, and seeing if other cases like multiple scintillation models have the same runtime profile.

HERE BE DRAGONS: the following are wacky ideas I've brainstormed and pursued just enough to see if they're feasible. They all need to be heavily refined before being implemented, but I'd rather have some expert input on them before spending too much time on them.

2. Combine all WLS into one model

Keeping with the model-action identification, we can just put all the WLS models into a single total one. It's summed XS will compete like the individual model XS, and if selected it can then do a sampling of the XS of the submodels to choose which parameters to pass to the executor. This has the small added benefit of only calculating the WLS submodel grids when WLS has been selected.

The main issue with this is how to calculate the MFPs of the total model. Because we're doing linear interpolations of MFPs, they need to be inverted to get the XS which we sum over and use to select models. When the MFP tables for a model have a non-trivial overlap, there's no direct way of making a summed MFP grid without introducing some error into the calculation (a very rough calculation gave me ~3% error of the XS on some sample data). This is dependent on the specifics of the grid. A different calculator might be useful, but then we'd need to code special behavior for the MFP calculation step based on the model type.

3. Multiple non-action models with a single model-action

The WLS models won't be actions; instead a common WLSModelAction is registered in the action registry and it's action ID is used by each WLS model. The models now adopt the meaning of "discrete physical processes competing with their own MFP and data", whlie the model action is independent of that information and is only responsible for executing the step.

Because the models don't have action IDs, we'd need to store some extra information in the track state (e.g. selected model ID) or pass it some other way.

4. Multiple action models and a single composite action

Have multiple models built and work as normal, except their step method is empty. Instead, a single composite action is added and its ConditionalTrackExecutor is predicated on a range of ActionIds. When actions are executed, the composite action will filter for all the WLS action IDs and runs a single kernel on those. If the IDs are sorted, then the composite action can just modify the track ID offsets to the total range of the submodel offsets (this requires the WLS models to be contiguous in action IDs, which I think is a simple restriction). The action IDs from the simulation view (which are just the submodel action IDs) let's the composite action select the correct parameters to pass to the executor / interactor.

This could be integrated with choice 1, where if a composite action is present then the WLS models won't execute their step functions, but if no composite action has been built then each model executes their usual step independently.

[hhollenb]

@sethrj @whokion you might find this interesting

[sethrj]

I think a single model makes the most sense to me: we'd basically be treating it like sampling from multiple optical "elements" and the cross section for the model/process would be a sum over cross sections for each element. (In regular EM physics we effectively pre-calculate all this.) I think we're close to the point where we want to enable more generic cross section calculations (and indeed we haven't yet implemented step length calculation for photons) so this is a good thing to keep in mind. Indeed we already basically do this with the MacroXsCalculator and LivermorePEMicroXsCalculator and LivermorePEExecutor sampling...

[amandalund]

I think both 1 and 2 are good options. I'm assuming in most typical cases there'd be only one or two wavelength shifters in a material, but I don't really know... are there applications that might use more? One other small thing to keep in mind for 2 if we're sampling the WLS submodel in the executor is if we run out of secondary storage and want to reallocate dynamically while maintaining reproducibility, we'll probably need an extra piece of state data to cache the sampled property (analogous to storing the sampled element in EM).

[whokion]

Thanks @hhollenb and @amandalund for several options and preferred approach which can be discussed further later even though I kind of support the single model approach as @sethrj is pushing as the first choice. As we all understand, WLS is different 1) from the scintillation process as scintillation does not have xs (i.e., not competing with other processes, but is always forced if the step is in a scintillation material) and its associated action is to generate secondary photons (which in some sense, equivalent or behavior like an interactor in a different context) according to multiple scintillation structures (yield, spectrum, timing, etc), 2) from a typical EM process with multiple models (like msc or brems) since cross sections can be overlaid each other (i.e., xs can be additive or combined if the applicable range is overlapped), but their sampling structure is totally independent each other and cannot be combined (i.e., can not be looped over sampling all structures/data for the final interaction output similar to the scintillation photon generation) 3) from a general/combined process (such as photon or neutron general processes) since there is no other process to compete for the step even though xs can be combined, but models/interactor codes are completely different (i.e., not like the case 2 in which the interactor is same, but with different data).

Even though there may be technical details how to combine absorption lengths provided by users (with different data in bins size and ranges) as Hayden pointed out, the following is my alternative proposal to consider: 1) combine absorption length (or inverse mfp) from an arbitrary number of WLS processes (ex. WLS, WLS2 as in Geant4, but can be extendable) and pass it to the discrete step competition with other processes (i.e., absorption, Raleigh, boundary, etc), 2) build all associated device data including the individual cross sections or relative weights (not the combined one), 3) then, randomly select one WLS model inside the interactor based on the relative weight to access the corresponding data 4) the rest code is common for all WLS processes. The down-sides of this approach are 1) to build and pass additional data to device (i.e, relative xs weight or cdf) and 2) to add the additional model selector based on the relative weight at the given energy. Other than those, the implementation is as clean/same as others and needs no special treatment outside the WLS process or model and equivalent to the multiple independent WLS model approach. Of course, we can discuss for the better or more performant approach.

[sethrj]

Sounds good, but I think the main hurdle will be whether we can sacrifice some accuracy in the MFP/XS calculation. We've been trying to preserve the Geant4 behavior of linearly interpolating on mfp rather than the inverse (i.e., the cross section). We can't precalculate a combined (summed over WLS models) MFP because the interpolation would result in different values outside of the grid points. Since the linear interpolation is somewhat arbitrary, my gut instinct is that we could use MFPs for everything...

[whokion]

As long as the interpolation is linear, I guess that the geometry sum of individual mfp should be equivalent to the inverse of the total cross section - so no approximation for the combined mfp except some pure technical tweak how to combine mfps together if their bin sizes are different or not fully overlapped (unlikely as the users want to overlapped the absorption transition region intentionally even though it is a hard restriction). The outside the grid points, it recovers the correct value for the total inverse cross section as $\sigma^{-1}_{t} = l_o l_x /(l_o + l_x) \rightarrow l_x$ when $l_o \gg l_x$ (which should be the usually assumption (i.e., $l_o \rightarrow$ DBL_MAX) as we discussed earlier in the slack physics channel). Okay, our hardest bound may be importing Geant4-style wls data structure and mimicking its approach as we are all royal to Geant4^^.

[hhollenb]

Did some calculations to estimate the relative error of using a linear interpolation for the summed models $\tilde{\lambda}$ versus the (harmonic) sum of linearly interpolated models $\lambda$. Assuming the relative step size
$$\alpha_{mi} = \Delta \lambda_{mi} / \lambda_{mi}$$
for each model $m$ on the $i$-th interval is small, then a rough upper bound to 2nd order in $\alpha_{mi}$ is:
$$\frac{|\lambda - \tilde{\lambda}|}{\lambda} \leq \frac{1}{4} \left( \max{\alpha_{mi}^2} - \min{\alpha_{mi}^2} \right)$$
Before some of the simplifications, the relative step sizes are summed as
$$\sum_m \frac{\alpha_{mi}}{\lambda_{mi}}$$
so any models with large grid values $\lambda_{mi}$ will have suppressed contribution to the error term. Thus the major source of error in doing a linear interpolation of the summed model grid will be in the non-trivial overlap regions, and will be small if the relative step sizes of the models are not far apart.

[sethrj]

Diverting a bit for our other WLS discussion: instead of having a secondary stack as part of the actual photon, maybe have a WLS "distribution" array like we do with cherenkov/scintillation (TODO: generalize distributions so that we can more easily support user distributions), and then inside the WLS interaction:

1. absorb if no photons produced
2. "scatter" + time delay existing track for first photon emitted
3. save a distribution with (N - 1) photons to generate from
   After the interaction, we'd do the same thing as we do with our existing distributions (sum, allocate, generate).