fix: use clamped timestep in ALP last-pulse iteration

  Clamp dt = min(coverageTimeStep, pulseTime - time) and push it to the
  surface model via a new setTimeStep hook to non-const dt_ so the final iteration
  integrates over the actual remaining pulse time instead of overshooting
  with the full coverageTimeStep.

Author: philipphaslhofer

include/viennaps/models/psSingleParticleALD.hpp

@@ -12,7 +12,7 @@ template <typename NumericType>
 class SingleParticleALDSurfaceModel : public SurfaceModel<NumericType> {
   using SurfaceModel<NumericType>::coverages;
 
-  const NumericType dt_;
+  NumericType dt_;
   const NumericType gpc_;
   const NumericType s0_;
 
@@ -37,6 +37,8 @@ class SingleParticleALDSurfaceModel : public SurfaceModel<NumericType> {
     coverages->insertNextScalarData(cov, "Coverage");
   }
 
+  void setTimeStep(NumericType dt) override { dt_ = dt; }
+
   SmartPointer<std::vector<NumericType>>
   calculateVelocities(SmartPointer<viennals::PointData<NumericType>> rates,
                       const std::vector<Vec3D<NumericType>> &coordinates,

include/viennaps/process/psALPStrategy.hpp

@@ -170,6 +170,10 @@ class ALPStrategy final : public ProcessStrategy<NumericType, D> {
           return ProcessResult::USER_INTERRUPTED;
 #endif
 
+        // Clamp last step to land exactly on pulseTime
+        double dt = std::min(coverageTimeStep, pulseTime - time);
+        context.model->getSurfaceModel()->setTimeStep(dt);
+
         // Calculate fluxes
         auto fluxes = SmartPointer<viennals::PointData<NumericType>>::New();
         PROCESS_CHECK(fluxEngine_->calculateFluxes(context, fluxes));
@@ -179,7 +183,7 @@ class ALPStrategy final : public ProcessStrategy<NumericType, D> {
 
         outputIntermediateResults(context, fluxes, pulseIteration);
 
-        time += coverageTimeStep;
+        time += dt;
         pulseIteration++;
 
         if (Logger::hasInfo()) {

include/viennaps/process/psSurfaceModel.hpp

@@ -32,6 +32,8 @@ template <typename NumericType> class SurfaceModel {
     // if no surface data get initialized here, they won't be used at all
   }
 
+  virtual void setTimeStep(NumericType dt) {}
+
   virtual SmartPointer<std::vector<NumericType>>
   calculateVelocities(SmartPointer<viennals::PointData<NumericType>> fluxes,
                       const std::vector<Vec3D<NumericType>> &coordinates,