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, but has been modified. Copyright for
36 modification:
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40
41 Copyright of original file:
42 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
43 http://lammps.sandia.gov, Sandia National Laboratories
44 Steve Plimpton, sjplimp@sandia.gov
45
46 Copyright (2003) Sandia Corporation. Under the terms of Contract
47 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
48 certain rights in this software. This software is distributed under
49 the GNU General Public License.
50 ------------------------------------------------------------------------- */
51
52 #include <mpi.h>
53 #include "compute_erotate_sphere.h"
54 #include "atom.h"
55 #include "atom_vec.h"
56 #include "update.h"
57 #include "force.h"
58 #include "domain.h"
59 #include "modify.h"
60 #include "fix_multisphere.h"
61 #include "group.h"
62 #include "error.h"
63
64 using namespace LAMMPS_NS;
65
66 #define INERTIA 0.4 // moment of inertia prefactor for sphere
67
68 /* ---------------------------------------------------------------------- */
69
ComputeERotateSphere(LAMMPS * lmp,int & iarg,int narg,char ** arg)70 ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int &iarg, int narg, char **arg) :
71 Compute(lmp, iarg, narg, arg),
72 pfactor(1.0),
73 halfstep(false),
74 fix_ms(NULL)
75 {
76 if (narg < iarg || narg > iarg+1) error->all(FLERR,"Illegal compute erotate/sphere command");
77
78 if (narg == iarg+1)
79 {
80 if (strcmp(arg[iarg++], "halfstep") == 0)
81 halfstep = true;
82 else
83 error->all(FLERR,"Illegal compute erotate/sphere option");
84 }
85
86 scalar_flag = 1;
87 extscalar = 1;
88
89 // error check
90
91 if (!atom->sphere_flag)
92 error->all(FLERR,"Compute erotate/sphere requires atom style sphere");
93 }
94
95 /* ---------------------------------------------------------------------- */
96
init()97 void ComputeERotateSphere::init()
98 {
99 pfactor = 0.5 * force->mvv2e * INERTIA;
100
101 fix_ms = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
102 }
103
104 /* ---------------------------------------------------------------------- */
105
compute_scalar()106 double ComputeERotateSphere::compute_scalar()
107 {
108 if (invoked_scalar == update->ntimestep)
109 return scalar;
110 invoked_scalar = update->ntimestep;
111
112 double **omega = atom->omega;
113 double **torque = atom->torque;
114 double *radius = atom->radius;
115 double *rmass = atom->rmass;
116 int *mask = atom->mask;
117 int nlocal = atom->nlocal;
118
119 // sum rotational energy for each particle
120 // point particles will not contribute, due to radius = 0.0
121
122 double erotate = 0.0;
123 for (int i = 0; i < nlocal; i++)
124 {
125 if (mask[i] & groupbit && (!fix_ms || fix_ms->belongs_to(i) < 0))
126 {
127 const double mult = halfstep ? update->dt*0.5/(INERTIA*radius[i]*radius[i]*rmass[i]) : 0.0;
128 const double omega0 = omega[i][0] + mult*torque[i][0];
129 const double omega1 = omega[i][1] + mult*torque[i][1];
130 const double omega2 = omega[i][2] + mult*torque[i][2];
131 erotate += (omega0*omega0 + omega1*omega1 + omega2*omega2) * radius[i]*radius[i]*rmass[i];
132 }
133 }
134
135 MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);
136 scalar *= pfactor;
137
138 return scalar;
139 }
140