core: Do the gamma check only once per tiem step

src/core/grid_based_algorithms/lb_particle_coupling.cpp

@@ -128,23 +128,13 @@ static Thermostat::GammaType lb_handle_particle_anisotropy(Particle const &p) {
 #ifdef THERMOSTAT_PER_PARTICLE
   auto const &partcl_gamma = p.gamma();
 #ifdef PARTICLE_ANISOTROPY
-  /*
-  lb does (at the moment) not support rotational particle coupling.
-  Consequently, anisotropic particles are also not supported.
-  */
-  auto const aniso_flag = (partcl_gamma[0] != partcl_gamma[1]) ||
-                          (partcl_gamma[1] != partcl_gamma[2]);
-  if (aniso_flag) {
-    runtimeErrorMsg() << "anisotropic particle (id " << p.id()
-                      << ") coupled to LB.";
-  }
   auto const default_gamma = Thermostat::GammaType::broadcast(lb_gamma);
 #else
   auto const default_gamma = lb_gamma;
 #endif // PARTICLE_ANISOTROPY
   return handle_particle_gamma(partcl_gamma, default_gamma);
 #else
-  return Thermostat::GammaType{lb_gamma};
+  return lb_gamma;
 #endif // THERMOSTAT_PER_PARTICLE
 }
 
@@ -294,16 +284,40 @@ bool CouplingBookkeeping::is_ghost_for_local_particle(Particle const &p) const {
   return not ::cell_structure.get_local_particle(p.id())->is_ghost();
 }
 
+#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
+static bool lb_coupling_sanity_checks(Particle const &p) {
+  /*
+  lb does (at the moment) not support rotational particle coupling.
+  Consequently, anisotropic particles are also not supported.
+  */
+  auto const &p_gamma = p.gamma();
+  if (p_gamma[0] != p_gamma[1] or p_gamma[1] != p_gamma[2]) {
+    runtimeErrorMsg() << "anisotropic particle (id " << p.id()
+                      << ") coupled to LB.";
+    return true;
+  }
+  return false;
+}
+#endif
+
 void couple_particles(bool couple_virtual, ParticleRange const &real_particles,
                       ParticleRange const &ghost_particles, double time_step) {
   ESPRESSO_PROFILER_CXX_MARK_FUNCTION;
   if (lattice_switch == ActiveLB::WALBERLA_LB) {
     if (lb_particle_coupling.couple_to_md) {
       ParticleCoupling coupling{couple_virtual, time_step};
       CouplingBookkeeping bookkeeping{};
+#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
+      bool do_sanity_checks = true;
+#endif
       for (auto const &particle_range : {real_particles, ghost_particles}) {
         for (auto &p : particle_range) {
           if (bookkeeping.should_be_coupled(p)) {
+#if defined(THERMOSTAT_PER_PARTICLE) and defined(PARTICLE_ANISOTROPY)
+            if (do_sanity_checks and lb_coupling_sanity_checks(p)) {
+              do_sanity_checks = false;
+            }
+#endif
             coupling.kernel(p);
           }
         }

testsuite/python/lb_thermostat.py

@@ -96,10 +96,7 @@ def test_fluid(self):
 
     def check_partcl_temp(self, partcl_vel):
         partcl_vel_rms = np.sqrt(
-            np.mean(
-                np.linalg.norm(
-                    partcl_vel,
-                    axis=2)**2))
+            np.mean(np.linalg.norm(partcl_vel, axis=2)**2))
 
         np.testing.assert_allclose(
             np.mean(partcl_vel), 0, atol=0.05 * partcl_vel_rms)
@@ -109,7 +106,8 @@ def check_partcl_temp(self, partcl_vel):
         vel_range = 2 * partcl_vel_rms
         n_bins = 7
         self.check_velocity_distribution(
-            partcl_vel.reshape((-1, 3)), vel_range, n_bins, 0.02, KT, mass=self.partcl_mass)
+            partcl_vel.reshape((-1, 3)), vel_range, n_bins, 0.02, KT,
+            mass=self.partcl_mass)
 
     def test_temperature_with_particles(self):
         n_partcls_per_type = 50
@@ -155,8 +153,8 @@ def test_friction(self):
             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
+        # 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
@@ -199,6 +197,13 @@ def test_friction(self):
             self.partcl_gamma,
             rtol=0.05)
 
+    @utx.skipIfMissingFeatures(["PARTICLE_ANISOTROPY",
+                               "THERMOSTAT_PER_PARTICLE"])
+    def test_exceptions(self):
+        self.system.part.add(pos=[0., 0., 0.], gamma=[1., 2., 3.], id=2)
+        with self.assertRaisesRegex(Exception, r"ERROR: anisotropic particle \(id 2\) coupled to LB"):
+            self.system.integrator.run(1)
+
 
 @utx.skipIfMissingFeatures(["WALBERLA"])
 class LBWalberlaThermostat(LBThermostatCommon, ut.TestCase):