small performance optimization

src/mesa_pd/kernel/VBondModel/IntegratedVDWContact.h

@@ -53,11 +53,17 @@ public:
                      const size_t p_idx2,
                      Accessor &ac) const;
 
-   static constexpr double eps_ = 0.07124;
-   static constexpr double A_ = 0.0223;
-   static constexpr double B_ = 1.31;
-   static constexpr double alf_ = 9.5;
-   static constexpr double bet_ = 4.0;
+   static constexpr real_t R_CNT = 6.78_r; ///< CNT radius
+   static constexpr real_t T = 2_r * R_CNT; ///< Height of a cylindrical segment
+   static constexpr real_t eps = 0.07124_r;
+   static constexpr real_t A = 0.0223_r;
+   static constexpr real_t B = 1.31_r;
+   static constexpr real_t alpha = 9.5_r;
+   static constexpr real_t beta = 4.0_r;
+
+   // vdW adhesion + linear repulsion potential.
+   const real_t r0 = R_CNT * (2_r + std::pow((alpha * A / (beta * B)), 1_r / (alpha - beta)));
+   const real_t u0 = 4_r * eps * (A / std::pow(r0 / R_CNT - 2_r, (alpha)) - B / std::pow(r0 / R_CNT - 2_r, (beta)));
 };
 
 /**
@@ -74,19 +80,15 @@ real_t IntegratedVDWContact::operator()(const size_t p_idx1,
                                         const size_t p_idx2,
                                         Accessor &ac) const
 {
-   constexpr real_t K_n = 200.0; // Good for fabrics modeling
+   constexpr real_t K_n = 200.0_r; // Good for fabrics modeling
 
    // particle centers
    Vec3 O1 = ac.getPosition(p_idx1);
    Vec3 O2 = ac.getPosition(p_idx2);
 
    // axial vectors
-   Vec3 b1 = ac.getRotation(p_idx1) * Vec3(1.0, 0.0, 0.0);
-   Vec3 b2 = ac.getRotation(p_idx2) * Vec3(1.0, 0.0, 0.0);
-
-   // vdW adhesion + linear repulsion potential.
-   const real_t r0 = R_CNT * (2 + std::pow((alf_ * A_ / (bet_ * B_)), 1 / (alf_ - bet_)));
-   const real_t u0 = 4 * eps_ * (A_ / pow(r0 / R_CNT - 2, (alf_)) - B_ / pow(r0 / R_CNT - 2, (bet_)));
+   Vec3 b1 = ac.getRotation(p_idx1).getMatrix() * Vec3(1_r, 0_r, 0_r);
+   Vec3 b2 = ac.getRotation(p_idx2).getMatrix() * Vec3(1_r, 0_r, 0_r);
 
    real_t U = 0, dU = 0; ///< total potential
    Vec3 force12(0), dforce12(0);  // Force 12
@@ -95,37 +97,37 @@ real_t IntegratedVDWContact::operator()(const size_t p_idx1,
    Vec3 moment21(0), dmoment21(0); // Total torque 21
 
    constexpr int Np = 5; // Number of integration points over each axis
-   constexpr double Npinv = 1.0 / Np;
+   constexpr real_t Npinv = 1.0_r / real_t(Np);
    for (int i = 0; i < Np; ++i) // integral dl1
    {
       for (int j = 0; j < Np; ++j) // integral dl2
       {
          // Levers
-         Vec3 l1 = (-0.5 * T + (0.5 + i) * T * Npinv) * b1;
-         Vec3 l2 = (-0.5 * T + (0.5 + j) * T * Npinv) * b2;
+         Vec3 l1 = (-0.5_r * T + (0.5_r + real_t(i)) * T * Npinv) * b1;
+         Vec3 l2 = (-0.5_r * T + (0.5_r + real_t(j)) * T * Npinv) * b2;
 
          /// radius vector between dl1 and dl2
          Vec3 R12 = (O2 + l2) - (O1 + l1);
 
-         Vec3 n12 = R12.getNormalized();
          real_t r12 = R12.length();
+         Vec3 n12 = R12 * (1_r / r12);
 
-         if (r12 < r0) // elastic interaction
+         if (r12 < r0)
          {
+            // elastic interaction
             dU = u0 + K_n * (r12 - r0) * (r12 - r0) * Npinv * Npinv;
-            dforce12 = n12 * K_n * ((r12 - r0)) * Npinv * Npinv;
+            dforce12 = n12 * K_n * (r12 - r0) * Npinv * Npinv;
             dforce21 = -dforce12;
-         }
-         if (r12 >= r0) // vdW interaction
+         } else
          {
-            const double powAlpha1 = pow(r12 / R_CNT - 2.0, (alf_ + 1.0));
-            const double powBeta1 = pow(r12 / R_CNT - 2.0, (bet_ + 1));
-            dU = 4.0 * eps_ * (A_ / pow(r12 / R_CNT - 2.0, (alf_)) - B_ / pow(r12 / R_CNT - 2.0, (bet_))) * Npinv *
-                 Npinv;
-            dforce12 =
-                  n12 * 4.0 * eps_ * (-(alf_ * A_) / powAlpha1 + (bet_ * B_) / powBeta1) / R_CNT * Npinv * Npinv;
-            dforce21 =
-                  -n12 * 4.0 * eps_ * (-(alf_ * A_) / powAlpha1 + (bet_ * B_) / powBeta1) / R_CNT * Npinv * Npinv;
+            // vdW interaction
+            const real_t normDistance = r12 / R_CNT - 2_r;
+            const real_t powAlpha = std::pow(normDistance, alpha);
+            const real_t powBeta = std::pow(normDistance, beta);
+            dU = 4_r * eps * (A / powAlpha - B / powBeta) * Npinv * Npinv;
+            dforce12 = n12 * 4_r * eps / R_CNT * Npinv * Npinv *
+                  (-(alpha * A) / (powAlpha * normDistance) + (beta * B) / (powBeta * normDistance));
+            dforce21 = -dforce12;
          }
 
          dmoment12 = l2 % dforce12;
@@ -144,15 +146,6 @@ real_t IntegratedVDWContact::operator()(const size_t p_idx1,
    addTorqueAtomic(p_idx1, ac,  -moment21);
    addTorqueAtomic(p_idx2, ac,  -moment12);
 
-   // Viscous addition
-   constexpr real_t dampFactor = 1052.0; // Calibrated in accordance with 2014 JAM
-   real_t damp = beta_ * dampFactor;
-   Vec3 relVelocity = ac.getLinearVel(p_idx1) - ac.getLinearVel(p_idx2);
-   Vec3 visc_force = -damp * relVelocity;
-   //WALBERLA_LOG_DEVEL( "Visc force: = " << visc_force );
-   addForceAtomic(p_idx1, ac,  visc_force);
-   addForceAtomic(p_idx2, ac, -visc_force);
-
    return U; // vdW adhesion energy
 }