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 <cmath>
47 #include <stdlib.h>
48 #include <string.h>
49 #include "fix_viscous.h"
50 #include "atom.h"
51 #include "update.h"
52 #include "respa.h"
53 #include "error.h"
54 #include "force.h"
55
56 using namespace LAMMPS_NS;
57 using namespace FixConst;
58
59 /* ---------------------------------------------------------------------- */
60
FixViscous(LAMMPS * lmp,int narg,char ** arg)61 FixViscous::FixViscous(LAMMPS *lmp, int narg, char **arg) :
62 Fix(lmp, narg, arg)
63 {
64 if (narg < 4) error->all(FLERR,"Illegal fix viscous command");
65
66 double gamma_one = force->numeric(FLERR,arg[3]);
67 gamma = new double[atom->ntypes+1];
68 for (int i = 1; i <= atom->ntypes; i++) gamma[i] = gamma_one;
69
70 // optional args
71
72 int iarg = 4;
73 while (iarg < narg) {
74 if (strcmp(arg[iarg],"scale") == 0) {
75 if (iarg+3 > narg) error->all(FLERR,"Illegal fix viscous command");
76 int itype = force->inumeric(FLERR,arg[iarg+1]);
77 double scale = force->numeric(FLERR,arg[iarg+2]);
78 if (itype <= 0 || itype > atom->ntypes)
79 error->all(FLERR,"Illegal fix viscous command");
80 gamma[itype] = gamma_one * scale;
81 iarg += 3;
82 } else error->all(FLERR,"Illegal fix viscous command");
83 }
84 }
85
86 /* ---------------------------------------------------------------------- */
87
~FixViscous()88 FixViscous::~FixViscous()
89 {
90 delete [] gamma;
91 }
92
93 /* ---------------------------------------------------------------------- */
94
setmask()95 int FixViscous::setmask()
96 {
97 int mask = 0;
98 mask |= POST_FORCE;
99 mask |= POST_FORCE_RESPA;
100 mask |= MIN_POST_FORCE;
101 return mask;
102 }
103
104 /* ---------------------------------------------------------------------- */
105
init()106 void FixViscous::init()
107 {
108 if (strstr(update->integrate_style,"respa"))
109 nlevels_respa = ((Respa *) update->integrate)->nlevels;
110 }
111
112 /* ---------------------------------------------------------------------- */
113
setup(int vflag)114 void FixViscous::setup(int vflag)
115 {
116 if (strstr(update->integrate_style,"verlet"))
117 post_force(vflag);
118 else {
119 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
120 post_force_respa(vflag,nlevels_respa-1,0);
121 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
122 }
123
124 // error checks on coarsegraining
125 if(force->cg_active())
126 error->cg(FLERR,this->style);
127 }
128
129 /* ---------------------------------------------------------------------- */
130
min_setup(int vflag)131 void FixViscous::min_setup(int vflag)
132 {
133 post_force(vflag);
134 }
135
136 /* ---------------------------------------------------------------------- */
137
post_force(int vflag)138 void FixViscous::post_force(int vflag)
139 {
140 // apply drag force to atoms in group
141 // direction is opposed to velocity vector
142 // magnitude depends on atom type
143
144 double **v = atom->v;
145 double **f = atom->f;
146 int *mask = atom->mask;
147 int *type = atom->type;
148 int nlocal = atom->nlocal;
149
150 double drag;
151
152 for (int i = 0; i < nlocal; i++)
153 if (mask[i] & groupbit) {
154 drag = gamma[type[i]];
155 f[i][0] -= drag*v[i][0];
156 f[i][1] -= drag*v[i][1];
157 f[i][2] -= drag*v[i][2];
158 }
159 }
160
161 /* ---------------------------------------------------------------------- */
162
post_force_respa(int vflag,int ilevel,int iloop)163 void FixViscous::post_force_respa(int vflag, int ilevel, int iloop)
164 {
165 if (ilevel == nlevels_respa-1) post_force(vflag);
166 }
167
168 /* ---------------------------------------------------------------------- */
169
min_post_force(int vflag)170 void FixViscous::min_post_force(int vflag)
171 {
172 post_force(vflag);
173 }
174