Implement framework for optical physics models

[hhollenb]

This is an initial implementation of optical models and their associated data structures.

General overview:

• OpticalModel: acts like both a process (building cross sections / mean free paths) and as a model (action to run the interactor).
• ImportedOpticalMaterials: physics data for optical models is really just optical material properties. This serves as a list of imported optical material data that models will access directly to get their physics data.
• OpticalModelBuilder: constructs the ImportedOpticalMaterials and is responsible for building the appropriate models based on the user's selection.
• OpticalPhysicsParams: responsible for building each optical model, adding them to the action registry, and building their MFP grids in OpticalPhysicsParamsData.
• OpticalMfpBuilder: A helper class to have optical models build MFP grids on a per optical material basis.
• AbsorptionModel: An almost trivial optical model. Has no data and just kills the track when called.
• RayleighModel: A more complex optical model. Included to make sure API for building model-specific data is sensible.

This PR isn't merge ready yet, and is meant to be a discussion point for how we want to organize optical models. Some functions and classes are empty, and I'm in the middle of adding tests.

Thoughts / discussion questions:

• May be better to have something like this in ImportData:

std::unordered_map<GeoMatIdx, ImportOpticalMaterialId> geo_optical_map;
std::vector<ImportOpticalMaterial> optical_materials;

• Will we want an OpticalMaterialParams? Is there ever a point where we need to process ImportOpticalMaterials and store the data outside of optical model specific data collections?
• How much flexibility do we want to give to user defined optical models? Should we have an ImportOpticalModel structure? A lambda table in ImportOpticalModel would duplicate the mean free paths in ImportOpticalMaterial...
• OpticalMfpBuilder is a bit awkward. There should be a better way to do this that is also portable to Cerenkov and scintillation parameter building.

[sethrj]

Does this supersede #1162 ?

[hhollenb]

Does this supersede #1162 ?

Yes, I'll close that one. I reused any relevant code but had to make enough major changes I thought it was worthwhile to start from a fresh branch.

[sethrj]

Sorry it's taken some time to go through this: this looks really solid; I think you've got basically all the skeleton there! I've got a thought or two about simplifying some aspects to accelerate the time it takes for us to flesh this out... perhaps the first PR could be just the interactors for rayleigh and absorption; so that way we can unit test them and make sure the initial OpticalInteraction class is good enough. (You know, we might even be able to leave out code like the secondary stack until later, just to keep things simple.) But I think virtually all of this code will eventually end up in the final version almost unchanged 😄

[sethrj]

The ImportedProcesses stuff in phys is a little janky; it also is some of the oldest code that we use. When I wrote it, I was hoping to reduce data copying and lifetime and so forth; but in reality there's just not that much data for us to worry about, especially with optical physics. I think it'd just be simpler and more effective to build each optical model directly from optical import data, so we don't even need this (admittedly lightweight) class.

[hhollenb]

In that case I think I'd like to change the optical map in ImportData to be a vector of ImportOpticalMaterials and change the map to be a map between GeoMatIds and OpticalMaterialIds. We need to just iterate over a list of optical materials in several places, and having a simple vector to index would be convenient.

[sethrj]

@hhollenb I fully support that! It's definitely more consistent with how we treat the other properties. Perhaps instead of a map between geo and optical matIDs, what about having an unsigned int optical_material_id field in ImportGeoMaterial? (You could define a static inline constexpr unsigned int unspecified = -1 like we have with ImportVolume.)

[sethrj]

I think we're going to need an optical interaction type—since it includes changes to polarisation, and I assume we're not ionizing anything so we never make any non-photon secondaries.

[sethrj]

I think we can put these guys in the detail namespace.

[sethrj]

Suggested change

	* Build the mean free paths for the model.
	* Build the mean free paths for one material for the model.

[sethrj]

It would be good to add a few notes about how this class differs from the main EM physics:

• We interpolate linearly on MFP vs energy, rather than macroscopic XS vs log energy (or scaled log energy)
• There's a one-to-one model/process correlation
• Only a single incident particle type is supported
• Not all materials in the problem have optical data (?)
• ....

[hhollenb]

* We interpolate linearly on MFP vs energy, rather than macroscopic XS vs log energy (or scaled log energy)

A thought on linear interpolation: this is effectively just due to our choice of grid, right? AFAIK, the code in OpticalPhysicsParams and determining the discrete action is agnostic of the grid, other than the fact it will calculate some XS / MFP for a given energy value. It's definitely important whether we're using MFP vs macroscopic XS, but otherwise the code behaves mostly the same as in the EM case when dealing with these grids.

[hhollenb]

Sorry it's taken some time to go through this: this looks really solid; I think you've got basically all the skeleton there! I've got a thought or two about simplifying some aspects to accelerate the time it takes for us to flesh this out... perhaps the first PR could be just the interactors for rayleigh and absorption; so that way we can unit test them and make sure the initial OpticalInteraction class is good enough. (You know, we might even be able to leave out code like the secondary stack until later, just to keep things simple.) But I think virtually all of this code will eventually end up in the final version almost unchanged 😄

I believe @whokion has some preliminary interactors implemented. The absorption interactor is trivial (just returns an absorbed state -- might be fine with just an executor) and doesn't need to be super rigorously tested. Rayleigh is non-trivial, but hopefully I got enough code in here to draft what setting up its data and model would look like. I'll discuss getting another PR up for just the interactions and tests over slack.

[whokion]

I can add OpticalRayInteractor (which I already have) once we decide how to build required optical data for the interactor - I was using OpticalRaylieghParams class, which can be absorbed to OpticalPhysicsParams or can be added too if we decide to have a separated data Params class for each model. Let's have a consensus first how we support optical material properties.

[sethrj]

@whokion I think the optical rayleigh model should have its own data independent of the physics. I assume it's not needed by the other models?

[whokion]

@sethrj Yes, the Rayleigh data is independent from other models - actually, there is no required optical data for the OpticalRayleigh interactor if all mfp are taken care by the OpticalPhysicsParams class. It was not clear to me whether each model is responsible to import its own mfp along with other properties or OpticalPhysicsParams manages all mfp.

[sethrj]

The plan is for the models to provide the MFP at setup time (usually by just processing the imported optical data) for the "optical physics" to manage all together at runtime.

[hhollenb]

Made a PR for absorption and Rayleigh interactors here: #1317

[hhollenb]

These model/Optical... files should've been deleted, ignore them and I'll remove them.

[sethrj]

There's a lot of great work here! I think the scope of this is just too large for one person and one PR; I apologize for not partitioning the work among us. We'll definitely use all the work you've done here, but we need to break it into more bite-size chunks.

To start, let's drill down on just constructing and evaluating interaction lengths:

• Model
• AbsorptionModel with "not implemented" for the step action
• MfpBuilder
• The minimal subset of physics data (not params yet)

Let's talk about this monday?

[sethrj]

This is redundant with the ImportOpticalMaterial.{rayleigh,absorption,wls} data which already tells us the available models and stores absorption_length (mfp): should we embed model information into the material data, or have separate model classes?

/*!
 * Store optical material properties.
 */
struct ImportOpticalMaterial
{
    ImportPhysicsVector refractive_index;
    ImportScintData scintillation;
    bool cerenkov{false};

    //! Whether minimal useful data is stored
    explicit operator bool() const
    {
        return static_cast<bool>(refractive_index);
    }
};

struct ImportOpticalWLSModel
{
    ImportPhysicsVector absorption_length;  //!< Absorption length [MeV, len]

    double mean_num_photons{};  //!< Mean number of re-emitted photons
    double time_constant{};  //!< Time delay between absorption and re-emission
    ImportPhysicsVector component;  //!< Re-emission population [MeV, unitless]

    explicit operator bool() const
    {
        return mfp && mean_num_photons > 0 && time_constant > 0
               && static_cast<bool>(absorption_length)
               && static_cast<bool>(component)
               && absorption_length.vector_type == ImportPhysicsVectorType::free
               && component.vector_type == absorption_length.vector_type;
    }
};

struct ImportOpticalAbsorptionModel
{
    ImportPhysicsVector absorption_length;  //!< Absorption length [MeV, len]

    explicit operator bool() const
    {
        return static_cast<bool>(absorption_length);
    }
};

[sethrj]

This can be obtained from the size of properties already in the table:

Suggested change

	return num_materials_;
	return this->host_ref().refractive_index.size();

[sethrj]

Delegate the constructor:

Suggested change

	//! Construct with defaults
	Model(ActionId id, std::string const& label, std::string const& description)
	: ConcreteAction(id, label, description)
	{
	}
	// Construct with ID, label, description
	using ConcreteAction::ConcreteAction;

[sethrj]

Suggested change

	virtual void build_mfps(detail::MfpBuilder build) const = 0;
	virtual void build_mfps(detail::MfpBuilder& build) const = 0;

[sethrj]

Add documentation

[sethrj]

The grid inserter should be outside this loop so that it can deduplicate data across models. I think we can also move most of this function into the implementation of MfpBuilder, which could include all the builders: see celeritas/optical/detail/MatScintSpecInserter.hh.