1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14
15 /* ----------------------------------------------------------------------
16 Contributing authors: Frances Mackay, Santtu Ollila, Colin Denniston (UWO)
17 ------------------------------------------------------------------------- */
18
19 #include "fix_lb_viscous.h"
20 #include <cstring>
21 #include "atom.h"
22 #include "update.h"
23 #include "respa.h"
24 #include "error.h"
25 #include "fix_lb_fluid.h"
26 #include "modify.h"
27 #include "group.h"
28
29 using namespace LAMMPS_NS;
30 using namespace FixConst;
31
32 /* ---------------------------------------------------------------------- */
33
FixLbViscous(LAMMPS * lmp,int narg,char ** arg)34 FixLbViscous::FixLbViscous(LAMMPS *lmp, int narg, char **arg) :
35 Fix(lmp, narg, arg)
36 {
37 if (narg < 3) error->all(FLERR,"Illegal fix lb/viscous command");
38
39 int groupbit_lb_fluid = 0;
40
41 for (int ifix=0; ifix<modify->nfix; ifix++)
42 if (strcmp(modify->fix[ifix]->style,"lb/fluid")==0) {
43 fix_lb_fluid = (FixLbFluid *)modify->fix[ifix];
44 groupbit_lb_fluid = group->bitmask[modify->fix[ifix]->igroup];
45 }
46
47 if (groupbit_lb_fluid == 0)
48 error->all(FLERR,"the lb/fluid fix must also be used if using the lb/viscous fix");
49
50 int *mask = atom->mask;
51 int nlocal = atom->nlocal;
52 for (int j=0; j<nlocal; j++) {
53 if ((mask[j] & groupbit) && !(mask[j] & groupbit_lb_fluid))
54 error->one(FLERR,"to apply a fluid force onto an atom, the lb/fluid fix must be used for that atom.");
55 }
56
57
58 }
59
60 /* ---------------------------------------------------------------------- */
61
~FixLbViscous()62 FixLbViscous::~FixLbViscous()
63 {
64
65 }
66
67 /* ---------------------------------------------------------------------- */
68
setmask()69 int FixLbViscous::setmask()
70 {
71 int mask = 0;
72 mask |= POST_FORCE;
73 mask |= POST_FORCE_RESPA;
74 mask |= MIN_POST_FORCE;
75 return mask;
76 }
77
78 /* ---------------------------------------------------------------------- */
79
init()80 void FixLbViscous::init()
81 {
82
83 if (utils::strmatch(update->integrate_style,"^respa"))
84 nlevels_respa = ((Respa *) update->integrate)->nlevels;
85
86 }
87
88 /* ---------------------------------------------------------------------- */
89
setup(int vflag)90 void FixLbViscous::setup(int vflag)
91 {
92 if (utils::strmatch(update->integrate_style,"^verlet"))
93 post_force(vflag);
94 else {
95 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
96 post_force_respa(vflag,nlevels_respa-1,0);
97 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
98 }
99
100 }
101
102 /* ---------------------------------------------------------------------- */
103
min_setup(int vflag)104 void FixLbViscous::min_setup(int vflag)
105 {
106 post_force(vflag);
107 }
108
109 /* ---------------------------------------------------------------------- */
110
post_force(int)111 void FixLbViscous::post_force(int /*vflag*/)
112 {
113 // apply drag force to atoms in group
114 // direction is opposed to velocity vector
115 // magnitude depends on atom type
116
117 double **f = atom->f;
118 int *mask = atom->mask;
119 int nlocal = atom->nlocal;
120 double **hydroF = fix_lb_fluid->hydroF;
121
122 for (int i = 0; i < nlocal; i++)
123 if (mask[i] & groupbit) {
124 f[i][0] += hydroF[i][0];
125 f[i][1] += hydroF[i][1];
126 f[i][2] += hydroF[i][2];
127
128 }
129
130 }
131
132 /* ---------------------------------------------------------------------- */
133
post_force_respa(int vflag,int ilevel,int)134 void FixLbViscous::post_force_respa(int vflag, int ilevel, int /*iloop*/)
135 {
136 if (ilevel == nlevels_respa-1) post_force(vflag);
137 }
138
139 /* ---------------------------------------------------------------------- */
140
min_post_force(int vflag)141 void FixLbViscous::min_post_force(int vflag)
142 {
143 post_force(vflag);
144 }
145
146