Add option to specify max substeps in field propagator (#1236)

celeritas-project · May 15, 2024 · 567576d · 567576d
1 parent 9e06741
commit 567576d
Show file tree

Hide file tree

Showing 13 changed files with 64 additions and 24 deletions.
diff --git a/app/celer-g4/GlobalSetup.cc b/app/celer-g4/GlobalSetup.cc
@@ -168,6 +168,7 @@ void GlobalSetup::ReadInput(std::string const& filename)
         options_->action_times = input_.action_times;
         options_->default_stream = input_.default_stream;
         options_->auto_flush = input_.auto_flush;
+        options_->max_field_substeps = input_.field_options.max_substeps;
     }
     else if (ends_with(filename, ".mac"))
     {

diff --git a/app/celer-sim/Runner.cc b/app/celer-sim/Runner.cc
@@ -381,7 +381,8 @@ void Runner::build_core_params(RunnerInput const& inp,
     params.rng = std::make_shared<RngParams>(inp.seed);
 
     // Construct simulation params
-    params.sim = SimParams::from_import(imported, params.particle);
+    params.sim = SimParams::from_import(
+        imported, params.particle, inp.field_options.max_substeps);
 
     // Get the total number of events
     auto num_events = this->build_events(inp, params.particle);

diff --git a/src/accel/SetupOptions.hh b/src/accel/SetupOptions.hh
@@ -138,6 +138,11 @@ struct SetupOptions
     AlongStepFactory make_along_step;
     //!@}
 
+    //!@{
+    //! \name Field options
+    short int max_field_substeps{100};
+    //!@}
+
     //!@{
     //! \name Sensitive detector options
     SDSetupOptions sd;

diff --git a/src/accel/SetupOptionsMessenger.cc b/src/accel/SetupOptionsMessenger.cc
@@ -193,6 +193,9 @@ SetupOptionsMessenger::SetupOptionsMessenger(SetupOptions* options)
     add_cmd(&options->auto_flush,
             "autoFlush",
             "Number of tracks to buffer before offloading");
+    add_cmd(&options->max_field_substeps,
+            "maxFieldSubsteps",
+            "Limit on substeps in the field propagator");
 
     directories_.emplace_back(new CelerDirectory(
         "/celer/detector/", "Celeritas sensitive detector setup options"));

diff --git a/src/accel/SetupOptionsMessenger.hh b/src/accel/SetupOptionsMessenger.hh
@@ -34,6 +34,7 @@ struct SetupOptions;
   maxInitializers      | Maximum number of track initializers
   secondaryStackFactor | At least the average number of secondaries per track
   autoFlush            | Number of tracks to buffer before offloading
+  maxFieldSubsteps     | Limit on substeps in field propagator
 
  * The following option is exposed in the \c /celer/detector/ command
  * "directory":

diff --git a/src/accel/SharedParams.cc b/src/accel/SharedParams.cc
@@ -552,7 +552,8 @@ void SharedParams::initialize_core(SetupOptions const& options)
     params.rng = std::make_shared<RngParams>(CLHEP::HepRandom::getTheSeed());
 
     // Construct simulation params
-    params.sim = SimParams::from_import(*imported, params.particle);
+    params.sim = SimParams::from_import(
+        *imported, params.particle, options.max_field_substeps);
 
     // Construct track initialization params
     params.init = [&options] {

diff --git a/src/celeritas/field/FieldDriver.hh b/src/celeritas/field/FieldDriver.hh
@@ -73,6 +73,11 @@ class FieldDriver
 
     //// ACCESSORS ////
 
+    CELER_FUNCTION short int max_substeps() const
+    {
+        return options_.max_substeps;
+    }
+
     CELER_FUNCTION real_type minimum_step() const
     {
         return options_.minimum_step;

diff --git a/src/celeritas/field/FieldDriverOptions.cc b/src/celeritas/field/FieldDriverOptions.cc
@@ -44,6 +44,8 @@ void validate_input(FieldDriverOptions const& opts)
         << "invalid max_stepping_decrease " << opts.max_stepping_decrease);
     CELER_VALIDATE(opts.max_nsteps > 0,
                    << "invalid max_nsteps " << opts.max_nsteps);
+    CELER_VALIDATE(opts.max_substeps > 0,
+                   << "invalid max_substeps " << opts.max_substeps);
     CELER_ENSURE(opts);
 }
 

diff --git a/src/celeritas/field/FieldDriverOptions.hh b/src/celeritas/field/FieldDriverOptions.hh
@@ -61,6 +61,9 @@ struct FieldDriverOptions
     //! Maximum number of integrations (or trials)
     short int max_nsteps = 100;
 
+    //! Maximum number of substeps in the field propagator
+    short int max_substeps = 100;
+
     //! Initial step tolerance
     static constexpr inline real_type initial_step_tol = 1e-6;
 
@@ -84,7 +87,7 @@ struct FieldDriverOptions
 	       && (safety > 0 && safety < 1)
 	       && (max_stepping_increase > 1)
 	       && (max_stepping_decrease > 0 && max_stepping_decrease < 1)
-	       && (max_nsteps > 0);
+	       && (max_nsteps > 0) && (max_substeps > 0);
         // clang-format on
     }
 };
@@ -107,6 +110,7 @@ operator==(FieldDriverOptions const& a, FieldDriverOptions const& b)
            && a.max_stepping_increase == b.max_stepping_increase
            && a.max_stepping_decrease == b.max_stepping_decrease
            && a.max_nsteps == b.max_nsteps
+           && a.max_substeps == b.max_substeps
            && a.initial_step_tol == b.initial_step_tol
            && a.dchord_tol == b.dchord_tol
            && a.min_chord_shrink == b.min_chord_shrink;

diff --git a/src/celeritas/field/FieldDriverOptionsIO.json.cc b/src/celeritas/field/FieldDriverOptionsIO.json.cc
@@ -41,6 +41,7 @@ void from_json(nlohmann::json const& j, FieldDriverOptions& opts)
     FDO_INPUT(max_stepping_increase);
     FDO_INPUT(max_stepping_decrease);
     FDO_INPUT(max_nsteps);
+    FDO_INPUT(max_substeps);
 
 #undef FDO_INPUT
 }
@@ -64,6 +65,7 @@ void to_json(nlohmann::json& j, FieldDriverOptions const& opts)
         CELER_JSON_PAIR(opts, max_stepping_increase),
         CELER_JSON_PAIR(opts, max_stepping_decrease),
         CELER_JSON_PAIR(opts, max_nsteps),
+        CELER_JSON_PAIR(opts, max_substeps),
     };
 }
 

diff --git a/src/celeritas/field/FieldPropagator.hh b/src/celeritas/field/FieldPropagator.hh
@@ -21,16 +21,6 @@
 
 namespace celeritas
 {
-//---------------------------------------------------------------------------//
-/*!
- * Configuration options for the field propagator.
- */
-struct FieldPropagatorOptions
-{
-    //! Limit on substeps
-    static constexpr short int max_substeps = 100;
-};
-
 //---------------------------------------------------------------------------//
 /*!
  * Propagate a charged particle in a field.
@@ -74,10 +64,7 @@ class FieldPropagator
     static CELER_CONSTEXPR_FUNCTION bool tracks_can_loop() { return true; }
 
     //! Limit on substeps
-    static CELER_CONSTEXPR_FUNCTION short int max_substeps()
-    {
-        return FieldPropagatorOptions::max_substeps;
-    }
+    inline CELER_FUNCTION short int max_substeps() const;
 
     // Intersection tolerance
     inline CELER_FUNCTION real_type delta_intersection() const;
@@ -353,6 +340,16 @@ CELER_FUNCTION real_type FieldPropagator<DriverT, GTV>::delta_intersection() con
     return driver_.delta_intersection();
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Maximum number of substeps.
+ */
+template<class DriverT, class GTV>
+CELER_FUNCTION short int FieldPropagator<DriverT, GTV>::max_substeps() const
+{
+    return driver_.max_substeps();
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Distance to bump or to consider a "zero" movement.

diff --git a/src/celeritas/track/SimParams.cc b/src/celeritas/track/SimParams.cc
@@ -10,7 +10,7 @@
 #include "corecel/Assert.hh"
 #include "corecel/data/CollectionBuilder.hh"
 #include "corecel/math/Algorithms.hh"
-#include "celeritas/field/FieldPropagator.hh"
+#include "celeritas/field/FieldDriverOptions.hh"
 #include "celeritas/io/ImportData.hh"
 #include "celeritas/phys/ParticleParams.hh"
 
@@ -23,9 +23,12 @@ namespace celeritas
  * Construct with imported data.
  */
 std::shared_ptr<SimParams>
-SimParams::from_import(ImportData const& data, SPConstParticles particle_params)
+SimParams::from_import(ImportData const& data,
+                       SPConstParticles particle_params,
+                       short int max_field_substeps)
 {
     CELER_EXPECT(particle_params);
+    CELER_EXPECT(max_field_substeps > 0);
     CELER_EXPECT(data.trans_params);
     CELER_EXPECT(data.trans_params.looping.size() == particle_params->size());
 
@@ -35,10 +38,9 @@ SimParams::from_import(ImportData const& data, SPConstParticles particle_params)
     // Calculate the maximum number of steps a track below the threshold energy
     // can take while looping (ceil(max Geant4 field propagator substeps / max
     // Celeritas field propagator substeps))
-    size_type max_substeps = FieldPropagatorOptions::max_substeps;
-    size_type imported_max_substeps = data.trans_params.max_substeps;
-    CELER_ASSERT(imported_max_substeps >= max_substeps);
-    auto max_subthreshold_steps = ceil_div(imported_max_substeps, max_substeps);
+    CELER_ASSERT(data.trans_params.max_substeps >= max_field_substeps);
+    auto max_subthreshold_steps = ceil_div<size_type>(
+        data.trans_params.max_substeps, max_field_substeps);
 
     for (auto pid : range(ParticleId{input.particles->size()}))
     {
@@ -58,6 +60,17 @@ SimParams::from_import(ImportData const& data, SPConstParticles particle_params)
     return std::make_shared<SimParams>(input);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Construct with imported data and default max field substeps.
+ */
+std::shared_ptr<SimParams>
+SimParams::from_import(ImportData const& data, SPConstParticles particle_params)
+{
+    return SimParams::from_import(
+        data, particle_params, FieldDriverOptions{}.max_substeps);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Construct with simulation options.

diff --git a/src/celeritas/track/SimParams.hh b/src/celeritas/track/SimParams.hh
@@ -41,10 +41,15 @@ class SimParams final : public ParamsDataInterface<SimParamsData>
     };
 
   public:
-    // Construct with imported data
+    // Construct with imported data and default max field substeps
     static std::shared_ptr<SimParams>
     from_import(ImportData const&, SPConstParticles);
 
+    // Construct with imported data
+    static std::shared_ptr<SimParams> from_import(ImportData const&,
+                                                  SPConstParticles,
+                                                  short int max_field_substeps);
+
     // Construct with simulation input data
     explicit SimParams(Input const&);