From d5e5f590e83ab943c15058ebaae5eb30a7564d2c Mon Sep 17 00:00:00 2001 From: jtclemm Date: Tue, 22 Jun 2021 14:14:25 -0600 Subject: [PATCH] Updating documentation, adding compute scalar --- doc/src/compute_fabric.rst | 21 +++-- src/GRANULAR/compute_fabric.cpp | 134 +++++++++++++++++++++++++++----- src/GRANULAR/compute_fabric.h | 2 + 3 files changed, 132 insertions(+), 25 deletions(-) diff --git a/doc/src/compute_fabric.rst b/doc/src/compute_fabric.rst index 1f01b1cbc1..3ded3d7adc 100644 --- a/doc/src/compute_fabric.rst +++ b/doc/src/compute_fabric.rst @@ -48,9 +48,15 @@ Description """"""""""" Define a compute that calculates various fabric tensors for pairwise -interactions :ref:`(Ouadfel) `. The *type* and *radius* settings -are used to select whether interactions cutoffs are determined by atom types -or by the sum of atomic radii (atom style sphere), respectively. +interaction :ref:`(Ouadfel) `. Fabric tensors are commonly used +to quantify the anisotropy or orientation of granular contacts but can also +be used to characterize the direction of pairwise interactions in general +systems. The *type* and *radius* settings are used to select whether interactions +cutoffs are determined by atom types or by the sum of atomic radii (atom +style sphere), respectively. Calling this compute is roughly the cost of a +pair style invocation as it involves a loop over the neighbor list. If the +normal or tangential force tensors are requested, it will be more expensive +than a pair style invocation as it will also recalculate all pair forces. Four fabric tensors are available: the contact, branch, normal force, or tangential force tensor. The contact tensor is calculated as @@ -137,17 +143,18 @@ Multiple *type/include* keywords may be added. Output info """"""""""" -This compute calculates a local vector of doubles. The vector stores the +This compute calculates a local vector of doubles and a scalar. 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 final entry is the +followed by the same components for all subsequent tensors. The length of the vector +is therefore six times the number of requested tensors. The scalar output is the number of pairwise interactions included in the calculation of the fabric tensor. -The length of the vector is therefore six times the number of requested tensors plus one. Restrictions """""""""""" This fix is part of the GRANULAR package. It is only enabled if LAMMPS -was built with that package. See the :doc:`Build package ` doc page for more info. +was built with that package. See the :doc:`Build package ` +doc page for more info. Currently, compute *fabric* does not support pair styles with many-body interactions. It also does not diff --git a/src/GRANULAR/compute_fabric.cpp b/src/GRANULAR/compute_fabric.cpp index aeb887c795..548c29eb10 100644 --- a/src/GRANULAR/compute_fabric.cpp +++ b/src/GRANULAR/compute_fabric.cpp @@ -130,9 +130,12 @@ ComputeFabric::ComputeFabric(LAMMPS *lmp, int narg, char **arg) : } vector_flag = 1; - size_vector = 1 + ntensors * 6; + size_vector = ntensors * 6; extvector = 0; + scalar_flag = 1; + extscalar = 0; + vector = new double[size_vector]; } @@ -186,6 +189,8 @@ void ComputeFabric::init_list(int /*id*/, NeighList *ptr) void ComputeFabric::compute_vector() { + invoked_vector = update->ntimestep; + int i, j, ii, jj, inum, jnum, itype, jtype; tagint itag, jtag; double xtmp, ytmp, ztmp, delx, dely, delz; @@ -200,7 +205,6 @@ void ComputeFabric::compute_vector() double D_ij[6] = {0.0}; double Xfn_ij[6] = {0.0}; double Xft_ij[6] = {0.0}; - int nc, temp_int; double temp_dbl[6]; int *ilist, *jlist, *numneigh, **firstneigh; @@ -224,8 +228,17 @@ void ComputeFabric::compute_vector() Pair *pair = force->pair; double **cutsq = force->pair->cutsq; + // invoke compute_scalar() to update the number of contacts, if needed + nc = compute_scalar(); + + // If no contacts, everything will be zero + if (nc == 0) { + for (i = 0; i < size_vector; i++) vector[i] = 0.0; + return; + } + ncinv = 1.0 / nc; + // First loop through and calculate contact fabric tensor - nc = 0; for (ii = 0; ii < inum; ii++) { i = ilist[ii]; if (!(mask[i] & groupbit)) continue; @@ -278,8 +291,6 @@ void ComputeFabric::compute_vector() if (rsq >= radsum * radsum) continue; } - nc += 1; - r = sqrt(rsq); rinv = 1.0 / r; nx = delx * rinv; @@ -295,19 +306,6 @@ void ComputeFabric::compute_vector() } } - //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; @@ -495,3 +493,103 @@ void ComputeFabric::compute_vector() } } } + +/* ---------------------------------------------------------------------- */ + +double ComputeFabric::compute_scalar() +{ + // Skip if already calculated on this timestep + if (invoked_scalar == update->ntimestep) return nc; + + invoked_scalar = update->ntimestep; + + int i, j, ii, jj, inum, jnum, itype, jtype; + tagint itag, jtag; + double xtmp, ytmp, ztmp, delx, dely, delz; + double rsq, radsum, temp_dbl; + + 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; + + 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; + + // itag = jtag is possible for long cutoffs that include images of self + + if (newton_pair == 0 && j >= nlocal) { + jtag = tag[j]; + if (itag > jtag) { + if ((itag + jtag) % 2 == 0) continue; + } else if (itag < jtag) { + if ((itag + jtag) % 2 == 1) continue; + } else { + if (x[j][2] < ztmp) continue; + if (x[j][2] == ztmp) { + if (x[j][1] < ytmp) continue; + if (x[j][1] == ytmp && x[j][0] < xtmp) continue; + } + } + } + + jtype = type[j]; + + if (type_filter) + if (type_filter[itype][jtype] == 0) continue; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + rsq = delx * delx + dely * dely + delz * delz; + + if (cutstyle == TYPE) { + if (rsq >= cutsq[itype][jtype]) continue; + } else { + radsum = radius[i] + radius[j]; + if (rsq >= radsum * radsum) continue; + } + + nc += 1.0; + } + } + + //Count total contacts across processors + MPI_Allreduce(&nc, &temp_dbl, 1, MPI_DOUBLE, MPI_SUM, world); + nc = temp_dbl; + + return nc; +} + diff --git a/src/GRANULAR/compute_fabric.h b/src/GRANULAR/compute_fabric.h index 6a2d9f4eb3..5e20ea6562 100644 --- a/src/GRANULAR/compute_fabric.h +++ b/src/GRANULAR/compute_fabric.h @@ -31,9 +31,11 @@ class ComputeFabric : public Compute { void init(); void init_list(int, class NeighList *); void compute_vector(); + double compute_scalar(); private: int ntensors, pstyle, cutstyle; + double nc; int *tensor_style; int **type_filter; class NeighList *list;