1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #include <cmath>
43 #include <stdlib.h>
44 #include <string.h>
45 #include "fix_property_atom.h"
46 #include "atom.h"
47 #include "memory.h"
48 #include "error.h"
49
50 #include "pair_gran.h"
51 #include "atom_vec.h"
52 #include "force.h"
53 #include "update.h"
54 #include "comm.h"
55 #include "modify.h"
56 #include "group.h"
57 #include "timer.h"
58 #include "neighbor.h"
59
60 #include "mpi_liggghts.h"
61
62 using namespace LAMMPS_NS;
63 using namespace FixConst;
64
65 #define EPSILON 0.001
66
67 /* ---------------------------------------------------------------------- */
68
FixPropertyAtom(LAMMPS * lmp,int narg,char ** arg,bool parse)69 FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg, bool parse) :
70 Fix(lmp, narg, arg),
71 propertyname(0),
72 property(0),
73 internal(false)
74 {
75
76 if(parse) parse_args(narg,arg);
77 }
78
parse_args(int narg,char ** arg)79 void FixPropertyAtom::parse_args(int narg, char **arg)
80 {
81 // Check args
82 if (narg < 9) error->all(FLERR,"Illegal fix property/atom command, not enough arguments");
83 if (narg > 29) error->warning(FLERR,"Vector length in fix property/atom larger than 20. Are you sure you want that?");
84
85 // Read args
86
87 int n = strlen(arg[3]) + 1;
88 variablename = new char[n];
89 strcpy(variablename,arg[3]);
90
91 bool vector_with_one_entry = false;
92 if (strcmp(arg[4],"scalar") == 0) data_style = FIXPROPERTY_ATOM_SCALAR;
93 else if (strcmp(arg[4],"vector") == 0) data_style = FIXPROPERTY_ATOM_VECTOR;
94 // This vector style allows for a vector to have only one entry. Under normal circumstances the scalar style should be chosen instead.
95 else if (strcmp(arg[4],"vector_one_entry") == 0)
96 {
97 vector_with_one_entry = true;
98 data_style = FIXPROPERTY_ATOM_VECTOR;
99 }
100 else error->all(FLERR,"Unknown style for fix property/atom. Valid styles are 'scalar', 'vector' or 'vector_one_entry'");
101
102 if (strcmp(arg[5],"yes") == 0)
103 {
104 restart_peratom = 1;
105 restart_global = 1;
106 }
107 else if (strcmp(arg[5],"no") == 0)
108 {
109 restart_peratom = 0;
110 restart_global = 0;
111 }
112 else error->all(FLERR,"Unknown restart style for fix property/atom. Valid styles are 'yes' or 'no'");
113
114 if (strcmp(arg[6],"yes") == 0) commGhost = 1;
115 else if (strcmp(arg[6],"no") == 0) commGhost = 0;
116 else error->all(FLERR,"Unknown communicate_ghost style for fix property/atom. Valid styles are 'yes' or 'no'");
117
118 if (strcmp(arg[7],"yes") == 0) commGhostRev = 1;
119 else if (strcmp(arg[7],"no") == 0) commGhostRev = 0;
120 else error->all(FLERR,"Unknown communicate_reverse_ghost style for fix property/atom. Valid styles are 'yes' or 'no'");
121
122 nvalues = narg - 8;
123
124 if ((nvalues == 1) && !(data_style == FIXPROPERTY_ATOM_SCALAR || (vector_with_one_entry && data_style == FIXPROPERTY_ATOM_VECTOR)))
125 error->all(FLERR,"Error in fix property/atom: Number of default values provided not consistent with vector style. Provide more than 1 value or use style 'scalar'");
126
127 if ((nvalues >1) && (data_style != FIXPROPERTY_ATOM_VECTOR))
128 error->all(FLERR,"Error in fix property/atom: Number of default values provided not consistent with scalar style. Provide 1 value or use style 'vector'");
129
130 defaultvalues = new double[nvalues];
131
132 // fix handles properties that need to be initialized at particle creation
133 create_attribute = 1;
134
135 // special case: only implemented for scalar default value from scalar
136 // initialize property from another existing scalar, which is found via Property class
137 // since this fix must exist already, it is initialized before this fix property
138 propertyname = 0;
139 if(FIXPROPERTY_ATOM_SCALAR == data_style)
140 {
141 char *prop = arg[8];
142 int n = strlen(prop);
143 bool is_digit = false;
144 for (int i = 0; i < n; i++)
145 if (isdigit(prop[i])) is_digit = true;
146
147 if(!is_digit)
148 {
149 // scalar property found
150 int len1,len2;
151 if(atom->get_properties()->find_property(prop,"scalar-atom",len1,len2))
152 {
153 propertyname = new char[n+1];
154 strcpy(propertyname,prop);
155 }
156 }
157 }
158
159 // set default values
160 if(!propertyname)
161 {
162 for (int j = 0; j < nvalues; j++)
163 {
164 // if any of the values is none, this fix will not init properties
165 if(strcmp(arg[8+j],"none") == 0)
166 {
167 create_attribute = 0;
168 continue;
169 }
170 defaultvalues[j] = force->numeric(FLERR,arg[8+j]);
171 }
172 }
173
174 if (data_style) size_peratom_cols = nvalues;
175 else size_peratom_cols = 0;
176
177 peratom_flag=1;
178 peratom_freq=1;
179 extvector=0;
180
181 if (commGhost) comm_forward = nvalues;
182 if (commGhostRev) comm_reverse = nvalues;
183
184 // perform initial allocation of atom-based array
185 // register with Atom class
186 vector_atom = NULL; array_atom = NULL;
187 grow_arrays(atom->nmax);
188 atom->add_callback(0);
189 if (restart_peratom) atom->add_callback(1);
190
191 // init all arrays since dump may access it on timestep 0
192 // or a variable may access it before first run
193
194 int nlocal = atom->nlocal;
195 if(create_attribute)
196 {
197 if(propertyname)
198 {
199 // not implemented
200 if(data_style)
201 error->all(FLERR,"internal error");
202
203 pre_set_arrays();
204 for (int i = 0; i < nlocal; i++)
205 vector_atom[i] = property[i];
206 }
207 else
208 {
209 for (int i = 0; i < nlocal; i++)
210 {
211 if (data_style)
212 {
213 for (int m = 0; m < nvalues; m++)
214 array_atom[i][m] = defaultvalues[m];
215 }
216 else vector_atom[i] = defaultvalues[0];
217 }
218 }
219 }
220
221 // check if there is already a fix that tries to register a property with the same name
222
223 for (int ifix = 0; ifix < modify->nfix; ifix++)
224 if ((modify->fix[ifix]) && (strcmp(modify->fix[ifix]->style,style) == 0) && (strcmp(((FixPropertyAtom*)(modify->fix[ifix]))->variablename,variablename)==0) )
225 error->fix_error(FLERR,this,"there is already a fix that registers a variable of the same name");
226
227 // flags for vector output
228 //vector_flag = 1;
229 size_vector = nvalues;
230 global_freq = 1;
231 extvector = 1;
232 }
233
234 /* ---------------------------------------------------------------------- */
235
~FixPropertyAtom()236 FixPropertyAtom::~FixPropertyAtom()
237 {
238 // unregister callbacks to this fix from Atom class
239 atom->delete_callback(id,0);
240 if (restart_peratom) atom->delete_callback(id,1);
241
242 // delete locally stored arrays
243 delete []variablename;
244 delete []defaultvalues;
245 if(propertyname) delete []propertyname;
246
247 if (data_style) memory->destroy(array_atom);
248 else memory->destroy(vector_atom);
249 }
250
251 /* ---------------------------------------------------------------------- */
252
check_fix(const char * varname,const char * svmstyle,int len1,int len2,const char * caller,bool errflag)253 Fix* FixPropertyAtom::check_fix(const char *varname,const char *svmstyle,int len1,int len2,const char *caller,bool errflag)
254 {
255 char errmsg[400];
256
257 if(strcmp(varname,variablename) == 0)
258 {
259 if(strcmp(svmstyle,"scalar") == 0) len1 = 1;
260
261 // check variable style
262 if(
263 (strcmp(svmstyle,"scalar") == 0 && data_style != FIXPROPERTY_ATOM_SCALAR) ||
264 (strcmp(svmstyle,"vector") == 0 && data_style != FIXPROPERTY_ATOM_VECTOR) ||
265 (strcmp(svmstyle,"vector2D") == 0 && data_style != FIXPROPERTY_ATOM_VECTOR2D) ||
266 (strcmp(svmstyle,"quaternion") == 0 && data_style != FIXPROPERTY_ATOM_QUATERNION)
267 )
268 {
269 if(errflag)
270 {
271 sprintf(errmsg,"%s style required for fix property/atom variable %s for usage with caller %s",
272 svmstyle,varname,caller);
273 error->all(FLERR,errmsg);
274 }
275 else return NULL;
276 }
277
278 // check length
279 if(len1 > nvalues)
280 {
281 if(errflag)
282 {
283 sprintf(errmsg,"Fix property/atom variable %s has wrong length (length is %d but length %d expected) for usage with caller %s",
284 varname,nvalues,len1,caller);
285 error->all(FLERR,errmsg);
286 }
287 else return NULL;
288 }
289
290 // success
291 return static_cast<Fix*>(this);
292 }
293 return NULL;
294 }
295
296 /* ---------------------------------------------------------------------- */
297
setmask()298 int FixPropertyAtom::setmask()
299 {
300 int mask = 0;
301 mask |= PRE_EXCHANGE;
302 return mask;
303 }
304
305 /* ----------------------------------------------------------------------
306 forward and backward comm to be used by other fixes as needed
307 ------------------------------------------------------------------------- */
308
do_forward_comm()309 void FixPropertyAtom::do_forward_comm()
310 {
311 timer->stamp();
312 if (commGhost) comm->forward_comm_fix(this);
313 else error->all(FLERR,"FixPropertyAtom: Faulty implementation - forward_comm invoked, but not registered");
314 timer->stamp(TIME_COMM);
315 }
316
do_reverse_comm()317 void FixPropertyAtom::do_reverse_comm()
318 {
319 timer->stamp();
320 if (commGhostRev) comm->reverse_comm_fix(this);
321 else error->all(FLERR,"FixPropertyAtom: Faulty implementation - reverse_comm invoked, but not registered");
322 timer->stamp(TIME_COMM);
323 }
324
325 /* ----------------------------------------------------------------------
326 memory usage of local atom-based arrays
327 ------------------------------------------------------------------------- */
328
memory_usage()329 double FixPropertyAtom::memory_usage()
330 {
331 int nmax = atom->nmax;
332 double bytes = nmax * nvalues * sizeof(double);
333 return bytes;
334 }
335
336 /* ----------------------------------------------------------------------
337 allocate local atom-based arrays
338 ------------------------------------------------------------------------- */
339
grow_arrays(int nmax)340 void FixPropertyAtom::grow_arrays(int nmax)
341 {
342 if (data_style) memory->grow(array_atom,nmax,nvalues,"FixPropertyAtom:array_atom");
343 else memory->grow(vector_atom, nmax, "FixPropertyAtom:vector_atom");
344
345 }
346
347 /* ----------------------------------------------------------------------
348 copy values within local atom-based arrays
349 ------------------------------------------------------------------------- */
350
copy_arrays(int i,int j,int delflag)351 void FixPropertyAtom::copy_arrays(int i, int j, int delflag)
352 {
353 if (data_style) for(int k=0;k<nvalues;k++) array_atom[j][k] = array_atom[i][k];
354 else vector_atom[j]=vector_atom[i];
355 }
356
357 /* ----------------------------------------------------------------------
358 called before set_arrays is called for each atom
359 ------------------------------------------------------------------------- */
360
pre_set_arrays()361 void FixPropertyAtom::pre_set_arrays()
362 {
363
364 property = 0;
365 if(propertyname)
366 {
367
368 int len1,len2;
369
370 property = (double*) atom->get_properties()->find_property(propertyname,"scalar-atom",len1,len2);
371 if(!property)
372 {
373 char errstr[200];
374 sprintf(errstr,"Property %s not found",propertyname);
375
376 error->fix_error(FLERR,this,errstr);
377 }
378
379 }
380 }
381
382 /* ----------------------------------------------------------------------
383 initialize one atom's array values, called when atom is created
384 ------------------------------------------------------------------------- */
385
set_arrays(int i)386 void FixPropertyAtom::set_arrays(int i)
387 {
388
389 if (data_style)
390 for(int k=0;k<nvalues;k++)
391 array_atom[i][k] = defaultvalues[k];
392 else vector_atom[i] = property ? property[i] : defaultvalues[0];
393 }
394
395 /* ----------------------------------------------------------------------
396 set all atoms values
397 ------------------------------------------------------------------------- */
398
set_all(double value,bool ghost)399 void FixPropertyAtom::set_all(double value, bool ghost)
400 {
401
402 int nall;
403
404 if(!ghost)
405 nall = atom->nlocal;
406 else
407 nall = atom->nlocal + atom->nghost;
408 if (data_style)
409 {
410 for(int i = 0; i < nall; i++)
411 {
412 for(int k=0;k<nvalues;k++)
413 array_atom[i][k] = value;
414 }
415 }
416 else
417 {
418 for(int i = 0; i < nall; i++)
419 vector_atom[i] = value;
420 }
421 }
422
423 /* ----------------------------------------------------------------------
424 pack values in local atom-based arrays for exchange with another proc
425 ------------------------------------------------------------------------- */
426
pack_exchange(int i,double * buf)427 int FixPropertyAtom::pack_exchange(int i, double *buf)
428 {
429 if (data_style) for(int k=0;k<nvalues;k++) buf[k] = array_atom[i][k];
430 else buf[0] = vector_atom[i];
431 return nvalues;
432 }
433
434 /* ----------------------------------------------------------------------
435 unpack values into local atom-based arrays after exchange
436 ------------------------------------------------------------------------- */
437
unpack_exchange(int nlocal,double * buf)438 int FixPropertyAtom::unpack_exchange(int nlocal, double *buf)
439 {
440 if (data_style) for(int k=0;k<nvalues;k++) array_atom[nlocal][k] = buf[k];
441 else vector_atom[nlocal]=buf[0];
442 return nvalues;
443 }
444
445 /* ----------------------------------------------------------------------
446 pack values in local atom-based arrays for restart file
447 ------------------------------------------------------------------------- */
448
pack_restart(int i,double * buf)449 int FixPropertyAtom::pack_restart(int i, double *buf)
450 {
451 buf[0] = static_cast<double>(nvalues+1);
452 if (data_style) for(int k=0;k<nvalues;k++) buf[k+1] = array_atom[i][k];
453 else buf[1] = vector_atom[i];
454
455 return (nvalues+1);
456 }
457
458 /* ----------------------------------------------------------------------
459 unpack values from atom->extra array to restart the fix
460 ------------------------------------------------------------------------- */
461
unpack_restart(int nlocal,int nth)462 void FixPropertyAtom::unpack_restart(int nlocal, int nth)
463 {
464 double **extra = atom->extra;
465
466 // skip to Nth set of extra values
467
468 int m = 0;
469 for (int i = 0; i < nth; i++) m += static_cast<int> (extra[nlocal][m]);
470 m++;
471
472 if (data_style) for(int k=0;k<nvalues;k++) array_atom[nlocal][k] = extra[nlocal][m++];
473 else vector_atom[nlocal] = extra[nlocal][m++];
474 }
475
476 /* ----------------------------------------------------------------------
477 maxsize of any atom's restart data
478 ------------------------------------------------------------------------- */
479
maxsize_restart()480 int FixPropertyAtom::maxsize_restart()
481 {
482 return nvalues+1;
483 }
484
485 /* ----------------------------------------------------------------------
486 size of atom nlocal's restart data
487 ------------------------------------------------------------------------- */
488
size_restart(int nlocal)489 int FixPropertyAtom::size_restart(int nlocal)
490 {
491 return nvalues+1;
492 }
493
494 /* ---------------------------------------------------------------------- */
495
496 /* ---------------------------------------------------------------------- */
497
pack_comm(int n,int * list,double * buf,int pbc_flag,int * pbc)498 int FixPropertyAtom::pack_comm(int n, int *list, double *buf,
499 int pbc_flag, int *pbc)
500 {
501 int i,j;
502 //we dont need to account for pbc here
503 int m = 0;
504 for (i = 0; i < n; i++) {
505 j = list[i];
506 if (data_style) for(int k=0;k<nvalues;k++) buf[m++] = array_atom[j][k];
507 else buf[m++] = vector_atom[j];
508 }
509 return nvalues;
510 }
511
512 /* ---------------------------------------------------------------------- */
513
unpack_comm(int n,int first,double * buf)514 void FixPropertyAtom::unpack_comm(int n, int first, double *buf)
515 {
516 int i,m,last;
517 m = 0;
518 last = first + n;
519 for (i = first; i < last; i++) {
520 if (data_style) for(int k=0;k<nvalues;k++) array_atom[i][k]=buf[m++];
521 else vector_atom[i]=buf[m++];
522 }
523
524 }
525
526 /* ---------------------------------------------------------------------- */
527
pack_reverse_comm(int n,int first,double * buf)528 int FixPropertyAtom::pack_reverse_comm(int n, int first, double *buf)
529 {
530 int i,m,last;
531 m = 0;
532 last = first + n;
533 for (i = first; i < last; i++) {
534 if (data_style) for(int k=0;k<nvalues;k++) buf[m++] = array_atom[i][k];
535 else buf[m++] = vector_atom[i];
536 }
537 return nvalues;
538 }
539
540 /* ---------------------------------------------------------------------- */
541
unpack_reverse_comm(int n,int * list,double * buf)542 void FixPropertyAtom::unpack_reverse_comm(int n, int *list, double *buf)
543 {
544 int i,j,m = 0;
545 for (i = 0; i < n; i++) {
546 j = list[i];
547 if (data_style) for(int k=0;k<nvalues;k++) array_atom[j][k]+=buf[m++];
548 else vector_atom[j]+=buf[m++];
549 }
550 }
551
552 /* ---------------------------------------------------------------------- */
553
write_restart(FILE * fp)554 void FixPropertyAtom::write_restart(FILE *fp)
555 {
556 int n = 0;
557 double list[1];
558 list[n++] = static_cast<double>(nvalues);
559
560 if (comm->me == 0) {
561 int size = n * sizeof(double);
562 fwrite(&size,sizeof(int),1,fp);
563 fwrite(list,sizeof(double),n,fp);
564 }
565 }
566
567 /* ----------------------------------------------------------------------
568 use state info from restart file to restart the Fix
569 ------------------------------------------------------------------------- */
570
restart(char * buf)571 void FixPropertyAtom::restart(char *buf)
572 {
573 int n = 0;
574 double *list = (double *) buf;
575 int nvalues_re;
576
577 nvalues_re = static_cast<int> (list[n++]);
578
579 if(nvalues_re != nvalues)
580 error->fix_error(FLERR,this,"restarted fix has incompatible data size");
581 }
582
583 /* ----------------------------------------------------------------------
584 return components of property sum on fix group, n = 0..nvalues-1
585 ------------------------------------------------------------------------- */
586
compute_vector(int n)587 double FixPropertyAtom::compute_vector(int n)
588 {
589 int nlocal = atom->nlocal;
590 int *mask = atom->mask;
591
592 double value = 0.;
593
594 for (int i = 0; i < nlocal; i++)
595 {
596 if (mask[i] & groupbit)
597 {
598 if (data_style) value += array_atom[i][n];
599 else value += vector_atom[i];
600 }
601 }
602
603 MPI_Sum_Scalar(value,world);
604 return value;
605 }
606