ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

clean-up and test example

Browse Source
This commit is contained in:
vyas
2024-11-18 04:42:17 -06:00
parent 5a826c67f6
commit db2f5dc407
8 changed files with 62 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

20
examples/granular/3ps.dat Normal file
View File

@ -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

30
examples/granular/in.sync_verlet Normal file
View File

@ -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

3
src/GRANULAR/gran_sub_mod.cpp
View File

@ -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];

3
src/GRANULAR/gran_sub_mod.h
View File

@ -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;

28
src/GRANULAR/gran_sub_mod_rolling.cpp
View File

@ -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);

88
src/GRANULAR/gran_sub_mod_tangential.cpp
View File

@ -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);
}

1
src/GRANULAR/gran_sub_mod_tangential.h
View File

@ -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; }

3
src/GRANULAR/granular_model.cpp
View File

@ -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];