Add builtin surface roughness models

[hhollenb]

Adds support for built-in roughness models: polished, uniform smear, and Gaussian roughness.

The BuiltinSurfaceModelBuilder builds surface models from supported input data in the form of std::map<PhysSurfaceId, InputT>. For reflectivity and interaction models (not yet supported), FakeSurfaceModels are built instead. BuiltinSurfaceModel provides a common interface for accessing surfaces from the input data and wrapping executors with an appropriate applier (cf. InteractionApplier in volumetric physics models).

The SurfaceModelView is a wrapper view to directly access physics relevant for a surface model. This encapsulates a lot of the bookkeeping of the surface traversal and represents a single physical surface crossing. This also hides some of the previous complications of SurfacePhysicsView, which now primarily handles subsurface traversal.

Should be direct to add reflectivity and interaction models once we get calculators for them implemented.

[github-actions]

Test summary

 4 564 files   7 293 suites   13m 37s ⏱️
 1 983 tests  1 967 ✅ 16 💤 0 ❌
25 335 runs  25 240 ✅ 95 💤 0 ❌

Results for commit c6aa9e1.

♻️ This comment has been updated with latest results.

[sethrj]

I'm still working through this but it's looking fantastic! I've got a few initial comments to suggest while I wrap up the vecgeom surface normal work.

[hhollenb]

@sethrj I've been seeing some errors with finding sincos in the GaussianRoughnessSampler when trying to build some tests. Not sure what the regression was, maybe a missing include?

[sethrj]

@hhollenb Can you post the errors? The only ones I'm seeing in the CI (about missing operators) is due to not including corecel/math/ArrayOperators.hh in surface/detail/InitBoundaryExecutor.hh.

[hhollenb]

Sorry I was looking at the build-docker tests:

 [429/1146] Building CUDA object src/celeritas/CMakeFiles/celeritas.dir/optical/surface/model/GaussianRoughnessModel.cu.o
FAILED: src/celeritas/CMakeFiles/celeritas.dir/optical/surface/model/GaussianRoughnessModel.cu.o 
/usr/local/cuda/bin/nvcc -forward-unknown-to-host-compiler -Dceleritas_EXPORTS -I/__w/celeritas/celeritas/src -I/__w/celeritas/celeritas/build/include -isystem /usr/local/cuda/targets/x86_64-linux/include -isystem /opt/view/include -isystem /opt/software/linux-rocky9-x86_64_v3/gcc-11.5.0/openmpi-5.0.6-qzolp62fnhdcym2hkvnbvxujqzgeqtf2/include -Werror all-warnings -g -std=c++17 "--generate-code=arch=compute_70,code=[compute_70,sm_70]" -Xcompiler=-fPIC -MD -MT src/celeritas/CMakeFiles/celeritas.dir/optical/surface/model/GaussianRoughnessModel.cu.o -MF src/celeritas/CMakeFiles/celeritas.dir/optical/surface/model/GaussianRoughnessModel.cu.o.d -x cu -c /__w/celeritas/celeritas/src/celeritas/optical/surface/model/GaussianRoughnessModel.cu -o src/celeritas/CMakeFiles/celeritas.dir/optical/surface/model/GaussianRoughnessModel.cu.o
/__w/celeritas/celeritas/src/celeritas/optical/surface/GaussianRoughnessSampler.hh(106): error: no instance of overloaded function "celeritas::sincos" matches the argument list
            argument types are: (celeritas::real_type, celeritas::real_type *, celeritas::real_type *)
          sincos(alpha, &sin_alpha, &cos_alpha);
          ^
/__w/celeritas/celeritas/src/corecel/math/Turn.hh(124): note #3326-D: function "celeritas::sincos(celeritas::IntQuarterTurn, int *, int *)" does not match because argument #1 does not match parameter
                          void sincos(IntQuarterTurn r, int* sinv, int* cosv)
                               ^
/__w/celeritas/celeritas/src/corecel/math/Turn.hh(98): note #3327-D: candidate function template "celeritas::sincos(celeritas::Turn_t<T>, T *, T *)" failed deduction
                            void sincos(Turn_t<T> r, T* sinv, T* cosv)
                                 ^

I'll fix the InitBoundaryExecutor missing include as well

[sethrj]

Aha, I see what's happening... the error is when building with NVCC, and there's a celeritas::sincos that shadows the builtin nvidia sincos 🤔 we need to get the std/global/nvidia namespace issues worked out. Maybe we just add celeritas functions that call the appropriate global/std namespace functions for host/nvidia? This issue is more complicated than it sounds due to functions like https://en.cppreference.com/w/cpp/numeric/math/fabs being for a single type in the global (C) namespace but overloaded in the std (C++) namespace.

[amandalund]

we need to get the std/global/nvidia namespace issues worked out

Not sure we made an issue for this, but a related reminder from the CUDA math API:

Note also that due to implementation constraints, certain math functions from std:: namespace may
be callable in device code even via explicitly qualified std:: names. However, such use is discouraged,
since this capability is unsupported, unverified, undocumented, not portable, and may change without
notice

[sethrj]

Yep I do remember that @amandalund . Created #1950 (which also explains why @hhollenb 's latest build is failing the CPU build...)

[sethrj]

Still working my way through this, and it's looking great! I'm really impressed.

[sethrj]

It might be good to clarify whether this is "internal" as in "implementation" or "sub-surface" (I can't remember what we ended up deciding or postponing...)

EDIT: I remember, originally this was a type alias inside SurfaceModel , and it's just a local index for the data. I wonder if it'd be better to call it ModelSurfaceId or something after all? :\

[hhollenb]

Maybe PerModelSurfaceId?

[sethrj]

It's generally not a good idea to put anonymous namespace/static functions in a header file because the code will be duplicated across multiple translation units (rather than eliminated by the linker). Since only operator+ is defined (no other math operators, and only with a mixed type) I recommend putting the function in SurfaceUtils and giving it a descriptive name like advance_subsurface_along(pos, dir)

[sethrj]

To improve the readability of the class synopsis,

Suggested change

	surface_model(Real3 const&, SurfacePhysicsOrder) const;
	surface_model(Real3 const& normal, SurfacePhysicsOrder) const;

[sethrj]

I guess it's direction not normal? Regardless I think the traversal direction should be an enum (boolean) value on the state, not a property of the geometry direction.

[hhollenb]

Because the traversal direction changes whenever a surface interaction occurs, we always call the function with

auto surface_physics = track.surface_physics();
auto surface_model = surface_physics.surface_model(surface_physics.traversal_direction(track.geometry().dir()), step);

so I collapsed it down to

auto surface_model = track.surface_physics().surface_model(track.geometry().dir(), step);

but I can go back and expand it again for readability.

[sethrj]

I guess it's direction not normal? Regardless I think the traversal direction should be an enum (boolean) value on the state, not a property of the geometry direction.

[hhollenb]

Added some changes to simplify the applier code, use maps as inputs to models, and cache the track traversal direction. This should remove most of the ambiguities about whether vectors refer to the surface normal or track direction in function arguments, and also remove the need to query geometry in the track slot executors. The traversal direction will need to be updated after a reflection / refraction in the interaction step which will be implemented later.

[sethrj]

@hhollenb This really is great work, but to make it easier to review we need to try to split these into smaller requests. I think there's a quadratic scaling with PR review time...

Still working on finishing this up but we need to decide about this Builtin thing first. (Or maybe we should just change that in a follow-on...)

[sethrj]

I think this is too much complexity for incomplete functionality (i.e., we have to inherit from it in the Roughness model, so we end up with an extra model class that doesn't do anything on its own: it's not an interface, it's an implementation).

We should choose to:

1. Continue like we've done with the other actions in the code and have a bunch of similar-looking code (~50% looks the same but it's spread among many different lines)
2. Commit fully to a template+traits approach: each traits instance will have surface physics order, input type, "internal surface" record type, input-to-record builder class, (likely another params-type type for "backend storage"), and executor.

For this PR, is it at all possible to do the first one, including only a single model, and using the current workaround for the "fake" processes? Then we can generalize in the next PR.

[sethrj]

@hhollenb @amandalund Notes from today, which should be incorporated into documentation once stable:

Physics surface crossing state

Logical directions:

• traversal direction: forward and reverse, based on pre -> post volume crossing
• orientation: whether we're going from low -> high subsurface IDs, or high ->
  low

Physical directions:

• geometry direction (updated during tracking at end of substep)
• polarization (also updated at end of substep)
• geometry normal is constant throguhout the surface traversal and cached, with
  sign normalized to point from post into pre
• Local

Changing during traversal:

• surface position
• facet normal

Cached:

• geometry normal
• pre and post volume material

Executor/applier/etc diagram

This is the best I could come up with since there's a mix of functional programming ($f(x; p) \to y$) and procedural programming ($f(), g(), h()$) that gets mixed together at the Applier level, as well as a mix of objects with different levels of "temporariness" (executor/applier/trackexecutor):

The red bold arrows are return values, dashed arrows are "creates during call operator", black arrows are data passing. Construction arguments are to the left, call arguments to the right. I'll create an equivalent for the EM interactor for comparison... maybe that'll help inform our design choices.

[sethrj]

Should this be done when entering a surface, reflecting, or here? It doesn't seem related to the roughness per se.

[sethrj]

For comparison, here's a typical EM interactor (bethe-heitler):

In the EM case, there's a clean break between the procedural components and functional: the executor doesn't return an interactor, it executes the interaction. I actually think philosophically the choice @hhollenb is nicer. Instead of building the executor and, embarrassingly, calling it only once:

    BetheHeitlerInteractor interact(
        params, particle, dir, allocate_secondaries, material, element);

    auto rng = track.rng();
    return interact(rng);

you would construct the interactor and let someone else handle it:

    return BetheHeitlerInteractor{...};

However, this isn't an executor since it's not really executing (though in the EM case it could also be sampling the element!). The "Executor" misnomer is one issue causing confusion in this PR.

I think the other complication is that forwarding and combining arguments can be confusing, especially if nested. That's done a few mote times here than in EM: EnteringSurfaceNormalSampler, BSM::make_executor in particular. And the templates on BSM make that even more confusing.

Finally, the roughness "executor" takes a few select track states rather than the whole track view. Even though this is the "correct" thing to do (pass around only the data you need, not all the data you could possibly want, in order to make the data flow and operations clearer) I think it adds another layer of cognitive burden.

So I think this design is good at heart, but to improve clarity we should maybe:

• Rename the Executor to RoughnessSampler{Builder/Factory/...} if its job is to return a functor that samples given an RNG.
• Consider having RoughnessApplier do more, or maybe eliminating it/combining it as a helper in the executor rather than the other way around.