From 9aa6c6c5a2f219604ddceb76f92e93515cb8b170 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Mon, 29 Mar 2010 23:59:18 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@3934
f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
src/compute_heat_flux.cpp | 234 ++++++++++++++------------------------
src/compute_heat_flux.h | 8 +-
src/fix_shake.cpp | 4 -
3 files changed, 87 insertions(+), 159 deletions(-)
diff --git a/src/compute_heat_flux.cpp b/src/compute_heat_flux.cpp
index cef6f017e0..28be7c5942 100644
--- a/src/compute_heat_flux.cpp
+++ b/src/compute_heat_flux.cpp
@@ -12,25 +12,18 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
- Contributing authors: Reese Jones (Sandia)
- Philip Howell (Siemens)
- Vikas Varsney (Air Force Research Laboratory)
+ Contributing authors: German Samolyuk (ORNL) and
+ Mario Pinto (Computational Research Lab, Pune, India)
------------------------------------------------------------------------- */
#include "math.h"
#include "string.h"
#include "compute_heat_flux.h"
#include "atom.h"
-#include "atom_vec.h"
#include "update.h"
-#include "force.h"
-#include "comm.h"
-#include "pair.h"
#include "modify.h"
+#include "force.h"
#include "group.h"
-#include "neighbor.h"
-#include "neigh_list.h"
-#include "neigh_request.h"
#include "error.h"
using namespace LAMMPS_NS;
@@ -42,23 +35,38 @@ using namespace LAMMPS_NS;
ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
- if (narg != 4) error->all("Illegal compute heat/flux command");
+ if (narg != 6) error->all("Illegal compute heat/flux command");
vector_flag = 1;
size_vector = 6;
extvector = 1;
- // store pe/atom ID used by heat flux computation
- // insure it is valid for pe/atom computation
+ // store ke/atom, pe/atom, stress/atom IDs used by heat flux computation
+ // insure they are valid for these computations
int n = strlen(arg[3]) + 1;
- id_atomPE = new char[n];
- strcpy(id_atomPE,arg[3]);
-
- int icompute = modify->find_compute(id_atomPE);
- if (icompute < 0) error->all("Could not find compute heat/flux compute ID");
- if (modify->compute[icompute]->peatomflag == 0)
+ id_ke = new char[n];
+ strcpy(id_ke,arg[3]);
+
+ n = strlen(arg[4]) + 1;
+ id_pe = new char[n];
+ strcpy(id_pe,arg[4]);
+
+ n = strlen(arg[5]) + 1;
+ id_stress = new char[n];
+ strcpy(id_stress,arg[5]);
+
+ int ike = modify->find_compute(id_ke);
+ int ipe = modify->find_compute(id_pe);
+ int istress = modify->find_compute(id_stress);
+ if (ike < 0 || ipe < 0 || istress < 0)
+ error->all("Could not find compute heat/flux compute ID");
+ if (strcmp(modify->compute[ike]->style,"ke/atom") != 0)
+ error->all("Compute heat/flux compute ID does not compute ke/atom");
+ if (modify->compute[ipe]->peatomflag == 0)
error->all("Compute heat/flux compute ID does not compute pe/atom");
+ if (modify->compute[istress]->pressatomflag == 0)
+ error->all("Compute heat/flux compute ID does not compute stress/atom");
vector = new double[6];
}
@@ -67,7 +75,9 @@ ComputeHeatFlux::ComputeHeatFlux(LAMMPS *lmp, int narg, char **arg) :
ComputeHeatFlux::~ComputeHeatFlux()
{
- delete [] id_atomPE;
+ delete [] id_ke;
+ delete [] id_pe;
+ delete [] id_stress;
delete [] vector;
}
@@ -77,155 +87,81 @@ void ComputeHeatFlux::init()
{
// error checks
- if (comm->ghost_velocity == 0)
- error->all("Compute heat/flux requires ghost atoms store velocity");
-
- if (force->pair == NULL || force->pair->single_enable == 0)
- error->all("Pair style does not support compute heat/flux");
-
- int icompute = modify->find_compute(id_atomPE);
- if (icompute < 0)
- error->all("Compute ID for compute heat/flux does not exist");
- atomPE = modify->compute[icompute];
-
- pair = force->pair;
- cutsq = force->pair->cutsq;
+ int ike = modify->find_compute(id_ke);
+ int ipe = modify->find_compute(id_pe);
+ int istress = modify->find_compute(id_stress);
+ if (ike < 0 || ipe < 0 || istress < 0)
+ error->all("Could not find compute heat/flux compute ID");
- // need an occasional half neighbor list
-
- int irequest = neighbor->request((void *) this);
- neighbor->requests[irequest]->pair = 0;
- neighbor->requests[irequest]->compute = 1;
- neighbor->requests[irequest]->occasional = 1;
-}
-
-/* ---------------------------------------------------------------------- */
-
-void ComputeHeatFlux::init_list(int id, NeighList *ptr)
-{
- list = ptr;
+ c_ke = modify->compute[ike];
+ c_pe = modify->compute[ipe];
+ c_stress = modify->compute[istress];
}
/* ---------------------------------------------------------------------- */
void ComputeHeatFlux::compute_vector()
{
- int i,j,ii,jj,inum,jnum,itype,jtype;
- double xtmp,ytmp,ztmp,delx,dely,delz;
- double rsq,eng,fpair,factor_coul,factor_lj,factor;
- double fdotv,massone,ke,pe;
- int *ilist,*jlist,*numneigh,**firstneigh;
-
invoked_vector = update->ntimestep;
- double **x = atom->x;
- double **v = atom->v;
- int *type = atom->type;
- int *mask = atom->mask;
- int nlocal = atom->nlocal;
- int nall = nlocal + atom->nghost;
- double *special_coul = force->special_coul;
- double *special_lj = force->special_lj;
- int newton_pair = force->newton_pair;
-
- // invoke half neighbor list (will copy or build if necessary)
-
- neighbor->build_one(list->index);
-
- // invoke ghost comm to insure ghost vels are up-to-date
-
- comm->forward_comm();
-
- inum = list->inum;
- ilist = list->ilist;
- numneigh = list->numneigh;
- firstneigh = list->firstneigh;
-
- // heat flux J = \sum_i e_i v_i + \sum_{i<j} (f_ij . v_j) x_ij
-
- // virial-like contribution
- // loop over neighbors of my atoms
- // require either i or j be in compute group
-
- double Jv[3] = {0.0,0.0,0.0};
-
- for (ii = 0; ii < inum; ii++) {
- i = ilist[ii];
- xtmp = x[i][0];
- ytmp = x[i][1];
- ztmp = x[i][2];
- itype = type[i];
- jlist = firstneigh[i];
- jnum = numneigh[i];
-
- for (jj = 0; jj < jnum; jj++) {
- j = jlist[jj];
-
- if (j < nall) factor_coul = factor_lj = 1.0;
- else {
- factor_coul = special_coul[j/nall];
- factor_lj = special_lj[j/nall];
- j %= nall;
- }
-
- if (!(mask[i] & groupbit) && !(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 (rsq < cutsq[itype][jtype]) {
- eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);
-
- if (newton_pair || j < nlocal) factor = 1.0;
- else factor = 0.5;
-
- // symmetrize velocities
-
- double vx = 0.5*(v[i][0]+v[j][0]);
- double vy = 0.5*(v[i][1]+v[j][1]);
- double vz = 0.5*(v[i][2]+v[j][2]);
- fdotv = factor * fpair * (delx*vx + dely*vy + delz*vz);
-
- Jv[0] += fdotv*delx;
- Jv[1] += fdotv*dely;
- Jv[2] += fdotv*delz;
- }
- }
+ // invoke 3 computes if they haven't been already
+
+ if (!(c_ke->invoked_flag & INVOKED_PERATOM)) {
+ c_ke->compute_peratom();
+ c_ke->invoked_flag |= INVOKED_PERATOM;
+ }
+ if (!(c_pe->invoked_flag & INVOKED_PERATOM)) {
+ c_pe->compute_peratom();
+ c_pe->invoked_flag |= INVOKED_PERATOM;
+ }
+ if (!(c_stress->invoked_flag & INVOKED_PERATOM)) {
+ c_stress->compute_peratom();
+ c_stress->invoked_flag |= INVOKED_PERATOM;
}
- // energy convection contribution
- // uses per-atom potential energy
+ // heat flux vector = jc[3] + jv[3]
+ // jc[3] = convective portion of heat flux = sum_i (ke_i + pe_i) v_i[3]
+ // jv[3] = virial portion of heat flux = sum_i (stress_tensor_i . v_i[3])
+ // normalization by volume is not included
- if (!(atomPE->invoked_flag & INVOKED_PERATOM)) {
- atomPE->compute_peratom();
- atomPE->invoked_flag |= INVOKED_PERATOM;
- }
+ double *ke = c_ke->vector_atom;
+ double *pe = c_pe->vector_atom;
+ double **stress = c_stress->array_atom;
- double *mass = atom->mass;
- double *rmass = atom->rmass;
- double mvv2e = force->mvv2e;
+ double **v = atom->v;
+ int *mask = atom->mask;
+ int nlocal = atom->nlocal;
- double Jc[3] = {0.0,0.0,0.0};
+ double jc[3] = {0.0,0.0,0.0};
+ double jv[3] = {0.0,0.0,0.0};
+ double eng;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
- massone = (rmass) ? rmass[i] : mass[type[i]];
- ke = mvv2e * 0.5 * massone *
- (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
- pe = atomPE->vector_atom[i];
- eng = pe + ke;
- Jc[0] += v[i][0]*eng;
- Jc[1] += v[i][1]*eng;
- Jc[2] += v[i][2]*eng;
+ eng = pe[i] + ke[i];
+ jc[0] += eng*v[i][0];
+ jc[1] += eng*v[i][1];
+ jc[2] += eng*v[i][2];
+ jv[0] -= stress[i][0]*v[i][0] + stress[i][3]*v[i][1] +
+ stress[i][4]*v[i][2];
+ jv[1] -= stress[i][3]*v[i][0] + stress[i][1]*v[i][1] +
+ stress[i][5]*v[i][2];
+ jv[2] -= stress[i][4]*v[i][0] + stress[i][5]*v[i][1] +
+ stress[i][2]*v[i][2];
}
}
- // total flux
+ // convert jv from stress*volume to energy units via nktv2p factor
+
+ double nktv2p = force->nktv2p;
+ jv[0] /= nktv2p;
+ jv[1] /= nktv2p;
+ jv[2] /= nktv2p;
- double data[6] = {Jv[0]+Jc[0],Jv[1]+Jc[1],Jv[2]+Jc[2],
- Jc[0],Jc[1],Jc[2]};
+ // sum across all procs
+ // 1st 3 terms are total heat flux
+ // 2nd 3 terms are just conductive portion
+
+ double data[6] = {jc[0]+jv[0],jc[1]+jv[1],jc[2]+jv[2],jc[0],jc[1],jc[2]};
MPI_Allreduce(data,vector,6,MPI_DOUBLE,MPI_SUM,world);
}
diff --git a/src/compute_heat_flux.h b/src/compute_heat_flux.h
index f332c79afa..2994b1e9a9 100644
--- a/src/compute_heat_flux.h
+++ b/src/compute_heat_flux.h
@@ -29,15 +29,11 @@ class ComputeHeatFlux : public Compute {
ComputeHeatFlux(class LAMMPS *, int, char **);
~ComputeHeatFlux();
void init();
- void init_list(int, class NeighList *);
void compute_vector();
private:
- double **cutsq;
- class Pair *pair;
- class NeighList *list;
- class Compute *atomPE;
- char *id_atomPE;
+ char *id_ke,*id_pe,*id_stress;
+ class Compute *c_ke,*c_pe,*c_stress;
};
}
diff --git a/src/fix_shake.cpp b/src/fix_shake.cpp
index 20ef80fe35..8f1ff82833 100644
--- a/src/fix_shake.cpp
+++ b/src/fix_shake.cpp
@@ -2031,10 +2031,6 @@ void FixShake::shake3angle(int m)
v_tally(nlist,list,3.0,v);
}
-
- if (i0 < 20)
- printf("AAA %d %d %d %g %g %g: %g\n",
- i0,i1,i2,lamda01,lamda02,lamda12,v[0]);
}
/* ----------------------------------------------------------------------
--
GitLab