From db2f5dc4076ffc329ec31e703e570d8b7b40fdff Mon Sep 17 00:00:00 2001 From: vyas Date: Mon, 18 Nov 2024 04:42:17 -0600 Subject: [PATCH] clean-up and test example --- examples/granular/3ps.dat | 20 ++++++ examples/granular/in.sync_verlet | 30 ++++++++ src/GRANULAR/gran_sub_mod.cpp | 3 +- src/GRANULAR/gran_sub_mod.h | 3 +- src/GRANULAR/gran_sub_mod_rolling.cpp | 28 -------- src/GRANULAR/gran_sub_mod_tangential.cpp | 88 ++---------------------- src/GRANULAR/gran_sub_mod_tangential.h | 1 + src/GRANULAR/granular_model.cpp | 3 + 8 files changed, 62 insertions(+), 114 deletions(-) create mode 100644 examples/granular/3ps.dat create mode 100644 examples/granular/in.sync_verlet diff --git a/examples/granular/3ps.dat b/examples/granular/3ps.dat new file mode 100644 index 0000000000..c8d8ff7f4d --- /dev/null +++ b/examples/granular/3ps.dat @@ -0,0 +1,20 @@ +LAMMPS data file + +3 atoms +2 atom types + +-0.03 0.03 xlo xhi +-0.02 0.02 ylo yhi +-0.02 0.02 zlo zhi + +Atoms # sphere + +1 2 0.0005 2500 0.005 2.7937481665024794E-05 0.0016005372526493155 0 0 0 +2 1 0.01 2500 0.01 0.0 0.0 0 0 0 +3 1 0.01 2500 0.0 0.0 0.0 0 0 0 + +Velocities + +1 0 0 0 0 0 0 +2 0 0 0 0 0 0 +3 0 0 0 0 0 0 \ No newline at end of file diff --git a/examples/granular/in.sync_verlet b/examples/granular/in.sync_verlet new file mode 100644 index 0000000000..7e2b6ca51a --- /dev/null +++ b/examples/granular/in.sync_verlet @@ -0,0 +1,30 @@ +# Example problem demonstrating the use of synchronized_verlet +# Refer https://arxiv.org/abs/2410.14798 for further details + +units si +atom_style sphere +newton off + +boundary p p f +read_data 3ps.dat + +group ps type 1 +group fine type 2 + +pair_style granular +pair_coeff * * hooke 1e4 0.5 tangential linear_history 8235 0 0.5 damping coeff_restitution #synchronized_verlet +# pair_coeff * * hertz/material 5e8 0.5 0.3 tangential mindlin NULL 0.0 0.5 damping coeff_restitution #synchronized_verlet + +timestep 1e-6 +fix frz ps freeze +fix g fine gravity 9.81 vector 0 0 -1 +fix 1 fine nve/sphere + +dump 1 all custom 3000 op.dump id x y z vx vy vz type diameter + +comm_modify vel yes + +thermo 3000 +thermo_style custom step ke +run_style verlet +run 300000 \ No newline at end of file diff --git a/src/GRANULAR/gran_sub_mod.cpp b/src/GRANULAR/gran_sub_mod.cpp index 0ac7923cba..323527d505 100644 --- a/src/GRANULAR/gran_sub_mod.cpp +++ b/src/GRANULAR/gran_sub_mod.cpp @@ -44,6 +44,7 @@ GranSubMod::GranSubMod(class GranularModel *gm, LAMMPS *lmp) : Pointers(lmp) beyond_contact = 0; num_coeffs = 0; contact_radius_flag = 0; + allow_synchronization = 0; nondefault_history_transfer = 0; transfer_history_factor = nullptr; @@ -132,7 +133,7 @@ double GranSubMod::mix_mean(double val1, double val2) return 0.5 * (val1 + val2); } -void GranSubMod::rotate_rescale_vec(double *v, double *n, double *ans) +void GranSubMod::rotate_rescale_vec(double *v, double *n, double *ans) { double rsht, shrmag, prjmag, temp_dbl, temp_array[3]; diff --git a/src/GRANULAR/gran_sub_mod.h b/src/GRANULAR/gran_sub_mod.h index 5e8101376d..fbd8c5bb26 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, double *ans); void allocate_coeffs(); std::string name; @@ -44,6 +44,7 @@ namespace Granular_NS { int beyond_contact; // If the sub model contact extends beyond overlap int allow_cohesion; // If the sub model works with a cohesive normal force int contact_radius_flag; // If the sub model requires contact radius + int allow_synchronization; // If the sub model works with synchronized verlet GranularModel *gm; diff --git a/src/GRANULAR/gran_sub_mod_rolling.cpp b/src/GRANULAR/gran_sub_mod_rolling.cpp index 413839a0ca..674716f7df 100644 --- a/src/GRANULAR/gran_sub_mod_rolling.cpp +++ b/src/GRANULAR/gran_sub_mod_rolling.cpp @@ -85,20 +85,6 @@ void GranSubModRollingSDS::calculate_forces() frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit); if (frameupdate) rotate_rescale_vec(hist_temp, gm->nxuse, hist_temp); - // { // rotate into tangential plane at half-step - // rollmag = len3(hist_temp); - // // projection - // scale3(rolldotn, gm->nxuse, temp_array); - // sub3(hist_temp, temp_array, hist_temp); - - // // also rescale to preserve magnitude - // prjmag = len3(hist_temp); - // if (prjmag > 0) - // scalefac = rollmag / prjmag; - // else - // scalefac = 0; - // scale3(scalefac, hist_temp); - // } // update history at half-step scale3(gm->dt, gm->vrl, temp_array); @@ -108,20 +94,6 @@ void GranSubModRollingSDS::calculate_forces() rolldotn = dot3(hist_temp, gm->nx); frameupdate = (fabs(rolldotn) * k) > (EPSILON * Frcrit); if (frameupdate) rotate_rescale_vec(hist_temp, gm->nx, hist_temp); - // { // rotate into tangential plane at full-step - // rollmag = len3(hist_temp); - // // projection - // scale3(rolldotn, gm->nx, temp_array); - // sub3(hist_temp, temp_array, hist_temp); - - // // also rescale to preserve magnitude - // prjmag = len3(hist_temp); - // if (prjmag > 0) - // scalefac = rollmag / prjmag; - // else - // scalefac = 0; - // scale3(scalefac, hist_temp); - // } } scaleadd3(-k, hist_temp, -gamma, gm->vrl, gm->fr); diff --git a/src/GRANULAR/gran_sub_mod_tangential.cpp b/src/GRANULAR/gran_sub_mod_tangential.cpp index 7e7cdebbb5..2740bc9c7e 100644 --- a/src/GRANULAR/gran_sub_mod_tangential.cpp +++ b/src/GRANULAR/gran_sub_mod_tangential.cpp @@ -95,6 +95,7 @@ GranSubModTangentialLinearHistory::GranSubModTangentialLinearHistory(GranularMod { num_coeffs = 3; size_history = 3; + allow_synchronization = 1; } /* ---------------------------------------------------------------------- */ @@ -128,20 +129,6 @@ void GranSubModTangentialLinearHistory::calculate_forces() frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit); if (frame_update) rotate_rescale_vec(history, gm->nxuse, history); - // shrmag = len3(history); - - // // first projection to half step normal - // scale3(rsht, gm->nxuse, temp_array); - // sub3(history, temp_array, history); - - // // also rescale to preserve magnitude - // prjmag = len3(history); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, history); - // } // update history, tangential force using velocities at half step // see e.g. eq. 18 of Thornton et al, Pow. Tech. 2013, v223,p30-46 @@ -150,47 +137,17 @@ void GranSubModTangentialLinearHistory::calculate_forces() if(gm->synchronized_verlet == 1) { rsht = dot3(history, gm->nx); - frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit); + frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit); //Second projection to nx (t+\Delta t) if (frame_update) rotate_rescale_vec(history, gm->nx, history); } - // rsht = dot3(history, gm->nx); - // frame_update = (fabs(rsht) * k) > (EPSILON * Fscrit); - - // if (frame_update) { - // shrmag = len3(history); - - // // second projection to full step normal - // scale3(rsht, gm->nx, temp_array); - // sub3(history, temp_array, history); - - // // also rescale to preserve magnitude - // prjmag = len3(history); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, history); - // } } // tangential forces = history + tangential velocity damping scale3(-k, history, gm->fs); //Rotating vtr for damping term in nx direction - // rsht = dot3(gm->vtr, gm->nx); if (frame_update && gm->synchronized_verlet == 1) { rotate_rescale_vec(gm->vtr, gm->nx, vtr2); - // shrmag = len3(gm->vtr); - // scale3(rsht, gm->nx, temp_array); - // sub3(gm->vtr, temp_array, vtr2); - - // // also rescale to preserve magnitude - // prjmag = len3(vtr2); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, vtr2); } else { copy3(gm->vtr, vtr2); } @@ -299,6 +256,7 @@ GranSubModTangentialMindlin::GranSubModTangentialMindlin(GranularModel *gm, LAMM mindlin_force = 0; mindlin_rescale = 0; contact_radius_flag = 1; + allow_synchronization = 1; } /* ---------------------------------------------------------------------- */ @@ -371,19 +329,6 @@ void GranSubModTangentialMindlin::calculate_forces() } if (frame_update) rotate_rescale_vec(history, gm->nxuse, history); - // { - // shrmag = len3(history); - // // projection - // scale3(rsht, gm->nxuse, temp_array); - // sub3(history, temp_array, history); - // // also rescale to preserve magnitude - // prjmag = len3(history); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, history); - // } // update history if (mindlin_force) { @@ -406,38 +351,13 @@ void GranSubModTangentialMindlin::calculate_forces() } if (frame_update) rotate_rescale_vec(history, gm->nx, history); - // { - // shrmag = len3(history); - // // projection - // scale3(rsht, gm->nx, temp_array); - // sub3(history, temp_array, history); - // // also rescale to preserve magnitude - // prjmag = len3(history); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, history); - // } } // tangential forces = history + tangential velocity damping //Rotating vtr for damping term in nx direction - // rsht = dot3(gm->vtr, gm->nx); - if (frame_update && gm->synchronized_verlet) + if (frame_update && gm->synchronized_verlet) { rotate_rescale_vec(gm->vtr, gm->nx, vtr2); - // shrmag = len3(gm->vtr); - // scale3(rsht, gm->nx, temp_array); - // sub3(gm->vtr, temp_array, vtr2); - - // // also rescale to preserve magnitude - // prjmag = len3(vtr2); - // if (prjmag > 0) - // temp_dbl = shrmag / prjmag; - // else - // temp_dbl = 0; - // scale3(temp_dbl, vtr2); } else { copy3(gm->vtr, vtr2); } diff --git a/src/GRANULAR/gran_sub_mod_tangential.h b/src/GRANULAR/gran_sub_mod_tangential.h index ca699fbf1e..0c807138f2 100644 --- a/src/GRANULAR/gran_sub_mod_tangential.h +++ b/src/GRANULAR/gran_sub_mod_tangential.h @@ -38,6 +38,7 @@ namespace Granular_NS { GranSubModTangential(class GranularModel *, class LAMMPS *); virtual void calculate_forces() = 0; + int get_synchronization_flag() const { return allow_synchronization; } double get_k() const { return k; } double get_damp() const { return damp; } double get_mu() const { return mu; } diff --git a/src/GRANULAR/granular_model.cpp b/src/GRANULAR/granular_model.cpp index 9d75f8769b..a0eb0d26df 100644 --- a/src/GRANULAR/granular_model.cpp +++ b/src/GRANULAR/granular_model.cpp @@ -271,6 +271,9 @@ void GranularModel::init() if (limit_damping && normal_model->get_cohesive_flag()) error->all(FLERR,"Cannot limit damping with a cohesive normal model, {}", normal_model->name); + if ((synchronized_verlet == 1) && (tangential_model->get_synchronization_flag()==0)) + error->all(FLERR,"Cannot use synchronized verlet with a non-synchronized tangential model, {}", tangential_model->name); + if (nondefault_history_transfer) { transfer_history_factor = new double[size_history];