documents abstract_particle_rhs

cpp/particles/abstract_particle_rhs.hpp

@@ -28,7 +28,17 @@
 
 #include "particles/interpolation/abstract_particle_set.hpp"
 
-// children of this class must provide a method to compute the rhs for a given particle set
+/** Particle model definition, i.e. provides right-hand-side of ODEs
+ *
+ * This abstract class prescribes the way that particle models shall be
+ * implemented in TurTLE.
+ * Children of this class must implement a "function call" operator,
+ * where the right-hand-side of the corresponding ODE shall be computed,
+ * as well as an "imposeModelConstraints" method where miscellaneous properties
+ * can be enforced as needed (an example would be orientation vectors being
+ * reset to unit lenght).
+ *
+ */
 class abstract_particle_rhs
 {
     protected:
@@ -37,21 +47,54 @@ class abstract_particle_rhs
         // destructor
         virtual ~abstract_particle_rhs(){}
 
-        // a way to probe the dimension of the ODE
+        /** Probe the dimension of the ODE
+         */
         virtual const int getStateSize() const = 0;
 
 
-        // important bit
-        // may resort particles, therefore it takes the
-        // parameter `additional_data` to possibly resort other arrays
+        /** Compute right-hand-side of the ODE
+         *
+         * Notes
+         * -----
+         *  1. This method may shuffle particle data in-memory, but
+         *  it may not redistribute particles between MPI processes.
+         *  The `additional_data` parameter allows for other particle data to be
+         *  shuffled in the same fashion.
+         *  In particular, you may use `additional_data` when computing substeps
+         *  for Runge Kutta methods to keep the different temporary arrays in
+         *  the same order.
+         *  Please see implementation of the Heun method for an example.
+         *
+         *  For completeness: shuffling of data will at least take place during
+         *  the computation of particle-to-particle interaction terms.
+         *
+         *  2. This method may not modify the current state of the particles.
+         *  I cannot mark the parameter `abstract_particle_set &pset` as
+         *  `const` because the arrays may need to be shuffled, but separate
+         *  parts of the code rely on the fact that this method does not change
+         *  the state of the particle set.
+         *  If your particle model has strict constraints that affect the
+         *  computation of the right-hand side, you must ensure that those
+         *  constraints are respected within this method *without* modifying
+         *  the current particle state.
+         *  Otherwise the ODE may not be integrated correctly.
+         */
         virtual int operator()(
                 double t,
                 abstract_particle_set &pset,
                 particle_rnumber *result,
-                std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> &additional_data) = 0;
-
-        // for post-timestep actions such as
-        // normalization of orientation vector
+                std::vector<std::unique_ptr<
+                    abstract_particle_set::particle_rnumber[]>> &additional_data) = 0;
+        /** Enforce model constraints on a set of particles
+         *
+         * Note:
+         * This is meant to be used as a safe-guard against growth of small
+         * errors.
+         * As an example, shaped particles for which orientation must be tracked
+         * through time may slowly develop non-unit-sized orientation vectors.
+         * This method will change the particle states such that any such
+         * constraints are respected.
+         */
         virtual int imposeModelConstraints(
                 abstract_particle_set &pset) const = 0;
 };