Add optical model importer and refactor imported optical materials

[hhollenb]

Refactoring imported optical materials

Added ImportedMaterials to act as a common storage for material data used by optical models (Rayleigh and WLS). Like with ImportedModels it helps prevent unnecessary copies of std::vector<ImportOpticalRayleigh> and std::vector<ImportWavelengthShift>.

RayleighModel was also updated to use an Input struct to concisely manage material dependencies for building MFP tables.

Model Importer

Similar to phys/ProcessBuilder, I added the ModelImporter to create models from imported data, as well as provide user build functionality like warn-and-ignore. The ModelBuilder concrete base class is meant to serve a similar purpose to phys/Process, and is just meant to build optical models with an action ID.

I'm trying to maintain the behavior of phys/PhysicsParams where importing processes/models is separate from being built and registered in the action registry, so a layer between ModelImporter and optical/PhysicsParams is necessary. The input for optical::PhysicsParams will look something like:

Input
{
    ActionRegistry* action_reg;
    std::vector<std::shared_ptr<ModelBuilder>> model_builders;
};

Simple builders are added to the model .cpp files and made accessible through static factory methods. There's some freedom in
determining where the builders live (or even using some type-erased lambdas if we want to be fancy).

N.B.: This loosely depends on #1519 so I'm leaving this as a draft PR, but can be more or less independently reviewed.

[github-actions]

Test summary

 4 478 files   6 904 suites   15m 34s ⏱️
 1 692 tests  1 686 ✅  6 💤 0 ❌
24 212 runs  24 128 ✅ 84 💤 0 ❌

Results for commit 01087c3.

♻️ This comment has been updated with latest results.

[hhollenb]

After discussing issue #1538 it might be worthwhile to call optical::ModelBuilder just optical::Process since they do roughly the same thing, and we can hide any multiple model details of WLS behind an overall WLS process.

Since #1519 was merged, this can now be reviewed.

[amandalund]

This looks great @hhollenb! I still have to finish up reviewing, but one thing I'm not totally sure about (and maybe @sethrj can comment) regarding the ModelBuilders is whether we still need to maintain the strict ordering in action ID when constructing the models/actions in the optical physics params. In the core physics we need the ordering because of how we access the various action IDs in the physics data (which seems a bit fragile anyway). We won't have msc, range, or integral rejection actions in optical, and discrete selection is a "pre-post" action, so we probably don't have to worry about registering the actions in a specific order and can just directly store the couple action IDs we need.

[sethrj]

I haven't had a chance to look at this yet, but I think it would be nice to preserve the contiguous mapping of action IDs from optical processes/models; that lets us do an indirection-free lookup of the post-step action from a sampled process. In other words, in the pre-step we're calculating cross sections into an array of P processes, and an interaction selects an index p from that list. If the processes map to sequential action IDs, then we can just add a constant start_action_id. Otherwise we'd need a separate array mapping p -> action id.

Given all the other work that has to be done in sampling and evaluating, that's probably a small cost, but I think it will be difficult to refactor to using the "less indirection" method if we choose the "easier to implement" method now.

[amandalund]

Right, we would still have contiguous model action IDs either way; I'm wondering whether we need the extra layer of the ModelBuilder to postpone registering the models so the model IDs are e.g. between the discrete select and failure action IDs.

[hhollenb]

My reasoning for using the ModelBuilders is more to try and separate out the user selection and building of the physics list, and the initialization of the physics list in PhysicsParams. We could pass in pre-registered models to the params, but then we'd need to check that their action IDs are ordered, and we'd have to coordinate physics options outside the construction of params.

[sethrj]

I feel like we're really in danger of overthinking/overengineering the import/process/model/mfps stuff. Let's spend our meeting today talking about the PR and the import stuff in general.

[hhollenb]

Some brainstorming thoughts on our discussion about importing data today:

1. I think we all readily agree that if a Params has been built, we should use it instead of passing its imported data to another Params (e.g. passing the core MaterialParams for the RayleighMfpBuilder instead of copying the material's temperature from import data).
2. It would be nice if the C++ API for the Params classes had type-safe and unit-safe inputs.
3. We need to be able to support different input modes (Geant4, ROOT, Json, etc.) as well as Celeritas being interfaced as a C++ library.
4. We must have lifetime safety for imported data. It would be nice to be able to free it after initialization, but we can't have any dangling references / pointers to freed memory.
5. It would be nice to never unnecessarily copy input data.

I'm personally partial towards a data-oriented design with strict separation of logic and data layers. Something like:

where all the logic on what data to import is done in the diamond importing classes (or specified in the options) and the rectangle classes just do consistency checks and store the data. The Celeritas input data struct can just be a loose collection of input structs that the params use. Params shouldn't rely upon previously built input data, so the input data can just be moved (red arrows).

Inputs consumed by param constructors can free or store what they need as necessary.

After the importers we can possibly have a user hook to modify / generate more imported data. I like the idea that there's a fairly static input that we don't do internal hidden decisions on, so users can be confident that what they put into Celeritas won't get modified under the hood.

(gotta head out real quick but I'll see if I can expand on the idea a bit more tomorrow)

[sethrj]

@hhollenb Sorry I didn't see your comment before the break, but I agree 100%. Let's discuss at our afternoon meeting if you're available?

[hhollenb]

A quick draft of what an input data consuming factory could look like:

// Concrete class
class ModelFactory final
{
  public:
    using ImportPhysicsTable = std::vector<ImportPhysicsVector>;

    ModelFactory(ImportPhysicsTable mfp_table)
        : mfp_table_(std::move(mfp_table))
    {
        // check table validity
    }

    void build_mfp_table(MfpBuilder&)
    {
        // copy the mfp table into the builder
    }

  private:
    ImportPhysicsTable mfp_table_;
};


// Diagnostic / logging data about the process
struct ProcessRecord
{
    std::string name;
    std::string description;
    ImportProcessClass process_class;
};

// Abstract base class
class ProcessFactory
{
  public:
    using ModelFactoryPair = std::tuple<std::shared_ptr<Model>, ModelFactory>;

    // consume process factory list to create model list
    static std::tuple<std::vector<std::shared_ptr<Model>>, std::vector<ModelFactory>>
        create_model_factory_list(ActionIdIter& iter,
                                  std::vector<std::unique_ptr<ProcessFactory>>&& processes)
    {
        std::vector<std::shared_ptr<Model>> models;
        std::vector<ModelFactory> factories;

        for (auto&& proc : processes)
        {
            auto model_factory_pairs = proc->create_models(iter);

            for (auto&& [model, factory] : model_factory_pairs)
            {
                models.push_back(std::move(model));
                factories.push_back(std::move(factory));
            }
        }

        return std::make_pair(std::move(models), std::move(factories));
    }

    // create a diagnostic record for the process
    virtual ProcessRecord initialize_record() const = 0;

    // disable copy constructor to ensure we only ever move factories
    ProcessFactory(ProcessFactory const&) = delete;
    ProcessFactory& operator=(ProcessFactory const&) = delete;

  protected:
    // override by specific process to create models and their corresponding factories
    virtual std::vector<ModelFactoryPair> create_models(ActionRegistry& reg) = 0;
};

class RayleighProcess : public ProcessFactory
{
  public:

    using SPConstMaterials = std::shared_ptr<MaterialParams const>;
    using SPConstCoreMaterials = std::shared_ptr<::celeritas::MaterialParams const>;


    struct Input
    {
        SPConstMaterials materials;
        SPConstCoreMaterials core_materials;
        std::vector<ImportOpticalRayleigh> imported_rayleigh;
    };


    RayleighProcess(std::vector<ImportPhysicsVector> mfp_table, Input input)
        : mfp_table_(std::move(mfp_table)), input_(std::move(input))
    {}

    ProcessRecord initialize_record() const final
    {
        return ProcessRecord{"rayleigh", "rayleigh process desc", ImportProcessClass::rayleigh};
    }

  protected:
    virtual std::vector<ModelFactoryPair> create_models(ActionIdIter& iter) final
    {
        auto rayleigh_model = std::make_shared<RayleighModel>(*iter++);

        /*
         * Use Rayleigh MFP calculator to fill missing entries in mfp_table_
         * Same logic as in RayleighModel::build_mfp_table
         */

        return std::vector<ModelFactoryPair>{
            std::make_tuple(std::move(rayleigh_model), ModelFactory{std::move(mfp_table_)})
        };
    }

  private:
    std::vector<ImportPhysicsVector> mfp_table_;
    Input input_;
};

The ProcessFactory mimics what Processes do in core physics, i.e. correspond to physical processes and create implementations which are the Model classes. I wanted to try and do a monadic approach of list of processes -> list of models where the processes get consumed and can move any of their imported data to models or model factories. This is done in the static method create_model_factory_list, and subclasses are responsible for overriding the create_models function. Making the create_models virtual function protected, and having the create_model_factory_list consume a list of unique pointers, means that create_models can move and invalidate the internal factory data, and it cannot be called multiple times on the same factory classes.

The ModelFactory is just a concrete data class to initialize the builders, and I associate them with their models via the tuples. They can be extended for the more complicated core physics tables by adding extra fields. Should only need to be written once - special cases like in Rayleigh can be handled in the process factory.

The ProcessRecord is just mock for data we might want to use for diagnostics. As far as I can tell, the processes in core physics never do any actual computation or logic, but it might still be handy to refer to a process' name or description.

Constructing PhysicsParams would be simply passing an appropriate list of unique_ptr<ProcessFactory>. Constructing that list can be handled by a ProcessBuilder or by the user or mocked. Models are completely separated from building their physics tables, and access is done solely through PhysicsTrackView. After initialization of PhysicsParams is done, then all of the factories can be freed and the data cleaned up.

As a test, the ImportData could be decomposed into something like

struct ImportOpticalData
{
    // Model data
    std::vector<ImportOpticalModel> models;

    // Material data
    std::vector<ImportOpticalProperty> properties;
    std::vector<ImportScintData> scintillation;
    std::vector<ImportOpticalRayleigh> rayleigh;
    std::vector<ImportWavelengthShift> wls;
};

Fields can then be moved into their respective factories in the optical::ModelImporter / phys::ProcessBuilder methods as necessary.

[sethrj]

Sorry it took me forever to get back to this. It also looks like there are some previous comments that weren't addressed?

[sethrj]

Thanks @hhollenb ! Sorry it took so long to get through. Maybe we could plan a little work half-day next week to update this to use inp?