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

Adding compute fabric

Browse Source
This commit is contained in:
jtclemm
2021-06-02 14:17:15 -06:00
parent 96ac2dc9f6
commit 479622e57d
4 changed files with 634 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

149
doc/src/compute_fabric.rst Normal file
View File

@ -0,0 +1,149 @@
.. index:: compute fabric
compute fabric command
==============================
Syntax
""""""
.. parsed-literal::
compute ID group-ID fabric cutoff attribute1 attribute2 ...
* ID, group-ID are documented in :doc:`compute <compute>` command
* fabric = style name of this compute command
* cutoff = *type* or *radius*
.. parsed-literal::
*type* = cutoffs determined for the types of the two atoms
*radius* = cutoffs determined based on finite size of atoms
* one or more attributes may be appended
.. parsed-literal::
*contact* = contact tensor
*branch* = branch tensor
*force/normal* = normal force tensor
*force/tangential* = tangential force tensor
Examples
""""""""
.. code-block:: LAMMPS
compute 1 all fabric type contact force/normal
compute 1 all fabric radius force/normal force/tangential
Description
"""""""""""
Define a computation that calculates various fabric tensors for pairwise
interactions. The *type* and *radius* settings are used to select whether
the compute will be used used with regular interactions with cutoffs
determined by atom types or with granular interactions with interaction
lengths determined by particle radii, respectively.
Four fabric tensors can be computed: a contact, branch, normal force, or
tangential force tensor. The contact tensor is calculated as
.. math::
C_{ab} = \frac{15}{2} (\phi_{ab} - \mathrm{tr}(\phi) \delta_{ab})
where :math:`a` and :math:`b` are the :math:`x`, :math:`y`, :math:`z`
directions, :math:`\delta_{ab}` is the Kronecker delta function, and
the tensor :math:`\phi` is defined as
.. math::
\phi_{ab} = \sum_{n = 1}^{N_p} \frac{r_{a} r_{b}}{r^2}
where :math:`N_p` is the number of pair interactions in the simulation,
:math:`r_{a}` is the :math:`a` component of the radial vector between the
two particles, and :math:`r` is the magnitude of the radial vector.
The branch tensor is calculated as
.. math::
B_{ab} = \frac{15}{6 tr(D)} (D_{ab} - tr(D) \delta_{ab})
where the tensor :math:`D` is defined as
.. math::
D_{ab} = \sum_{n = 1}^{N_p}
\frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
\frac{r_{a} r_{b}}{r}
where :math:`N_c` is the total number of contacts in the system and the subscripts
:math:`c` and :math:`d` are summed according to Einstein notation.
The normal force fabric tensor is calculated as
.. math::
F^n_{ab} = \frac{15}{6 tr(N)} (N_{ab} - tr(N) \delta_{ab})
where the tensor :math:`N` is defined as
.. math::
N_{ab} = \sum_{n = 1}^{N_p}
\frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
\frac{r_{a} r_{b}}{r^2} f_n
where :math:`f_n` is the magnitude of the normal, central-body force between the two atoms.
Finally, the tangential force fabric tensor is only defined for pair styles that
apply tangential forces to particles, namely granular pair styles. It is calculated
as
.. math::
F^t_{ab} = \frac{15}{9 tr(N)} (T_{ab} - tr(T) \delta_{ab})
where the tensor :math:`T` is defined as
.. math::
T_{ab} = \sum_{n = 1}^{N_p}
\frac{1}{N_c (r^2 + C_{cd} r_c r_d)}
\frac{r_{a} r_{b}}{r^2} f_t
where :math:`f_t` is the magnitude of the tagential force between the two atoms.
Output info
"""""""""""
This compute calculates a local vector of doubles. The vector stores the
unique components of the first requested tensor in the order
xx, yy, zz, xy, xz, yz followed by the same components for all
subsequent tensors. The length of the vector is therefore six times
the number of requested tensors.
Restrictions
""""""""""""
Currently, compute *fabric* does not support pair styles
with many-body interactions. It also does not
support models with long-range Coulombic or dispersion forces,
i.e. the kspace_style command in LAMMPS. It also does not support the
following fixes which add rigid-body constraints: :doc:`fix shake
<fix_shake>`, :doc:`fix rattle <fix_shake>`, :doc:`fix rigid
<fix_rigid>`, :doc:`fix rigid/small <fix_rigid>`. It does not support
granular pair styles that extend beyond the contact (e.g. JKR and DMT).
Related commands
""""""""""""""""
none
Default
"""""""
none

2
src/.gitignore vendored
View File

@ -379,6 +379,8 @@
/compute_erotate_rigid.h
/compute_event_displace.cpp
/compute_event_displace.h
/compute_fabric.cpp
/compute_fabric.h
/compute_fep.cpp
/compute_fep.h
/compute_force_tally.cpp

407
src/GRANULAR/compute_fabric.cpp Normal file
View File

@ -0,0 +1,407 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "compute_fabric.h"
#include <cmath>
#include <cstring>
#include <mpi.h>
#include "atom.h"
#include "update.h"
#include "force.h"
#include "pair.h"
#include "modify.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "neigh_list.h"
#include "error.h"
using namespace LAMMPS_NS;
enum{OTHER,GRANULAR};
enum{TYPE,RADIUS};
enum{CN,BR,FN,FT};
/* ---------------------------------------------------------------------- */
ComputeFabric::ComputeFabric(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg),
tensor_style(NULL)
{
if (narg < 3) error->all(FLERR,"Illegal compute fabric command");
if (strcmp(arg[3],"type") == 0) cutstyle = TYPE;
else if (strcmp(arg[3],"radius") == 0) cutstyle = RADIUS;
else error->all(FLERR,"Illegal compute fabric command");
if (cutstyle == RADIUS && !atom->radius_flag)
error->all(FLERR,"Compute fabric radius style requires atom attribute radius");
ntensors = narg - 4;
tensor_style = new int[ntensors];
cn_flag = 0;
br_flag = 0;
fn_flag = 0;
ft_flag = 0;
ntensors = 0;
int iarg = 4;
while (iarg < narg) {
if (strcmp(arg[iarg],"contact") == 0) {
cn_flag = 1;
tensor_style[ntensors++] = CN;
} else if (strcmp(arg[iarg],"branch") == 0) {
br_flag = 1;
tensor_style[ntensors++] = BR;
} else if (strcmp(arg[iarg],"force/normal") == 0) {
fn_flag = 1;
tensor_style[ntensors++] = FN;
} else if (strcmp(arg[iarg],"force/tangential") == 0) {
ft_flag = 1;
tensor_style[ntensors++] = FT;
} else error->all(FLERR,"Illegal compute fabric command");
iarg++;
}
vector_flag = 1;
size_vector = 1+ntensors*6;
extvector = 0;
vector = new double[size_vector];
}
/* ---------------------------------------------------------------------- */
ComputeFabric::~ComputeFabric()
{
delete [] vector;
delete [] tensor_style;
}
/* ---------------------------------------------------------------------- */
void ComputeFabric::init()
{
if (force->pair == NULL)
error->all(FLERR,"No pair style is defined for compute fabric");
if (force->pair->single_enable == 0 && (fn_flag || ft_flag))
error->all(FLERR,"Pair style does not support compute fabric normal or tangential force");
// Find if granular or gran
pstyle = OTHER;
if (force->pair_match("granular",0) ||
force->pair_match("gran/hooke",0) ||
force->pair_match("gran/hertz",0) ||
force->pair_match("gran/hooke/history",0) ||
force->pair_match("gran/hertz/history",0)) pstyle = GRANULAR;
if (pstyle != GRANULAR && ft_flag)
error->all(FLERR,"Pair style does not have tangential forces for compute fabric");
if(force->pair->beyond_contact)
error->all(FLERR, "Compute fabric does not support pair styles that extend beyond contact");
// need an occasional half neighbor list
// set size to same value as request made by force->pair
// this should enable it to always be a copy list (e.g. for granular pstyle)
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->compute = 1;
neighbor->requests[irequest]->occasional = 1;
NeighRequest *pairrequest = neighbor->find_request((void *) force->pair);
if (pairrequest) neighbor->requests[irequest]->size = pairrequest->size;
}
/* ---------------------------------------------------------------------- */
void ComputeFabric::init_list(int /*id*/, NeighList *ptr)
{
list = ptr;
}
/* ---------------------------------------------------------------------- */
void ComputeFabric::compute_vector()
{
int i,j,ii,jj,inum,jnum,itype,jtype;
tagint itag,jtag;
double xtmp,ytmp,ztmp,delx,dely,delz;
double r,rinv,rsq,radsum,eng,fpair;
double nx, ny, nz, fx, fy ,fz;
double ncinv, denom, fn, ft, prefactor;
double cn_tensor[6], br_tensor[6], ft_tensor[6], fn_tensor[6];
double trace_phi, trace_D, trace_Xfn, trace_Xft;
double phi_ij[6] = {0.0};
double Ac_ij[6] = {0.0};
double D_ij[6] = {0.0};
double Xfn_ij[6] = {0.0};
double Xft_ij[6] = {0.0};
int nc;
int temp_int;
double temp_dbl[6];
int *ilist,*jlist,*numneigh,**firstneigh;
double **x = atom->x;
double *radius = atom->radius;
tagint *tag = atom->tag;
int *type = atom->type;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list);
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
Pair *pair = force->pair;
double **cutsq = force->pair->cutsq;
// First loop through and calculate contact fabric tensor
nc = 0;
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itag = tag[i];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (cutstyle == TYPE) {
if (rsq >= cutsq[itype][jtype]) continue;
} else {
radsum = radius[i] + radius[j];
if (rsq >= radsum*radsum) continue;
}
nc += 1;
r = sqrt(rsq);
rinv = 1.0/r;
nx = delx*rinv;
ny = dely*rinv;
nz = delz*rinv;
phi_ij[0] += nx*nx;
phi_ij[1] += ny*ny;
phi_ij[2] += nz*nz;
phi_ij[3] += nx*ny;
phi_ij[4] += nx*nz;
phi_ij[5] += ny*nz;
}
}
//Count total contacts across processors
MPI_Allreduce(&nc,&temp_int,1,MPI_INT,MPI_SUM,world);
nc = temp_int;
// If no contacts, everything will be zero
if (nc == 0) {
for (i = 0; i < size_vector; i ++) vector[i] = 0.0;
return;
}
vector[ntensors*6] = nc;
ncinv = 1.0/nc;
//Sum phi across processors
MPI_Allreduce(phi_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) phi_ij[i] = temp_dbl[i]*ncinv;
trace_phi = (1.0/3.0)*(phi_ij[0]+phi_ij[1]+phi_ij[2]);
Ac_ij[0] = (15.0/2.0)*(phi_ij[0] - trace_phi);
Ac_ij[1] = (15.0/2.0)*(phi_ij[1] - trace_phi);
Ac_ij[2] = (15.0/2.0)*(phi_ij[2] - trace_phi);
Ac_ij[3] = (15.0/2.0)*(phi_ij[3]);
Ac_ij[4] = (15.0/2.0)*(phi_ij[4]);
Ac_ij[5] = (15.0/2.0)*(phi_ij[5]);
// If needed, loop through and calculate other fabric tensors
if(br_flag || fn_flag || ft_flag) {
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
if (!(mask[i] & groupbit)) continue;
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
itag = tag[i];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
if (!(mask[j] & groupbit)) continue;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
jtype = type[j];
if (cutstyle == TYPE) {
if (rsq >= cutsq[itype][jtype]) continue;
} else {
radsum = radius[i] + radius[j];
if (rsq >= radsum*radsum) continue;
}
if (fn_flag || ft_flag)
eng = pair->single(i,j,itype,jtype,rsq,1.0,1.0,fpair);
r = sqrt(rsq);
rinv = 1.0/r;
nx = delx*rinv;
ny = dely*rinv;
nz = delz*rinv;
denom = 1 + Ac_ij[0]*nx*nx + Ac_ij[1]*ny*ny + Ac_ij[2]*nz*nz;
denom += 2*Ac_ij[3]*nx*ny + 2*Ac_ij[4]*nx*nz + 2*Ac_ij[5]*ny*nz;
prefactor = ncinv/denom;
if (br_flag) {
D_ij[0] += prefactor*nx*nx*r;
D_ij[1] += prefactor*ny*ny*r;
D_ij[2] += prefactor*nz*nz*r;
D_ij[3] += prefactor*nx*ny*r;
D_ij[4] += prefactor*nx*nz*r;
D_ij[5] += prefactor*ny*nz*r;
}
if (fn_flag || ft_flag) {
fn = r*fpair;
Xfn_ij[0] += prefactor*nx*nx*fn;
Xfn_ij[1] += prefactor*ny*ny*fn;
Xfn_ij[2] += prefactor*nz*nz*fn;
Xfn_ij[3] += prefactor*nx*ny*fn;
Xfn_ij[4] += prefactor*nx*nz*fn;
Xfn_ij[5] += prefactor*ny*nz*fn;
if (ft_flag) {
ft = force->pair->svector[3];
Xft_ij[0] += prefactor*nx*nx*ft;
Xft_ij[1] += prefactor*ny*ny*ft;
Xft_ij[2] += prefactor*nz*nz*ft;
Xft_ij[3] += prefactor*nx*ny*ft;
Xft_ij[4] += prefactor*nx*nz*ft;
Xft_ij[5] += prefactor*ny*nz*ft;
}
}
}
}
}
// Output results
if (cn_flag) {
for (i = 0; i < ntensors; i ++) {
if (tensor_style[i] == CN) {
for (j = 0; j < 6; j++) vector[6*i + j] = Ac_ij[j];
}
}
}
if (br_flag) {
MPI_Allreduce(D_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) D_ij[i] = temp_dbl[i];
trace_D = (1.0/3.0)*(D_ij[0]+D_ij[1]+D_ij[2]);
br_tensor[0] = (15.0/(6.0*trace_D))*(D_ij[0] - trace_D);
br_tensor[1] = (15.0/(6.0*trace_D))*(D_ij[1] - trace_D);
br_tensor[2] = (15.0/(6.0*trace_D))*(D_ij[2] - trace_D);
br_tensor[3] = (15.0/(6.0*trace_D))*(D_ij[3]);
br_tensor[4] = (15.0/(6.0*trace_D))*(D_ij[4]);
br_tensor[5] = (15.0/(6.0*trace_D))*(D_ij[5]);
for (i = 0; i < ntensors; i ++) {
if (tensor_style[i] == BR) {
for (j = 0; j < 6; j++) vector[6*i + j] = br_tensor[j];
}
}
}
if (fn_flag || ft_flag) {
MPI_Allreduce(Xfn_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) Xfn_ij[i] = temp_dbl[i];
trace_Xfn = (1.0/3.0)*(Xfn_ij[0]+Xfn_ij[1]+Xfn_ij[2]);
}
if (fn_flag) {
fn_tensor[0] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[0] - trace_Xfn);
fn_tensor[1] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[1] - trace_Xfn);
fn_tensor[2] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[2] - trace_Xfn);
fn_tensor[3] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[3]);
fn_tensor[4] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[4]);
fn_tensor[5] = (15.0/(6.0*trace_Xfn))*(Xfn_ij[5]);
for (i = 0; i < ntensors; i ++) {
if (tensor_style[i] == FN) {
for (j = 0; j < 6; j++) vector[6*i + j] = fn_tensor[j];
}
}
}
if (ft_flag) {
MPI_Allreduce(Xft_ij,temp_dbl,6,MPI_DOUBLE,MPI_SUM,world);
for (i = 0; i < 6; i++) Xft_ij[i] = temp_dbl[i];
trace_Xft = (1.0/3.0)*(Xft_ij[0]+Xft_ij[1]+Xft_ij[2]);
ft_tensor[0] = (15.0/(9.0*trace_Xfn))*(Xft_ij[0] - trace_Xft);
ft_tensor[1] = (15.0/(9.0*trace_Xfn))*(Xft_ij[1] - trace_Xft);
ft_tensor[2] = (15.0/(9.0*trace_Xfn))*(Xft_ij[2] - trace_Xft);
ft_tensor[3] = (15.0/(9.0*trace_Xfn))*(Xft_ij[3]);
ft_tensor[4] = (15.0/(9.0*trace_Xfn))*(Xft_ij[4]);
ft_tensor[5] = (15.0/(9.0*trace_Xfn))*(Xft_ij[5]);
for (i = 0; i < ntensors; i ++) {
if (tensor_style[i] == FT) {
for (j = 0; j < 6; j++) vector[6*i + j] = ft_tensor[j];
}
}
}
}

76
src/GRANULAR/compute_fabric.h Normal file
View File

@ -0,0 +1,76 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
ComputeStyle(fabric,ComputeFabric)
#else
#ifndef LMP_COMPUTE_FABRIC_H
#define LMP_COMPUTE_FABRIC_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeFabric : public Compute {
public:
ComputeFabric(class LAMMPS *, int, char **);
~ComputeFabric();
void init();
void init_list(int, class NeighList *);
void compute_vector();
private:
int ntensors, pstyle, cutstyle;
int *tensor_style;
class NeighList *list;
int cn_flag, br_flag, fn_flag, ft_flag;
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute fabric radius style requires atom attribute radius
Self-explanatory.
E: No pair style is defined for compute fabric
Self-explanatory.
E: Pair style does not support compute fabric normal or tangential force
Pair style must be single enabled to calculate the normal or tangential force tensors
E: Pair style does not have tangential forces for compute fabric
The tangential force tensor can only be calculated for granular pair styles
E: Compute fabric does not support pair styles that extend beyond contact
Granular pair styles that extend beyond contact such as JKR or DMT are not supported
*/