Lb: per particle gamma

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -24,6 +24,7 @@
 #include "errorhandling.hpp"
 #include "grid.hpp"
 #include "random.hpp"
+#include "thermostat.hpp"
 
 #include "grid_based_algorithms/lb_interface.hpp"
 #include "grid_based_algorithms/lb_interpolation.hpp"
@@ -122,6 +123,23 @@ Utils::Vector3d lb_particle_coupling_drift_vel_offset(const Particle &p) {
   return vel_offset;
 }
 
+double get_particle_lb_gamma(Particle const &p) {
+  auto gamma = lb_lbcoupling_get_gamma();
+
+#ifdef THERMOSTAT_PER_PARTICLE
+  // override default if particle-specific gamma
+#ifdef PARTICLE_ANISOTROPY
+  // lb gamma not anisotropic
+  auto const particle_gamma = p.gamma()[0];
+#else
+  auto const particle_gamma = p.gamma();
+#endif
+  gamma = (p.gamma() >= Thermostat::GammaType{}) ? particle_gamma : gamma;
+#endif
+
+  return gamma;
+}
+
 Utils::Vector3d lb_drag_force(Particle const &p,
                               Utils::Vector3d const &shifted_pos,
                               Utils::Vector3d const &vel_offset) {
@@ -131,9 +149,11 @@ Utils::Vector3d lb_drag_force(Particle const &p,
       lb_lbinterpolation_get_interpolated_velocity(shifted_pos) *
       LB::get_lattice_speed();
 
+  auto const gamma = get_particle_lb_gamma(p);
+
   Utils::Vector3d v_drift = interpolated_u + vel_offset;
   /* calculate viscous force (eq. (9) @cite ahlrichs99a) */
-  return -lb_lbcoupling_get_gamma() * (p.v() - v_drift);
+  return -gamma * (p.v() - v_drift);
 }
 
 template <class T, std::size_t N>
@@ -319,28 +339,33 @@ void lb_lbcoupling_calc_particle_lattice_ia(bool couple_virtual,
       }
       using Utils::sqr;
       auto const kT = LB::get_kT() * sqr(LB::get_lattice_speed());
+
       /* Eq. (16) @cite ahlrichs99a.
        * The factor 12 comes from the fact that we use random numbers
        * from -0.5 to 0.5 (equally distributed) which have variance 1/12.
        * time_step comes from the discretization.
        */
-      auto const noise_amplitude =
-          (kT > 0.)
-              ? std::sqrt(12. * 2. * lb_lbcoupling_get_gamma() * kT / time_step)
-              : 0.0;
 
       std::unordered_set<int> coupled_ghost_particles;
 
       /* Couple particles ranges */
       for (auto &p : particles) {
         if (should_be_coupled(p, coupled_ghost_particles)) {
+          auto const noise_amplitude =
+              (kT > 0.) ? std::sqrt(12. * 2. * get_particle_lb_gamma(p) * kT /
+                                    time_step)
+                        : 0.0;
           couple_particle(p, couple_virtual, noise_amplitude,
                           lb_particle_coupling.rng_counter_coupling, time_step);
         }
       }
 
       for (auto &p : more_particles) {
         if (should_be_coupled(p, coupled_ghost_particles)) {
+          auto const noise_amplitude =
+              (kT > 0.) ? std::sqrt(12. * 2. * get_particle_lb_gamma(p) * kT /
+                                    time_step)
+                        : 0.0;
           couple_particle(p, couple_virtual, noise_amplitude,
                           lb_particle_coupling.rng_counter_coupling, time_step);
         }

testsuite/python/lb_thermostat.py

@@ -29,82 +29,165 @@
 distribution.
 """
 
-KT = 0.9
+KT = 0.2
 AGRID = 0.8
 node_volume = AGRID**3
-VISC = 6
-DENS = 1.7
+KIN_VISC = 0.3
+DENS = 0.8
 TIME_STEP = 0.005
-GAMMA = 2
+
 LB_PARAMS = {'agrid': AGRID,
              'density': DENS,
-             'viscosity': VISC,
+             'viscosity': KIN_VISC,
              'tau': TIME_STEP,
              'kT': KT,
              'seed': 123}
 
 
 class LBThermostatCommon(thermostats_common.ThermostatsCommon):
 
-    """Check the LB thermostat."""
+    """
+    Check the LB thermostat.
+    All temperature checks rely on https://en.wikipedia.org/wiki/Thermal_velocity
+    with the root-mean-squared-velocity.
+    The temperature check together with the friction test 
+    ensure that the noise amplitude is correct.
+    """
 
     system = espressomd.System(box_l=[AGRID * 12] * 3)
     system.time_step = TIME_STEP
     system.cell_system.skin = 0.4 * AGRID
+    # use "small" particles such that radius=hydrodynamic radius
+    global_gamma = 0.1 * AGRID * 6 * np.pi * KIN_VISC * DENS
+    partcl_gamma = 0.05 * AGRID * 6 * np.pi * KIN_VISC * DENS
+    # relaxation time 0.2 sim time units
+    partcl_mass = 0.2 * partcl_gamma 
+
     np.random.seed(41)
 
     def setUp(self):
         self.lbf = self.lb_class(**LB_PARAMS, **self.lb_params)
         self.system.actors.add(self.lbf)
-        self.system.thermostat.set_lb(LB_fluid=self.lbf, seed=5, gamma=GAMMA)
+        self.system.thermostat.set_lb(LB_fluid=self.lbf, seed=5, gamma=self.global_gamma)
 
     def tearDown(self):
         self.system.actors.clear()
         self.system.thermostat.turn_off()
+        self.system.part.clear()
+
+    def get_lb_kT(self, lbf):
+        nodes_mass = lbf[:, :, :].density * node_volume
+        nodes_vel_sq = np.sum(np.square(lbf[:, :, :].velocity), axis=3)    
+        return np.mean(nodes_mass * nodes_vel_sq) / 3.
 
-    def get_lb_momentum(self, lbf):
-        nodes_den = lbf[:, :, :].density
-        nodes_vel = np.sum(np.square(lbf[:, :, :].velocity), axis=3)
-        return np.multiply(nodes_den, nodes_vel)
+    def get_lb_velocity(self, lbf):
+        return np.mean(lbf[:, :, :].velocity, axis=(0, 1, 2))
 
     def test_fluid(self):
         self.system.integrator.run(100)
-        fluid_temps = []
+        fluid_kTs = []
         for _ in range(100):
-            lb_momentum = self.get_lb_momentum(self.lbf)
-            fluid_temps.append(np.average(lb_momentum) * node_volume)
+            fluid_kTs.append(self.get_lb_kT(self.lbf))
             self.system.integrator.run(3)
 
-        fluid_temp = np.average(fluid_temps) / 3
-        self.assertAlmostEqual(fluid_temp, KT, delta=0.05)
+        fluid_kT = np.mean(fluid_kTs)
+        np.testing.assert_allclose(fluid_kT, KT, rtol=0.05)
 
-    def test_with_particles(self):
-        self.system.part.add(
-            pos=np.random.random((100, 3)) * self.system.box_l)
-        self.system.integrator.run(120)
-        N = len(self.system.part)
-        loops = 250
-        v_particles = np.zeros((loops, N, 3))
-        fluid_temps = []
+    def check_partcl_temp(self, partcl_vel):
+        partcl_vel_rms = np.sqrt(np.mean(np.linalg.norm(partcl_vel, axis=2)**2))
 
-        for i in range(loops):
-            self.system.integrator.run(3)
-            if i % 10 == 0:
-                lb_momentum = self.get_lb_momentum(self.lbf)
-                fluid_temps.append(np.average(lb_momentum) * node_volume)
-            v_particles[i] = self.system.part.all().v
-        fluid_temp = np.average(fluid_temps) / 3.
+        np.testing.assert_allclose(np.mean(partcl_vel), 0, atol=0.05 * partcl_vel_rms)
+        np.testing.assert_allclose(partcl_vel_rms**2 * self.partcl_mass / 3., 
+                                   KT, rtol=0.05)
 
-        np.testing.assert_allclose(np.average(v_particles), 0, atol=0.033)
-        np.testing.assert_allclose(np.var(v_particles), KT, atol=0.033)
-
-        minmax = 3
+        vel_range = 2 * partcl_vel_rms
         n_bins = 7
-        error_tol = 0.016
         self.check_velocity_distribution(
-            v_particles.reshape((-1, 3)), minmax, n_bins, error_tol, KT)
-
-        np.testing.assert_allclose(fluid_temp, KT, atol=5e-3)
+            partcl_vel.reshape((-1, 3)), vel_range, n_bins, 0.016, KT, mass=self.partcl_mass)
+
+    def test_temperature_with_particles(self):
+        n_partcls_per_type = 50
+        partcls_global_gamma = self.system.part.add(
+            pos=np.random.random((n_partcls_per_type, 3)) * self.system.box_l,
+            mass=n_partcls_per_type * [self.partcl_mass])
+        partcls_per_part_gamma = self.system.part.add(
+            pos=np.random.random((n_partcls_per_type, 3)) * self.system.box_l,
+            mass=n_partcls_per_type * [self.partcl_mass],
+            gamma=n_partcls_per_type * [self.partcl_gamma])
+        self.system.integrator.run(50)
+        n_samples = 250
+        vel_global_gamma = np.zeros((n_samples, n_partcls_per_type, 3))
+        vel_per_part_gamma = np.zeros_like(vel_global_gamma)
+        fluid_kTs = []
+
+        for i in range(n_samples):
+            self.system.integrator.run(3)
+            if i % 10 == 0:
+                fluid_kTs.append(self.get_lb_kT(self.lbf))
+            vel_global_gamma[i] = partcls_global_gamma.v
+            vel_per_part_gamma[i] = partcls_per_part_gamma.v
+        fluid_kT = np.mean(fluid_kTs)
+
+        self.check_partcl_temp(vel_global_gamma)
+        self.check_partcl_temp(vel_per_part_gamma)
+
+        np.testing.assert_allclose(fluid_kT, KT, rtol=0.02)
+
+    def test_friction(self):
+        """apply force and measure if the average velocity matches expectation"""
+
+        # large force to get better signal  to noise ratio
+        ext_force = np.array([1.2, 2.1, 1.1])
+        n_partcls_each_type = 50
+        partcls_global_gamma = self.system.part.add(
+            pos=np.random.random((n_partcls_each_type, 3)) * self.system.box_l,
+            ext_force=n_partcls_each_type * [ext_force],
+            mass=n_partcls_each_type * [self.partcl_mass])
+        partcls_per_partcl_gamma = self.system.part.add(
+            pos=np.random.random((n_partcls_each_type, 3)) * self.system.box_l,
+            ext_force=n_partcls_each_type * [ext_force],
+            mass=n_partcls_each_type * [self.partcl_mass],
+            gamma=n_partcls_each_type * [self.partcl_gamma])
+
+        # add counterforce to fluid such that velocities cannot increase indefinitely
+        # and we get a steady state relative velocity
+        total_force_partcls = ext_force * 2 * n_partcls_each_type
+        lb_volume = np.prod(self.system.box_l)
+        self.lbf.ext_force_density = -total_force_partcls / lb_volume
+
+        # let particles and fluid accelerate towards steady state before measuring
+        self.system.integrator.run(int(0.8 / self.system.time_step))
+        startpos_global_gamma = np.copy(partcls_global_gamma.pos)
+        startpos_per_p_gamma = np.copy(partcls_per_partcl_gamma.pos)     
+
+        run_time = 0.3
+        n_steps = int(run_time / self.system.time_step)
+        steps_per_sample = 10
+        lb_vels = []
+        for _ in range(int(n_steps / steps_per_sample)):
+            self.system.integrator.run(steps_per_sample)
+            lb_vels.append(self.get_lb_velocity(self.lbf))
+        lb_vel = np.mean(lb_vels, axis=0)
+
+        vel_global_gamma = np.mean(
+            partcls_global_gamma.pos - startpos_global_gamma,
+            axis=0) / run_time
+        vel_per_part_gamma = np.mean(
+            partcls_per_partcl_gamma.pos - startpos_per_p_gamma,
+            axis=0) / run_time
+        # get velocity relative to the fluid
+        vel_global_gamma -= lb_vel
+        vel_per_part_gamma -= lb_vel
+
+        # average over 3 spatial directions (isotropic friction)
+        global_gamma_measured = np.mean(ext_force / vel_global_gamma)
+        per_partcl_gamma_measured = np.mean(ext_force / vel_per_part_gamma)
+
+        np.testing.assert_allclose(global_gamma_measured, self.global_gamma, rtol=0.05)
+        np.testing.assert_allclose(
+            per_partcl_gamma_measured,
+            self.partcl_gamma,
+            rtol=0.05)
 
 
 @utx.skipIfMissingFeatures(["WALBERLA"])

testsuite/python/thermostats_common.py

@@ -23,31 +23,31 @@
 import espressomd.observables
 
 
-def single_component_maxwell(x1, x2, kT):
+def single_component_maxwell(x1, x2, kT, mass=1.):
     """Integrate the probability density from x1 to x2 using the trapezoidal rule"""
     x = np.linspace(x1, x2, 1000)
-    return np.trapz(np.exp(-x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT)
+    maxwell_distr = np.exp(- mass * x**2 / (2. * kT)) * np.sqrt(mass / (2. * np.pi * kT))
+    return np.trapz(maxwell_distr, x=x)
 
 
 class ThermostatsCommon:
 
     """Tests the velocity distribution created by a thermostat."""
 
-    def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT):
+    def check_velocity_distribution(self, vel, minmax, n_bins, error_tol, kT, mass=1.):
         """Check the recorded particle distributions in velocity against a
            histogram with n_bins bins. Drop velocities outside minmax. Check
            individual histogram bins up to an accuracy of error_tol against
            the analytical result for kT."""
         for i in range(3):
-            hist = np.histogram(
+            hist, bin_edges = np.histogram(
                 vel[:, i], range=(-minmax, minmax), bins=n_bins, density=False)
-            data = hist[0] / float(vel.shape[0])
-            bins = hist[1]
+            hist = hist / len(vel)
             expected = []
             for j in range(n_bins):
                 expected.append(single_component_maxwell(
-                    bins[j], bins[j + 1], kT))
-            np.testing.assert_allclose(data[:n_bins], expected, atol=error_tol)
+                    bin_edges[j], bin_edges[j + 1], kT, mass=mass))
+            np.testing.assert_allclose(hist, expected, atol=error_tol)
 
     def test_00_verify_single_component_maxwell(self):
         """Verifies the normalization of the analytical expression."""