Merge pull request #69 from stochasticHydroTools/deterministic_drift

Deterministic thermal drift tests

Author: RaulPPelaez

include/MobilityInterface/random_finite_differences.h

@@ -31,8 +31,7 @@ template <typename Mdot>
 void random_finite_differences(Mdot mdot, device_span<const real> positions,
                                device_span<real> ilinear,
                                device_span<real> iangular, uint seed,
-                               real prefactor = 1) {
-  constexpr real delta = 1e-3;
+                               const real delta, real prefactor = 1) {
   const uint N = ilinear.size() / 3;
   using device_vector =
       thrust::device_vector<real, allocator::thrust_cached_allocator<real>>;

solvers/DPStokes/extra/uammd_interface.h

@@ -31,6 +31,7 @@ struct PyParameters {
   real zmin, zmax;
   // Tolerance will be ignored in DP mode, TP will use only tolerance and nxy/nz
   real tolerance = 1e-5;
+  real delta = 1e-3; // RFD step size
   real w, w_d;
   real hydrodynamicRadius = -1;
   real3 beta = {-1.0, -1.0, -1.0};

solvers/DPStokes/mobility.h

@@ -222,8 +222,9 @@ class DPStokes : public libmobility::Mobility {
     uint seed = rng();
     device_vector thermal_drift_m(ilinear.size(), 0);
     device_vector thermal_drift_d(iangular.size(), 0);
-    libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m,
-                                           thermal_drift_d, seed, prefactor);
+    libmobility::random_finite_differences(
+        mdot, original_pos, thermal_drift_m, thermal_drift_d, seed,
+        this->dppar.delta * this->dppar.hydrodynamicRadius, prefactor);
     this->setPositions(original_pos);
     thrust::transform(thrust::cuda::par, thermal_drift_m.begin(),
                       thermal_drift_m.end(), linear.begin(), linear.begin(),

solvers/DPStokes/python_wrapper.cu

@@ -22,6 +22,8 @@ zmax : float
 		The maximum value of the z coordinate. This is the position of the top wall if the Z periodicity is `two_walls`.
 allowChangingBoxSize : bool
     Whether the periodic extents Lx & Ly can be modified during parameter selection. Default: false.
+delta : float
+    The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default is 1e-3.
 )pbdoc";
 
 static const char *docstring = R"pbdoc(
@@ -35,18 +37,20 @@ The periodicity must be set to `periodic` in the X and Y directions. The Z perio
 )pbdoc";
 
 MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(
-    DPStokes, solver.def(
-                  "setParameters",
-                  [](DPStokes &self, real Lx, real Ly, real zmin, real zmax,
-                     bool allowChangingBoxSize) {
-                    DPStokesParameters params;
-                    params.Lx = Lx;
-                    params.Ly = Ly;
-                    params.zmin = zmin;
-                    params.zmax = zmax;
-                    params.allowChangingBoxSize = allowChangingBoxSize;
-                    self.setParametersDPStokes(params);
-                  },
-                  "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a,
-                  "allowChangingBoxSize"_a = false, setparameters_docstring);
+    DPStokes,
+    solver.def(
+        "setParameters",
+        [](DPStokes &self, real Lx, real Ly, real zmin, real zmax,
+           bool allowChangingBoxSize, real delta) {
+          DPStokesParameters params;
+          params.Lx = Lx;
+          params.Ly = Ly;
+          params.zmin = zmin;
+          params.zmax = zmax;
+          params.allowChangingBoxSize = allowChangingBoxSize;
+          params.delta = delta;
+          self.setParametersDPStokes(params);
+        },
+        "Lx"_a, "Ly"_a, "zmin"_a, "zmax"_a, "allowChangingBoxSize"_a = false,
+        "delta"_a = 1e-3, setparameters_docstring);
     , docstring);

solvers/NBody/mobility.h

@@ -28,6 +28,7 @@ class NBody : public libmobility::Mobility {
   nbody_rpy::algorithm algorithm = nbody_rpy::algorithm::advise;
 
   real wallHeight; // location of the wall in z
+  real delta;      // finite difference step size for random finite differences
 
   // Batched functionality configuration
   int Nbatch = -1;
@@ -53,6 +54,7 @@ class NBody : public libmobility::Mobility {
     int Nbatch = -1;
     int NperBatch = -1;
     std::optional<real> wallHeight = std::nullopt;
+    real delta = 1e-3; // in units of particle radius
   };
   /**
    * @brief Sets the parameters for the N-body computation
@@ -94,6 +96,7 @@ class NBody : public libmobility::Mobility {
           "you want to use a wall, set periodicityZ to single_wall in the "
           "configuration.");
     }
+    this->delta = par.delta;
   }
 
   void initialize(Parameters ipar) override {
@@ -195,8 +198,9 @@ class NBody : public libmobility::Mobility {
     device_vector thermal_drift_m(ilinear.size(), 0);
     device_vector thermal_drift_d(iangular.size(), 0);
 
-    libmobility::random_finite_differences(mdot, original_pos, thermal_drift_m,
-                                           thermal_drift_d, seed, prefactor);
+    libmobility::random_finite_differences(
+        mdot, original_pos, thermal_drift_m, thermal_drift_d, seed,
+        this->delta * this->hydrodynamicRadius, prefactor);
     device_adapter<real> linear(ilinear, device::cuda);
     this->setPositions(original_pos);
     thrust::transform(thrust::cuda::par, thermal_drift_m.begin(),

solvers/NBody/python_wrapper.cu

@@ -5,19 +5,21 @@
 #include <nanobind/stl/optional.h>
 
 static const char *docstringSetParameters = R"pbdoc(
-        Set the parameters for the NBody solver.
+Set the parameters for the NBody solver.
 
-        Parameters
-        ----------
-        algorithm : str
-                The algorithm to use. Options are "naive", "fast", "block" and "advise". Default is "advise".
-        NBatch : int
-                The number of batches to use. If -1 (default), the number of batches is automatically determined.
-        NperBatch : int
-                The number of particles per batch. If -1 (default), the number of particles per batch is automatically determined.
-        wallHeight : float
-                The height of the wall. Only valid if periodicityZ is single_wall.
-        )pbdoc";
+Parameters
+----------
+algorithm : str
+        The algorithm to use. Options are "naive", "fast", "block" and "advise". Default is "advise".
+NBatch : int
+        The number of batches to use. If -1 (default), the number of batches is automatically determined.
+NperBatch : int
+        The number of particles per batch. If -1 (default), the number of particles per batch is automatically determined.
+wallHeight : float
+        The height of the wall. Only valid if periodicityZ is single_wall.
+delta : float
+        The finite difference step size for random finite differences. Specified in units of hydrodynamicRadius. Default is 1e-3.
+)pbdoc";
 
 static const char *docstring = R"pbdoc(
 This module computes hydrodynamic interactions using an :math:`O(N^2)` algorithm.
@@ -50,10 +52,10 @@ MOBILITY_PYTHONIFY_WITH_EXTRA_CODE(
     solver.def(
         "setParameters",
         [](NBody &myself, std::string algo, int NBatch, int NperBatch,
-           std::optional<real> wallHeight) {
+           std::optional<real> wallHeight, real delta) {
           myself.setParametersNBody({nbody_rpy::string2NBodyAlgorithm(algo),
-                                     NBatch, NperBatch, wallHeight});
+                                     NBatch, NperBatch, wallHeight, delta});
         },
         docstringSetParameters, "algorithm"_a = "advise", "Nbatch"_a = -1,
-        "NperBatch"_a = -1, "wallHeight"_a = std::nullopt);
+        "NperBatch"_a = -1, "wallHeight"_a = std::nullopt, "delta"_a = 1e-3);
     , docstring);