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 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2014- DCS Computing GmbH, Linz
39 ------------------------------------------------------------------------- */
40
41 #include <cmath>
42 #include <stdlib.h>
43 #include <string.h>
44 #include "fix_buoyancy.h"
45 #include "atom.h"
46 #include "update.h"
47 #include "force.h"
48 #include "respa.h"
49 #include "domain.h"
50 #include "modify.h"
51 #include "fix_gravity.h"
52 #include "region.h"
53 #include "vector_liggghts.h"
54 #include "math_extra_liggghts.h"
55 #include "mpi_liggghts.h"
56 #include "error.h"
57
58 using namespace LAMMPS_NS;
59 using namespace FixConst;
60
61 /* ---------------------------------------------------------------------- */
62
FixBuoyancy(LAMMPS * lmp,int narg,char ** arg)63 FixBuoyancy::FixBuoyancy(LAMMPS *lmp, int narg, char **arg) :
64 FixBaseLiggghts(lmp, narg, arg),
65 density_(0.),
66 dim_(-1),
67 direction_(1.),
68 fluid_level_(0.),
69 fix_gravity_(0),
70 force_flag_(false)
71 {
72 vectorZeroize3D(buyoancy_force_total_);
73
74 do_support_respa();
75
76 if (narg < 3) error->fix_error(FLERR,this,"not enough arguments");
77
78 vector_flag = 1;
79 size_vector = 3;
80 global_freq = 1;
81 extvector = 1;
82
83 // args
84
85 int iarg = 3;
86 while (iarg < narg) {
87 if (strcmp(arg[iarg],"density") == 0) {
88 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments for 'density'");
89 density_ = atof(arg[iarg+1]);
90 if (density_ <= 0.)
91 error->fix_error(FLERR,this,"'density' > 0 required");
92 iarg += 2;
93 } else if (strcmp(arg[iarg],"level") == 0) {
94 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments for 'level'");
95 fluid_level_ = atof(arg[iarg+1]);
96 if (fluid_level_ <= 0.)
97 error->fix_error(FLERR,this,"'level' > 0 required");
98 iarg += 2;
99 } else if (strcmp(arg[iarg],"dim") == 0) {
100 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments for 'dim'");
101 if(0 == strcmp(arg[iarg+1],"x"))
102 dim_ = 0;
103 else if(0 == strcmp(arg[iarg+1],"y"))
104 dim_ = 1;
105 else if(0 == strcmp(arg[iarg+1],"z"))
106 dim_ = 2;
107 else
108 error->fix_error(FLERR,this,"expecting 'x' or 'y' or 'z' after 'dim'");
109 iarg += 2;
110 } else if (strcmp(arg[iarg],"region") == 0) {
111 if (iarg+2 > narg) error->fix_error(FLERR,this,"not enough arguments for 'region'");
112 process_region(arg[iarg+1]);
113 iarg += 2;
114 } else error->fix_error(FLERR,this," expecting 'density' or 'region'");
115 }
116
117 if(-1 == dim_)
118 error->fix_error(FLERR,this," you have to specify 'dim'");
119
120 if(MathExtraLiggghts::compDouble(density_,0.,1e-6))
121 error->fix_error(FLERR,this," you have to specify 'density'");
122
123 }
124
125 /* ---------------------------------------------------------------------- */
126
~FixBuoyancy()127 FixBuoyancy::~FixBuoyancy()
128 {
129 }
130
131 /* ---------------------------------------------------------------------- */
132
setmask()133 int FixBuoyancy::setmask()
134 {
135 int mask = 0;
136 mask |= POST_FORCE;
137 mask |= POST_FORCE_RESPA;
138 mask |= MIN_POST_FORCE;
139 return mask;
140 }
141
142 /* ---------------------------------------------------------------------- */
143
init()144 void FixBuoyancy::init()
145 {
146 FixBaseLiggghts::init();
147
148 if(!atom->density_flag)
149 error->fix_error(FLERR,this,"requires per-atom density");
150
151 if(1 != modify->n_fixes_style_strict("gravity"))
152 error->fix_error(FLERR,this,"need exactly one fix gravity");
153
154 fix_gravity_ = static_cast<FixGravity*>(modify->find_fix_style_strict("gravity",0));
155
156 // test if dim of fix gravity and this fix coincide
157 test_direction();
158 }
159
160 /* ---------------------------------------------------------------------- */
161
test_direction()162 void FixBuoyancy::test_direction()
163 {
164 double gravity[3];
165 fix_gravity_->get_gravity(gravity);
166
167 // test if dim of fix gravity and this fix coincide
168 double test[3] = {1.,1.,1.};
169 test[dim_] = 0.;
170 if(MathExtraLiggghts::abs(vectorDot3D(gravity,test)) > 0.)
171 error->fix_error(FLERR,this," fix gravity direction is not equal to 'dim'");
172
173 vectorZeroize3D(test);
174 test[dim_] = 1.;
175
176 direction_ = vectorDot3D(gravity,test);
177 direction_ = direction_ / MathExtraLiggghts::abs(direction_);
178
179 }
180
181 /* ---------------------------------------------------------------------- */
182
setup(int vflag)183 void FixBuoyancy::setup(int vflag)
184 {
185 FixBaseLiggghts::setup(vflag);
186
187 // error checks on coarsegraining
188 if(force->cg_active())
189 error->cg(FLERR,this->style);
190 }
191
192 /* ---------------------------------------------------------------------- */
193
min_setup(int vflag)194 void FixBuoyancy::min_setup(int vflag)
195 {
196 post_force(vflag);
197 }
198
199 /* ---------------------------------------------------------------------- */
200
post_force(int vflag)201 void FixBuoyancy::post_force(int vflag)
202 {
203 // apply drag force to atoms in group
204 // direction is opposed to velocity vector
205 // magnitude depends on atom type
206
207 double **x = atom->x;
208 double **f = atom->f;
209 double *radius = atom->radius;
210 int *mask = atom->mask;
211 int nlocal = atom->nlocal;
212
213 double buyoancy_force;
214 vectorZeroize3D(buyoancy_force_total_);
215 force_flag_ = false;
216
217 double gravity[3],gravity_mag;
218 fix_gravity_->get_gravity(gravity);
219 gravity_mag = vectorMag3D(gravity);
220
221 // test if dim of fix gravity and this fix coincide
222 // need to re-do this since orientation of gravity might change
223 test_direction();
224
225 double c = 4.*M_PI/3.;
226 double pi6 = M_PI/6.;
227
228 for (int i = 0; i < nlocal; i++)
229 if (mask[i] & groupbit) {
230
231 if (region_ && !region_->match(x[i][0],x[i][1],x[i][2]))
232 continue;
233
234 double h = (x[i][dim_] - fluid_level_)*direction_;
235 double hmag = MathExtraLiggghts::abs(h);
236 double Vimmersed = 0.;
237
238 // particle not immersed
239 if(hmag > radius[i] && h < 0)
240 {
241
242 continue;
243 }
244 // particle fully immersed
245 else if(hmag > radius[i] && h > 0)
246 {
247
248 Vimmersed = c*radius[i]*radius[i]*radius[i];
249 }
250 // larger part of sphere immersed
251 else if(h > 0)
252 {
253
254 h = radius[i] - h;
255 double a = sqrt(2.*radius[i]*h - h*h);
256 double Vsegment = h*pi6 * (3.*a*a + h*h);
257 Vimmersed = (c*radius[i]*radius[i]*radius[i] - Vsegment);
258 }
259 // smaller part of sphere immersed
260 else
261 {
262
263 h = -h;
264 h = radius[i] - h;
265 double a = sqrt(2.*radius[i]*h - h*h);
266 Vimmersed = h * pi6 * (3.*a*a + h*h);
267 }
268
269 buyoancy_force = - direction_*gravity_mag*Vimmersed*(density_);
270 f[i][dim_] += buyoancy_force;
271 buyoancy_force_total_[dim_] += buyoancy_force;
272
273 }
274 }
275
276 /* ---------------------------------------------------------------------- */
277
post_force_respa(int vflag,int ilevel,int iloop)278 void FixBuoyancy::post_force_respa(int vflag, int ilevel, int iloop)
279 {
280 if (ilevel == nlevels_respa-1) post_force(vflag);
281 }
282
283 /* ---------------------------------------------------------------------- */
284
min_post_force(int vflag)285 void FixBuoyancy::min_post_force(int vflag)
286 {
287 post_force(vflag);
288 }
289
290 /* ----------------------------------------------------------------------
291 return components of total force on fix group before force was changed
292 ------------------------------------------------------------------------- */
293
compute_vector(int n)294 double FixBuoyancy::compute_vector(int n)
295 {
296 // only sum across procs one time
297
298 if (!force_flag_) {
299 MPI_Sum_Vector(buyoancy_force_total_,3,world);
300 force_flag_ = true;
301 }
302 return buyoancy_force_total_[n];
303 }
304