From 4e9e279dd9bf13481ebd24fb440217bafab73b79 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Mon, 18 May 2015 18:04:29 -0400 Subject: [PATCH] mostly functional example code --- src/compute_tally_stress.cpp | 161 +++++++++++++++++++++++++++++++++-- src/compute_tally_stress.h | 16 +++- 2 files changed, 166 insertions(+), 11 deletions(-) diff --git a/src/compute_tally_stress.cpp b/src/compute_tally_stress.cpp index a033e689e0..c5fd61a0c3 100644 --- a/src/compute_tally_stress.cpp +++ b/src/compute_tally_stress.cpp @@ -16,6 +16,7 @@ #include "atom.h" #include "pair.h" #include "update.h" +#include "memory.h" #include "error.h" #include "force.h" @@ -28,9 +29,18 @@ ComputeTallyStress::ComputeTallyStress(LAMMPS *lmp, int narg, char **arg) : { if (narg < 3) error->all(FLERR,"Illegal compute tally/stress command"); scalar_flag = 1; - extscalar = 0; - peflag = 1; // we need Pair::ev_tally() to be run + vector_flag = 1; + peratom_flag = 1; timeflag = 1; + + comm_reverse = size_peratom_cols = 6; + extscalar = 0; + extvector = 0; + size_vector = 6; + peflag = 1; // we need Pair::ev_tally() to be run + nmax = -1; + stress = NULL; + vector = new double[size_vector]; } /* ---------------------------------------------------------------------- */ @@ -38,6 +48,7 @@ ComputeTallyStress::ComputeTallyStress(LAMMPS *lmp, int narg, char **arg) : ComputeTallyStress::~ComputeTallyStress() { if (force->pair) force->pair->del_tally_callback(this); + delete[] vector; } /* ---------------------------------------------------------------------- */ @@ -49,14 +60,150 @@ void ComputeTallyStress::init() else force->pair->add_tally_callback(this); - did_compute = 0; + did_compute = -1; } /* ---------------------------------------------------------------------- */ -void ComputeTallyStress::pair_tally_callback(int i, int j, int, int, - double, double, double, - double, double, double) +void ComputeTallyStress::pair_tally_callback(int i, int j, int nlocal, int newton, + double, double, double fpair, + double dx, double dy, double dz) { - did_compute = 1; + const int ntotal = atom->nlocal + atom->nghost; + + // do setup work that needs to be done only once per timestep + + if (did_compute != update->ntimestep) { + did_compute = update->ntimestep; + + // grow local stress array if necessary + // needs to be atom->nmax in length + + if (atom->nmax > nmax) { + memory->destroy(stress); + nmax = atom->nmax; + memory->create(stress,nmax,size_peratom_cols,"tally/stress:stress"); + array_atom = stress; + } + + // clear storage as needed + + if (newton) + memset(&stress[0][0],sizeof(double)*size_peratom_cols*nmax,0); + else + memset(&stress[0][0],sizeof(double)*size_peratom_cols*nlocal,0); + + memset(virial,sizeof(double)*6,0); + } + + fpair *= 0.5; + const double v0 = dx*dx*fpair; + const double v1 = dy*dy*fpair; + const double v2 = dz*dz*fpair; + const double v3 = dx*dy*fpair; + const double v4 = dx*dz*fpair; + const double v5 = dy*dz*fpair; + + + if (newton || i < nlocal) { + virial[0] += v0; stress[i][0] += v0; + virial[1] += v1; stress[i][1] += v1; + virial[2] += v2; stress[i][2] += v2; + virial[3] += v3; stress[i][3] += v3; + virial[4] += v4; stress[i][4] += v4; + virial[5] += v5; stress[i][5] += v5; + } + if (newton || j < nlocal) { + virial[0] += v0; stress[j][0] += v0; + virial[1] += v1; stress[j][1] += v1; + virial[2] += v2; stress[j][2] += v2; + virial[3] += v3; stress[j][3] += v3; + virial[4] += v4; stress[j][4] += v4; + virial[5] += v5; stress[j][5] += v5; + } } + +/* ---------------------------------------------------------------------- */ + +int ComputeTallyStress::pack_reverse_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + buf[m++] = stress[i][0]; + buf[m++] = stress[i][1]; + buf[m++] = stress[i][2]; + buf[m++] = stress[i][3]; + buf[m++] = stress[i][4]; + buf[m++] = stress[i][5]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTallyStress::unpack_reverse_comm(int n, int *list, double *buf) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + stress[j][0] += buf[m++]; + stress[j][1] += buf[m++]; + stress[j][2] += buf[m++]; + stress[j][3] += buf[m++]; + stress[j][4] += buf[m++]; + stress[j][5] += buf[m++]; + } +} + +/* ---------------------------------------------------------------------- */ + +double ComputeTallyStress::compute_scalar() +{ + invoked_scalar = update->ntimestep; + if (update->eflag_global != invoked_scalar) + error->all(FLERR,"Energy was not tallied on needed timestep"); + + scalar = (virial[0]+virial[1]+virial[2])/3.0; + return scalar; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTallyStress::compute_vector() +{ + invoked_vector = update->ntimestep; + if (update->eflag_global != invoked_vector) + error->all(FLERR,"Energy was not tallied on needed timestep"); + + for (int i=0; i < 6; ++i) + vector[i] = virial[i]; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeTallyStress::compute_peratom() +{ + invoked_peratom = update->ntimestep; + if (update->eflag_global != invoked_peratom) + error->all(FLERR,"Energy was not tallied on needed timestep"); + + if (force->newton_pair) + comm->reverse_comm_compute(this); +} + + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeTallyStress::memory_usage() +{ + double bytes = nmax*6 * sizeof(double); + return bytes; +} + diff --git a/src/compute_tally_stress.h b/src/compute_tally_stress.h index c3504b6aba..30618a6c97 100644 --- a/src/compute_tally_stress.h +++ b/src/compute_tally_stress.h @@ -25,6 +25,7 @@ ComputeStyle(tally/stress,ComputeTallyStress) namespace LAMMPS_NS { class ComputeTallyStress : public Compute { + public: ComputeTallyStress(class LAMMPS *, int, char **); virtual ~ComputeTallyStress(); @@ -32,16 +33,23 @@ class ComputeTallyStress : public Compute { void init(); void setup() { did_compute = -1;} - double compute_scalar() { return 0.0; } - void compute_vector() {} - void compute_peratom() {} - + double compute_scalar(); + void compute_vector(); + void compute_peratom(); + + int pack_reverse_comm(int, int, double *); + void unpack_reverse_comm(int, int *, double *); + double memory_usage(); + void pair_tally_callback(int, int, int, int, double, double, double, double, double, double); private: bigint did_compute; + int nmax; + double **stress; + double virial[6]; }; }