1 /* ----------------------------------------------------------------------
2     This is the
3 
4     ██╗     ██╗ ██████╗  ██████╗  ██████╗ ██╗  ██╗████████╗███████╗
5     ██║     ██║██╔════╝ ██╔════╝ ██╔════╝ ██║  ██║╚══██╔══╝██╔════╝
6     ██║     ██║██║  ███╗██║  ███╗██║  ███╗███████║   ██║   ███████╗
7     ██║     ██║██║   ██║██║   ██║██║   ██║██╔══██║   ██║   ╚════██║
8     ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║  ██║   ██║   ███████║
9     ╚══════╝╚═╝ ╚═════╝  ╚═════╝  ╚═════╝ ╚═╝  ╚═╝   ╚═╝   ╚══════╝®
10 
11     DEM simulation engine, released by
12     DCS Computing Gmbh, Linz, Austria
13     http://www.dcs-computing.com, office@dcs-computing.com
14 
15     LIGGGHTS® is part of CFDEM®project:
16     http://www.liggghts.com | http://www.cfdem.com
17 
18     Core developer and main author:
19     Christoph Kloss, christoph.kloss@dcs-computing.com
20 
21     LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22     License, version 2 or later. It is distributed in the hope that it will
23     be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24     of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25     received a copy of the GNU General Public License along with LIGGGHTS®.
26     If not, see http://www.gnu.org/licenses . See also top-level README
27     and LICENSE files.
28 
29     LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30     the producer of the LIGGGHTS® software and the CFDEM®coupling software
31     See http://www.cfdem.com/terms-trademark-policy for details.
32 
33 -------------------------------------------------------------------------
34     Contributing author and copyright for this file:
35     This file is from LAMMPS
36     LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37     http://lammps.sandia.gov, Sandia National Laboratories
38     Steve Plimpton, sjplimp@sandia.gov
39 
40     Copyright (2003) Sandia Corporation.  Under the terms of Contract
41     DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42     certain rights in this software.  This software is distributed under
43     the GNU General Public License.
44 ------------------------------------------------------------------------- */
45 
46 #include <string.h>
47 #include "compute_ke_atom.h"
48 #include "atom.h"
49 #include "update.h"
50 #include "modify.h"
51 #include "comm.h"
52 #include "fix_multisphere.h"
53 #include "force.h"
54 #include "memory.h"
55 #include "error.h"
56 
57 using namespace LAMMPS_NS;
58 
59 /* ---------------------------------------------------------------------- */
60 
ComputeKEAtom(LAMMPS * lmp,int & iarg,int narg,char ** arg)61 ComputeKEAtom::ComputeKEAtom(LAMMPS *lmp, int &iarg, int narg, char **arg) :
62   Compute(lmp, iarg, narg, arg)
63 {
64   if (narg != iarg) error->all(FLERR,"Illegal compute ke/atom command");
65 
66   peratom_flag = 1;
67   size_peratom_cols = 0;
68 
69   nmax = 0;
70   ke = NULL;
71   fix_ms = NULL;
72 }
73 
74 /* ---------------------------------------------------------------------- */
75 
~ComputeKEAtom()76 ComputeKEAtom::~ComputeKEAtom()
77 {
78   memory->destroy(ke);
79 }
80 
81 /* ---------------------------------------------------------------------- */
82 
init()83 void ComputeKEAtom::init()
84 {
85   int count = 0;
86   for (int i = 0; i < modify->ncompute; i++)
87     if (strcmp(modify->compute[i]->style,"ke/atom") == 0) count++;
88   if (count > 1 && comm->me == 0)
89     error->warning(FLERR,"More than one compute ke/atom");
90   fix_ms =  static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
91 }
92 
93 /* ---------------------------------------------------------------------- */
94 
compute_peratom()95 void ComputeKEAtom::compute_peratom()
96 {
97   invoked_peratom = update->ntimestep;
98 
99   // grow ke array if necessary
100 
101   if (atom->nlocal > nmax) {
102     memory->destroy(ke);
103     nmax = atom->nmax;
104     memory->create(ke,nmax,"ke/atom:ke");
105     vector_atom = ke;
106   }
107 
108   // compute kinetic energy for each atom in group
109 
110   double mvv2e = force->mvv2e;
111   double **v = atom->v;
112   double *mass = atom->mass;
113   double *rmass = atom->rmass;
114   int *mask = atom->mask;
115   int *type = atom->type;
116   int nlocal = atom->nlocal;
117 
118   if (rmass)
119     for (int i = 0; i < nlocal; i++) {
120       if (mask[i] & groupbit && (!fix_ms || fix_ms->belongs_to(i) < 0)) {
121         ke[i] = 0.5 * mvv2e * rmass[i] *
122           (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
123       } else ke[i] = 0.0;
124     }
125 
126   else
127     for (int i = 0; i < nlocal; i++) {
128       if (mask[i] & groupbit) {
129         ke[i] = 0.5 * mvv2e * mass[type[i]] *
130           (v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2]);
131       } else ke[i] = 0.0;
132     }
133 }
134 
135 /* ----------------------------------------------------------------------
136    memory usage of local atom-based array
137 ------------------------------------------------------------------------- */
138 
memory_usage()139 double ComputeKEAtom::memory_usage()
140 {
141   double bytes = nmax * sizeof(double);
142   return bytes;
143 }
144