1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/, Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 /* ----------------------------------------------------------------------
15 Contributing authors: Jeremy Lechman (SNL), Pieter in 't Veld (BASF)
16 ------------------------------------------------------------------------- */
17
18 #include "fix_srd.h"
19
20 #include "atom.h"
21 #include "atom_vec_ellipsoid.h"
22 #include "atom_vec_line.h"
23 #include "atom_vec_tri.h"
24 #include "citeme.h"
25 #include "comm.h"
26 #include "domain.h"
27 #include "error.h"
28 #include "fix_deform.h"
29 #include "fix_wall_srd.h"
30 #include "force.h"
31 #include "group.h"
32 #include "math_const.h"
33 #include "math_extra.h"
34 #include "memory.h"
35 #include "modify.h"
36 #include "neighbor.h"
37 #include "random_mars.h"
38 #include "random_park.h"
39 #include "update.h"
40
41 #include <cmath>
42 #include <cstring>
43
44 using namespace LAMMPS_NS;
45 using namespace FixConst;
46 using namespace MathConst;
47
48 enum { SLIP, NOSLIP };
49 enum { SPHERE, ELLIPSOID, LINE, TRIANGLE, WALL };
50 enum { INSIDE_ERROR, INSIDE_WARN, INSIDE_IGNORE };
51 enum { BIG_MOVE, SRD_MOVE, SRD_ROTATE };
52 enum { CUBIC_ERROR, CUBIC_WARN };
53 enum { SHIFT_NO, SHIFT_YES, SHIFT_POSSIBLE };
54
55 #define EINERTIA 0.2 // moment of inertia prefactor for ellipsoid
56
57 #define ATOMPERBIN 30
58 #define BIG 1.0e20
59 #define VBINSIZE 5
60 #define TOLERANCE 0.00001
61 #define MAXITER 20
62
63 static const char cite_fix_srd[] = "fix srd command:\n\n"
64 "@Article{Petersen10,\n"
65 " author = {M. K. Petersen, J. B. Lechman, S. J. Plimpton, G. "
66 "S. Grest, P. J. in 't Veld, P. R. Schunk},\n"
67 " title = {Mesoscale Hydrodynamics via Stochastic Rotation "
68 "Dynamics: Comparison with Lennard-Jones Fluid},"
69 " journal = {J.~Chem.~Phys.},\n"
70 " year = 2010,\n"
71 " volume = 132,\n"
72 " pages = {174106}\n"
73 "}\n\n";
74
75 //#define SRD_DEBUG 1
76 //#define SRD_DEBUG_ATOMID 58
77 //#define SRD_DEBUG_TIMESTEP 449
78
79 /* ---------------------------------------------------------------------- */
80
FixSRD(LAMMPS * lmp,int narg,char ** arg)81 FixSRD::FixSRD(LAMMPS *lmp, int narg, char **arg) :
82 Fix(lmp, narg, arg), wallfix(nullptr), wallwhich(nullptr), xwall(nullptr), xwallhold(nullptr),
83 vwall(nullptr), fwall(nullptr), avec_ellipsoid(nullptr), avec_line(nullptr), avec_tri(nullptr),
84 random(nullptr), randomshift(nullptr), flocal(nullptr), tlocal(nullptr), biglist(nullptr),
85 binhead(nullptr), binnext(nullptr), sbuf1(nullptr), sbuf2(nullptr), rbuf1(nullptr),
86 rbuf2(nullptr), nbinbig(nullptr), binbig(nullptr), binsrd(nullptr), stencil(nullptr)
87 {
88 if (lmp->citeme) lmp->citeme->add(cite_fix_srd);
89
90 if (narg < 8) error->all(FLERR, "Illegal fix srd command");
91
92 restart_pbc = 1;
93 vector_flag = 1;
94 size_vector = 12;
95 global_freq = 1;
96 extvector = 0;
97
98 nevery = utils::inumeric(FLERR, arg[3], false, lmp);
99
100 bigexist = 1;
101 if (strcmp(arg[4], "NULL") == 0)
102 bigexist = 0;
103 else
104 biggroup = group->find(arg[4]);
105
106 temperature_srd = utils::numeric(FLERR, arg[5], false, lmp);
107 gridsrd = utils::numeric(FLERR, arg[6], false, lmp);
108 int seed = utils::inumeric(FLERR, arg[7], false, lmp);
109
110 // parse options
111
112 lamdaflag = 0;
113 collidestyle = NOSLIP;
114 overlap = 0;
115 insideflag = INSIDE_ERROR;
116 exactflag = 1;
117 radfactor = 1.0;
118 maxbounceallow = 0;
119 gridsearch = gridsrd;
120 cubicflag = CUBIC_ERROR;
121 cubictol = 0.01;
122 shiftuser = SHIFT_NO;
123 shiftseed = 0;
124 tstat = 0;
125 rescale_rotate = rescale_collide = 1;
126
127 int iarg = 8;
128 while (iarg < narg) {
129 if (strcmp(arg[iarg], "lamda") == 0) {
130 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
131 lamda = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
132 lamdaflag = 1;
133 iarg += 2;
134 } else if (strcmp(arg[iarg], "collision") == 0) {
135 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
136 if (strcmp(arg[iarg + 1], "slip") == 0)
137 collidestyle = SLIP;
138 else if (strcmp(arg[iarg + 1], "noslip") == 0)
139 collidestyle = NOSLIP;
140 else
141 error->all(FLERR, "Illegal fix srd command");
142 iarg += 2;
143 } else if (strcmp(arg[iarg], "overlap") == 0) {
144 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
145 if (strcmp(arg[iarg + 1], "yes") == 0)
146 overlap = 1;
147 else if (strcmp(arg[iarg + 1], "no") == 0)
148 overlap = 0;
149 else
150 error->all(FLERR, "Illegal fix srd command");
151 iarg += 2;
152 } else if (strcmp(arg[iarg], "inside") == 0) {
153 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
154 if (strcmp(arg[iarg + 1], "error") == 0)
155 insideflag = INSIDE_ERROR;
156 else if (strcmp(arg[iarg + 1], "warn") == 0)
157 insideflag = INSIDE_WARN;
158 else if (strcmp(arg[iarg + 1], "ignore") == 0)
159 insideflag = INSIDE_IGNORE;
160 else
161 error->all(FLERR, "Illegal fix srd command");
162 iarg += 2;
163 } else if (strcmp(arg[iarg], "exact") == 0) {
164 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
165 if (strcmp(arg[iarg + 1], "yes") == 0)
166 exactflag = 1;
167 else if (strcmp(arg[iarg + 1], "no") == 0)
168 exactflag = 0;
169 else
170 error->all(FLERR, "Illegal fix srd command");
171 iarg += 2;
172 } else if (strcmp(arg[iarg], "radius") == 0) {
173 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
174 radfactor = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
175 iarg += 2;
176 } else if (strcmp(arg[iarg], "bounce") == 0) {
177 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
178 maxbounceallow = utils::inumeric(FLERR, arg[iarg + 1], false, lmp);
179 iarg += 2;
180 } else if (strcmp(arg[iarg], "search") == 0) {
181 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
182 gridsearch = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
183 iarg += 2;
184 } else if (strcmp(arg[iarg], "cubic") == 0) {
185 if (iarg + 3 > narg) error->all(FLERR, "Illegal fix srd command");
186 if (strcmp(arg[iarg + 1], "error") == 0)
187 cubicflag = CUBIC_ERROR;
188 else if (strcmp(arg[iarg + 1], "warn") == 0)
189 cubicflag = CUBIC_WARN;
190 else
191 error->all(FLERR, "Illegal fix srd command");
192 cubictol = utils::numeric(FLERR, arg[iarg + 2], false, lmp);
193 iarg += 3;
194 } else if (strcmp(arg[iarg], "shift") == 0) {
195 if (iarg + 3 > narg)
196 error->all(FLERR, "Illegal fix srd command");
197 else if (strcmp(arg[iarg + 1], "no") == 0)
198 shiftuser = SHIFT_NO;
199 else if (strcmp(arg[iarg + 1], "yes") == 0)
200 shiftuser = SHIFT_YES;
201 else if (strcmp(arg[iarg + 1], "possible") == 0)
202 shiftuser = SHIFT_POSSIBLE;
203 else
204 error->all(FLERR, "Illegal fix srd command");
205 shiftseed = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
206 iarg += 3;
207 } else if (strcmp(arg[iarg], "tstat") == 0) {
208 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
209 if (strcmp(arg[iarg + 1], "no") == 0)
210 tstat = 0;
211 else if (strcmp(arg[iarg + 1], "yes") == 0)
212 tstat = 1;
213 else
214 error->all(FLERR, "Illegal fix srd command");
215 iarg += 2;
216 } else if (strcmp(arg[iarg], "rescale") == 0) {
217 if (iarg + 2 > narg) error->all(FLERR, "Illegal fix srd command");
218 if (strcmp(arg[iarg + 1], "no") == 0)
219 rescale_rotate = rescale_collide = 0;
220 else if (strcmp(arg[iarg + 1], "yes") == 0)
221 rescale_rotate = rescale_collide = 1;
222 else if (strcmp(arg[iarg + 1], "rotate") == 0) {
223 rescale_rotate = 1;
224 rescale_collide = 0;
225 } else if (strcmp(arg[iarg + 1], "collide") == 0) {
226 rescale_rotate = 0;
227 rescale_collide = 1;
228 } else
229 error->all(FLERR, "Illegal fix srd command");
230 iarg += 2;
231 } else
232 error->all(FLERR, "Illegal fix srd command");
233 }
234
235 // error check
236
237 if (nevery <= 0) error->all(FLERR, "Illegal fix srd command");
238 if (bigexist && biggroup < 0) error->all(FLERR, "Could not find fix srd group ID");
239 if (gridsrd <= 0.0) error->all(FLERR, "Illegal fix srd command");
240 if (temperature_srd <= 0.0) error->all(FLERR, "Illegal fix srd command");
241 if (seed <= 0) error->all(FLERR, "Illegal fix srd command");
242 if (radfactor <= 0.0) error->all(FLERR, "Illegal fix srd command");
243 if (maxbounceallow < 0) error->all(FLERR, "Illegal fix srd command");
244 if (lamdaflag && lamda <= 0.0) error->all(FLERR, "Illegal fix srd command");
245 if (gridsearch <= 0.0) error->all(FLERR, "Illegal fix srd command");
246 if (cubictol < 0.0 || cubictol > 1.0) error->all(FLERR, "Illegal fix srd command");
247 if ((shiftuser == SHIFT_YES || shiftuser == SHIFT_POSSIBLE) && shiftseed <= 0)
248 error->all(FLERR, "Illegal fix srd command");
249
250 // initialize Marsaglia RNG with processor-unique seed
251
252 me = comm->me;
253 nprocs = comm->nprocs;
254
255 random = new RanMars(lmp, seed + me);
256
257 // if requested, initialize shift RNG, same on every proc
258
259 if (shiftuser == SHIFT_YES || shiftuser == SHIFT_POSSIBLE)
260 randomshift = new RanPark(lmp, shiftseed);
261 else
262 randomshift = nullptr;
263
264 // initialize data structs and flags
265
266 if (bigexist)
267 biggroupbit = group->bitmask[biggroup];
268 else
269 biggroupbit = 0;
270
271 nmax = 0;
272 binhead = nullptr;
273 maxbin1 = 0;
274 binnext = nullptr;
275 maxbuf = 0;
276 sbuf1 = sbuf2 = rbuf1 = rbuf2 = nullptr;
277
278 shifts[0].maxvbin = shifts[1].maxvbin = 0;
279 shifts[0].vbin = shifts[1].vbin = nullptr;
280
281 shifts[0].maxbinsq = shifts[1].maxbinsq = 0;
282 for (int ishift = 0; ishift < 2; ishift++)
283 for (int iswap = 0; iswap < 6; iswap++)
284 shifts[ishift].bcomm[iswap].sendlist = shifts[ishift].bcomm[iswap].recvlist = nullptr;
285
286 maxbin2 = 0;
287 nbinbig = nullptr;
288 binbig = nullptr;
289 binsrd = nullptr;
290
291 nstencil = maxstencil = 0;
292 stencil = nullptr;
293
294 maxbig = 0;
295 biglist = nullptr;
296
297 stats_flag = 1;
298 for (int i = 0; i < size_vector; i++) stats_all[i] = 0.0;
299
300 initflag = 0;
301
302 srd_bin_temp = 0.0;
303 srd_bin_count = 0;
304
305 // atom style pointers to particles that store bonus info
306
307 avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
308 avec_line = (AtomVecLine *) atom->style_match("line");
309 avec_tri = (AtomVecTri *) atom->style_match("tri");
310
311 // fix parameters
312
313 if (collidestyle == SLIP)
314 comm_reverse = 3;
315 else
316 comm_reverse = 6;
317 force_reneighbor = 1;
318 }
319
320 /* ---------------------------------------------------------------------- */
321
~FixSRD()322 FixSRD::~FixSRD()
323 {
324 delete random;
325 delete randomshift;
326
327 memory->destroy(binhead);
328 memory->destroy(binnext);
329 memory->destroy(sbuf1);
330 memory->destroy(sbuf2);
331 memory->destroy(rbuf1);
332 memory->destroy(rbuf2);
333
334 memory->sfree(shifts[0].vbin);
335 memory->sfree(shifts[1].vbin);
336 for (int ishift = 0; ishift < 2; ishift++)
337 for (int iswap = 0; iswap < 6; iswap++) {
338 memory->destroy(shifts[ishift].bcomm[iswap].sendlist);
339 memory->destroy(shifts[ishift].bcomm[iswap].recvlist);
340 }
341
342 memory->destroy(nbinbig);
343 memory->destroy(binbig);
344 memory->destroy(binsrd);
345 memory->destroy(stencil);
346 memory->sfree(biglist);
347 }
348
349 /* ---------------------------------------------------------------------- */
350
setmask()351 int FixSRD::setmask()
352 {
353 int mask = 0;
354 mask |= PRE_NEIGHBOR;
355 mask |= POST_FORCE;
356 return mask;
357 }
358
359 /* ---------------------------------------------------------------------- */
360
init()361 void FixSRD::init()
362 {
363 // error checks
364
365 if (force->newton_pair == 0) error->all(FLERR, "Fix srd requires newton pair on");
366 if (bigexist && comm->ghost_velocity == 0)
367 error->all(FLERR, "Fix srd requires ghost atoms store velocity");
368 if (bigexist && collidestyle == NOSLIP && !atom->torque_flag)
369 error->all(FLERR, "Fix srd no-slip requires atom attribute torque");
370 if (initflag && update->dt != dt_big)
371 error->all(FLERR, "Cannot change timestep once fix srd is setup");
372 if (comm->style != 0)
373 error->universe_all(FLERR, "Fix srd can only currently be used with comm_style brick");
374
375 // orthogonal vs triclinic simulation box
376 // could be static or shearing box
377
378 triclinic = domain->triclinic;
379
380 // wallexist = 1 if SRD wall(s) are defined
381
382 wallexist = 0;
383 for (int m = 0; m < modify->nfix; m++) {
384 if (strcmp(modify->fix[m]->style, "wall/srd") == 0) {
385 if (wallexist) error->all(FLERR, "Cannot use fix wall/srd more than once");
386 wallexist = 1;
387 wallfix = (FixWallSRD *) modify->fix[m];
388 nwall = wallfix->nwall;
389 wallvarflag = wallfix->varflag;
390 wallwhich = wallfix->wallwhich;
391 xwall = wallfix->xwall;
392 xwallhold = wallfix->xwallhold;
393 vwall = wallfix->vwall;
394 fwall = wallfix->fwall;
395 walltrigger = 0.5 * neighbor->skin;
396 if (wallfix->overlap && overlap == 0 && me == 0)
397 error->warning(FLERR, "Fix SRD walls overlap but fix srd overlap not set");
398 }
399 }
400
401 // set change_flags if box size or shape changes
402
403 change_size = change_shape = deformflag = 0;
404 if (domain->nonperiodic == 2) change_size = 1;
405
406 Fix **fixes = modify->fix;
407 for (int i = 0; i < modify->nfix; i++) {
408 if (fixes[i]->box_change & BOX_CHANGE_SIZE) change_size = 1;
409 if (fixes[i]->box_change & BOX_CHANGE_SHAPE) change_shape = 1;
410 if (strcmp(fixes[i]->style, "deform") == 0) {
411 deformflag = 1;
412 FixDeform *deform = (FixDeform *) modify->fix[i];
413 if ((deform->box_change & BOX_CHANGE_SHAPE) && deform->remapflag != Domain::V_REMAP)
414 error->all(FLERR, "Using fix srd with inconsistent fix deform remap option");
415 }
416 }
417
418 if (deformflag && tstat == 0 && me == 0)
419 error->warning(FLERR, "Using fix srd with box deformation but no SRD thermostat");
420
421 // parameterize based on current box volume
422
423 dimension = domain->dimension;
424 parameterize();
425
426 // limit initial SRD velocities if necessary
427
428 double **v = atom->v;
429 int *mask = atom->mask;
430 int nlocal = atom->nlocal;
431
432 double vsq;
433 nrescale = 0;
434
435 for (int i = 0; i < nlocal; i++)
436 if (mask[i] & groupbit) {
437 vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
438 if (vsq > vmaxsq) {
439 nrescale++;
440 MathExtra::scale3(vmax / sqrt(vsq), v[i]);
441 }
442 }
443
444 int all;
445 MPI_Allreduce(&nrescale, &all, 1, MPI_INT, MPI_SUM, world);
446 if (me == 0) utils::logmesg(lmp, " # of rescaled SRD velocities = {}\n", all);
447
448 velocity_stats(igroup);
449 if (bigexist) velocity_stats(biggroup);
450
451 // zero per-run stats
452
453 nrescale = 0;
454 bouncemaxnum = 0;
455 bouncemax = 0;
456 reneighcount = 0;
457 initflag = 1;
458
459 next_reneighbor = -1;
460 }
461
462 /* ---------------------------------------------------------------------- */
463
setup(int)464 void FixSRD::setup(int /*vflag*/)
465 {
466 setup_bounds();
467
468 if (dist_srd_reneigh < nevery * dt_big * vmax && me == 0)
469 error->warning(FLERR, "Fix srd SRD moves may trigger frequent reneighboring");
470
471 // setup search bins and search stencil based on these distances
472
473 if (bigexist || wallexist) {
474 setup_search_bins();
475 setup_search_stencil();
476 } else
477 nbins2 = 0;
478
479 // perform first binning of SRD and big particles and walls
480 // set reneighflag to turn off SRD rotation
481 // don't do SRD rotation in setup, only during timestepping
482
483 reneighflag = BIG_MOVE;
484 pre_neighbor();
485 }
486
487 /* ----------------------------------------------------------------------
488 assign SRD particles to bins
489 assign big particles to all bins they overlap
490 ------------------------------------------------------------------------- */
491
pre_neighbor()492 void FixSRD::pre_neighbor()
493 {
494 int i, j, m, ix, iy, iz, jx, jy, jz, ibin, jbin, lo, hi;
495 double rsq, cutbinsq;
496
497 // grow SRD per-atom bin arrays if necessary
498
499 if (atom->nmax > nmax) {
500 nmax = atom->nmax;
501 memory->destroy(binsrd);
502 memory->destroy(binnext);
503 memory->create(binsrd, nmax, "fix/srd:binsrd");
504 memory->create(binnext, nmax, "fix/srd:binnext");
505 }
506
507 // setup and grow BIG info list if necessary
508 // set index ptrs to BIG particles and to WALLS
509 // big_static() adds static properties to info list
510
511 if (bigexist || wallexist) {
512 if (bigexist) {
513 if (biggroup == atom->firstgroup)
514 nbig = atom->nfirst + atom->nghost;
515 else {
516 int *mask = atom->mask;
517 int nlocal = atom->nlocal;
518 nbig = atom->nghost;
519 for (i = 0; i < nlocal; i++)
520 if (mask[i] & biggroupbit) nbig++;
521 }
522 } else
523 nbig = 0;
524
525 int ninfo = nbig;
526 if (wallexist) ninfo += nwall;
527
528 if (ninfo > maxbig) {
529 maxbig = ninfo;
530 memory->destroy(biglist);
531 biglist = (Big *) memory->smalloc(maxbig * sizeof(Big), "fix/srd:biglist");
532 }
533
534 if (bigexist) {
535 int *mask = atom->mask;
536 int nlocal = atom->nlocal;
537 if (biggroup == atom->firstgroup) nlocal = atom->nfirst;
538 nbig = 0;
539 for (i = 0; i < nlocal; i++)
540 if (mask[i] & biggroupbit) biglist[nbig++].index = i;
541 int nall = atom->nlocal + atom->nghost;
542 for (i = atom->nlocal; i < nall; i++)
543 if (mask[i] & biggroupbit) biglist[nbig++].index = i;
544 big_static();
545 }
546
547 if (wallexist) {
548 for (m = 0; m < nwall; m++) {
549 biglist[nbig + m].index = m;
550 biglist[nbig + m].type = WALL;
551 }
552 wallfix->wall_params(1);
553 }
554 }
555
556 // if simulation box size changes, reset velocity bins
557 // if big particles exist, reset search bins if box size or shape changes,
558 // b/c finite-size particles will overlap different bins as the box tilts
559
560 if (change_size) setup_bounds();
561 if (change_size) setup_velocity_bins();
562 if ((change_size || change_shape) && (bigexist || wallexist)) {
563 setup_search_bins();
564 setup_search_stencil();
565 }
566
567 // map each owned & ghost big particle to search bins it overlaps
568 // zero out bin counters for big particles
569 // if firstgroup is defined, only loop over first and ghost particles
570 // for each big particle: loop over stencil to find overlap bins
571
572 int *mask = atom->mask;
573 double **x = atom->x;
574 int nlocal = atom->nlocal;
575 int nall = nlocal + atom->nghost;
576 int nfirst = nlocal;
577 if (bigexist && biggroup == atom->firstgroup) nfirst = atom->nfirst;
578
579 if (bigexist || wallexist)
580 for (i = 0; i < nbins2; i++) nbinbig[i] = 0;
581
582 if (bigexist) {
583 i = nbig = 0;
584 while (i < nall) {
585 if (mask[i] & biggroupbit) {
586 ix = static_cast<int>((x[i][0] - xblo2) * bininv2x);
587 iy = static_cast<int>((x[i][1] - yblo2) * bininv2y);
588 iz = static_cast<int>((x[i][2] - zblo2) * bininv2z);
589 ibin = iz * nbin2y * nbin2x + iy * nbin2x + ix;
590
591 if (ix < 0 || ix >= nbin2x || iy < 0 || iy >= nbin2y || iz < 0 || iz >= nbin2z)
592 error->one(FLERR, "Fix SRD: bad search bin assignment");
593
594 cutbinsq = biglist[nbig].cutbinsq;
595 for (j = 0; j < nstencil; j++) {
596 jx = ix + stencil[j][0];
597 jy = iy + stencil[j][1];
598 jz = iz + stencil[j][2];
599
600 if (jx < 0 || jx >= nbin2x || jy < 0 || jy >= nbin2y || jz < 0 || jz >= nbin2z)
601 error->one(FLERR,
602 "Fix SRD: bad stencil bin for big particle\n"
603 "Big particle {} {} {:.8} {:.8} {:.8}\n"
604 "Bin indices: {} {} {}, {} {} {}, {} {} {}\n",
605 atom->tag[i], i, x[i][0], x[i][1], x[i][2], ix, iy, iz, jx, jy, jz, nbin2x,
606 nbin2y, nbin2z);
607 rsq = point_bin_distance(x[i], jx, jy, jz);
608 if (rsq < cutbinsq) {
609 jbin = ibin + stencil[j][3];
610 if (nbinbig[jbin] == ATOMPERBIN)
611 error->one(FLERR, "Fix SRD: too many big particles in bin");
612 binbig[jbin][nbinbig[jbin]++] = nbig;
613 }
614 }
615 nbig++;
616 }
617
618 i++;
619 if (i == nfirst) i = nlocal;
620 }
621 }
622
623 // map each wall to search bins it covers, up to non-periodic boundary
624 // if wall moves, add walltrigger to its position
625 // this insures it is added to all search bins it may move into
626 // may not overlap any of my search bins
627
628 if (wallexist) {
629 double delta = 0.0;
630 if (wallvarflag) delta = walltrigger;
631
632 for (m = 0; m < nwall; m++) {
633 int dim = wallwhich[m] / 2;
634 int side = wallwhich[m] % 2;
635
636 if (dim == 0) {
637 if (side == 0) {
638 hi = static_cast<int>((xwall[m] + delta - xblo2) * bininv2x);
639 if (hi < 0) continue;
640 if (hi >= nbin2x) error->all(FLERR, "Fix SRD: bad search bin assignment");
641 lo = 0;
642 } else {
643 lo = static_cast<int>((xwall[m] - delta - xblo2) * bininv2x);
644 if (lo >= nbin2x) continue;
645 if (lo < 0) error->all(FLERR, "Fix SRD: bad search bin assignment");
646 hi = nbin2x - 1;
647 }
648
649 for (ix = lo; ix <= hi; ix++)
650 for (iy = 0; iy < nbin2y; iy++)
651 for (iz = 0; iz < nbin2z; iz++) {
652 ibin = iz * nbin2y * nbin2x + iy * nbin2x + ix;
653 if (nbinbig[ibin] == ATOMPERBIN) error->all(FLERR, "Fix SRD: too many walls in bin");
654 binbig[ibin][nbinbig[ibin]++] = nbig + m;
655 }
656
657 } else if (dim == 1) {
658 if (side == 0) {
659 hi = static_cast<int>((xwall[m] + delta - yblo2) * bininv2y);
660 if (hi < 0) continue;
661 if (hi >= nbin2y) error->all(FLERR, "Fix SRD: bad search bin assignment");
662 lo = 0;
663 } else {
664 lo = static_cast<int>((xwall[m] - delta - yblo2) * bininv2y);
665 if (lo >= nbin2y) continue;
666 if (lo < 0) error->all(FLERR, "Fix SRD: bad search bin assignment");
667 hi = nbin2y - 1;
668 }
669
670 for (iy = lo; iy <= hi; iy++)
671 for (ix = 0; ix < nbin2x; ix++)
672 for (iz = 0; iz < nbin2z; iz++) {
673 ibin = iz * nbin2y * nbin2x + iy * nbin2x + ix;
674 if (nbinbig[ibin] == ATOMPERBIN) error->all(FLERR, "Fix SRD: too many walls in bin");
675 binbig[ibin][nbinbig[ibin]++] = nbig + m;
676 }
677
678 } else if (dim == 2) {
679 if (side == 0) {
680 hi = static_cast<int>((xwall[m] + delta - zblo2) * bininv2z);
681 if (hi < 0) continue;
682 if (hi >= nbin2z) error->all(FLERR, "Fix SRD: bad search bin assignment");
683 lo = 0;
684 } else {
685 lo = static_cast<int>((xwall[m] - delta - zblo2) * bininv2z);
686 if (lo >= nbin2z) continue;
687 if (lo < 0) error->all(FLERR, "Fix SRD: bad search bin assignment");
688 hi = nbin2z - 1;
689 }
690
691 for (iz = lo; iz < hi; iz++)
692 for (ix = 0; ix < nbin2x; ix++)
693 for (iy = 0; iy < nbin2y; iy++) {
694 ibin = iz * nbin2y * nbin2x + iy * nbin2x + ix;
695 if (nbinbig[ibin] == ATOMPERBIN) error->all(FLERR, "Fix SRD: too many walls in bin");
696 binbig[ibin][nbinbig[ibin]++] = nbig + m;
697 }
698 }
699 }
700 }
701
702 // rotate SRD velocities on SRD timestep
703 // done now since all SRDs are currently inside my sub-domain
704
705 if (reneighflag == SRD_ROTATE) reset_velocities();
706
707 // log stats if reneighboring occurred b/c SRDs moved too far
708
709 if (reneighflag == SRD_MOVE) reneighcount++;
710 reneighflag = BIG_MOVE;
711 }
712
713 /* ----------------------------------------------------------------------
714 advect SRD particles and detect collisions between SRD and BIG particles
715 when collision occurs, change x,v of SRD, force,torque of BIG particle
716 ------------------------------------------------------------------------- */
717
post_force(int)718 void FixSRD::post_force(int /*vflag*/)
719 {
720 int i, m, ix, iy, iz;
721
722 // zero per-timestep stats
723
724 stats_flag = 0;
725 ncheck = ncollide = nbounce = ninside = 0;
726
727 // zero ghost forces & torques on BIG particles
728
729 double **f = atom->f;
730 double **torque = atom->torque;
731 int nlocal = atom->nlocal;
732 int nall = nlocal + atom->nghost;
733 if (bigexist == 0) nall = 0;
734
735 for (i = nlocal; i < nall; i++) f[i][0] = f[i][1] = f[i][2] = 0.0;
736
737 if (collidestyle == NOSLIP)
738 for (i = nlocal; i < nall; i++) torque[i][0] = torque[i][1] = torque[i][2] = 0.0;
739
740 // advect SRD particles
741 // assign to search bins if big particles or walls exist
742
743 int *mask = atom->mask;
744 double **x = atom->x;
745 double **v = atom->v;
746
747 if (bigexist || wallexist) {
748 for (i = 0; i < nlocal; i++)
749 if (mask[i] & groupbit) {
750 x[i][0] += dt_big * v[i][0];
751 x[i][1] += dt_big * v[i][1];
752 x[i][2] += dt_big * v[i][2];
753
754 ix = static_cast<int>((x[i][0] - xblo2) * bininv2x);
755 iy = static_cast<int>((x[i][1] - yblo2) * bininv2y);
756 iz = static_cast<int>((x[i][2] - zblo2) * bininv2z);
757 binsrd[i] = iz * nbin2y * nbin2x + iy * nbin2x + ix;
758
759 if (ix < 0 || ix >= nbin2x || iy < 0 || iy >= nbin2y || iz < 0 || iz >= nbin2z)
760 error->one(FLERR,
761 "Fix SRD: bad bin assignment for SRD advection\n"
762 "SRD particle {} on step {}\n"
763 "v = {:.8} {:.8} {:.8}\nx = {:.8} {:.8} {:.8}\n"
764 "ix,iy,iz nx,ny,nz = {} {} {} {} {} {}\n",
765 atom->tag[i], update->ntimestep, v[i][0], v[i][1], v[i][2], x[i][0], x[i][1],
766 x[i][2], ix, iy, iz, nbin2x, nbin2y, nbin2z);
767 }
768 } else {
769 for (i = 0; i < nlocal; i++)
770 if (mask[i] & groupbit) {
771 x[i][0] += dt_big * v[i][0];
772 x[i][1] += dt_big * v[i][1];
773 x[i][2] += dt_big * v[i][2];
774 }
775 }
776
777 // detect collision of SRDs with BIG particles or walls
778
779 if (bigexist || wallexist) {
780 if (bigexist) big_dynamic();
781 if (wallexist) wallfix->wall_params(0);
782 if (overlap)
783 collisions_multi();
784 else
785 collisions_single();
786 }
787
788 // reverse communicate forces & torques on BIG particles
789
790 if (bigexist) {
791 flocal = f;
792 tlocal = torque;
793 comm->reverse_comm_fix(this);
794 }
795
796 // if any SRD particle has moved too far, trigger reneigh on next step
797 // for triclinic, perform check in lamda units
798
799 int flag = 0;
800
801 if (triclinic) domain->x2lamda(nlocal);
802 for (i = 0; i < nlocal; i++)
803 if (mask[i] & groupbit) {
804 if (x[i][0] < srdlo_reneigh[0] || x[i][0] > srdhi_reneigh[0] || x[i][1] < srdlo_reneigh[1] ||
805 x[i][1] > srdhi_reneigh[1] || x[i][2] < srdlo_reneigh[2] || x[i][2] > srdhi_reneigh[2])
806 flag = 1;
807 }
808 if (triclinic) domain->lamda2x(nlocal);
809
810 int flagall;
811 MPI_Allreduce(&flag, &flagall, 1, MPI_INT, MPI_SUM, world);
812 if (flagall) {
813 next_reneighbor = update->ntimestep + 1;
814 reneighflag = SRD_MOVE;
815 }
816
817 // if wall has moved too far, trigger reneigh on next step
818 // analogous to neighbor check for big particle moving 1/2 of skin distance
819
820 if (wallexist) {
821 for (m = 0; m < nwall; m++)
822 if (fabs(xwall[m] - xwallhold[m]) > walltrigger) next_reneighbor = update->ntimestep + 1;
823 }
824
825 // if next timestep is SRD timestep, trigger reneigh
826
827 if ((update->ntimestep + 1) % nevery == 0) {
828 next_reneighbor = update->ntimestep + 1;
829 reneighflag = SRD_ROTATE;
830 }
831 }
832
833 /* ----------------------------------------------------------------------
834 reset SRD velocities
835 may perform random shifting by up to 1/2 bin in each dimension
836 called at pre-neighbor stage when all SRDs are now inside my sub-domain
837 if tstat, then thermostat SRD particles as well, including streaming effects
838 ------------------------------------------------------------------------- */
839
reset_velocities()840 void FixSRD::reset_velocities()
841 {
842 int i, j, n, ix, iy, iz, ibin, axis, sign, irandom;
843 double u[3], vsum[3];
844 double vsq, tbin, scale;
845 double *vave, *xlamda;
846 double vstream[3];
847
848 // if requested, perform a dynamic shift of bin positions
849
850 if (shiftflag) {
851 double *boxlo;
852 if (triclinic == 0)
853 boxlo = domain->boxlo;
854 else
855 boxlo = domain->boxlo_lamda;
856 shifts[1].corner[0] = boxlo[0] - binsize1x * randomshift->uniform();
857 shifts[1].corner[1] = boxlo[1] - binsize1y * randomshift->uniform();
858 if (dimension == 3)
859 shifts[1].corner[2] = boxlo[2] - binsize1z * randomshift->uniform();
860 else
861 shifts[1].corner[2] = boxlo[2];
862 setup_velocity_shift(1, 1);
863 }
864
865 double *corner = shifts[shiftflag].corner;
866 int *binlo = shifts[shiftflag].binlo;
867 int *binhi = shifts[shiftflag].binhi;
868 int nbins = shifts[shiftflag].nbins;
869 int nbinx = shifts[shiftflag].nbinx;
870 int nbiny = shifts[shiftflag].nbiny;
871 BinAve *vbin = shifts[shiftflag].vbin;
872
873 // binhead = 1st SRD particle in each bin
874 // binnext = index of next particle in bin
875 // bin assignment is done in lamda units for triclinic
876
877 int *mask = atom->mask;
878 double **x = atom->x;
879 double **v = atom->v;
880 int nlocal = atom->nlocal;
881
882 if (triclinic) domain->x2lamda(nlocal);
883
884 for (i = 0; i < nbins; i++) binhead[i] = -1;
885
886 for (i = 0; i < nlocal; i++)
887 if (mask[i] & groupbit) {
888 ix = static_cast<int>((x[i][0] - corner[0]) * bininv1x);
889 ix = MAX(ix, binlo[0]);
890 ix = MIN(ix, binhi[0]);
891 iy = static_cast<int>((x[i][1] - corner[1]) * bininv1y);
892 iy = MAX(iy, binlo[1]);
893 iy = MIN(iy, binhi[1]);
894 iz = static_cast<int>((x[i][2] - corner[2]) * bininv1z);
895 iz = MAX(iz, binlo[2]);
896 iz = MIN(iz, binhi[2]);
897
898 ibin = (iz - binlo[2]) * nbiny * nbinx + (iy - binlo[1]) * nbinx + (ix - binlo[0]);
899 binnext[i] = binhead[ibin];
900 binhead[ibin] = i;
901 }
902
903 if (triclinic) domain->lamda2x(nlocal);
904
905 // for each bin I have particles contributing to:
906 // compute summed v of particles in that bin
907 // if I own the bin, set its random value, else set to 0.0
908
909 for (i = 0; i < nbins; i++) {
910 n = 0;
911 vsum[0] = vsum[1] = vsum[2] = 0.0;
912 for (j = binhead[i]; j >= 0; j = binnext[j]) {
913 vsum[0] += v[j][0];
914 vsum[1] += v[j][1];
915 vsum[2] += v[j][2];
916 n++;
917 }
918
919 vbin[i].vsum[0] = vsum[0];
920 vbin[i].vsum[1] = vsum[1];
921 vbin[i].vsum[2] = vsum[2];
922 vbin[i].n = n;
923 if (vbin[i].owner)
924 vbin[i].random = random->uniform();
925 else
926 vbin[i].random = 0.0;
927 }
928
929 // communicate bin info for bins which more than 1 proc contribute to
930
931 if (shifts[shiftflag].commflag) vbin_comm(shiftflag);
932
933 // for each bin I have particles contributing to:
934 // compute vave over particles in bin
935 // thermal velocity of each particle = v - vave
936 // rotate thermal vel of each particle around one of 6 random axes
937 // add vave back to each particle
938 // thermostat if requested:
939 // if no deformation, rescale thermal vel to temperature
940 // if deformation, rescale thermal vel and change vave to vstream
941 // these are settings for extra dof_temp, dof_tstat to subtract
942 // (not sure why these settings work best)
943 // no deformation, no tstat: dof_temp = 1
944 // yes deformation, no tstat: doesn't matter, system will not equilibrate
945 // no deformation, yes tstat: dof_temp = dof_tstat = 1
946 // yes deformation, yes tstat: dof_temp = dof_tstat = 0
947 // accumulate final T_srd for each bin I own
948
949 double tfactor = force->mvv2e * mass_srd / (dimension * force->boltz);
950 int dof_temp = 1;
951 int dof_tstat;
952 if (tstat) {
953 if (deformflag)
954 dof_tstat = dof_temp = 0;
955 else
956 dof_tstat = 1;
957 }
958
959 srd_bin_temp = 0.0;
960 srd_bin_count = 0;
961
962 if (dimension == 2) axis = 2;
963 double *h_rate = domain->h_rate;
964 double *h_ratelo = domain->h_ratelo;
965
966 for (i = 0; i < nbins; i++) {
967 vbin[i].value[0] = 0.0;
968 n = vbin[i].n;
969 if (n == 0) continue;
970 vave = vbin[i].vsum;
971 vave[0] /= n;
972 vave[1] /= n;
973 vave[2] /= n;
974
975 irandom = static_cast<int>(6.0 * vbin[i].random);
976 sign = irandom % 2;
977 if (dimension == 3) axis = irandom / 2;
978
979 vsq = 0.0;
980 for (j = binhead[i]; j >= 0; j = binnext[j]) {
981 if (axis == 0) {
982 u[0] = v[j][0] - vave[0];
983 u[1] = sign ? v[j][2] - vave[2] : vave[2] - v[j][2];
984 u[2] = sign ? vave[1] - v[j][1] : v[j][1] - vave[1];
985 } else if (axis == 1) {
986 u[1] = v[j][1] - vave[1];
987 u[0] = sign ? v[j][2] - vave[2] : vave[2] - v[j][2];
988 u[2] = sign ? vave[0] - v[j][0] : v[j][0] - vave[0];
989 } else {
990 u[2] = v[j][2] - vave[2];
991 u[1] = sign ? v[j][0] - vave[0] : vave[0] - v[j][0];
992 u[0] = sign ? vave[1] - v[j][1] : v[j][1] - vave[1];
993 }
994 vsq += u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
995 v[j][0] = u[0] + vave[0];
996 v[j][1] = u[1] + vave[1];
997 v[j][2] = u[2] + vave[2];
998 }
999
1000 if (n > 1) vbin[i].value[0] = vsq;
1001 }
1002
1003 if (shifts[shiftflag].commflag) xbin_comm(shiftflag, 1);
1004
1005 if (tstat) {
1006 for (i = 0; i < nbins; i++) {
1007 n = vbin[i].n;
1008 if (n <= 1) continue;
1009
1010 // vsum is already average velocity
1011
1012 vave = vbin[i].vsum;
1013
1014 if (deformflag) {
1015 xlamda = vbin[i].xctr;
1016 vstream[0] =
1017 h_rate[0] * xlamda[0] + h_rate[5] * xlamda[1] + h_rate[4] * xlamda[2] + h_ratelo[0];
1018 vstream[1] = h_rate[1] * xlamda[1] + h_rate[3] * xlamda[2] + h_ratelo[1];
1019 vstream[2] = h_rate[2] * xlamda[2] + h_ratelo[2];
1020 } else {
1021 vstream[0] = vave[0];
1022 vstream[1] = vave[1];
1023 vstream[2] = vave[2];
1024 }
1025
1026 // tbin = thermal temperature of particles in bin
1027 // scale = scale factor for thermal velocity
1028
1029 tbin = vbin[i].value[0] / (n - dof_tstat) * tfactor;
1030 scale = sqrt(temperature_srd / tbin);
1031 vsq = 0.0;
1032 for (j = binhead[i]; j >= 0; j = binnext[j]) {
1033 u[0] = (v[j][0] - vave[0]) * scale;
1034 u[1] = (v[j][1] - vave[1]) * scale;
1035 u[2] = (v[j][2] - vave[2]) * scale;
1036 vsq += u[0] * u[0] + u[1] * u[1] + u[2] * u[2];
1037 v[j][0] = u[0] + vstream[0];
1038 v[j][1] = u[1] + vstream[1];
1039 v[j][2] = u[2] + vstream[2];
1040 }
1041 vbin[i].value[0] = vsq;
1042 }
1043
1044 if (shifts[shiftflag].commflag) xbin_comm(shiftflag, 1);
1045 }
1046
1047 for (i = 0; i < nbins; i++) {
1048 if (vbin[i].owner) {
1049 if (vbin[i].n > 1) {
1050 srd_bin_temp += vbin[i].value[0] / (vbin[i].n - dof_temp);
1051 srd_bin_count++;
1052 }
1053 }
1054 }
1055
1056 srd_bin_temp *= tfactor;
1057
1058 // rescale any too-large velocities
1059
1060 if (rescale_rotate) {
1061 for (i = 0; i < nlocal; i++)
1062 if (mask[i] & groupbit) {
1063 vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
1064 if (vsq > vmaxsq) {
1065 nrescale++;
1066 MathExtra::scale3(vmax / sqrt(vsq), v[i]);
1067 }
1068 }
1069 }
1070 }
1071
1072 /* ----------------------------------------------------------------------
1073 communicate summed particle info for bins that overlap 1 or more procs
1074 ------------------------------------------------------------------------- */
1075
vbin_comm(int ishift)1076 void FixSRD::vbin_comm(int ishift)
1077 {
1078 BinComm *bcomm1, *bcomm2;
1079 MPI_Request request1, request2;
1080
1081 // send/recv bins in both directions in each dimension
1082 // don't send if nsend = 0
1083 // due to static bins aligning with proc boundary
1084 // due to dynamic bins across non-periodic global boundary
1085 // copy to self if sendproc = me
1086 // MPI send to another proc if sendproc != me
1087 // don't recv if nrecv = 0
1088 // copy from self if recvproc = me
1089 // MPI recv from another proc if recvproc != me
1090
1091 BinAve *vbin = shifts[ishift].vbin;
1092 int *procgrid = comm->procgrid;
1093
1094 int iswap = 0;
1095 for (int idim = 0; idim < dimension; idim++) {
1096 bcomm1 = &shifts[ishift].bcomm[iswap++];
1097 bcomm2 = &shifts[ishift].bcomm[iswap++];
1098
1099 if (procgrid[idim] == 1) {
1100 if (bcomm1->nsend) vbin_pack(vbin, bcomm1->nsend, bcomm1->sendlist, sbuf1);
1101 if (bcomm2->nsend) vbin_pack(vbin, bcomm2->nsend, bcomm2->sendlist, sbuf2);
1102 if (bcomm1->nrecv) vbin_unpack(sbuf1, vbin, bcomm1->nrecv, bcomm1->recvlist);
1103 if (bcomm2->nrecv) vbin_unpack(sbuf2, vbin, bcomm2->nrecv, bcomm2->recvlist);
1104
1105 } else {
1106 if (bcomm1->nrecv)
1107 MPI_Irecv(rbuf1, bcomm1->nrecv * VBINSIZE, MPI_DOUBLE, bcomm1->recvproc, 0, world,
1108 &request1);
1109 if (bcomm2->nrecv)
1110 MPI_Irecv(rbuf2, bcomm2->nrecv * VBINSIZE, MPI_DOUBLE, bcomm2->recvproc, 0, world,
1111 &request2);
1112 if (bcomm1->nsend) {
1113 vbin_pack(vbin, bcomm1->nsend, bcomm1->sendlist, sbuf1);
1114 MPI_Send(sbuf1, bcomm1->nsend * VBINSIZE, MPI_DOUBLE, bcomm1->sendproc, 0, world);
1115 }
1116 if (bcomm2->nsend) {
1117 vbin_pack(vbin, bcomm2->nsend, bcomm2->sendlist, sbuf2);
1118 MPI_Send(sbuf2, bcomm2->nsend * VBINSIZE, MPI_DOUBLE, bcomm2->sendproc, 0, world);
1119 }
1120 if (bcomm1->nrecv) {
1121 MPI_Wait(&request1, MPI_STATUS_IGNORE);
1122 vbin_unpack(rbuf1, vbin, bcomm1->nrecv, bcomm1->recvlist);
1123 }
1124 if (bcomm2->nrecv) {
1125 MPI_Wait(&request2, MPI_STATUS_IGNORE);
1126 vbin_unpack(rbuf2, vbin, bcomm2->nrecv, bcomm2->recvlist);
1127 }
1128 }
1129 }
1130 }
1131
1132 /* ----------------------------------------------------------------------
1133 pack velocity bin data into a message buffer for sending
1134 ------------------------------------------------------------------------- */
1135
vbin_pack(BinAve * vbin,int n,int * list,double * buf)1136 void FixSRD::vbin_pack(BinAve *vbin, int n, int *list, double *buf)
1137 {
1138 int j;
1139 int m = 0;
1140 for (int i = 0; i < n; i++) {
1141 j = list[i];
1142 buf[m++] = vbin[j].n;
1143 buf[m++] = vbin[j].vsum[0];
1144 buf[m++] = vbin[j].vsum[1];
1145 buf[m++] = vbin[j].vsum[2];
1146 buf[m++] = vbin[j].random;
1147 }
1148 }
1149
1150 /* ----------------------------------------------------------------------
1151 unpack velocity bin data from a message buffer and sum values to my bins
1152 ------------------------------------------------------------------------- */
1153
vbin_unpack(double * buf,BinAve * vbin,int n,int * list)1154 void FixSRD::vbin_unpack(double *buf, BinAve *vbin, int n, int *list)
1155 {
1156 int j;
1157 int m = 0;
1158 for (int i = 0; i < n; i++) {
1159 j = list[i];
1160 vbin[j].n += static_cast<int>(buf[m++]);
1161 vbin[j].vsum[0] += buf[m++];
1162 vbin[j].vsum[1] += buf[m++];
1163 vbin[j].vsum[2] += buf[m++];
1164 vbin[j].random += buf[m++];
1165 }
1166 }
1167
1168 /* ----------------------------------------------------------------------
1169 communicate summed particle vsq info for bins that overlap 1 or more procs
1170 ------------------------------------------------------------------------- */
1171
xbin_comm(int ishift,int nval)1172 void FixSRD::xbin_comm(int ishift, int nval)
1173 {
1174 BinComm *bcomm1, *bcomm2;
1175 MPI_Request request1, request2;
1176
1177 // send/recv bins in both directions in each dimension
1178 // don't send if nsend = 0
1179 // due to static bins aligning with proc boundary
1180 // due to dynamic bins across non-periodic global boundary
1181 // copy to self if sendproc = me
1182 // MPI send to another proc if sendproc != me
1183 // don't recv if nrecv = 0
1184 // copy from self if recvproc = me
1185 // MPI recv from another proc if recvproc != me
1186
1187 BinAve *vbin = shifts[ishift].vbin;
1188 int *procgrid = comm->procgrid;
1189
1190 int iswap = 0;
1191 for (int idim = 0; idim < dimension; idim++) {
1192 bcomm1 = &shifts[ishift].bcomm[iswap++];
1193 bcomm2 = &shifts[ishift].bcomm[iswap++];
1194
1195 if (procgrid[idim] == 1) {
1196 if (bcomm1->nsend) xbin_pack(vbin, bcomm1->nsend, bcomm1->sendlist, sbuf1, nval);
1197 if (bcomm2->nsend) xbin_pack(vbin, bcomm2->nsend, bcomm2->sendlist, sbuf2, nval);
1198 if (bcomm1->nrecv) xbin_unpack(sbuf1, vbin, bcomm1->nrecv, bcomm1->recvlist, nval);
1199 if (bcomm2->nrecv) xbin_unpack(sbuf2, vbin, bcomm2->nrecv, bcomm2->recvlist, nval);
1200
1201 } else {
1202 if (bcomm1->nrecv)
1203 MPI_Irecv(rbuf1, bcomm1->nrecv * nval, MPI_DOUBLE, bcomm1->recvproc, 0, world, &request1);
1204 if (bcomm2->nrecv)
1205 MPI_Irecv(rbuf2, bcomm2->nrecv * nval, MPI_DOUBLE, bcomm2->recvproc, 0, world, &request2);
1206 if (bcomm1->nsend) {
1207 xbin_pack(vbin, bcomm1->nsend, bcomm1->sendlist, sbuf1, nval);
1208 MPI_Send(sbuf1, bcomm1->nsend * nval, MPI_DOUBLE, bcomm1->sendproc, 0, world);
1209 }
1210 if (bcomm2->nsend) {
1211 xbin_pack(vbin, bcomm2->nsend, bcomm2->sendlist, sbuf2, nval);
1212 MPI_Send(sbuf2, bcomm2->nsend * nval, MPI_DOUBLE, bcomm2->sendproc, 0, world);
1213 }
1214 if (bcomm1->nrecv) {
1215 MPI_Wait(&request1, MPI_STATUS_IGNORE);
1216 xbin_unpack(rbuf1, vbin, bcomm1->nrecv, bcomm1->recvlist, nval);
1217 }
1218 if (bcomm2->nrecv) {
1219 MPI_Wait(&request2, MPI_STATUS_IGNORE);
1220 xbin_unpack(rbuf2, vbin, bcomm2->nrecv, bcomm2->recvlist, nval);
1221 }
1222 }
1223 }
1224 }
1225
1226 /* ----------------------------------------------------------------------
1227 pack velocity bin data into a message buffer for sending
1228 ------------------------------------------------------------------------- */
1229
xbin_pack(BinAve * vbin,int n,int * list,double * buf,int nval)1230 void FixSRD::xbin_pack(BinAve *vbin, int n, int *list, double *buf, int nval)
1231 {
1232 int j, k;
1233 int m = 0;
1234 for (int i = 0; i < n; i++) {
1235 j = list[i];
1236 for (k = 0; k < nval; k++) buf[m++] = vbin[j].value[k];
1237 }
1238 }
1239
1240 /* ----------------------------------------------------------------------
1241 unpack velocity bin data from a message buffer and sum values to my bins
1242 ------------------------------------------------------------------------- */
1243
xbin_unpack(double * buf,BinAve * vbin,int n,int * list,int nval)1244 void FixSRD::xbin_unpack(double *buf, BinAve *vbin, int n, int *list, int nval)
1245 {
1246 int j, k;
1247 int m = 0;
1248 for (int i = 0; i < n; i++) {
1249 j = list[i];
1250 for (k = 0; k < nval; k++) vbin[j].value[k] += buf[m++];
1251 }
1252 }
1253
1254 /* ----------------------------------------------------------------------
1255 detect all collisions between SRD and BIG particles or WALLS
1256 assume SRD can be inside at most one BIG particle or WALL at a time
1257 unoverlap SRDs for each collision
1258 ------------------------------------------------------------------------- */
1259
collisions_single()1260 void FixSRD::collisions_single()
1261 {
1262 int i, j, k, m, type, nbig, ibin, ibounce, inside, collide_flag;
1263 double dt, t_remain;
1264 double norm[3], xscoll[3], xbcoll[3], vsnew[3];
1265 Big *big;
1266
1267 // outer loop over SRD particles
1268 // inner loop over BIG particles or WALLS that overlap SRD particle bin
1269 // if overlap between SRD and BIG particle or wall:
1270 // for exact, compute collision pt in time
1271 // for inexact, push SRD to surf of BIG particle or WALL
1272 // update x,v of SRD and f,torque on BIG particle
1273 // re-bin SRD particle after collision
1274 // iterate until the SRD particle has no overlaps with BIG particles or WALLS
1275
1276 double **x = atom->x;
1277 double **v = atom->v;
1278 double **f = atom->f;
1279 double **torque = atom->torque;
1280 tagint *tag = atom->tag;
1281 int *mask = atom->mask;
1282 int nlocal = atom->nlocal;
1283
1284 for (i = 0; i < nlocal; i++) {
1285 if (!(mask[i] & groupbit)) continue;
1286
1287 ibin = binsrd[i];
1288 if (nbinbig[ibin] == 0) continue;
1289
1290 ibounce = 0;
1291 collide_flag = 1;
1292 dt = dt_big;
1293
1294 while (collide_flag) {
1295 nbig = nbinbig[ibin];
1296 if (ibounce == 0) ncheck += nbig;
1297
1298 collide_flag = 0;
1299 for (m = 0; m < nbig; m++) {
1300 k = binbig[ibin][m];
1301 big = &biglist[k];
1302 j = big->index;
1303 type = big->type;
1304
1305 if (type == SPHERE)
1306 inside = inside_sphere(x[i], x[j], big);
1307 else if (type == ELLIPSOID)
1308 inside = inside_ellipsoid(x[i], x[j], big);
1309 else
1310 inside = inside_wall(x[i], j);
1311
1312 if (inside) {
1313 if (exactflag) {
1314 if (type == SPHERE)
1315 t_remain = collision_sphere_exact(x[i], x[j], v[i], v[j], big, xscoll, xbcoll, norm);
1316 else if (type == ELLIPSOID)
1317 t_remain =
1318 collision_ellipsoid_exact(x[i], x[j], v[i], v[j], big, xscoll, xbcoll, norm);
1319 else
1320 t_remain = collision_wall_exact(x[i], j, v[i], xscoll, xbcoll, norm);
1321
1322 } else {
1323 t_remain = 0.5 * dt;
1324 if (type == SPHERE)
1325 collision_sphere_inexact(x[i], x[j], big, xscoll, xbcoll, norm);
1326 else if (type == ELLIPSOID)
1327 collision_ellipsoid_inexact(x[i], x[j], big, xscoll, xbcoll, norm);
1328 else
1329 collision_wall_inexact(x[i], j, xscoll, xbcoll, norm);
1330 }
1331
1332 #ifdef SRD_DEBUG
1333 if (update->ntimestep == SRD_DEBUG_TIMESTEP && tag[i] == SRD_DEBUG_ATOMID)
1334 print_collision(i, j, ibounce, t_remain, dt, xscoll, xbcoll, norm, type);
1335 #endif
1336
1337 if (t_remain > dt) {
1338 ninside++;
1339 if (insideflag == INSIDE_ERROR || insideflag == INSIDE_WARN) {
1340 std::string mesg;
1341 if (type != WALL)
1342 mesg = fmt::format("SRD particle {} started inside big particle {} on step {} "
1343 " bounce {}",
1344 tag[i], tag[j], update->ntimestep, ibounce + 1);
1345 else
1346 mesg = fmt::format("SRD particle {} started inside wall {} on step {} "
1347 "bounce {}",
1348 tag[i], j, update->ntimestep, ibounce + 1);
1349
1350 if (insideflag == INSIDE_ERROR)
1351 error->one(FLERR, mesg);
1352 else
1353 error->warning(FLERR, mesg);
1354 }
1355 break;
1356 }
1357
1358 if (collidestyle == SLIP) {
1359 if (type != WALL)
1360 slip(v[i], v[j], x[j], big, xscoll, norm, vsnew);
1361 else
1362 slip_wall(v[i], j, norm, vsnew);
1363 } else {
1364 if (type != WALL)
1365 noslip(v[i], v[j], x[j], big, -1, xscoll, norm, vsnew);
1366 else
1367 noslip(v[i], nullptr, x[j], big, j, xscoll, norm, vsnew);
1368 }
1369
1370 if (dimension == 2) vsnew[2] = 0.0;
1371
1372 // check on rescaling of vsnew
1373
1374 if (rescale_collide) {
1375 double vsq = vsnew[0] * vsnew[0] + vsnew[1] * vsnew[1] + vsnew[2] * vsnew[2];
1376 if (vsq > vmaxsq) {
1377 nrescale++;
1378 MathExtra::scale3(vmax / sqrt(vsq), vsnew);
1379 }
1380 }
1381
1382 // update BIG particle and WALL and SRD
1383 // BIG particle is not torqued if sphere and SLIP collision
1384
1385 if (collidestyle == SLIP && type == SPHERE)
1386 force_torque(v[i], vsnew, xscoll, xbcoll, f[j], nullptr);
1387 else if (type != WALL)
1388 force_torque(v[i], vsnew, xscoll, xbcoll, f[j], torque[j]);
1389 else if (type == WALL)
1390 force_wall(v[i], vsnew, j);
1391
1392 ibin = binsrd[i] = update_srd(i, t_remain, xscoll, vsnew, x[i], v[i]);
1393
1394 if (ibounce == 0) ncollide++;
1395 ibounce++;
1396 if (ibounce < maxbounceallow || maxbounceallow == 0) collide_flag = 1;
1397 dt = t_remain;
1398 break;
1399 }
1400 }
1401 }
1402
1403 nbounce += ibounce;
1404 if (maxbounceallow && ibounce >= maxbounceallow) bouncemaxnum++;
1405 if (ibounce > bouncemax) bouncemax = ibounce;
1406 }
1407 }
1408
1409 /* ----------------------------------------------------------------------
1410 detect all collisions between SRD and big particles
1411 an SRD can be inside more than one big particle at a time
1412 requires finding which big particle SRD collided with first
1413 unoverlap SRDs for each collision
1414 ------------------------------------------------------------------------- */
1415
collisions_multi()1416 void FixSRD::collisions_multi()
1417 {
1418 int i, j, k, m, type, nbig, ibin, ibounce, inside, jfirst, typefirst, jlast;
1419 double dt, t_remain, t_first;
1420 double norm[3], xscoll[3], xbcoll[3], vsnew[3];
1421 double normfirst[3], xscollfirst[3], xbcollfirst[3];
1422 Big *big;
1423
1424 // outer loop over SRD particles
1425 // inner loop over BIG particles or WALLS that overlap SRD particle bin
1426 // loop over all BIG and WALLS to find which one SRD collided with first
1427 // if overlap between SRD and BIG particle or wall:
1428 // compute collision pt in time
1429 // update x,v of SRD and f,torque on BIG particle
1430 // re-bin SRD particle after collision
1431 // iterate until the SRD particle has no overlaps with BIG particles or WALLS
1432
1433 double **x = atom->x;
1434 double **v = atom->v;
1435 double **f = atom->f;
1436 double **torque = atom->torque;
1437 tagint *tag = atom->tag;
1438 int *mask = atom->mask;
1439 int nlocal = atom->nlocal;
1440
1441 for (i = 0; i < nlocal; i++) {
1442 if (!(mask[i] & groupbit)) continue;
1443
1444 ibin = binsrd[i];
1445 if (nbinbig[ibin] == 0) continue;
1446
1447 ibounce = 0;
1448 jlast = -1;
1449 dt = dt_big;
1450
1451 while (1) {
1452 nbig = nbinbig[ibin];
1453 if (ibounce == 0) ncheck += nbig;
1454
1455 t_first = 0.0;
1456 for (m = 0; m < nbig; m++) {
1457 k = binbig[ibin][m];
1458 big = &biglist[k];
1459 j = big->index;
1460 if (j == jlast) continue;
1461 type = big->type;
1462
1463 if (type == SPHERE)
1464 inside = inside_sphere(x[i], x[j], big);
1465 else if (type == ELLIPSOID)
1466 inside = inside_ellipsoid(x[i], x[j], big);
1467 else if (type == LINE)
1468 inside = inside_line(x[i], x[j], v[i], v[j], big, dt);
1469 else if (type == TRIANGLE)
1470 inside = inside_tri(x[i], x[j], v[i], v[j], big, dt);
1471 else
1472 inside = inside_wall(x[i], j);
1473
1474 if (inside) {
1475 if (type == SPHERE)
1476 t_remain = collision_sphere_exact(x[i], x[j], v[i], v[j], big, xscoll, xbcoll, norm);
1477 else if (type == ELLIPSOID)
1478 t_remain = collision_ellipsoid_exact(x[i], x[j], v[i], v[j], big, xscoll, xbcoll, norm);
1479 else if (type == LINE)
1480 t_remain = collision_line_exact(x[i], x[j], v[i], v[j], big, dt, xscoll, xbcoll, norm);
1481 else if (type == TRIANGLE)
1482 t_remain = collision_tri_exact(x[i], x[j], v[i], v[j], big, dt, xscoll, xbcoll, norm);
1483 else
1484 t_remain = collision_wall_exact(x[i], j, v[i], xscoll, xbcoll, norm);
1485
1486 #ifdef SRD_DEBUG
1487 if (update->ntimestep == SRD_DEBUG_TIMESTEP && tag[i] == SRD_DEBUG_ATOMID)
1488 print_collision(i, j, ibounce, t_remain, dt, xscoll, xbcoll, norm, type);
1489 #endif
1490
1491 if (t_remain > dt || t_remain < 0.0) {
1492 ninside++;
1493 if (insideflag == INSIDE_ERROR || insideflag == INSIDE_WARN) {
1494 std::string mesg;
1495 if (type != WALL)
1496 mesg = fmt::format("SRD particle {} started inside big particle {} on step {} "
1497 "bounce {}",
1498 tag[i], tag[j], update->ntimestep, ibounce + 1);
1499 else
1500 mesg = fmt::format("SRD particle {} started inside wall {} on step {} "
1501 "bounce {}",
1502 tag[i], j, update->ntimestep, ibounce + 1);
1503
1504 if (insideflag == INSIDE_ERROR) error->one(FLERR, mesg);
1505 error->warning(FLERR, mesg);
1506 }
1507 t_first = 0.0;
1508 break;
1509 }
1510
1511 if (t_remain > t_first) {
1512 t_first = t_remain;
1513 jfirst = j;
1514 typefirst = type;
1515 xscollfirst[0] = xscoll[0];
1516 xscollfirst[1] = xscoll[1];
1517 xscollfirst[2] = xscoll[2];
1518 xbcollfirst[0] = xbcoll[0];
1519 xbcollfirst[1] = xbcoll[1];
1520 xbcollfirst[2] = xbcoll[2];
1521 normfirst[0] = norm[0];
1522 normfirst[1] = norm[1];
1523 normfirst[2] = norm[2];
1524 }
1525 }
1526 }
1527
1528 if (t_first == 0.0) break;
1529 j = jlast = jfirst;
1530 type = typefirst;
1531 xscoll[0] = xscollfirst[0];
1532 xscoll[1] = xscollfirst[1];
1533 xscoll[2] = xscollfirst[2];
1534 xbcoll[0] = xbcollfirst[0];
1535 xbcoll[1] = xbcollfirst[1];
1536 xbcoll[2] = xbcollfirst[2];
1537 norm[0] = normfirst[0];
1538 norm[1] = normfirst[1];
1539 norm[2] = normfirst[2];
1540
1541 if (collidestyle == SLIP) {
1542 if (type != WALL)
1543 slip(v[i], v[j], x[j], big, xscoll, norm, vsnew);
1544 else
1545 slip_wall(v[i], j, norm, vsnew);
1546 } else {
1547 if (type != WALL)
1548 noslip(v[i], v[j], x[j], big, -1, xscoll, norm, vsnew);
1549 else
1550 noslip(v[i], nullptr, x[j], big, j, xscoll, norm, vsnew);
1551 }
1552
1553 if (dimension == 2) vsnew[2] = 0.0;
1554
1555 // check on rescaling of vsnew
1556
1557 if (rescale_collide) {
1558 double vsq = vsnew[0] * vsnew[0] + vsnew[1] * vsnew[1] + vsnew[2] * vsnew[2];
1559 if (vsq > vmaxsq) {
1560 nrescale++;
1561 MathExtra::scale3(vmax / sqrt(vsq), vsnew);
1562 }
1563 }
1564
1565 // update BIG particle and WALL and SRD
1566 // BIG particle is not torqued if sphere and SLIP collision
1567
1568 if (collidestyle == SLIP && type == SPHERE)
1569 force_torque(v[i], vsnew, xscoll, xbcoll, f[j], nullptr);
1570 else if (type != WALL)
1571 force_torque(v[i], vsnew, xscoll, xbcoll, f[j], torque[j]);
1572 else if (type == WALL)
1573 force_wall(v[i], vsnew, j);
1574
1575 ibin = binsrd[i] = update_srd(i, t_first, xscoll, vsnew, x[i], v[i]);
1576
1577 if (ibounce == 0) ncollide++;
1578 ibounce++;
1579 if (ibounce == maxbounceallow) break;
1580 dt = t_first;
1581 }
1582
1583 nbounce += ibounce;
1584 if (maxbounceallow && ibounce >= maxbounceallow) bouncemaxnum++;
1585 if (ibounce > bouncemax) bouncemax = ibounce;
1586 }
1587 }
1588
1589 /* ----------------------------------------------------------------------
1590 check if SRD particle S is inside spherical big particle B
1591 ------------------------------------------------------------------------- */
1592
inside_sphere(double * xs,double * xb,Big * big)1593 int FixSRD::inside_sphere(double *xs, double *xb, Big *big)
1594 {
1595 double dx, dy, dz;
1596
1597 dx = xs[0] - xb[0];
1598 dy = xs[1] - xb[1];
1599 dz = xs[2] - xb[2];
1600
1601 if (dx * dx + dy * dy + dz * dz <= big->radsq) return 1;
1602 return 0;
1603 }
1604
1605 /* ----------------------------------------------------------------------
1606 check if SRD particle S is inside ellipsoidal big particle B
1607 ------------------------------------------------------------------------- */
1608
inside_ellipsoid(double * xs,double * xb,Big * big)1609 int FixSRD::inside_ellipsoid(double *xs, double *xb, Big *big)
1610 {
1611 double x, y, z;
1612
1613 double *ex = big->ex;
1614 double *ey = big->ey;
1615 double *ez = big->ez;
1616
1617 double xs_xb[3];
1618 xs_xb[0] = xs[0] - xb[0];
1619 xs_xb[1] = xs[1] - xb[1];
1620 xs_xb[2] = xs[2] - xb[2];
1621
1622 x = xs_xb[0] * ex[0] + xs_xb[1] * ex[1] + xs_xb[2] * ex[2];
1623 y = xs_xb[0] * ey[0] + xs_xb[1] * ey[1] + xs_xb[2] * ey[2];
1624 z = xs_xb[0] * ez[0] + xs_xb[1] * ez[1] + xs_xb[2] * ez[2];
1625
1626 if (x * x * big->aradsqinv + y * y * big->bradsqinv + z * z * big->cradsqinv <= 1.0) return 1;
1627 return 0;
1628 }
1629
1630 /* ----------------------------------------------------------------------
1631 check if SRD particle S is inside line big particle B
1632 collision only possible if:
1633 S starts on positive side of infinite line,
1634 which means it will collide with outside of rigid body made of lines
1635 since line segments have outward normals,
1636 when vector from first to last point is crossed with +z axis
1637 S ends on negative side of infinite line
1638 unlike most other inside() routines, then calculate exact collision:
1639 solve for collision pt along infinite line
1640 collision if pt is within endpoints of B
1641 ------------------------------------------------------------------------- */
1642
inside_line(double * xs,double * xb,double * vs,double * vb,Big * big,double dt_step)1643 int FixSRD::inside_line(double *xs, double *xb, double *vs, double *vb, Big *big, double dt_step)
1644 {
1645 double pmc0[2], pmc1[2], n0[2], n1[2];
1646
1647 // 1 and 2 = start and end of timestep
1648 // pmc = P - C, where P = position of S, C = position of B
1649 // n = normal to line = [-sin(theta),cos(theta)], theta = orientation of B
1650 // (P-C) dot N = side of line that S is on
1651 // side0 = -1,0,1 for which side of line B that S is on at start of step
1652 // side1 = -1,0,1 for which side of line B that S is on at end of step
1653
1654 xs1[0] = xs[0];
1655 xs1[1] = xs[1];
1656 xb1[0] = xb[0];
1657 xb1[1] = xb[1];
1658
1659 xs0[0] = xs1[0] - dt_step * vs[0];
1660 xs0[1] = xs1[1] - dt_step * vs[1];
1661 xb0[0] = xb1[0] - dt_step * vb[0];
1662 xb0[1] = xb1[1] - dt_step * vb[1];
1663
1664 theta1 = big->theta;
1665 theta0 = theta1 - dt_step * big->omega[2];
1666
1667 pmc0[0] = xs0[0] - xb0[0];
1668 pmc0[1] = xs0[1] - xb0[1];
1669 n0[0] = sin(theta0);
1670 n0[1] = -cos(theta0);
1671
1672 pmc1[0] = xs1[0] - xb1[0];
1673 pmc1[1] = xs1[1] - xb1[1];
1674 n1[0] = sin(theta1);
1675 n1[1] = -cos(theta1);
1676
1677 double side0 = pmc0[0] * n0[0] + pmc0[1] * n0[1];
1678 double side1 = pmc1[0] * n1[0] + pmc1[1] * n1[1];
1679
1680 if (side0 <= 0.0 || side1 >= 0.0) return 0;
1681
1682 // solve for time t (0 to 1) at which moving particle
1683 // crosses infinite moving/rotating line
1684
1685 // Newton-Raphson solve of full non-linear parametric equation
1686
1687 tfraction = newton_raphson(0.0, 1.0);
1688
1689 // quadratic equation solve of approximate parametric equation
1690
1691 /*
1692 n1_n0[0] = n1[0]-n0[0]; n1_n0[1] = n1[1]-n0[1];
1693 pmc1_pmc0[0] = pmc1[0]-pmc0[0]; pmc1_pmc0[1] = pmc1[1]-pmc0[1];
1694
1695 double a = pmc1_pmc0[0]*n1_n0[0] + pmc1_pmc0[1]*n1_n0[1];
1696 double b = pmc1_pmc0[0]*n0[0] + pmc1_pmc0[1]*n0[1] +
1697 n1_n0[0]*pmc0[0] + n1_n0[1]*pmc0[1];
1698 double c = pmc0[0]*n0[0] + pmc0[1]*n0[1];
1699
1700 if (a == 0.0) {
1701 double dot0 = pmc0[0]*n0[0] + pmc0[1]*n0[1];
1702 double dot1 = pmc1[0]*n0[0] + pmc1[1]*n0[1];
1703 double root = -dot0 / (dot1 - dot0);
1704 //printf("Linear root: %g %g\n",root,tfraction);
1705 tfraction = root;
1706
1707 } else {
1708
1709 double term = sqrt(b*b - 4.0*a*c);
1710 double root1 = (-b + term) / (2.0*a);
1711 double root2 = (-b - term) / (2.0*a);
1712
1713 //printf("ABC vecs: %g %g: %g %g\n",
1714 // pmc1_pmc0[0],pmc1_pmc0[1],n1_n0[0],n1_n0[1]);
1715 //printf("ABC vecs: %g %g: %g %g: %g %g %g\n",
1716 // n0[0],n0[1],n1[0],n1[1],theta0,theta1,big->omega[2]);
1717 //printf("ABC root: %g %g %g: %g %g %g\n",a,b,c,root1,root2,tfraction);
1718
1719 if (0.0 <= root1 && root1 <= 1.0) tfraction = root1;
1720 else if (0.0 <= root2 && root2 <= 1.0) tfraction = root2;
1721 else error->one(FLERR,"Bad quadratic solve for particle/line collision");
1722 }
1723 */
1724
1725 // check if collision pt is within line segment at collision time
1726
1727 xsc[0] = xs0[0] + tfraction * (xs1[0] - xs0[0]);
1728 xsc[1] = xs0[1] + tfraction * (xs1[1] - xs0[1]);
1729 xbc[0] = xb0[0] + tfraction * (xb1[0] - xb0[0]);
1730 xbc[1] = xb0[1] + tfraction * (xb1[1] - xb0[1]);
1731 double delx = xsc[0] - xbc[0];
1732 double dely = xsc[1] - xbc[1];
1733 double rsq = delx * delx + dely * dely;
1734 if (rsq > 0.25 * big->length * big->length) return 0;
1735
1736 //nbc[0] = n0[0] + tfraction*(n1[0]-n0[0]);
1737 //nbc[1] = n0[1] + tfraction*(n1[1]-n0[1]);
1738
1739 nbc[0] = sin(theta0 + tfraction * (theta1 - theta0));
1740 nbc[1] = -cos(theta0 + tfraction * (theta1 - theta0));
1741
1742 return 1;
1743 }
1744
1745 /* ----------------------------------------------------------------------
1746 check if SRD particle S is inside triangle big particle B
1747 collision only possible if:
1748 S starts on positive side of triangle plane,
1749 which means it will collide with outside of rigid body made of tris
1750 since triangles have outward normals,
1751 S ends on negative side of triangle plane
1752 unlike most other inside() routines, then calculate exact collision:
1753 solve for collision pt on triangle plane
1754 collision if pt is inside triangle B
1755 ------------------------------------------------------------------------- */
1756
inside_tri(double * xs,double * xb,double * vs,double * vb,Big * big,double dt_step)1757 int FixSRD::inside_tri(double *xs, double *xb, double *vs, double *vb, Big *big, double dt_step)
1758 {
1759 double pmc0[3], pmc1[3], n0[3];
1760 double n1_n0[3], pmc1_pmc0[3];
1761 double excoll[3], eycoll[3], ezcoll[3];
1762 double dc1[3], dc2[3], dc3[3];
1763 double c1[3], c2[3], c3[3];
1764 double c2mc1[3], c3mc2[3], c1mc3[3];
1765 double pvec[3], xproduct[3];
1766
1767 // 1 and 2 = start and end of timestep
1768 // pmc = P - C, where P = position of S, C = position of B
1769 // n = normal to triangle
1770 // (P-C) dot N = side of tri that S is on
1771 // side0 = -1,0,1 for which side of tri B that S is on at start of step
1772 // side1 = -1,0,1 for which side of tri B that S is on at end of step
1773
1774 double *omega = big->omega;
1775 double *n1 = big->norm;
1776
1777 n0[0] = n1[0] - dt_step * (omega[1] * n1[2] - omega[2] * n1[1]);
1778 n0[1] = n1[1] - dt_step * (omega[2] * n1[0] - omega[0] * n1[2]);
1779 n0[2] = n1[2] - dt_step * (omega[0] * n1[1] - omega[1] * n1[0]);
1780
1781 pmc0[0] = xs[0] - dt_step * vs[0] - xb[0] + dt_step * vb[0];
1782 pmc0[1] = xs[1] - dt_step * vs[1] - xb[1] + dt_step * vb[1];
1783 pmc0[2] = xs[2] - dt_step * vs[2] - xb[2] + dt_step * vb[2];
1784 pmc1[0] = xs[0] - xb[0];
1785 pmc1[1] = xs[1] - xb[1];
1786 pmc1[2] = xs[2] - xb[2];
1787
1788 double side0 = MathExtra::dot3(pmc0, n0);
1789 double side1 = MathExtra::dot3(pmc1, n1);
1790
1791 if (side0 <= 0.0 || side1 >= 0.0) return 0;
1792
1793 // solve for time t (0 to 1) at which moving particle
1794 // crosses moving/rotating tri
1795 // quadratic equation solve of approximate parametric equation
1796
1797 n1_n0[0] = n1[0] - n0[0];
1798 n1_n0[1] = n1[1] - n0[1];
1799 n1_n0[2] = n1[2] - n0[2];
1800 pmc1_pmc0[0] = pmc1[0] - pmc0[0];
1801 pmc1_pmc0[1] = pmc1[1] - pmc0[1];
1802 pmc1_pmc0[2] = pmc1[2] - pmc0[2];
1803
1804 double a = MathExtra::dot3(pmc1_pmc0, n1_n0);
1805 double b = MathExtra::dot3(pmc1_pmc0, n0) + MathExtra::dot3(n1_n0, pmc0);
1806 double c = MathExtra::dot3(pmc0, n0);
1807
1808 if (a == 0.0) {
1809 double dot0 = MathExtra::dot3(pmc0, n0);
1810 double dot1 = MathExtra::dot3(pmc1, n0);
1811 double root = -dot0 / (dot1 - dot0);
1812 tfraction = root;
1813 } else {
1814 double term = sqrt(b * b - 4.0 * a * c);
1815 double root1 = (-b + term) / (2.0 * a);
1816 double root2 = (-b - term) / (2.0 * a);
1817 if (0.0 <= root1 && root1 <= 1.0)
1818 tfraction = root1;
1819 else if (0.0 <= root2 && root2 <= 1.0)
1820 tfraction = root2;
1821 else
1822 error->one(FLERR, "Bad quadratic solve for particle/tri collision");
1823 }
1824
1825 // calculate position/orientation of S and B at collision time
1826 // dt = time previous to now at which collision occurs
1827 // point = S position in plane of triangle at collision time
1828 // Excoll,Eycoll,Ezcoll = orientation of tri at collision time
1829 // c1,c2,c3 = corner points of triangle at collision time
1830 // nbc = normal to plane of triangle at collision time
1831
1832 AtomVecTri::Bonus *tbonus;
1833 tbonus = avec_tri->bonus;
1834
1835 double *ex = big->ex;
1836 double *ey = big->ey;
1837 double *ez = big->ez;
1838 int index = atom->tri[big->index];
1839 double *c1body = tbonus[index].c1;
1840 double *c2body = tbonus[index].c2;
1841 double *c3body = tbonus[index].c3;
1842
1843 double dt = (1.0 - tfraction) * dt_step;
1844
1845 xsc[0] = xs[0] - dt * vs[0];
1846 xsc[1] = xs[1] - dt * vs[1];
1847 xsc[2] = xs[2] - dt * vs[2];
1848 xbc[0] = xb[0] - dt * vb[0];
1849 xbc[1] = xb[1] - dt * vb[1];
1850 xbc[2] = xb[2] - dt * vb[2];
1851
1852 excoll[0] = ex[0] - dt * (omega[1] * ex[2] - omega[2] * ex[1]);
1853 excoll[1] = ex[1] - dt * (omega[2] * ex[0] - omega[0] * ex[2]);
1854 excoll[2] = ex[2] - dt * (omega[0] * ex[1] - omega[1] * ex[0]);
1855
1856 eycoll[0] = ey[0] - dt * (omega[1] * ey[2] - omega[2] * ey[1]);
1857 eycoll[1] = ey[1] - dt * (omega[2] * ey[0] - omega[0] * ey[2]);
1858 eycoll[2] = ey[2] - dt * (omega[0] * ey[1] - omega[1] * ey[0]);
1859
1860 ezcoll[0] = ez[0] - dt * (omega[1] * ez[2] - omega[2] * ez[1]);
1861 ezcoll[1] = ez[1] - dt * (omega[2] * ez[0] - omega[0] * ez[2]);
1862 ezcoll[2] = ez[2] - dt * (omega[0] * ez[1] - omega[1] * ez[0]);
1863
1864 MathExtra::matvec(excoll, eycoll, ezcoll, c1body, dc1);
1865 MathExtra::matvec(excoll, eycoll, ezcoll, c2body, dc2);
1866 MathExtra::matvec(excoll, eycoll, ezcoll, c3body, dc3);
1867
1868 MathExtra::add3(xbc, dc1, c1);
1869 MathExtra::add3(xbc, dc2, c2);
1870 MathExtra::add3(xbc, dc3, c3);
1871
1872 MathExtra::sub3(c2, c1, c2mc1);
1873 MathExtra::sub3(c3, c2, c3mc2);
1874 MathExtra::sub3(c1, c3, c1mc3);
1875
1876 MathExtra::cross3(c2mc1, c3mc2, nbc);
1877 MathExtra::norm3(nbc);
1878
1879 // check if collision pt is within triangle
1880 // pvec = vector from tri vertex to intersection point
1881 // xproduct = cross product of edge vec with pvec
1882 // if dot product of xproduct with nbc < 0.0 for any of 3 edges,
1883 // intersection point is outside tri
1884
1885 MathExtra::sub3(xsc, c1, pvec);
1886 MathExtra::cross3(c2mc1, pvec, xproduct);
1887 if (MathExtra::dot3(xproduct, nbc) < 0.0) return 0;
1888
1889 MathExtra::sub3(xsc, c2, pvec);
1890 MathExtra::cross3(c3mc2, pvec, xproduct);
1891 if (MathExtra::dot3(xproduct, nbc) < 0.0) return 0;
1892
1893 MathExtra::sub3(xsc, c3, pvec);
1894 MathExtra::cross3(c1mc3, pvec, xproduct);
1895 if (MathExtra::dot3(xproduct, nbc) < 0.0) return 0;
1896
1897 return 1;
1898 }
1899
1900 /* ----------------------------------------------------------------------
1901 check if SRD particle S is inside wall IWALL
1902 ------------------------------------------------------------------------- */
1903
inside_wall(double * xs,int iwall)1904 int FixSRD::inside_wall(double *xs, int iwall)
1905 {
1906 int dim = wallwhich[iwall] / 2;
1907 int side = wallwhich[iwall] % 2;
1908
1909 if (side == 0 && xs[dim] < xwall[iwall]) return 1;
1910 if (side && xs[dim] > xwall[iwall]) return 1;
1911 return 0;
1912 }
1913
1914 /* ----------------------------------------------------------------------
1915 collision of SRD particle S with surface of spherical big particle B
1916 exact because compute time of collision
1917 dt = time previous to now at which collision occurs
1918 xscoll = collision pt = position of SRD at time of collision
1919 xbcoll = position of big particle at time of collision
1920 norm = surface normal of collision pt at time of collision
1921 ------------------------------------------------------------------------- */
1922
collision_sphere_exact(double * xs,double * xb,double * vs,double * vb,Big * big,double * xscoll,double * xbcoll,double * norm)1923 double FixSRD::collision_sphere_exact(double *xs, double *xb, double *vs, double *vb, Big *big,
1924 double *xscoll, double *xbcoll, double *norm)
1925 {
1926 double vs_dot_vs, vb_dot_vb, vs_dot_vb;
1927 double vs_dot_xb, vb_dot_xs, vs_dot_xs, vb_dot_xb;
1928 double xs_dot_xs, xb_dot_xb, xs_dot_xb;
1929 double a, b, c, scale;
1930
1931 vs_dot_vs = vs[0] * vs[0] + vs[1] * vs[1] + vs[2] * vs[2];
1932 vb_dot_vb = vb[0] * vb[0] + vb[1] * vb[1] + vb[2] * vb[2];
1933 vs_dot_vb = vs[0] * vb[0] + vs[1] * vb[1] + vs[2] * vb[2];
1934
1935 vs_dot_xb = vs[0] * xb[0] + vs[1] * xb[1] + vs[2] * xb[2];
1936 vb_dot_xs = vb[0] * xs[0] + vb[1] * xs[1] + vb[2] * xs[2];
1937 vs_dot_xs = vs[0] * xs[0] + vs[1] * xs[1] + vs[2] * xs[2];
1938 vb_dot_xb = vb[0] * xb[0] + vb[1] * xb[1] + vb[2] * xb[2];
1939
1940 xs_dot_xs = xs[0] * xs[0] + xs[1] * xs[1] + xs[2] * xs[2];
1941 xb_dot_xb = xb[0] * xb[0] + xb[1] * xb[1] + xb[2] * xb[2];
1942 xs_dot_xb = xs[0] * xb[0] + xs[1] * xb[1] + xs[2] * xb[2];
1943
1944 a = vs_dot_vs + vb_dot_vb - 2.0 * vs_dot_vb;
1945 b = 2.0 * (vs_dot_xb + vb_dot_xs - vs_dot_xs - vb_dot_xb);
1946 c = xs_dot_xs + xb_dot_xb - 2.0 * xs_dot_xb - big->radsq;
1947
1948 double dt = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
1949
1950 xscoll[0] = xs[0] - dt * vs[0];
1951 xscoll[1] = xs[1] - dt * vs[1];
1952 xscoll[2] = xs[2] - dt * vs[2];
1953
1954 xbcoll[0] = xb[0] - dt * vb[0];
1955 xbcoll[1] = xb[1] - dt * vb[1];
1956 xbcoll[2] = xb[2] - dt * vb[2];
1957
1958 norm[0] = xscoll[0] - xbcoll[0];
1959 norm[1] = xscoll[1] - xbcoll[1];
1960 norm[2] = xscoll[2] - xbcoll[2];
1961 scale = 1.0 / sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
1962 norm[0] *= scale;
1963 norm[1] *= scale;
1964 norm[2] *= scale;
1965
1966 return dt;
1967 }
1968
1969 /* ----------------------------------------------------------------------
1970 collision of SRD particle S with surface of spherical big particle B
1971 inexact because just push SRD to surface of big particle at end of step
1972 time of collision = end of step
1973 xscoll = collision pt = position of SRD at time of collision
1974 xbcoll = xb = position of big particle at time of collision
1975 norm = surface normal of collision pt at time of collision
1976 ------------------------------------------------------------------------- */
1977
collision_sphere_inexact(double * xs,double * xb,Big * big,double * xscoll,double * xbcoll,double * norm)1978 void FixSRD::collision_sphere_inexact(double *xs, double *xb, Big *big, double *xscoll,
1979 double *xbcoll, double *norm)
1980 {
1981 double scale;
1982
1983 norm[0] = xs[0] - xb[0];
1984 norm[1] = xs[1] - xb[1];
1985 norm[2] = xs[2] - xb[2];
1986 scale = 1.0 / sqrt(norm[0] * norm[0] + norm[1] * norm[1] + norm[2] * norm[2]);
1987 norm[0] *= scale;
1988 norm[1] *= scale;
1989 norm[2] *= scale;
1990
1991 xscoll[0] = xb[0] + big->radius * norm[0];
1992 xscoll[1] = xb[1] + big->radius * norm[1];
1993 xscoll[2] = xb[2] + big->radius * norm[2];
1994
1995 xbcoll[0] = xb[0];
1996 xbcoll[1] = xb[1];
1997 xbcoll[2] = xb[2];
1998 }
1999
2000 /* ----------------------------------------------------------------------
2001 collision of SRD particle S with surface of ellipsoidal big particle B
2002 exact because compute time of collision
2003 dt = time previous to now at which collision occurs
2004 xscoll = collision pt = position of SRD at time of collision
2005 xbcoll = position of big particle at time of collision
2006 norm = surface normal of collision pt at time of collision
2007 ------------------------------------------------------------------------- */
2008
collision_ellipsoid_exact(double * xs,double * xb,double * vs,double * vb,Big * big,double * xscoll,double * xbcoll,double * norm)2009 double FixSRD::collision_ellipsoid_exact(double *xs, double *xb, double *vs, double *vb, Big *big,
2010 double *xscoll, double *xbcoll, double *norm)
2011 {
2012 double vs_vb[3], xs_xb[3], omega_ex[3], omega_ey[3], omega_ez[3];
2013 double excoll[3], eycoll[3], ezcoll[3], delta[3], xbody[3], nbody[3];
2014 double ax, bx, cx, ay, by, cy, az, bz, cz;
2015 double a, b, c;
2016
2017 double *omega = big->omega;
2018 double *ex = big->ex;
2019 double *ey = big->ey;
2020 double *ez = big->ez;
2021
2022 vs_vb[0] = vs[0] - vb[0];
2023 vs_vb[1] = vs[1] - vb[1];
2024 vs_vb[2] = vs[2] - vb[2];
2025 xs_xb[0] = xs[0] - xb[0];
2026 xs_xb[1] = xs[1] - xb[1];
2027 xs_xb[2] = xs[2] - xb[2];
2028
2029 MathExtra::cross3(omega, ex, omega_ex);
2030 MathExtra::cross3(omega, ey, omega_ey);
2031 MathExtra::cross3(omega, ez, omega_ez);
2032
2033 ax = vs_vb[0] * omega_ex[0] + vs_vb[1] * omega_ex[1] + vs_vb[2] * omega_ex[2];
2034 bx = -(vs_vb[0] * ex[0] + vs_vb[1] * ex[1] + vs_vb[2] * ex[2]);
2035 bx -= xs_xb[0] * omega_ex[0] + xs_xb[1] * omega_ex[1] + xs_xb[2] * omega_ex[2];
2036 cx = xs_xb[0] * ex[0] + xs_xb[1] * ex[1] + xs_xb[2] * ex[2];
2037
2038 ay = vs_vb[0] * omega_ey[0] + vs_vb[1] * omega_ey[1] + vs_vb[2] * omega_ey[2];
2039 by = -(vs_vb[0] * ey[0] + vs_vb[1] * ey[1] + vs_vb[2] * ey[2]);
2040 by -= xs_xb[0] * omega_ey[0] + xs_xb[1] * omega_ey[1] + xs_xb[2] * omega_ey[2];
2041 cy = xs_xb[0] * ey[0] + xs_xb[1] * ey[1] + xs_xb[2] * ey[2];
2042
2043 az = vs_vb[0] * omega_ez[0] + vs_vb[1] * omega_ez[1] + vs_vb[2] * omega_ez[2];
2044 bz = -(vs_vb[0] * ez[0] + vs_vb[1] * ez[1] + vs_vb[2] * ez[2]);
2045 bz -= xs_xb[0] * omega_ez[0] + xs_xb[1] * omega_ez[1] + xs_xb[2] * omega_ez[2];
2046 cz = xs_xb[0] * ez[0] + xs_xb[1] * ez[1] + xs_xb[2] * ez[2];
2047
2048 a = (bx * bx + 2.0 * ax * cx) * big->aradsqinv + (by * by + 2.0 * ay * cy) * big->bradsqinv +
2049 (bz * bz + 2.0 * az * cz) * big->cradsqinv;
2050 b = 2.0 * (bx * cx * big->aradsqinv + by * cy * big->bradsqinv + bz * cz * big->cradsqinv);
2051 c = cx * cx * big->aradsqinv + cy * cy * big->bradsqinv + cz * cz * big->cradsqinv - 1.0;
2052
2053 double dt = (-b + sqrt(b * b - 4.0 * a * c)) / (2.0 * a);
2054
2055 xscoll[0] = xs[0] - dt * vs[0];
2056 xscoll[1] = xs[1] - dt * vs[1];
2057 xscoll[2] = xs[2] - dt * vs[2];
2058
2059 xbcoll[0] = xb[0] - dt * vb[0];
2060 xbcoll[1] = xb[1] - dt * vb[1];
2061 xbcoll[2] = xb[2] - dt * vb[2];
2062
2063 // calculate normal to ellipsoid at collision pt
2064 // Excoll,Eycoll,Ezcoll = orientation of ellipsoid at collision time
2065 // nbody = normal in body frame of ellipsoid (Excoll,Eycoll,Ezcoll)
2066 // norm = normal in space frame
2067 // only worry about normalizing final norm vector
2068
2069 excoll[0] = ex[0] - dt * (omega[1] * ex[2] - omega[2] * ex[1]);
2070 excoll[1] = ex[1] - dt * (omega[2] * ex[0] - omega[0] * ex[2]);
2071 excoll[2] = ex[2] - dt * (omega[0] * ex[1] - omega[1] * ex[0]);
2072
2073 eycoll[0] = ey[0] - dt * (omega[1] * ey[2] - omega[2] * ey[1]);
2074 eycoll[1] = ey[1] - dt * (omega[2] * ey[0] - omega[0] * ey[2]);
2075 eycoll[2] = ey[2] - dt * (omega[0] * ey[1] - omega[1] * ey[0]);
2076
2077 ezcoll[0] = ez[0] - dt * (omega[1] * ez[2] - omega[2] * ez[1]);
2078 ezcoll[1] = ez[1] - dt * (omega[2] * ez[0] - omega[0] * ez[2]);
2079 ezcoll[2] = ez[2] - dt * (omega[0] * ez[1] - omega[1] * ez[0]);
2080
2081 MathExtra::sub3(xscoll, xbcoll, delta);
2082 MathExtra::transpose_matvec(excoll, eycoll, ezcoll, delta, xbody);
2083
2084 nbody[0] = xbody[0] * big->aradsqinv;
2085 nbody[1] = xbody[1] * big->bradsqinv;
2086 nbody[2] = xbody[2] * big->cradsqinv;
2087
2088 MathExtra::matvec(excoll, eycoll, ezcoll, nbody, norm);
2089 MathExtra::norm3(norm);
2090
2091 return dt;
2092 }
2093
2094 /* ----------------------------------------------------------------------
2095 collision of SRD particle S with surface of ellipsoidal big particle B
2096 inexact because just push SRD to surface of big particle at end of step
2097 time of collision = end of step
2098 xscoll = collision pt = position of SRD at time of collision
2099 xbcoll = xb = position of big particle at time of collision
2100 norm = surface normal of collision pt at time of collision
2101 ------------------------------------------------------------------------- */
2102
collision_ellipsoid_inexact(double * xs,double * xb,Big * big,double * xscoll,double * xbcoll,double * norm)2103 void FixSRD::collision_ellipsoid_inexact(double *xs, double *xb, Big *big, double *xscoll,
2104 double *xbcoll, double *norm)
2105 {
2106 double xs_xb[3], delta[3], xbody[3], nbody[3];
2107
2108 double *ex = big->ex;
2109 double *ey = big->ey;
2110 double *ez = big->ez;
2111
2112 MathExtra::sub3(xs, xb, xs_xb);
2113 double x = MathExtra::dot3(xs_xb, ex);
2114 double y = MathExtra::dot3(xs_xb, ey);
2115 double z = MathExtra::dot3(xs_xb, ez);
2116
2117 double scale =
2118 1.0 / sqrt(x * x * big->aradsqinv + y * y * big->bradsqinv + z * z * big->cradsqinv);
2119 x *= scale;
2120 y *= scale;
2121 z *= scale;
2122
2123 xscoll[0] = x * ex[0] + y * ey[0] + z * ez[0] + xb[0];
2124 xscoll[1] = x * ex[1] + y * ey[1] + z * ez[1] + xb[1];
2125 xscoll[2] = x * ex[2] + y * ey[2] + z * ez[2] + xb[2];
2126
2127 xbcoll[0] = xb[0];
2128 xbcoll[1] = xb[1];
2129 xbcoll[2] = xb[2];
2130
2131 // calculate normal to ellipsoid at collision pt
2132 // nbody = normal in body frame of ellipsoid
2133 // norm = normal in space frame
2134 // only worry about normalizing final norm vector
2135
2136 MathExtra::sub3(xscoll, xbcoll, delta);
2137 MathExtra::transpose_matvec(ex, ey, ez, delta, xbody);
2138
2139 nbody[0] = xbody[0] * big->aradsqinv;
2140 nbody[1] = xbody[1] * big->bradsqinv;
2141 nbody[2] = xbody[2] * big->cradsqinv;
2142
2143 MathExtra::matvec(ex, ey, ez, nbody, norm);
2144 MathExtra::norm3(norm);
2145 }
2146
2147 /* ----------------------------------------------------------------------
2148 collision of SRD particle S with surface of line big particle B
2149 exact because compute time of collision
2150 dt = time previous to now at which collision occurs
2151 xscoll = collision pt = position of SRD at time of collision
2152 xbcoll = position of big particle at time of collision
2153 norm = surface normal of collision pt at time of collision
2154 ------------------------------------------------------------------------- */
2155
collision_line_exact(double *,double *,double *,double *,Big *,double dt_step,double * xscoll,double * xbcoll,double * norm)2156 double FixSRD::collision_line_exact(double * /*xs*/, double * /*xb*/, double * /*vs*/,
2157 double * /*vb*/, Big * /*big*/, double dt_step, double *xscoll,
2158 double *xbcoll, double *norm)
2159 {
2160 xscoll[0] = xsc[0];
2161 xscoll[1] = xsc[1];
2162 xscoll[2] = 0.0;
2163 xbcoll[0] = xbc[0];
2164 xbcoll[1] = xbc[1];
2165 xbcoll[2] = 0.0;
2166
2167 norm[0] = nbc[0];
2168 norm[1] = nbc[1];
2169 norm[2] = 0.0;
2170
2171 return (1.0 - tfraction) * dt_step;
2172 }
2173
2174 /* ----------------------------------------------------------------------
2175 collision of SRD particle S with surface of triangle big particle B
2176 exact because compute time of collision
2177 dt = time previous to now at which collision occurs
2178 xscoll = collision pt = position of SRD at time of collision
2179 xbcoll = position of big particle at time of collision
2180 norm = surface normal of collision pt at time of collision
2181 ------------------------------------------------------------------------- */
2182
collision_tri_exact(double *,double *,double *,double *,Big *,double dt_step,double * xscoll,double * xbcoll,double * norm)2183 double FixSRD::collision_tri_exact(double * /*xs*/, double * /*xb*/, double * /*vs*/,
2184 double * /*vb*/, Big * /*big*/, double dt_step, double *xscoll,
2185 double *xbcoll, double *norm)
2186 {
2187 xscoll[0] = xsc[0];
2188 xscoll[1] = xsc[1];
2189 xscoll[2] = xsc[2];
2190 xbcoll[0] = xbc[0];
2191 xbcoll[1] = xbc[1];
2192 xbcoll[2] = xbc[2];
2193
2194 norm[0] = nbc[0];
2195 norm[1] = nbc[1];
2196 norm[2] = nbc[2];
2197
2198 return (1.0 - tfraction) * dt_step;
2199 }
2200
2201 /* ----------------------------------------------------------------------
2202 collision of SRD particle S with wall IWALL
2203 exact because compute time of collision
2204 dt = time previous to now at which collision occurs
2205 xscoll = collision pt = position of SRD at time of collision
2206 xbcoll = position of wall at time of collision
2207 norm = surface normal of collision pt at time of collision
2208 ------------------------------------------------------------------------- */
2209
collision_wall_exact(double * xs,int iwall,double * vs,double * xscoll,double * xbcoll,double * norm)2210 double FixSRD::collision_wall_exact(double *xs, int iwall, double *vs, double *xscoll,
2211 double *xbcoll, double *norm)
2212 {
2213 int dim = wallwhich[iwall] / 2;
2214
2215 double dt = (xs[dim] - xwall[iwall]) / (vs[dim] - vwall[iwall]);
2216 xscoll[0] = xs[0] - dt * vs[0];
2217 xscoll[1] = xs[1] - dt * vs[1];
2218 xscoll[2] = xs[2] - dt * vs[2];
2219
2220 xbcoll[0] = xbcoll[1] = xbcoll[2] = 0.0;
2221 xbcoll[dim] = xwall[iwall] - dt * vwall[iwall];
2222
2223 int side = wallwhich[iwall] % 2;
2224 norm[0] = norm[1] = norm[2] = 0.0;
2225 if (side == 0)
2226 norm[dim] = 1.0;
2227 else
2228 norm[dim] = -1.0;
2229
2230 return dt;
2231 }
2232
2233 /* ----------------------------------------------------------------------
2234 collision of SRD particle S with wall IWALL
2235 inexact because just push SRD to surface of wall at end of step
2236 time of collision = end of step
2237 xscoll = collision pt = position of SRD at time of collision
2238 xbcoll = position of wall at time of collision
2239 norm = surface normal of collision pt at time of collision
2240 ------------------------------------------------------------------------- */
2241
collision_wall_inexact(double * xs,int iwall,double * xscoll,double * xbcoll,double * norm)2242 void FixSRD::collision_wall_inexact(double *xs, int iwall, double *xscoll, double *xbcoll,
2243 double *norm)
2244 {
2245 int dim = wallwhich[iwall] / 2;
2246
2247 xscoll[0] = xs[0];
2248 xscoll[1] = xs[1];
2249 xscoll[2] = xs[2];
2250 xscoll[dim] = xwall[iwall];
2251
2252 xbcoll[0] = xbcoll[1] = xbcoll[2] = 0.0;
2253 xbcoll[dim] = xwall[iwall];
2254
2255 int side = wallwhich[iwall] % 2;
2256 norm[0] = norm[1] = norm[2] = 0.0;
2257 if (side == 0)
2258 norm[dim] = 1.0;
2259 else
2260 norm[dim] = -1.0;
2261 }
2262
2263 /* ----------------------------------------------------------------------
2264 SLIP collision with BIG particle with omega
2265 vs = velocity of SRD, vb = velocity of BIG
2266 xb = position of BIG, omega = rotation of BIG
2267 xsurf = collision pt on surf of BIG
2268 norm = unit normal from surface of BIG at collision pt
2269 v of BIG particle in direction of surf normal is added to v of SRD
2270 includes component due to rotation of BIG
2271 return vsnew of SRD
2272 ------------------------------------------------------------------------- */
2273
slip(double * vs,double * vb,double * xb,Big * big,double * xsurf,double * norm,double * vsnew)2274 void FixSRD::slip(double *vs, double *vb, double *xb, Big *big, double *xsurf, double *norm,
2275 double *vsnew)
2276 {
2277 double r1, r2, vnmag, vs_dot_n, vsurf_dot_n;
2278 double tangent[3], vsurf[3];
2279 double *omega = big->omega;
2280
2281 while (1) {
2282 r1 = sigma * random->gaussian();
2283 r2 = sigma * random->gaussian();
2284 vnmag = sqrt(r1 * r1 + r2 * r2);
2285 if (vnmag * vnmag <= vmaxsq) break;
2286 }
2287
2288 vs_dot_n = vs[0] * norm[0] + vs[1] * norm[1] + vs[2] * norm[2];
2289
2290 tangent[0] = vs[0] - vs_dot_n * norm[0];
2291 tangent[1] = vs[1] - vs_dot_n * norm[1];
2292 tangent[2] = vs[2] - vs_dot_n * norm[2];
2293
2294 // vsurf = velocity of collision pt = translation/rotation of BIG particle
2295 // NOTE: for sphere could just use vsurf = vb, since w x (xsurf-xb)
2296 // is orthogonal to norm and thus doesn't contribute to vsurf_dot_n
2297
2298 vsurf[0] = vb[0] + omega[1] * (xsurf[2] - xb[2]) - omega[2] * (xsurf[1] - xb[1]);
2299 vsurf[1] = vb[1] + omega[2] * (xsurf[0] - xb[0]) - omega[0] * (xsurf[2] - xb[2]);
2300 vsurf[2] = vb[2] + omega[0] * (xsurf[1] - xb[1]) - omega[1] * (xsurf[0] - xb[0]);
2301
2302 vsurf_dot_n = vsurf[0] * norm[0] + vsurf[1] * norm[1] + vsurf[2] * norm[2];
2303
2304 vsnew[0] = (vnmag + vsurf_dot_n) * norm[0] + tangent[0];
2305 vsnew[1] = (vnmag + vsurf_dot_n) * norm[1] + tangent[1];
2306 vsnew[2] = (vnmag + vsurf_dot_n) * norm[2] + tangent[2];
2307 }
2308
2309 /* ----------------------------------------------------------------------
2310 SLIP collision with wall IWALL
2311 vs = velocity of SRD
2312 norm = unit normal from WALL at collision pt
2313 v of WALL in direction of surf normal is added to v of SRD
2314 return vsnew of SRD
2315 ------------------------------------------------------------------------- */
2316
slip_wall(double * vs,int iwall,double * norm,double * vsnew)2317 void FixSRD::slip_wall(double *vs, int iwall, double *norm, double *vsnew)
2318 {
2319 double vs_dot_n, scale, r1, r2, vnmag, vtmag1, vtmag2;
2320 double tangent1[3], tangent2[3];
2321
2322 vs_dot_n = vs[0] * norm[0] + vs[1] * norm[1] + vs[2] * norm[2];
2323
2324 tangent1[0] = vs[0] - vs_dot_n * norm[0];
2325 tangent1[1] = vs[1] - vs_dot_n * norm[1];
2326 tangent1[2] = vs[2] - vs_dot_n * norm[2];
2327 scale =
2328 1.0 / sqrt(tangent1[0] * tangent1[0] + tangent1[1] * tangent1[1] + tangent1[2] * tangent1[2]);
2329 tangent1[0] *= scale;
2330 tangent1[1] *= scale;
2331 tangent1[2] *= scale;
2332
2333 tangent2[0] = norm[1] * tangent1[2] - norm[2] * tangent1[1];
2334 tangent2[1] = norm[2] * tangent1[0] - norm[0] * tangent1[2];
2335 tangent2[2] = norm[0] * tangent1[1] - norm[1] * tangent1[0];
2336
2337 while (1) {
2338 r1 = sigma * random->gaussian();
2339 r2 = sigma * random->gaussian();
2340 vnmag = sqrt(r1 * r1 + r2 * r2);
2341 vtmag1 = sigma * random->gaussian();
2342 vtmag2 = sigma * random->gaussian();
2343 if (vnmag * vnmag + vtmag1 * vtmag1 + vtmag2 * vtmag2 <= vmaxsq) break;
2344 }
2345
2346 vsnew[0] = vnmag * norm[0] + vtmag1 * tangent1[0] + vtmag2 * tangent2[0];
2347 vsnew[1] = vnmag * norm[1] + vtmag1 * tangent1[1] + vtmag2 * tangent2[1];
2348 vsnew[2] = vnmag * norm[2] + vtmag1 * tangent1[2] + vtmag2 * tangent2[2];
2349
2350 // add in velocity of collision pt = velocity of wall
2351
2352 int dim = wallwhich[iwall] / 2;
2353 vsnew[dim] += vwall[iwall];
2354 }
2355
2356 /* ----------------------------------------------------------------------
2357 NO-SLIP collision with BIG particle including WALL
2358 vs = velocity of SRD, vb = velocity of BIG
2359 xb = position of BIG, omega = rotation of BIG
2360 xsurf = collision pt on surf of BIG
2361 norm = unit normal from surface of BIG at collision pt
2362 v of collision pt is added to v of SRD
2363 includes component due to rotation of BIG
2364 return vsnew of SRD
2365 ------------------------------------------------------------------------- */
2366
noslip(double * vs,double * vb,double * xb,Big * big,int iwall,double * xsurf,double * norm,double * vsnew)2367 void FixSRD::noslip(double *vs, double *vb, double *xb, Big *big, int iwall, double *xsurf,
2368 double *norm, double *vsnew)
2369 {
2370 double vs_dot_n, scale, r1, r2, vnmag, vtmag1, vtmag2;
2371 double tangent1[3], tangent2[3];
2372
2373 vs_dot_n = vs[0] * norm[0] + vs[1] * norm[1] + vs[2] * norm[2];
2374
2375 tangent1[0] = vs[0] - vs_dot_n * norm[0];
2376 tangent1[1] = vs[1] - vs_dot_n * norm[1];
2377 tangent1[2] = vs[2] - vs_dot_n * norm[2];
2378 scale =
2379 1.0 / sqrt(tangent1[0] * tangent1[0] + tangent1[1] * tangent1[1] + tangent1[2] * tangent1[2]);
2380 tangent1[0] *= scale;
2381 tangent1[1] *= scale;
2382 tangent1[2] *= scale;
2383
2384 tangent2[0] = norm[1] * tangent1[2] - norm[2] * tangent1[1];
2385 tangent2[1] = norm[2] * tangent1[0] - norm[0] * tangent1[2];
2386 tangent2[2] = norm[0] * tangent1[1] - norm[1] * tangent1[0];
2387
2388 while (1) {
2389 r1 = sigma * random->gaussian();
2390 r2 = sigma * random->gaussian();
2391 vnmag = sqrt(r1 * r1 + r2 * r2);
2392 vtmag1 = sigma * random->gaussian();
2393 vtmag2 = sigma * random->gaussian();
2394 if (vnmag * vnmag + vtmag1 * vtmag1 + vtmag2 * vtmag2 <= vmaxsq) break;
2395 }
2396
2397 vsnew[0] = vnmag * norm[0] + vtmag1 * tangent1[0] + vtmag2 * tangent2[0];
2398 vsnew[1] = vnmag * norm[1] + vtmag1 * tangent1[1] + vtmag2 * tangent2[1];
2399 vsnew[2] = vnmag * norm[2] + vtmag1 * tangent1[2] + vtmag2 * tangent2[2];
2400
2401 // add in velocity of collision pt
2402 // for WALL: velocity of wall in one dim
2403 // else: translation/rotation of BIG particle
2404
2405 if (big->type == WALL) {
2406 int dim = wallwhich[iwall] / 2;
2407 vsnew[dim] += vwall[iwall];
2408
2409 } else {
2410 double *omega = big->omega;
2411 vsnew[0] += vb[0] + omega[1] * (xsurf[2] - xb[2]) - omega[2] * (xsurf[1] - xb[1]);
2412 vsnew[1] += vb[1] + omega[2] * (xsurf[0] - xb[0]) - omega[0] * (xsurf[2] - xb[2]);
2413 vsnew[2] += vb[2] + omega[0] * (xsurf[1] - xb[1]) - omega[1] * (xsurf[0] - xb[0]);
2414 }
2415 }
2416
2417 /* ----------------------------------------------------------------------
2418 impart force and torque to BIG particle
2419 force on BIG particle = -dp/dt of SRD particle
2420 torque on BIG particle = r cross (-dp/dt)
2421 ------------------------------------------------------------------------- */
2422
force_torque(double * vsold,double * vsnew,double * xs,double * xb,double * fb,double * tb)2423 void FixSRD::force_torque(double *vsold, double *vsnew, double *xs, double *xb, double *fb,
2424 double *tb)
2425 {
2426 double dpdt[3], xs_xb[3];
2427
2428 double factor = mass_srd / dt_big / force->ftm2v;
2429 dpdt[0] = factor * (vsnew[0] - vsold[0]);
2430 dpdt[1] = factor * (vsnew[1] - vsold[1]);
2431 dpdt[2] = factor * (vsnew[2] - vsold[2]);
2432
2433 fb[0] -= dpdt[0];
2434 fb[1] -= dpdt[1];
2435 fb[2] -= dpdt[2];
2436
2437 // no torque if SLIP collision and BIG is a sphere
2438
2439 if (tb) {
2440 xs_xb[0] = xs[0] - xb[0];
2441 xs_xb[1] = xs[1] - xb[1];
2442 xs_xb[2] = xs[2] - xb[2];
2443 tb[0] -= xs_xb[1] * dpdt[2] - xs_xb[2] * dpdt[1];
2444 tb[1] -= xs_xb[2] * dpdt[0] - xs_xb[0] * dpdt[2];
2445 tb[2] -= xs_xb[0] * dpdt[1] - xs_xb[1] * dpdt[0];
2446 }
2447 }
2448
2449 /* ----------------------------------------------------------------------
2450 impart force to WALL
2451 force on WALL = -dp/dt of SRD particle
2452 ------------------------------------------------------------------------- */
2453
force_wall(double * vsold,double * vsnew,int iwall)2454 void FixSRD::force_wall(double *vsold, double *vsnew, int iwall)
2455
2456 {
2457 double dpdt[3];
2458
2459 double factor = mass_srd / dt_big / force->ftm2v;
2460 dpdt[0] = factor * (vsnew[0] - vsold[0]);
2461 dpdt[1] = factor * (vsnew[1] - vsold[1]);
2462 dpdt[2] = factor * (vsnew[2] - vsold[2]);
2463
2464 fwall[iwall][0] -= dpdt[0];
2465 fwall[iwall][1] -= dpdt[1];
2466 fwall[iwall][2] -= dpdt[2];
2467 }
2468
2469 /* ----------------------------------------------------------------------
2470 update SRD particle position & velocity & search bin assignment
2471 check if SRD moved outside of valid region
2472 if so, may overlap off-processor BIG particle
2473 ------------------------------------------------------------------------- */
2474
update_srd(int i,double dt,double * xscoll,double * vsnew,double * xs,double * vs)2475 int FixSRD::update_srd(int i, double dt, double *xscoll, double *vsnew, double *xs, double *vs)
2476 {
2477 int ix, iy, iz;
2478
2479 vs[0] = vsnew[0];
2480 vs[1] = vsnew[1];
2481 vs[2] = vsnew[2];
2482
2483 xs[0] = xscoll[0] + dt * vsnew[0];
2484 xs[1] = xscoll[1] + dt * vsnew[1];
2485 xs[2] = xscoll[2] + dt * vsnew[2];
2486
2487 if (triclinic) domain->x2lamda(xs, xs);
2488
2489 if (xs[0] < srdlo[0] || xs[0] > srdhi[0] || xs[1] < srdlo[1] || xs[1] > srdhi[1] ||
2490 xs[2] < srdlo[2] || xs[2] > srdhi[2]) {
2491 if (screen)
2492 error->warning(FLERR,
2493 "Fix srd particle moved outside valid domain\n"
2494 " particle {} on proc {} at timestep {}\n"
2495 " xnew {:.8} {:.8} {:.8}\n"
2496 " srdlo/hi x {:.8} {:.8}\n"
2497 " srdlo/hi y {:.8} {:.8}\n"
2498 " srdlo/hi z {:.8} {:.8}\n",
2499 atom->tag[i], me, update->ntimestep, xs[0], xs[1], xs[2], srdlo[0], srdhi[0],
2500 srdlo[1], srdhi[1], srdlo[2], srdhi[2]);
2501 }
2502
2503 if (triclinic) domain->lamda2x(xs, xs);
2504
2505 ix = static_cast<int>((xs[0] - xblo2) * bininv2x);
2506 iy = static_cast<int>((xs[1] - yblo2) * bininv2y);
2507 iz = static_cast<int>((xs[2] - zblo2) * bininv2z);
2508 return iz * nbin2y * nbin2x + iy * nbin2x + ix;
2509 }
2510
2511 /* ----------------------------------------------------------------------
2512 setup all SRD parameters with big particles
2513 ------------------------------------------------------------------------- */
2514
parameterize()2515 void FixSRD::parameterize()
2516 {
2517 // timesteps
2518
2519 dt_big = update->dt;
2520 dt_srd = nevery * update->dt;
2521
2522 // maxbigdiam,minbigdiam = max/min diameter of any big particle
2523 // big particle must either have radius > 0 or shape > 0 defined
2524 // apply radfactor at end
2525
2526 AtomVecEllipsoid::Bonus *ebonus;
2527 if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
2528 AtomVecLine::Bonus *lbonus;
2529 if (avec_line) lbonus = avec_line->bonus;
2530 AtomVecTri::Bonus *tbonus;
2531 if (avec_tri) tbonus = avec_tri->bonus;
2532 double *radius = atom->radius;
2533 int *ellipsoid = atom->ellipsoid;
2534 int *line = atom->line;
2535 int *tri = atom->tri;
2536 int *mask = atom->mask;
2537 int nlocal = atom->nlocal;
2538
2539 int any_ellipsoids = 0;
2540 int any_lines = 0;
2541 int any_tris = 0;
2542 maxbigdiam = 0.0;
2543 minbigdiam = BIG;
2544
2545 for (int i = 0; i < nlocal; i++)
2546 if (mask[i] & biggroupbit) {
2547 if (radius && radius[i] > 0.0) {
2548 maxbigdiam = MAX(maxbigdiam, 2.0 * radius[i]);
2549 minbigdiam = MIN(minbigdiam, 2.0 * radius[i]);
2550 } else if (ellipsoid && ellipsoid[i] >= 0) {
2551 any_ellipsoids = 1;
2552 double *shape = ebonus[ellipsoid[i]].shape;
2553 maxbigdiam = MAX(maxbigdiam, 2.0 * shape[0]);
2554 maxbigdiam = MAX(maxbigdiam, 2.0 * shape[1]);
2555 maxbigdiam = MAX(maxbigdiam, 2.0 * shape[2]);
2556 minbigdiam = MIN(minbigdiam, 2.0 * shape[0]);
2557 minbigdiam = MIN(minbigdiam, 2.0 * shape[1]);
2558 minbigdiam = MIN(minbigdiam, 2.0 * shape[2]);
2559 } else if (line && line[i] >= 0) {
2560 any_lines = 1;
2561 double length = lbonus[line[i]].length;
2562 maxbigdiam = MAX(maxbigdiam, length);
2563 minbigdiam = MIN(minbigdiam, length);
2564 } else if (tri && tri[i] >= 0) {
2565 any_tris = 1;
2566 double length1 = MathExtra::len3(tbonus[tri[i]].c1);
2567 double length2 = MathExtra::len3(tbonus[tri[i]].c2);
2568 double length3 = MathExtra::len3(tbonus[tri[i]].c3);
2569 double length = MAX(length1, length2);
2570 length = MAX(length, length3);
2571 maxbigdiam = MAX(maxbigdiam, length);
2572 minbigdiam = MIN(minbigdiam, length);
2573 } else
2574 error->one(FLERR, "Big particle in fix srd cannot be point particle");
2575 }
2576
2577 double tmp = maxbigdiam;
2578 MPI_Allreduce(&tmp, &maxbigdiam, 1, MPI_DOUBLE, MPI_MAX, world);
2579 tmp = minbigdiam;
2580 MPI_Allreduce(&tmp, &minbigdiam, 1, MPI_DOUBLE, MPI_MIN, world);
2581
2582 maxbigdiam *= radfactor;
2583 minbigdiam *= radfactor;
2584
2585 int itmp = any_ellipsoids;
2586 MPI_Allreduce(&itmp, &any_ellipsoids, 1, MPI_INT, MPI_MAX, world);
2587 itmp = any_lines;
2588 MPI_Allreduce(&itmp, &any_lines, 1, MPI_INT, MPI_MAX, world);
2589 itmp = any_tris;
2590 MPI_Allreduce(&itmp, &any_tris, 1, MPI_INT, MPI_MAX, world);
2591
2592 if (any_lines && overlap == 0)
2593 error->all(FLERR, "Cannot use lines with fix srd unless overlap is set");
2594 if (any_tris && overlap == 0)
2595 error->all(FLERR, "Cannot use tris with fix srd unless overlap is set");
2596
2597 // big particles are only torqued if ellipsoids/lines/tris or NOSLIP
2598
2599 if (any_ellipsoids == 0 && any_lines == 0 && any_tris == 0 && collidestyle == SLIP)
2600 torqueflag = 0;
2601 else
2602 torqueflag = 1;
2603
2604 // mass of SRD particles, require monodispersity
2605
2606 double *rmass = atom->rmass;
2607 double *mass = atom->mass;
2608 int *type = atom->type;
2609
2610 int flag = 0;
2611 mass_srd = 0.0;
2612
2613 for (int i = 0; i < nlocal; i++)
2614 if (mask[i] & groupbit) {
2615 if (rmass) {
2616 if (mass_srd == 0.0)
2617 mass_srd = rmass[i];
2618 else if (rmass[i] != mass_srd)
2619 flag = 1;
2620 } else {
2621 if (mass_srd == 0.0)
2622 mass_srd = mass[type[i]];
2623 else if (mass[type[i]] != mass_srd)
2624 flag = 1;
2625 }
2626 }
2627
2628 int flagall;
2629 MPI_Allreduce(&flag, &flagall, 1, MPI_INT, MPI_MAX, world);
2630 if (flagall) error->all(FLERR, "Fix srd requires SRD particles all have same mass");
2631
2632 // set temperature and lamda of SRD particles from each other
2633 // lamda = dt_srd * sqrt(boltz * temperature_srd / mass_srd);
2634
2635 if (lamdaflag == 0)
2636 lamda = dt_srd * sqrt(force->boltz * temperature_srd / mass_srd / force->mvv2e);
2637 else
2638 temperature_srd = force->mvv2e * (lamda / dt_srd) * (lamda / dt_srd) * mass_srd / force->boltz;
2639
2640 // vmax = maximum velocity of an SRD particle
2641 // dmax = maximum distance an SRD can move = 4*lamda = vmax * dt_srd
2642
2643 sigma = lamda / dt_srd;
2644 dmax = 4.0 * lamda;
2645 vmax = dmax / dt_srd;
2646 vmaxsq = vmax * vmax;
2647
2648 // volbig = total volume of all big particles
2649 // LINE/TRI particles have no volume
2650 // incorrect if part of rigid particles, so add fudge factor with WIDTH
2651 // apply radfactor to reduce final volume
2652
2653 double volbig = 0.0;
2654 double WIDTH = 1.0;
2655
2656 if (dimension == 3) {
2657 for (int i = 0; i < nlocal; i++)
2658 if (mask[i] & biggroupbit) {
2659 if (radius && radius[i] > 0.0) {
2660 double r = radfactor * radius[i];
2661 volbig += 4.0 / 3.0 * MY_PI * r * r * r;
2662 ;
2663 } else if (ellipsoid && ellipsoid[i] >= 0) {
2664 double *shape = ebonus[ellipsoid[i]].shape;
2665 volbig += 4.0 / 3.0 * MY_PI * shape[0] * shape[1] * shape[2] * radfactor * radfactor *
2666 radfactor;
2667 } else if (tri && tri[i] >= 0) {
2668 double *c1 = tbonus[tri[i]].c1;
2669 double *c2 = tbonus[tri[i]].c2;
2670 double *c3 = tbonus[tri[i]].c3;
2671 double c2mc1[3], c3mc1[3], cross[3];
2672 MathExtra::sub3(c2, c1, c2mc1);
2673 MathExtra::sub3(c3, c1, c3mc1);
2674 MathExtra::cross3(c2mc1, c3mc1, cross);
2675 volbig += 0.5 * MathExtra::len3(cross);
2676 }
2677 }
2678 } else {
2679 for (int i = 0; i < nlocal; i++)
2680 if (mask[i] & biggroupbit) {
2681 if (radius && radius[i] > 0.0) {
2682 double r = radfactor * radius[i];
2683 volbig += MY_PI * r * r;
2684 } else if (ellipsoid && ellipsoid[i] >= 0) {
2685 double *shape = ebonus[ellipsoid[i]].shape;
2686 volbig += MY_PI * shape[0] * shape[1] * radfactor * radfactor;
2687 } else if (line && line[i] >= 0) {
2688 double length = lbonus[line[i]].length;
2689 volbig += length * WIDTH;
2690 }
2691 }
2692 }
2693
2694 tmp = volbig;
2695 MPI_Allreduce(&tmp, &volbig, 1, MPI_DOUBLE, MPI_SUM, world);
2696
2697 // particle counts
2698
2699 bigint mbig = 0;
2700 if (bigexist) mbig = group->count(biggroup);
2701 bigint nsrd = group->count(igroup);
2702
2703 // mass_big = total mass of all big particles
2704
2705 mass_big = 0.0;
2706 for (int i = 0; i < nlocal; i++)
2707 if (mask[i] & biggroupbit) {
2708 if (rmass)
2709 mass_big += rmass[i];
2710 else
2711 mass_big += mass[type[i]];
2712 }
2713
2714 tmp = mass_big;
2715 MPI_Allreduce(&tmp, &mass_big, 1, MPI_DOUBLE, MPI_SUM, world);
2716
2717 // mass density ratio = big / SRD
2718
2719 double density_big = 0.0;
2720 if (bigexist) density_big = mass_big / volbig;
2721
2722 double volsrd, density_srd;
2723 if (dimension == 3) {
2724 volsrd = (domain->xprd * domain->yprd * domain->zprd) - volbig;
2725 density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd * domain->zprd - volbig);
2726 } else {
2727 volsrd = (domain->xprd * domain->yprd) - volbig;
2728 density_srd = nsrd * mass_srd / (domain->xprd * domain->yprd - volbig);
2729 }
2730
2731 double mdratio = density_big / density_srd;
2732
2733 // create grid for binning/rotating SRD particles from gridsrd
2734
2735 setup_velocity_bins();
2736
2737 // binsize3 = binsize1 in box units (not lamda units for triclinic)
2738
2739 double binsize3x = binsize1x;
2740 double binsize3y = binsize1y;
2741 double binsize3z = binsize1z;
2742 if (triclinic) {
2743 binsize3x = binsize1x * domain->xprd;
2744 binsize3y = binsize1y * domain->yprd;
2745 binsize3z = binsize1z * domain->zprd;
2746 }
2747
2748 // srd_per_grid = # of SRD particles per SRD grid cell
2749
2750 double ncell;
2751 if (dimension == 3)
2752 ncell = volsrd / (binsize3x * binsize3y * binsize3z);
2753 else
2754 ncell = volsrd / (binsize3x * binsize3y);
2755
2756 srd_per_cell = (double) nsrd / ncell;
2757
2758 // kinematic viscosity of SRD fluid
2759 // output in cm^2/sec units, converted by xxt2kmu
2760
2761 double viscosity;
2762 if (dimension == 3)
2763 viscosity =
2764 gridsrd * gridsrd / (18.0 * dt_srd) * (1.0 - (1.0 - exp(-srd_per_cell)) / srd_per_cell) +
2765 (force->boltz * temperature_srd * dt_srd / (4.0 * mass_srd * force->mvv2e)) *
2766 ((srd_per_cell + 2.0) / (srd_per_cell - 1.0));
2767 else
2768 viscosity = (force->boltz * temperature_srd * dt_srd / (2.0 * mass_srd * force->mvv2e)) *
2769 (srd_per_cell / (srd_per_cell - 1.0 + exp(-srd_per_cell)) - 1.0) +
2770 (gridsrd * gridsrd) / (12.0 * dt_srd) *
2771 ((srd_per_cell - 1.0 + exp(-srd_per_cell)) / srd_per_cell);
2772 viscosity *= force->xxt2kmu;
2773
2774 // print SRD parameters
2775
2776 if (me == 0) {
2777 std::string mesg = "SRD info:\n";
2778 mesg += fmt::format(" SRD/big particles = {} {}\n", nsrd, mbig);
2779 mesg += fmt::format(" big particle diameter max/min = {:.8} {:.8}\n", maxbigdiam, minbigdiam);
2780 mesg += fmt::format(" SRD temperature & lamda = {:.8} {:.8}\n", temperature_srd, lamda);
2781 mesg += fmt::format(" SRD max distance & max velocity = {:.8} {:.8}\n", dmax, vmax);
2782 mesg += fmt::format(" SRD grid counts: {} {} {}\n", nbin1x, nbin1y, nbin1z);
2783 mesg += fmt::format(" SRD grid size: request, actual (xyz) = {:.8}, {:.8} {:.8} {:.8}\n",
2784 gridsrd, binsize3x, binsize3y, binsize3z);
2785 mesg += fmt::format(" SRD per actual grid cell = {:.8}\n", srd_per_cell);
2786 mesg += fmt::format(" SRD viscosity = {:.8}\n", viscosity);
2787 mesg += fmt::format(" big/SRD mass density ratio = {:.8}\n", mdratio);
2788 utils::logmesg(lmp, mesg);
2789 }
2790
2791 // error if less than 1 SRD bin per processor in some dim
2792
2793 if (nbin1x < comm->procgrid[0] || nbin1y < comm->procgrid[1] || nbin1z < comm->procgrid[2])
2794 error->all(FLERR, "Fewer SRD bins than processors in some dimension");
2795
2796 // check if SRD bins are within tolerance for shape and size
2797
2798 int tolflag = 0;
2799 if (binsize3y / binsize3x > 1.0 + cubictol || binsize3x / binsize3y > 1.0 + cubictol) tolflag = 1;
2800 if (dimension == 3) {
2801 if (binsize3z / binsize3x > 1.0 + cubictol || binsize3x / binsize3z > 1.0 + cubictol)
2802 tolflag = 1;
2803 }
2804
2805 if (tolflag) {
2806 if (cubicflag == CUBIC_ERROR) error->all(FLERR, "SRD bins for fix srd are not cubic enough");
2807 if (me == 0) error->warning(FLERR, "SRD bins for fix srd are not cubic enough");
2808 }
2809
2810 tolflag = 0;
2811 if (binsize3x / gridsrd > 1.0 + cubictol || gridsrd / binsize3x > 1.0 + cubictol) tolflag = 1;
2812 if (binsize3y / gridsrd > 1.0 + cubictol || gridsrd / binsize3y > 1.0 + cubictol) tolflag = 1;
2813 if (dimension == 3) {
2814 if (binsize3z / gridsrd > 1.0 + cubictol || gridsrd / binsize3z > 1.0 + cubictol) tolflag = 1;
2815 }
2816
2817 if (tolflag) {
2818 if (cubicflag == CUBIC_ERROR)
2819 error->all(FLERR, "SRD bin size for fix srd differs from user request");
2820 if (me == 0) error->warning(FLERR, "SRD bin size for fix srd differs from user request");
2821 }
2822
2823 // error if lamda < 0.6 of SRD grid size and no shifting allowed
2824 // turn on shifting in this case if allowed
2825
2826 double maxgridsrd = MAX(binsize3x, binsize3y);
2827 if (dimension == 3) maxgridsrd = MAX(maxgridsrd, binsize3z);
2828
2829 shiftflag = 0;
2830 if (lamda < 0.6 * maxgridsrd && shiftuser == SHIFT_NO)
2831 error->all(FLERR, "Fix srd lamda must be >= 0.6 of SRD grid size");
2832 else if (lamda < 0.6 * maxgridsrd && shiftuser == SHIFT_POSSIBLE) {
2833 shiftflag = 1;
2834 if (me == 0) error->warning(FLERR, "SRD bin shifting turned on due to small lamda");
2835 } else if (shiftuser == SHIFT_YES)
2836 shiftflag = 1;
2837
2838 // warnings
2839
2840 if (bigexist && maxgridsrd > 0.25 * minbigdiam && me == 0)
2841 error->warning(FLERR, "Fix srd grid size > 1/4 of big particle diameter");
2842 if (viscosity < 0.0 && me == 0)
2843 error->warning(FLERR, "Fix srd viscosity < 0.0 due to low SRD density");
2844 if (bigexist && dt_big * vmax > minbigdiam && me == 0)
2845 error->warning(FLERR, "Fix srd particles may move > big particle diameter");
2846 }
2847
2848 /* ----------------------------------------------------------------------
2849 set static parameters of each big particle, owned and ghost
2850 called each reneighboring
2851 use radfactor in distance parameters as appropriate
2852 ------------------------------------------------------------------------- */
2853
big_static()2854 void FixSRD::big_static()
2855 {
2856 int i;
2857 double rad, arad, brad, crad, length, length1, length2, length3;
2858 double *shape, *c1, *c2, *c3;
2859 double c2mc1[3], c3mc1[3];
2860
2861 AtomVecEllipsoid::Bonus *ebonus;
2862 if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
2863 AtomVecLine::Bonus *lbonus;
2864 if (avec_line) lbonus = avec_line->bonus;
2865 AtomVecTri::Bonus *tbonus;
2866 if (avec_tri) tbonus = avec_tri->bonus;
2867 double *radius = atom->radius;
2868 int *ellipsoid = atom->ellipsoid;
2869 int *line = atom->line;
2870 int *tri = atom->tri;
2871
2872 double skinhalf = 0.5 * neighbor->skin;
2873
2874 for (int k = 0; k < nbig; k++) {
2875 i = biglist[k].index;
2876
2877 // sphere
2878 // set radius and radsq and cutoff based on radius
2879
2880 if (radius && radius[i] > 0.0) {
2881 biglist[k].type = SPHERE;
2882 rad = radfactor * radius[i];
2883 biglist[k].radius = rad;
2884 biglist[k].radsq = rad * rad;
2885 biglist[k].cutbinsq = (rad + skinhalf) * (rad + skinhalf);
2886
2887 // ellipsoid
2888 // set abc radsqinv and cutoff based on max radius
2889
2890 } else if (ellipsoid && ellipsoid[i] >= 0) {
2891 shape = ebonus[ellipsoid[i]].shape;
2892 biglist[k].type = ELLIPSOID;
2893 arad = radfactor * shape[0];
2894 brad = radfactor * shape[1];
2895 crad = radfactor * shape[2];
2896 biglist[k].aradsqinv = 1.0 / (arad * arad);
2897 biglist[k].bradsqinv = 1.0 / (brad * brad);
2898 biglist[k].cradsqinv = 1.0 / (crad * crad);
2899 rad = MAX(arad, brad);
2900 rad = MAX(rad, crad);
2901 biglist[k].cutbinsq = (rad + skinhalf) * (rad + skinhalf);
2902
2903 // line
2904 // set length and cutoff based on 1/2 length
2905
2906 } else if (line && line[i] >= 0) {
2907 length = lbonus[line[i]].length;
2908 biglist[k].type = LINE;
2909 biglist[k].length = length;
2910 rad = 0.5 * length;
2911 biglist[k].cutbinsq = (rad + skinhalf) * (rad + skinhalf);
2912
2913 // tri
2914 // set normbody based on c1,c2,c3
2915 // set cutoff based on point furthest from centroid
2916
2917 } else if (tri && tri[i] >= 0) {
2918 biglist[k].type = TRIANGLE;
2919 c1 = tbonus[tri[i]].c1;
2920 c2 = tbonus[tri[i]].c2;
2921 c3 = tbonus[tri[i]].c3;
2922 MathExtra::sub3(c2, c1, c2mc1);
2923 MathExtra::sub3(c3, c1, c3mc1);
2924 MathExtra::cross3(c2mc1, c3mc1, biglist[k].normbody);
2925 length1 = MathExtra::len3(c1);
2926 length2 = MathExtra::len3(c2);
2927 length3 = MathExtra::len3(c3);
2928 rad = MAX(length1, length2);
2929 rad = MAX(rad, length3);
2930 biglist[k].cutbinsq = (rad + skinhalf) * (rad + skinhalf);
2931 }
2932 }
2933 }
2934
2935 /* ----------------------------------------------------------------------
2936 set dynamic parameters of each big particle, owned and ghost
2937 called each timestep
2938 ------------------------------------------------------------------------- */
2939
big_dynamic()2940 void FixSRD::big_dynamic()
2941 {
2942 int i;
2943 double *shape, *quat, *inertia;
2944 double inertiaone[3];
2945
2946 AtomVecEllipsoid::Bonus *ebonus;
2947 if (avec_ellipsoid) ebonus = avec_ellipsoid->bonus;
2948 AtomVecLine::Bonus *lbonus;
2949 if (avec_line) lbonus = avec_line->bonus;
2950 AtomVecTri::Bonus *tbonus;
2951 if (avec_tri) tbonus = avec_tri->bonus;
2952 double **omega = atom->omega;
2953 double **angmom = atom->angmom;
2954 double *rmass = atom->rmass;
2955 int *ellipsoid = atom->ellipsoid;
2956 int *line = atom->line;
2957 int *tri = atom->tri;
2958
2959 for (int k = 0; k < nbig; k++) {
2960 i = biglist[k].index;
2961
2962 // sphere
2963 // set omega from atom->omega directly
2964
2965 if (biglist[k].type == SPHERE) {
2966 biglist[k].omega[0] = omega[i][0];
2967 biglist[k].omega[1] = omega[i][1];
2968 biglist[k].omega[2] = omega[i][2];
2969
2970 // ellipsoid
2971 // set ex,ey,ez from quaternion
2972 // set omega from angmom & ex,ey,ez
2973
2974 } else if (biglist[k].type == ELLIPSOID) {
2975 quat = ebonus[ellipsoid[i]].quat;
2976 MathExtra::q_to_exyz(quat, biglist[k].ex, biglist[k].ey, biglist[k].ez);
2977 shape = ebonus[ellipsoid[i]].shape;
2978 inertiaone[0] = EINERTIA * rmass[i] * (shape[1] * shape[1] + shape[2] * shape[2]);
2979 inertiaone[1] = EINERTIA * rmass[i] * (shape[0] * shape[0] + shape[2] * shape[2]);
2980 inertiaone[2] = EINERTIA * rmass[i] * (shape[0] * shape[0] + shape[1] * shape[1]);
2981 MathExtra::angmom_to_omega(angmom[i], biglist[k].ex, biglist[k].ey, biglist[k].ez, inertiaone,
2982 biglist[k].omega);
2983
2984 // line
2985 // set omega from atom->omega directly
2986
2987 } else if (biglist[k].type == LINE) {
2988 biglist[k].theta = lbonus[line[i]].theta;
2989 biglist[k].omega[0] = omega[i][0];
2990 biglist[k].omega[1] = omega[i][1];
2991 biglist[k].omega[2] = omega[i][2];
2992
2993 // tri
2994 // set ex,ey,ez from quaternion
2995 // set omega from angmom & ex,ey,ez
2996 // set unit space-frame norm from body-frame norm
2997
2998 } else if (biglist[k].type == TRIANGLE) {
2999 quat = tbonus[tri[i]].quat;
3000 MathExtra::q_to_exyz(quat, biglist[k].ex, biglist[k].ey, biglist[k].ez);
3001 inertia = tbonus[tri[i]].inertia;
3002 MathExtra::angmom_to_omega(angmom[i], biglist[k].ex, biglist[k].ey, biglist[k].ez, inertia,
3003 biglist[k].omega);
3004 MathExtra::matvec(biglist[k].ex, biglist[k].ey, biglist[k].ez, biglist[k].normbody,
3005 biglist[k].norm);
3006 MathExtra::norm3(biglist[k].norm);
3007 }
3008 }
3009 }
3010
3011 /* ----------------------------------------------------------------------
3012 set bounds for big and SRD particle movement
3013 called at setup() and when box size changes (but not shape)
3014 ------------------------------------------------------------------------- */
3015
setup_bounds()3016 void FixSRD::setup_bounds()
3017 {
3018 // triclinic scale factors
3019 // convert a real distance (perpendicular to box face) to a lamda distance
3020
3021 double length0, length1, length2;
3022 if (triclinic) {
3023 double *h_inv = domain->h_inv;
3024 length0 = sqrt(h_inv[0] * h_inv[0] + h_inv[5] * h_inv[5] + h_inv[4] * h_inv[4]);
3025 length1 = sqrt(h_inv[1] * h_inv[1] + h_inv[3] * h_inv[3]);
3026 length2 = h_inv[2];
3027 }
3028
3029 // collision object = CO = big particle or wall
3030 // big particles can be owned or ghost or unknown, walls are all owned
3031 // dist_ghost = distance from sub-domain (SD) that
3032 // owned/ghost CO may move to before reneigh,
3033 // used to bound search bins in setup_search_bins()
3034 // dist_srd = distance from SD at which SRD could collide with unknown CO,
3035 // used to error check bounds of SRD movement after collisions via srdlo/hi
3036 // dist_srd_reneigh = distance from SD at which an SRD should trigger
3037 // a reneigh, b/c next SRD move might overlap with unknown CO,
3038 // used for SRD triggering of reneighboring via srdlo/hi_reneigh
3039 // onemove = max distance an SRD can move in one step
3040 // if big_particles (and possibly walls):
3041 // dist_ghost = cut + 1/2 skin due to moving away before reneigh
3042 // dist_srd = cut - 1/2 skin - 1/2 diam due to ghost CO moving towards
3043 // dist_reneigh = dist_srd - onemove
3044 // if walls and no big particles:
3045 // dist_ghost = 0.0, since not used
3046 // if no big particles or walls:
3047 // dist_ghost and dist_srd = 0.0, since not used since no search bins
3048 // dist_srd_reneigh = subsize - onemove =
3049 // max distance to move without being lost during comm->exchange()
3050 // subsize = perp distance between sub-domain faces (orthog or triclinic)
3051
3052 double cut = MAX(neighbor->cutneighmax, comm->cutghostuser);
3053 double onemove = dt_big * vmax;
3054
3055 if (bigexist) {
3056 dist_ghost = cut + 0.5 * neighbor->skin;
3057 dist_srd = cut - 0.5 * neighbor->skin - 0.5 * maxbigdiam;
3058 dist_srd_reneigh = dist_srd - onemove;
3059 } else if (wallexist) {
3060 dist_ghost = 4 * onemove;
3061 dist_srd = 4 * onemove;
3062 dist_srd_reneigh = 4 * onemove - onemove;
3063 } else {
3064 dist_ghost = dist_srd = 0.0;
3065 double subsize;
3066 if (triclinic == 0) {
3067 subsize = domain->prd[0] / comm->procgrid[0];
3068 subsize = MIN(subsize, domain->prd[1] / comm->procgrid[1]);
3069 if (dimension == 3) subsize = MIN(subsize, domain->prd[2] / comm->procgrid[2]);
3070 } else {
3071 subsize = 1.0 / comm->procgrid[0] / length0;
3072 subsize = MIN(subsize, 1.0 / comm->procgrid[1] / length1);
3073 if (dimension == 3) subsize = MIN(subsize, 1.0 / comm->procgrid[2] / length2);
3074 }
3075 dist_srd_reneigh = subsize - onemove;
3076 }
3077
3078 // lo/hi = bbox on this proc which SRD particles must stay inside
3079 // lo/hi reneigh = bbox on this proc outside of which SRDs trigger a reneigh
3080 // for triclinic, these bbox are in lamda units
3081
3082 if (triclinic == 0) {
3083 srdlo[0] = domain->sublo[0] - dist_srd;
3084 srdhi[0] = domain->subhi[0] + dist_srd;
3085 srdlo[1] = domain->sublo[1] - dist_srd;
3086 srdhi[1] = domain->subhi[1] + dist_srd;
3087 srdlo[2] = domain->sublo[2] - dist_srd;
3088 srdhi[2] = domain->subhi[2] + dist_srd;
3089
3090 srdlo_reneigh[0] = domain->sublo[0] - dist_srd_reneigh;
3091 srdhi_reneigh[0] = domain->subhi[0] + dist_srd_reneigh;
3092 srdlo_reneigh[1] = domain->sublo[1] - dist_srd_reneigh;
3093 srdhi_reneigh[1] = domain->subhi[1] + dist_srd_reneigh;
3094 srdlo_reneigh[2] = domain->sublo[2] - dist_srd_reneigh;
3095 srdhi_reneigh[2] = domain->subhi[2] + dist_srd_reneigh;
3096
3097 } else {
3098 srdlo[0] = domain->sublo_lamda[0] - dist_srd * length0;
3099 srdhi[0] = domain->subhi_lamda[0] + dist_srd * length0;
3100 srdlo[1] = domain->sublo_lamda[1] - dist_srd * length1;
3101 srdhi[1] = domain->subhi_lamda[1] + dist_srd * length1;
3102 srdlo[2] = domain->sublo_lamda[2] - dist_srd * length2;
3103 srdhi[2] = domain->subhi_lamda[2] + dist_srd * length2;
3104
3105 srdlo_reneigh[0] = domain->sublo_lamda[0] - dist_srd_reneigh * length0;
3106 srdhi_reneigh[0] = domain->subhi_lamda[0] + dist_srd_reneigh * length0;
3107 srdlo_reneigh[1] = domain->sublo_lamda[1] - dist_srd_reneigh * length1;
3108 srdhi_reneigh[1] = domain->subhi_lamda[1] + dist_srd_reneigh * length1;
3109 srdlo_reneigh[2] = domain->sublo_lamda[2] - dist_srd_reneigh * length2;
3110 srdhi_reneigh[2] = domain->subhi_lamda[2] + dist_srd_reneigh * length2;
3111 }
3112 }
3113
3114 /* ----------------------------------------------------------------------
3115 setup bins used for binning SRD particles for velocity reset
3116 gridsrd = desired bin size
3117 also setup bin shifting parameters
3118 also setup comm of bins that straddle processor boundaries
3119 called at beginning of each run
3120 called every reneighbor if box size changes, but not if box shape changes
3121 ------------------------------------------------------------------------- */
3122
setup_velocity_bins()3123 void FixSRD::setup_velocity_bins()
3124 {
3125 // require integer # of bins across global domain
3126
3127 nbin1x = static_cast<int>(domain->xprd / gridsrd + 0.5);
3128 nbin1y = static_cast<int>(domain->yprd / gridsrd + 0.5);
3129 nbin1z = static_cast<int>(domain->zprd / gridsrd + 0.5);
3130 if (dimension == 2) nbin1z = 1;
3131
3132 if (nbin1x == 0) nbin1x = 1;
3133 if (nbin1y == 0) nbin1y = 1;
3134 if (nbin1z == 0) nbin1z = 1;
3135
3136 if (triclinic == 0) {
3137 binsize1x = domain->xprd / nbin1x;
3138 binsize1y = domain->yprd / nbin1y;
3139 binsize1z = domain->zprd / nbin1z;
3140 bininv1x = 1.0 / binsize1x;
3141 bininv1y = 1.0 / binsize1y;
3142 bininv1z = 1.0 / binsize1z;
3143 } else {
3144 binsize1x = 1.0 / nbin1x;
3145 binsize1y = 1.0 / nbin1y;
3146 binsize1z = 1.0 / nbin1z;
3147 bininv1x = nbin1x;
3148 bininv1y = nbin1y;
3149 bininv1z = nbin1z;
3150 }
3151
3152 nbins1 = nbin1x * nbin1y * nbin1z;
3153
3154 // setup two shifts, 0 = no shift, 1 = shift
3155 // initialize no shift case since static
3156 // shift case is dynamic, has to be initialized each time shift occurs
3157 // setup_velocity_shift allocates memory for vbin and sendlist/recvlist
3158
3159 double *boxlo;
3160 if (triclinic == 0)
3161 boxlo = domain->boxlo;
3162 else
3163 boxlo = domain->boxlo_lamda;
3164
3165 shifts[0].corner[0] = boxlo[0];
3166 shifts[0].corner[1] = boxlo[1];
3167 shifts[0].corner[2] = boxlo[2];
3168 setup_velocity_shift(0, 0);
3169
3170 shifts[1].corner[0] = boxlo[0];
3171 shifts[1].corner[1] = boxlo[1];
3172 shifts[1].corner[2] = boxlo[2];
3173 setup_velocity_shift(1, 0);
3174
3175 // allocate binhead based on max # of bins in either shift
3176
3177 int max = shifts[0].nbins;
3178 max = MAX(max, shifts[1].nbins);
3179
3180 if (max > maxbin1) {
3181 memory->destroy(binhead);
3182 maxbin1 = max;
3183 memory->create(binhead, max, "fix/srd:binhead");
3184 }
3185
3186 // allocate sbuf,rbuf based on biggest bin message
3187
3188 max = 0;
3189 for (int ishift = 0; ishift < 2; ishift++)
3190 for (int iswap = 0; iswap < 2 * dimension; iswap++) {
3191 max = MAX(max, shifts[ishift].bcomm[iswap].nsend);
3192 max = MAX(max, shifts[ishift].bcomm[iswap].nrecv);
3193 }
3194
3195 if (max > maxbuf) {
3196 memory->destroy(sbuf1);
3197 memory->destroy(sbuf2);
3198 memory->destroy(rbuf1);
3199 memory->destroy(rbuf2);
3200 maxbuf = max;
3201 memory->create(sbuf1, max * VBINSIZE, "fix/srd:sbuf");
3202 memory->create(sbuf2, max * VBINSIZE, "fix/srd:sbuf");
3203 memory->create(rbuf1, max * VBINSIZE, "fix/srd:rbuf");
3204 memory->create(rbuf2, max * VBINSIZE, "fix/srd:rbuf");
3205 }
3206
3207 // commflag = 1 if any comm required due to bins overlapping proc boundaries
3208
3209 shifts[0].commflag = 0;
3210 if (nbin1x % comm->procgrid[0]) shifts[0].commflag = 1;
3211 if (nbin1y % comm->procgrid[1]) shifts[0].commflag = 1;
3212 if (nbin1z % comm->procgrid[2]) shifts[0].commflag = 1;
3213 shifts[1].commflag = 1;
3214 }
3215
3216 /* ----------------------------------------------------------------------
3217 setup velocity shift parameters
3218 set binlo[]/binhi[] and nbins,nbinx,nbiny,nbinz for this proc
3219 set bcomm[6] params based on bin overlaps with proc boundaries
3220 no comm of bins across non-periodic global boundaries
3221 set vbin owner flags for bins I am owner of
3222 ishift = 0, dynamic = 0:
3223 set all settings since are static
3224 allocate and set bcomm params and vbins
3225 do not comm bins that align with proc boundaries
3226 ishift = 1, dynamic = 0:
3227 set max bounds on bin counts and message sizes
3228 allocate and set bcomm params and vbins based on max bounds
3229 other settings will later change dynamically
3230 ishift = 1, dynamic = 1:
3231 set actual bin bounds and counts for specific shift
3232 set bcomm params and vbins (already allocated)
3233 called by setup_velocity_bins() and reset_velocities()
3234 ------------------------------------------------------------------------- */
3235
setup_velocity_shift(int ishift,int dynamic)3236 void FixSRD::setup_velocity_shift(int ishift, int dynamic)
3237 {
3238 int i, j, k, m, id, nsend;
3239 int *sendlist;
3240 BinComm *first, *second;
3241 BinAve *vbin;
3242
3243 double *sublo, *subhi;
3244 if (triclinic == 0) {
3245 sublo = domain->sublo;
3246 subhi = domain->subhi;
3247 } else {
3248 sublo = domain->sublo_lamda;
3249 subhi = domain->subhi_lamda;
3250 }
3251
3252 int *binlo = shifts[ishift].binlo;
3253 int *binhi = shifts[ishift].binhi;
3254 double *corner = shifts[ishift].corner;
3255 int *procgrid = comm->procgrid;
3256 int *myloc = comm->myloc;
3257
3258 binlo[0] = static_cast<int>((sublo[0] - corner[0]) * bininv1x);
3259 binlo[1] = static_cast<int>((sublo[1] - corner[1]) * bininv1y);
3260 binlo[2] = static_cast<int>((sublo[2] - corner[2]) * bininv1z);
3261 if (dimension == 2) shifts[ishift].binlo[2] = 0;
3262
3263 binhi[0] = static_cast<int>((subhi[0] - corner[0]) * bininv1x);
3264 binhi[1] = static_cast<int>((subhi[1] - corner[1]) * bininv1y);
3265 binhi[2] = static_cast<int>((subhi[2] - corner[2]) * bininv1z);
3266 if (dimension == 2) shifts[ishift].binhi[2] = 0;
3267
3268 if (ishift == 0) {
3269 if (myloc[0] * nbin1x % procgrid[0] == 0) binlo[0] = myloc[0] * nbin1x / procgrid[0];
3270 if (myloc[1] * nbin1y % procgrid[1] == 0) binlo[1] = myloc[1] * nbin1y / procgrid[1];
3271 if (myloc[2] * nbin1z % procgrid[2] == 0) binlo[2] = myloc[2] * nbin1z / procgrid[2];
3272
3273 if ((myloc[0] + 1) * nbin1x % procgrid[0] == 0)
3274 binhi[0] = (myloc[0] + 1) * nbin1x / procgrid[0] - 1;
3275 if ((myloc[1] + 1) * nbin1y % procgrid[1] == 0)
3276 binhi[1] = (myloc[1] + 1) * nbin1y / procgrid[1] - 1;
3277 if ((myloc[2] + 1) * nbin1z % procgrid[2] == 0)
3278 binhi[2] = (myloc[2] + 1) * nbin1z / procgrid[2] - 1;
3279 }
3280
3281 int nbinx = binhi[0] - binlo[0] + 1;
3282 int nbiny = binhi[1] - binlo[1] + 1;
3283 int nbinz = binhi[2] - binlo[2] + 1;
3284
3285 // allow for one extra bin if shifting will occur
3286
3287 if (ishift == 1 && dynamic == 0) {
3288 nbinx++;
3289 nbiny++;
3290 if (dimension == 3) nbinz++;
3291 }
3292
3293 int nbins = nbinx * nbiny * nbinz;
3294 int nbinxy = nbinx * nbiny;
3295 int nbinsq = nbinx * nbiny;
3296 nbinsq = MAX(nbiny * nbinz, nbinsq);
3297 nbinsq = MAX(nbinx * nbinz, nbinsq);
3298
3299 shifts[ishift].nbins = nbins;
3300 shifts[ishift].nbinx = nbinx;
3301 shifts[ishift].nbiny = nbiny;
3302 shifts[ishift].nbinz = nbinz;
3303
3304 int reallocflag = 0;
3305 if (dynamic == 0 && nbinsq > shifts[ishift].maxbinsq) {
3306 shifts[ishift].maxbinsq = nbinsq;
3307 reallocflag = 1;
3308 }
3309
3310 // bcomm neighbors
3311 // first = send in lo direction, recv from hi direction
3312 // second = send in hi direction, recv from lo direction
3313
3314 if (dynamic == 0) {
3315 shifts[ishift].bcomm[0].sendproc = comm->procneigh[0][0];
3316 shifts[ishift].bcomm[0].recvproc = comm->procneigh[0][1];
3317 shifts[ishift].bcomm[1].sendproc = comm->procneigh[0][1];
3318 shifts[ishift].bcomm[1].recvproc = comm->procneigh[0][0];
3319 shifts[ishift].bcomm[2].sendproc = comm->procneigh[1][0];
3320 shifts[ishift].bcomm[2].recvproc = comm->procneigh[1][1];
3321 shifts[ishift].bcomm[3].sendproc = comm->procneigh[1][1];
3322 shifts[ishift].bcomm[3].recvproc = comm->procneigh[1][0];
3323 shifts[ishift].bcomm[4].sendproc = comm->procneigh[2][0];
3324 shifts[ishift].bcomm[4].recvproc = comm->procneigh[2][1];
3325 shifts[ishift].bcomm[5].sendproc = comm->procneigh[2][1];
3326 shifts[ishift].bcomm[5].recvproc = comm->procneigh[2][0];
3327 }
3328
3329 // set nsend,nrecv and sendlist,recvlist for each swap in x,y,z
3330 // set nsend,nrecv = 0 if static bins align with proc boundary
3331 // or to prevent dynamic bin swapping across non-periodic global boundary
3332 // allocate sendlist,recvlist only for dynamic = 0
3333
3334 first = &shifts[ishift].bcomm[0];
3335 second = &shifts[ishift].bcomm[1];
3336
3337 first->nsend = first->nrecv = second->nsend = second->nrecv = nbiny * nbinz;
3338 if (ishift == 0) {
3339 if (myloc[0] * nbin1x % procgrid[0] == 0) first->nsend = second->nrecv = 0;
3340 if ((myloc[0] + 1) * nbin1x % procgrid[0] == 0) second->nsend = first->nrecv = 0;
3341 } else {
3342 if (domain->xperiodic == 0) {
3343 if (myloc[0] == 0) first->nsend = second->nrecv = 0;
3344 if (myloc[0] == procgrid[0] - 1) second->nsend = first->nrecv = 0;
3345 }
3346 }
3347
3348 if (reallocflag) {
3349 memory->destroy(first->sendlist);
3350 memory->destroy(first->recvlist);
3351 memory->destroy(second->sendlist);
3352 memory->destroy(second->recvlist);
3353 memory->create(first->sendlist, nbinsq, "fix/srd:sendlist");
3354 memory->create(first->recvlist, nbinsq, "fix/srd:sendlist");
3355 memory->create(second->sendlist, nbinsq, "fix/srd:sendlist");
3356 memory->create(second->recvlist, nbinsq, "fix/srd:sendlist");
3357 }
3358
3359 m = 0;
3360 i = 0;
3361 for (j = 0; j < nbiny; j++)
3362 for (k = 0; k < nbinz; k++) {
3363 id = k * nbinxy + j * nbinx + i;
3364 first->sendlist[m] = second->recvlist[m] = id;
3365 m++;
3366 }
3367 m = 0;
3368 i = nbinx - 1;
3369 for (j = 0; j < nbiny; j++)
3370 for (k = 0; k < nbinz; k++) {
3371 id = k * nbinxy + j * nbinx + i;
3372 second->sendlist[m] = first->recvlist[m] = id;
3373 m++;
3374 }
3375
3376 first = &shifts[ishift].bcomm[2];
3377 second = &shifts[ishift].bcomm[3];
3378
3379 first->nsend = first->nrecv = second->nsend = second->nrecv = nbinx * nbinz;
3380 if (ishift == 0) {
3381 if (myloc[1] * nbin1y % procgrid[1] == 0) first->nsend = second->nrecv = 0;
3382 if ((myloc[1] + 1) * nbin1y % procgrid[1] == 0) second->nsend = first->nrecv = 0;
3383 } else {
3384 if (domain->yperiodic == 0) {
3385 if (myloc[1] == 0) first->nsend = second->nrecv = 0;
3386 if (myloc[1] == procgrid[1] - 1) second->nsend = first->nrecv = 0;
3387 }
3388 }
3389
3390 if (reallocflag) {
3391 memory->destroy(first->sendlist);
3392 memory->destroy(first->recvlist);
3393 memory->destroy(second->sendlist);
3394 memory->destroy(second->recvlist);
3395 memory->create(first->sendlist, nbinsq, "fix/srd:sendlist");
3396 memory->create(first->recvlist, nbinsq, "fix/srd:sendlist");
3397 memory->create(second->sendlist, nbinsq, "fix/srd:sendlist");
3398 memory->create(second->recvlist, nbinsq, "fix/srd:sendlist");
3399 }
3400
3401 m = 0;
3402 j = 0;
3403 for (i = 0; i < nbinx; i++)
3404 for (k = 0; k < nbinz; k++) {
3405 id = k * nbinxy + j * nbinx + i;
3406 first->sendlist[m] = second->recvlist[m] = id;
3407 m++;
3408 }
3409 m = 0;
3410 j = nbiny - 1;
3411 for (i = 0; i < nbinx; i++)
3412 for (k = 0; k < nbinz; k++) {
3413 id = k * nbinxy + j * nbinx + i;
3414 second->sendlist[m] = first->recvlist[m] = id;
3415 m++;
3416 }
3417
3418 if (dimension == 3) {
3419 first = &shifts[ishift].bcomm[4];
3420 second = &shifts[ishift].bcomm[5];
3421
3422 first->nsend = first->nrecv = second->nsend = second->nrecv = nbinx * nbiny;
3423 if (ishift == 0) {
3424 if (myloc[2] * nbin1z % procgrid[2] == 0) first->nsend = second->nrecv = 0;
3425 if ((myloc[2] + 1) * nbin1z % procgrid[2] == 0) second->nsend = first->nrecv = 0;
3426 } else {
3427 if (domain->zperiodic == 0) {
3428 if (myloc[2] == 0) first->nsend = second->nrecv = 0;
3429 if (myloc[2] == procgrid[2] - 1) second->nsend = first->nrecv = 0;
3430 }
3431 }
3432
3433 if (reallocflag) {
3434 memory->destroy(first->sendlist);
3435 memory->destroy(first->recvlist);
3436 memory->destroy(second->sendlist);
3437 memory->destroy(second->recvlist);
3438 memory->create(first->sendlist, nbinx * nbiny, "fix/srd:sendlist");
3439 memory->create(first->recvlist, nbinx * nbiny, "fix/srd:sendlist");
3440 memory->create(second->sendlist, nbinx * nbiny, "fix/srd:sendlist");
3441 memory->create(second->recvlist, nbinx * nbiny, "fix/srd:sendlist");
3442 }
3443
3444 m = 0;
3445 k = 0;
3446 for (i = 0; i < nbinx; i++)
3447 for (j = 0; j < nbiny; j++) {
3448 id = k * nbinxy + j * nbinx + i;
3449 first->sendlist[m] = second->recvlist[m] = id;
3450 m++;
3451 }
3452 m = 0;
3453 k = nbinz - 1;
3454 for (i = 0; i < nbinx; i++)
3455 for (j = 0; j < nbiny; j++) {
3456 id = k * nbinxy + j * nbinx + i;
3457 second->sendlist[m] = first->recvlist[m] = id;
3458 m++;
3459 }
3460 }
3461
3462 // allocate vbins, only for dynamic = 0
3463
3464 if (dynamic == 0 && nbins > shifts[ishift].maxvbin) {
3465 memory->destroy(shifts[ishift].vbin);
3466 shifts[ishift].maxvbin = nbins;
3467 shifts[ishift].vbin = (BinAve *) memory->smalloc(nbins * sizeof(BinAve), "fix/srd:vbin");
3468 }
3469
3470 // for vbins I own, set owner = 1
3471 // if bin never sent to anyone, I own it
3472 // if bin sent to lower numbered proc, I do not own it
3473 // if bin sent to self, I do not own it on even swap (avoids double counting)
3474
3475 vbin = shifts[ishift].vbin;
3476 for (i = 0; i < nbins; i++) vbin[i].owner = 1;
3477 for (int iswap = 0; iswap < 2 * dimension; iswap++) {
3478 if (shifts[ishift].bcomm[iswap].sendproc > me) continue;
3479 if (shifts[ishift].bcomm[iswap].sendproc == me && iswap % 2 == 0) continue;
3480 nsend = shifts[ishift].bcomm[iswap].nsend;
3481 sendlist = shifts[ishift].bcomm[iswap].sendlist;
3482 for (m = 0; m < nsend; m++) vbin[sendlist[m]].owner = 0;
3483 }
3484
3485 // if tstat and deformflag:
3486 // set xctr for all nbins in lamda units so can later compute vstream of bin
3487
3488 if (tstat && deformflag) {
3489 m = 0;
3490 for (k = 0; k < nbinz; k++)
3491 for (j = 0; j < nbiny; j++)
3492 for (i = 0; i < nbinx; i++) {
3493 vbin[m].xctr[0] = corner[0] + (i + binlo[0] + 0.5) / nbin1x;
3494 vbin[m].xctr[1] = corner[1] + (j + binlo[1] + 0.5) / nbin1y;
3495 vbin[m].xctr[2] = corner[2] + (k + binlo[2] + 0.5) / nbin1z;
3496 m++;
3497 }
3498 }
3499 }
3500
3501 /* ----------------------------------------------------------------------
3502 setup bins used for big and SRD particle searching
3503 gridsearch = desired bin size
3504 bins are orthogonal whether simulation box is orthogonal or triclinic
3505 for orthogonal boxes, called at each setup since cutghost may change
3506 for triclinic boxes, called at every reneigh, since tilt may change
3507 sets the following:
3508 nbin2 xyz = # of bins in each dim
3509 binsize2 and bininv2 xyz = size of bins in each dim
3510 xyz blo2 = origin of bins
3511 ------------------------------------------------------------------------- */
3512
setup_search_bins()3513 void FixSRD::setup_search_bins()
3514 {
3515 // subboxlo/hi = real space bbox which
3516 // owned/ghost big particles or walls can be in
3517 // start with bounding box for my sub-domain, add dist_ghost
3518 // for triclinic, need to:
3519 // a) convert dist_ghost to lamda space via length012
3520 // b) lo/hi = sub-domain big particle bbox in lamda space
3521 // c) convert lo/hi to real space bounding box via domain->bbox()
3522 // similar to neighbor::setup_bins() and comm::cutghost[] calculation
3523
3524 double subboxlo[3], subboxhi[3];
3525
3526 if (triclinic == 0) {
3527 subboxlo[0] = domain->sublo[0] - dist_ghost;
3528 subboxlo[1] = domain->sublo[1] - dist_ghost;
3529 subboxlo[2] = domain->sublo[2] - dist_ghost;
3530 subboxhi[0] = domain->subhi[0] + dist_ghost;
3531 subboxhi[1] = domain->subhi[1] + dist_ghost;
3532 subboxhi[2] = domain->subhi[2] + dist_ghost;
3533 } else {
3534 double *h_inv = domain->h_inv;
3535 double length0, length1, length2;
3536 length0 = sqrt(h_inv[0] * h_inv[0] + h_inv[5] * h_inv[5] + h_inv[4] * h_inv[4]);
3537 length1 = sqrt(h_inv[1] * h_inv[1] + h_inv[3] * h_inv[3]);
3538 length2 = h_inv[2];
3539 double lo[3], hi[3];
3540 lo[0] = domain->sublo_lamda[0] - dist_ghost * length0;
3541 lo[1] = domain->sublo_lamda[1] - dist_ghost * length1;
3542 lo[2] = domain->sublo_lamda[2] - dist_ghost * length2;
3543 hi[0] = domain->subhi_lamda[0] + dist_ghost * length0;
3544 hi[1] = domain->subhi_lamda[1] + dist_ghost * length1;
3545 hi[2] = domain->subhi_lamda[2] + dist_ghost * length2;
3546 domain->bbox(lo, hi, subboxlo, subboxhi);
3547 }
3548
3549 // require integer # of bins for that volume
3550
3551 nbin2x = static_cast<int>((subboxhi[0] - subboxlo[0]) / gridsearch);
3552 nbin2y = static_cast<int>((subboxhi[1] - subboxlo[1]) / gridsearch);
3553 nbin2z = static_cast<int>((subboxhi[2] - subboxlo[2]) / gridsearch);
3554 if (dimension == 2) nbin2z = 1;
3555
3556 if (nbin2x == 0) nbin2x = 1;
3557 if (nbin2y == 0) nbin2y = 1;
3558 if (nbin2z == 0) nbin2z = 1;
3559
3560 binsize2x = (subboxhi[0] - subboxlo[0]) / nbin2x;
3561 binsize2y = (subboxhi[1] - subboxlo[1]) / nbin2y;
3562 binsize2z = (subboxhi[2] - subboxlo[2]) / nbin2z;
3563 bininv2x = 1.0 / binsize2x;
3564 bininv2y = 1.0 / binsize2y;
3565 bininv2z = 1.0 / binsize2z;
3566
3567 // add bins on either end due to extent of big particles
3568 // radmax = max distance from central bin that biggest particle overlaps
3569 // includes skin movement
3570 // nx,ny,nz = max # of bins to search away from central bin
3571
3572 double radmax = 0.5 * maxbigdiam + 0.5 * neighbor->skin;
3573
3574 int nx = static_cast<int>(radmax / binsize2x) + 1;
3575 int ny = static_cast<int>(radmax / binsize2y) + 1;
3576 int nz = static_cast<int>(radmax / binsize2z) + 1;
3577 if (dimension == 2) nz = 0;
3578
3579 nbin2x += 2 * nx;
3580 nbin2y += 2 * ny;
3581 nbin2z += 2 * nz;
3582
3583 xblo2 = subboxlo[0] - nx * binsize2x;
3584 yblo2 = subboxlo[1] - ny * binsize2y;
3585 zblo2 = subboxlo[2] - nz * binsize2z;
3586 if (dimension == 2) zblo2 = domain->boxlo[2];
3587
3588 // allocate bins
3589 // first deallocate previous bins if necessary
3590
3591 nbins2 = nbin2x * nbin2y * nbin2z;
3592 if (nbins2 > maxbin2) {
3593 memory->destroy(nbinbig);
3594 memory->destroy(binbig);
3595 maxbin2 = nbins2;
3596 memory->create(nbinbig, nbins2, "fix/srd:nbinbig");
3597 memory->create(binbig, nbins2, ATOMPERBIN, "fix/srd:binbig");
3598 }
3599 }
3600
3601 /* ----------------------------------------------------------------------
3602 compute stencil of bin offsets for a big particle overlapping search bins
3603 ------------------------------------------------------------------------- */
3604
setup_search_stencil()3605 void FixSRD::setup_search_stencil()
3606 {
3607 // radmax = max distance from central bin that any big particle overlaps
3608 // includes skin movement
3609 // nx,ny,nz = max # of bins to search away from central bin
3610
3611 double radmax = 0.5 * maxbigdiam + 0.5 * neighbor->skin;
3612 double radsq = radmax * radmax;
3613
3614 int nx = static_cast<int>(radmax / binsize2x) + 1;
3615 int ny = static_cast<int>(radmax / binsize2y) + 1;
3616 int nz = static_cast<int>(radmax / binsize2z) + 1;
3617 if (dimension == 2) nz = 0;
3618
3619 // allocate stencil array
3620 // deallocate previous stencil if necessary
3621
3622 int max = (2 * nx + 1) * (2 * ny + 1) * (2 * nz + 1);
3623 if (max > maxstencil) {
3624 memory->destroy(stencil);
3625 maxstencil = max;
3626 memory->create(stencil, max, 4, "fix/srd:stencil");
3627 }
3628
3629 // loop over all bins
3630 // add bin to stencil:
3631 // if distance of closest corners of bin and central bin is within cutoff
3632
3633 nstencil = 0;
3634 for (int k = -nz; k <= nz; k++)
3635 for (int j = -ny; j <= ny; j++)
3636 for (int i = -nx; i <= nx; i++)
3637 if (bin_bin_distance(i, j, k) < radsq) {
3638 stencil[nstencil][0] = i;
3639 stencil[nstencil][1] = j;
3640 stencil[nstencil][2] = k;
3641 stencil[nstencil][3] = k * nbin2y * nbin2x + j * nbin2x + i;
3642 nstencil++;
3643 }
3644 }
3645
3646 /* ----------------------------------------------------------------------
3647 compute closest squared distance between point x and bin ibin
3648 ------------------------------------------------------------------------- */
3649
point_bin_distance(double * x,int i,int j,int k)3650 double FixSRD::point_bin_distance(double *x, int i, int j, int k)
3651 {
3652 double delx, dely, delz;
3653
3654 double xlo = xblo2 + i * binsize2x;
3655 double xhi = xlo + binsize2x;
3656 double ylo = yblo2 + j * binsize2y;
3657 double yhi = ylo + binsize2y;
3658 double zlo = zblo2 + k * binsize2z;
3659 double zhi = zlo + binsize2z;
3660
3661 if (x[0] < xlo)
3662 delx = xlo - x[0];
3663 else if (x[0] > xhi)
3664 delx = x[0] - xhi;
3665 else
3666 delx = 0.0;
3667
3668 if (x[1] < ylo)
3669 dely = ylo - x[1];
3670 else if (x[1] > yhi)
3671 dely = x[1] - yhi;
3672 else
3673 dely = 0.0;
3674
3675 if (x[2] < zlo)
3676 delz = zlo - x[2];
3677 else if (x[2] > zhi)
3678 delz = x[2] - zhi;
3679 else
3680 delz = 0.0;
3681
3682 return (delx * delx + dely * dely + delz * delz);
3683 }
3684
3685 /* ----------------------------------------------------------------------
3686 compute closest squared distance between 2 bins
3687 central bin (0,0,0) and bin (i,j,k)
3688 ------------------------------------------------------------------------- */
3689
bin_bin_distance(int i,int j,int k)3690 double FixSRD::bin_bin_distance(int i, int j, int k)
3691 {
3692 double delx, dely, delz;
3693
3694 if (i > 0)
3695 delx = (i - 1) * binsize2x;
3696 else if (i == 0)
3697 delx = 0.0;
3698 else
3699 delx = (i + 1) * binsize2x;
3700
3701 if (j > 0)
3702 dely = (j - 1) * binsize2y;
3703 else if (j == 0)
3704 dely = 0.0;
3705 else
3706 dely = (j + 1) * binsize2y;
3707
3708 if (k > 0)
3709 delz = (k - 1) * binsize2z;
3710 else if (k == 0)
3711 delz = 0.0;
3712 else
3713 delz = (k + 1) * binsize2z;
3714
3715 return (delx * delx + dely * dely + delz * delz);
3716 }
3717
3718 /* ---------------------------------------------------------------------- */
3719
pack_reverse_comm(int n,int first,double * buf)3720 int FixSRD::pack_reverse_comm(int n, int first, double *buf)
3721 {
3722 int i, m, last;
3723
3724 m = 0;
3725 last = first + n;
3726
3727 if (torqueflag) {
3728 for (i = first; i < last; i++) {
3729 buf[m++] = flocal[i][0];
3730 buf[m++] = flocal[i][1];
3731 buf[m++] = flocal[i][2];
3732 buf[m++] = tlocal[i][0];
3733 buf[m++] = tlocal[i][1];
3734 buf[m++] = tlocal[i][2];
3735 }
3736
3737 } else {
3738 for (i = first; i < last; i++) {
3739 buf[m++] = flocal[i][0];
3740 buf[m++] = flocal[i][1];
3741 buf[m++] = flocal[i][2];
3742 }
3743 }
3744
3745 return m;
3746 }
3747
3748 /* ---------------------------------------------------------------------- */
3749
unpack_reverse_comm(int n,int * list,double * buf)3750 void FixSRD::unpack_reverse_comm(int n, int *list, double *buf)
3751 {
3752 int i, j, m;
3753
3754 m = 0;
3755
3756 if (torqueflag) {
3757 for (i = 0; i < n; i++) {
3758 j = list[i];
3759 flocal[j][0] += buf[m++];
3760 flocal[j][1] += buf[m++];
3761 flocal[j][2] += buf[m++];
3762 tlocal[j][0] += buf[m++];
3763 tlocal[j][1] += buf[m++];
3764 tlocal[j][2] += buf[m++];
3765 }
3766
3767 } else {
3768 for (i = 0; i < n; i++) {
3769 j = list[i];
3770 flocal[j][0] += buf[m++];
3771 flocal[j][1] += buf[m++];
3772 flocal[j][2] += buf[m++];
3773 }
3774 }
3775 }
3776
3777 /* ----------------------------------------------------------------------
3778 SRD stats
3779 ------------------------------------------------------------------------- */
3780
compute_vector(int n)3781 double FixSRD::compute_vector(int n)
3782 {
3783 // only sum across procs one time
3784
3785 if (stats_flag == 0) {
3786 stats[0] = ncheck;
3787 stats[1] = ncollide;
3788 stats[2] = nbounce;
3789 stats[3] = ninside;
3790 stats[4] = nrescale;
3791 stats[5] = nbins2;
3792 stats[6] = nbins1;
3793 stats[7] = srd_bin_count;
3794 stats[8] = srd_bin_temp;
3795 stats[9] = bouncemaxnum;
3796 stats[10] = bouncemax;
3797 stats[11] = reneighcount;
3798
3799 MPI_Allreduce(stats, stats_all, 10, MPI_DOUBLE, MPI_SUM, world);
3800 MPI_Allreduce(&stats[10], &stats_all[10], 1, MPI_DOUBLE, MPI_MAX, world);
3801 if (stats_all[7] != 0.0) stats_all[8] /= stats_all[7];
3802 stats_all[6] /= nprocs;
3803
3804 stats_flag = 1;
3805 }
3806
3807 return stats_all[n];
3808 }
3809
3810 /* ---------------------------------------------------------------------- */
3811
velocity_stats(int groupnum)3812 void FixSRD::velocity_stats(int groupnum)
3813 {
3814 int bitmask = group->bitmask[groupnum];
3815
3816 double **v = atom->v;
3817 int *mask = atom->mask;
3818 int nlocal = atom->nlocal;
3819
3820 double vone;
3821 double vave = 0.0;
3822 double vmax = 0.0;
3823
3824 for (int i = 0; i < nlocal; i++)
3825 if (mask[i] & bitmask) {
3826 vone = sqrt(v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2]);
3827 vave += vone;
3828 if (vone > vmax) vmax = vone;
3829 }
3830
3831 double all;
3832 MPI_Allreduce(&vave, &all, 1, MPI_DOUBLE, MPI_SUM, world);
3833 double count = group->count(groupnum);
3834 if (count != 0.0)
3835 vave = all / count;
3836 else
3837 vave = 0.0;
3838
3839 MPI_Allreduce(&vmax, &all, 1, MPI_DOUBLE, MPI_MAX, world);
3840 vmax = all;
3841
3842 if (me == 0)
3843 utils::logmesg(lmp, " ave/max {} velocity = {:.8} {:.8}\n", group->names[groupnum], vave,
3844 vmax);
3845 }
3846
3847 /* ---------------------------------------------------------------------- */
3848
newton_raphson(double t1,double t2)3849 double FixSRD::newton_raphson(double t1, double t2)
3850 {
3851 double f1, df, tlo, thi;
3852 lineside(t1, f1, df);
3853 if (f1 < 0.0) {
3854 tlo = t1;
3855 thi = t2;
3856 } else {
3857 thi = t1;
3858 tlo = t2;
3859 }
3860
3861 double f;
3862 double t = 0.5 * (t1 + t2);
3863 double dtold = fabs(t2 - t1);
3864 double dt = dtold;
3865 lineside(t, f, df);
3866
3867 double temp;
3868 for (int i = 0; i < MAXITER; i++) {
3869 if ((((t - thi) * df - f) * ((t - tlo) * df - f) > 0.0) || (fabs(2.0 * f) > fabs(dtold * df))) {
3870 dtold = dt;
3871 dt = 0.5 * (thi - tlo);
3872 t = tlo + dt;
3873 if (tlo == t) return t;
3874 } else {
3875 dtold = dt;
3876 dt = f / df;
3877 temp = t;
3878 t -= dt;
3879 if (temp == t) return t;
3880 }
3881 if (fabs(dt) < TOLERANCE) return t;
3882 lineside(t, f, df);
3883 if (f < 0.0)
3884 tlo = t;
3885 else
3886 thi = t;
3887 }
3888
3889 return t;
3890 }
3891
3892 /* ---------------------------------------------------------------------- */
3893
lineside(double t,double & f,double & df)3894 void FixSRD::lineside(double t, double &f, double &df)
3895 {
3896 double p[2], c[2];
3897
3898 p[0] = xs0[0] + (xs1[0] - xs0[0]) * t;
3899 p[1] = xs0[1] + (xs1[1] - xs0[1]) * t;
3900 c[0] = xb0[0] + (xb1[0] - xb0[0]) * t;
3901 c[1] = xb0[1] + (xb1[1] - xb0[1]) * t;
3902 double dtheta = theta1 - theta0;
3903 double theta = theta0 + dtheta * t;
3904 double cosT = cos(theta);
3905 double sinT = sin(theta);
3906
3907 f = (p[1] - c[1]) * cosT - (p[0] - c[0]) * sinT;
3908 df = ((xs1[1] - xs0[1]) - (xb1[1] - xb0[1])) * cosT - (p[1] - c[1]) * sinT * dtheta -
3909 ((xs1[0] - xs0[0]) - (xb1[0] - xb0[0])) * sinT - (p[0] - c[0]) * cosT * dtheta;
3910 }
3911
3912 /* ---------------------------------------------------------------------- */
3913
triside(double t,double & f,double & df)3914 void FixSRD::triside(double t, double &f, double &df)
3915 {
3916 double p[2], c[2];
3917
3918 p[0] = xs0[0] + (xs1[0] - xs0[0]) * t;
3919 p[1] = xs0[1] + (xs1[1] - xs0[1]) * t;
3920 c[0] = xb0[0] + (xb1[0] - xb0[0]) * t;
3921 c[1] = xb0[1] + (xb1[1] - xb0[1]) * t;
3922 double dtheta = theta1 - theta0;
3923 double theta = theta0 + dtheta * t;
3924 double cosT = cos(theta);
3925 double sinT = sin(theta);
3926
3927 f = (p[1] - c[1]) * cosT - (p[0] - c[0]) * sinT;
3928 df = ((xs1[1] - xs0[1]) - (xb1[1] - xb0[1])) * cosT - (p[1] - c[1]) * sinT * dtheta -
3929 ((xs1[0] - xs0[0]) - (xb1[0] - xb0[0])) * sinT - (p[0] - c[0]) * cosT * dtheta;
3930 }
3931
3932 /* ---------------------------------------------------------------------- */
3933
memory_usage()3934 double FixSRD::memory_usage()
3935 {
3936 double bytes = 0.0;
3937 bytes += (double) (shifts[0].nbins + shifts[1].nbins) * sizeof(BinAve);
3938 bytes += (double) nmax * sizeof(int);
3939 if (bigexist) {
3940 bytes += (double) nbins2 * sizeof(int);
3941 bytes += (double) nbins2 * ATOMPERBIN * sizeof(int);
3942 }
3943 bytes += (double) nmax * sizeof(int);
3944 return bytes;
3945 }
3946
3947 /* ----------------------------------------------------------------------
3948 useful debugging functions
3949 ------------------------------------------------------------------------- */
3950
distance(int i,int j)3951 double FixSRD::distance(int i, int j)
3952 {
3953 double dx = atom->x[i][0] - atom->x[j][0];
3954 double dy = atom->x[i][1] - atom->x[j][1];
3955 double dz = atom->x[i][2] - atom->x[j][2];
3956 return sqrt(dx * dx + dy * dy + dz * dz);
3957 }
3958
3959 /* ---------------------------------------------------------------------- */
3960
print_collision(int i,int j,int ibounce,double t_remain,double dt,double * xscoll,double * xbcoll,double * norm,int type)3961 void FixSRD::print_collision(int i, int j, int ibounce, double t_remain, double dt, double *xscoll,
3962 double *xbcoll, double *norm, int type)
3963 {
3964 double xsstart[3], xbstart[3];
3965 double **x = atom->x;
3966 double **v = atom->v;
3967
3968 if (type != WALL) {
3969 printf("COLLISION between SRD " TAGINT_FORMAT " and BIG " TAGINT_FORMAT "\n", atom->tag[i],
3970 atom->tag[j]);
3971 printf(" bounce # = %d\n", ibounce + 1);
3972 printf(" local indices: %d %d\n", i, j);
3973 printf(" timestep = %g\n", dt);
3974 printf(" time remaining post-collision = %g\n", t_remain);
3975
3976 xsstart[0] = x[i][0] - dt * v[i][0];
3977 xsstart[1] = x[i][1] - dt * v[i][1];
3978 xsstart[2] = x[i][2] - dt * v[i][2];
3979 xbstart[0] = x[j][0] - dt * v[j][0];
3980 xbstart[1] = x[j][1] - dt * v[j][1];
3981 xbstart[2] = x[j][2] - dt * v[j][2];
3982
3983 printf(" SRD start position = %g %g %g\n", xsstart[0], xsstart[1], xsstart[2]);
3984 printf(" BIG start position = %g %g %g\n", xbstart[0], xbstart[1], xbstart[2]);
3985 printf(" SRD coll position = %g %g %g\n", xscoll[0], xscoll[1], xscoll[2]);
3986 printf(" BIG coll position = %g %g %g\n", xbcoll[0], xbcoll[1], xbcoll[2]);
3987 printf(" SRD end position = %g %g %g\n", x[i][0], x[i][1], x[i][2]);
3988 printf(" BIG end position = %g %g %g\n", x[j][0], x[j][1], x[j][2]);
3989
3990 printf(" SRD vel = %g %g %g\n", v[i][0], v[i][1], v[i][2]);
3991 printf(" BIG vel = %g %g %g\n", v[j][0], v[j][1], v[j][2]);
3992 printf(" surf norm = %g %g %g\n", norm[0], norm[1], norm[2]);
3993
3994 double rstart = sqrt((xsstart[0] - xbstart[0]) * (xsstart[0] - xbstart[0]) +
3995 (xsstart[1] - xbstart[1]) * (xsstart[1] - xbstart[1]) +
3996 (xsstart[2] - xbstart[2]) * (xsstart[2] - xbstart[2]));
3997 double rcoll = sqrt((xscoll[0] - xbcoll[0]) * (xscoll[0] - xbcoll[0]) +
3998 (xscoll[1] - xbcoll[1]) * (xscoll[1] - xbcoll[1]) +
3999 (xscoll[2] - xbcoll[2]) * (xscoll[2] - xbcoll[2]));
4000 double rend =
4001 sqrt((x[i][0] - x[j][0]) * (x[i][0] - x[j][0]) + (x[i][1] - x[j][1]) * (x[i][1] - x[j][1]) +
4002 (x[i][2] - x[j][2]) * (x[i][2] - x[j][2]));
4003
4004 printf(" separation at start = %g\n", rstart);
4005 printf(" separation at coll = %g\n", rcoll);
4006 printf(" separation at end = %g\n", rend);
4007
4008 } else {
4009 int dim = wallwhich[j] / 2;
4010
4011 printf("COLLISION between SRD " TAGINT_FORMAT " and WALL %d\n", atom->tag[i], j);
4012 printf(" bounce # = %d\n", ibounce + 1);
4013 printf(" local indices: %d %d\n", i, j);
4014 printf(" timestep = %g\n", dt);
4015 printf(" time remaining post-collision = %g\n", t_remain);
4016
4017 xsstart[0] = x[i][0] - dt * v[i][0];
4018 xsstart[1] = x[i][1] - dt * v[i][1];
4019 xsstart[2] = x[i][2] - dt * v[i][2];
4020 xbstart[0] = xbstart[1] = xbstart[2] = 0.0;
4021 xbstart[dim] = xwall[j] - dt * vwall[j];
4022
4023 printf(" SRD start position = %g %g %g\n", xsstart[0], xsstart[1], xsstart[2]);
4024 printf(" WALL start position = %g\n", xbstart[dim]);
4025 printf(" SRD coll position = %g %g %g\n", xscoll[0], xscoll[1], xscoll[2]);
4026 printf(" WALL coll position = %g\n", xbcoll[dim]);
4027 printf(" SRD end position = %g %g %g\n", x[i][0], x[i][1], x[i][2]);
4028 printf(" WALL end position = %g\n", xwall[j]);
4029
4030 printf(" SRD vel = %g %g %g\n", v[i][0], v[i][1], v[i][2]);
4031 printf(" WALL vel = %g\n", vwall[j]);
4032 printf(" surf norm = %g %g %g\n", norm[0], norm[1], norm[2]);
4033
4034 double rstart = xsstart[dim] - xbstart[dim];
4035 double rcoll = xscoll[dim] - xbcoll[dim];
4036 double rend = x[dim][0] - xwall[j];
4037
4038 printf(" separation at start = %g\n", rstart);
4039 printf(" separation at coll = %g\n", rcoll);
4040 printf(" separation at end = %g\n", rend);
4041 }
4042 }
4043