diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index 539c2620a9..025aa180fe 100644 --- a/doc/src/pair_granular.rst +++ b/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 `, 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 `, 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. ---------- diff --git a/src/GRANULAR/gran_sub_mod.cpp b/src/GRANULAR/gran_sub_mod.cpp index 734f8f7f3d..c312d11b48 100644 --- a/src/GRANULAR/gran_sub_mod.cpp +++ b/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); } diff --git a/src/GRANULAR/gran_sub_mod.h b/src/GRANULAR/gran_sub_mod.h index fbd8c5bb26..d93c7147cc 100644 --- a/src/GRANULAR/gran_sub_mod.h +++ b/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; diff --git a/src/GRANULAR/gran_sub_mod_rolling.cpp b/src/GRANULAR/gran_sub_mod_rolling.cpp index d39fb75744..867b69cf9d 100644 --- a/src/GRANULAR/gran_sub_mod_rolling.cpp +++ b/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); } } diff --git a/src/GRANULAR/gran_sub_mod_tangential.cpp b/src/GRANULAR/gran_sub_mod_tangential.cpp index fc1938928b..d4ea2b79fe 100644 --- a/src/GRANULAR/gran_sub_mod_tangential.cpp +++ b/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); }