Updated rotation and documentation

doc/src/pair_granular.rst

@@ -670,7 +670,7 @@ attractive force. This keyword cannot be used with the JKR or DMT models.
 
 ----------
 
-The standard velocity-Verlet integration scheme's half-step staggering of position and velocity can introduce inaccuracies in frictional tangential force calculations, resulting in unphysical kinematics in certain systems. These effects are particularly pronounced in polydisperse frictional flows characterized by large-to-small size ratios exceeding three. The *synchronized_verlet* flag implements an improved Velocity-Verlet integration scheme, as detailed in :ref:`Vyas et al <Vyas2025>`, that synchronizes position and velocity updates for force evaluation. By refining tangential force calculations, the *synchronized_verlet* method ensures physically consistent results without significantly impacting computational cost.
+The standard velocity-Verlet integration scheme's half-step staggering of position and velocity can introduce inaccuracies in frictional tangential force calculations, resulting in unphysical kinematics in certain systems. These effects are particularly pronounced in polydisperse frictional flows characterized by large-to-small size ratios exceeding three. The *synchronized_verlet* flag implements an alternate Velocity-Verlet integration scheme, as detailed in :ref:`Vyas et al <Vyas2025>`, that synchronizes position and velocity updates for force evaluation. By refining tangential force calculations, the *synchronized_verlet* method ensures physically consistent results without significantly impacting computational cost.
 
 ----------
 

src/GRANULAR/gran_sub_mod.cpp

@@ -137,7 +137,7 @@ double GranSubMod::mix_mean(double val1, double val2)
   rotate-rescale vector v so it is perpendicular to unit vector n
   and has the same magnitude as before
   ---------------------------------------------------------------------- */
-void GranSubMod::rotate_rescale_vec(double *v, double *n, double *rotated_v)
+void GranSubMod::rotate_rescale_vec(double *v, double *n)
 {
   double rsht, shrmag, prjmag, temp_dbl, temp_array[3];
 
@@ -154,5 +154,4 @@ void GranSubMod::rotate_rescale_vec(double *v, double *n, double *rotated_v)
   else
     temp_dbl = 0;
   scale3(temp_dbl, v);
-  copy3(v, rotated_v);
 }

src/GRANULAR/gran_sub_mod.h

@@ -31,7 +31,7 @@ namespace Granular_NS {
     virtual void mix_coeffs(double *, double *);
     virtual void coeffs_to_local(){};
     virtual void init(){};    // called after all sub models + coeffs defined
-    void   rotate_rescale_vec(double *hislocal, double *n, double *ans);
+    void   rotate_rescale_vec(double *hislocal, double *n);
 
     void allocate_coeffs();
     std::string name;

src/GRANULAR/gran_sub_mod_rolling.cpp

@@ -85,7 +85,7 @@ void GranSubModRollingSDS::calculate_forces()
     rolldotn = dot3(hist_temp, gm->nxuse);
 
     frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit);
-    if (frameupdate) rotate_rescale_vec(hist_temp, gm->nxuse, hist_temp);
+    if (frameupdate) rotate_rescale_vec(hist_temp, gm->nxuse);
 
     // update history at half-step
     scale3(gm->dt, gm->vrl, temp_array);
@@ -95,7 +95,7 @@ void GranSubModRollingSDS::calculate_forces()
     if (gm->synchronized_verlet == 1) {
       rolldotn = dot3(hist_temp, gm->nx);
       frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit);
-      if (frameupdate) rotate_rescale_vec(hist_temp, gm->nx, hist_temp);
+      if (frameupdate) rotate_rescale_vec(hist_temp, gm->nx);
     }
   }
 

src/GRANULAR/gran_sub_mod_tangential.cpp

@@ -128,7 +128,7 @@ void GranSubModTangentialLinearHistory::calculate_forces()
     rsht = dot3(history, gm->nxuse);
     frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit);
 
-    if (frame_update) rotate_rescale_vec(history, gm->nxuse, history);
+    if (frame_update) rotate_rescale_vec(history, gm->nxuse);
 
     // update history, tangential force using velocities at half step
     // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46
@@ -139,7 +139,7 @@ void GranSubModTangentialLinearHistory::calculate_forces()
       rsht = dot3(history, gm->nx);
       frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit);
       //Second projection to nx (t+\Delta t)
-      if (frame_update) rotate_rescale_vec(history, gm->nx, history);
+      if (frame_update) rotate_rescale_vec(history, gm->nx);
     }
   }
 
@@ -147,7 +147,8 @@ void GranSubModTangentialLinearHistory::calculate_forces()
   scale3(-k, history, gm->fs);
   //Rotating vtr for damping term in nx direction
   if (frame_update && gm->synchronized_verlet == 1) {
-    rotate_rescale_vec(gm->vtr, gm->nx, vtr2);
+    copy3(gm->vtr, vtr2);
+    rotate_rescale_vec(vtr2, gm->nx);
   } else {
     copy3(gm->vtr, vtr2);
   }
@@ -328,7 +329,7 @@ void GranSubModTangentialMindlin::calculate_forces()
       frame_update = (fabs(rsht) * k_scaled) > (EPSILON * Fscrit);
     }
 
-    if (frame_update) rotate_rescale_vec(history, gm->nxuse, history);
+    if (frame_update) rotate_rescale_vec(history, gm->nxuse);
 
     // update history
     if (mindlin_force) {
@@ -350,7 +351,7 @@ void GranSubModTangentialMindlin::calculate_forces()
       } else {
         frame_update = (fabs(rsht) * k_scaled) > (EPSILON * Fscrit);
       }
-      if (frame_update) rotate_rescale_vec(history, gm->nx, history);
+      if (frame_update) rotate_rescale_vec(history, gm->nx);
     }
   }
 
@@ -358,7 +359,8 @@ void GranSubModTangentialMindlin::calculate_forces()
   //Rotating vtr for damping term in nx direction
   if (frame_update && gm->synchronized_verlet)
   {
-    rotate_rescale_vec(gm->vtr, gm->nx, vtr2);
+    copy3(gm->vtr, vtr2);
+    rotate_rescale_vec(vtr2, gm->nx);
   } else {
     copy3(gm->vtr, vtr2);
   }