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 no contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2013- DCS Computing GmbH, Linz
39 ------------------------------------------------------------------------- */
40
41 #include <cmath>
42 #include <stdlib.h>
43 #include <string.h>
44 #include "atom.h"
45 #include "force.h"
46 #include "update.h"
47 #include "modify.h"
48 #include "memory.h"
49 #include "error.h"
50 #include "group.h"
51 #include "region.h"
52 #include "domain.h"
53 #include "mpi_liggghts.h"
54 #include "math_extra_liggghts.h"
55 #include "fix_property_atom_region_tracer.h"
56
57 using namespace LAMMPS_NS;
58 using namespace FixConst;
59
60 /* ---------------------------------------------------------------------- */
61
FixPropertyAtomRegionTracer(LAMMPS * lmp,int narg,char ** arg,bool parse)62 FixPropertyAtomRegionTracer::FixPropertyAtomRegionTracer(LAMMPS *lmp, int narg, char **arg,bool parse) :
63 FixPropertyAtom(lmp, narg, arg, false),
64 iarg_(3),
65 check_every_(10)
66 {
67
68 if (strcmp(style, "property/atom/timetracer") == 0)
69 error->fix_error(FLERR, this, "Style of fix property/atom/timetracer is now property/atom/regiontracer/time");
70 // do the derived class stuff
71
72 // parse args
73
74 bool hasargs = true;
75 while(iarg_ < narg && hasargs)
76 {
77 hasargs = false;
78
79 if(strcmp(arg[iarg_],"add_region") == 0) {
80 if(narg < iarg_+2)
81 error->fix_error(FLERR,this,"not enough arguments for 'add_region'");
82 iarg_++;
83 int n = strlen(arg[iarg_]) + 1;
84 char *idreg = new char[n];
85 strcpy(idreg,arg[iarg_]);
86 int ireg = domain->find_region(arg[iarg_++]);
87 if (ireg == -1)
88 error->fix_error(FLERR,this,"Region ID does not exist");
89 iregion_.push_back(ireg);
90 idregion_.push_back(idreg);
91 hasargs = true;
92 } else if(strcmp(arg[iarg_],"check_region_every") == 0) {
93 if(narg < iarg_+2)
94 error->fix_error(FLERR,this,"not enough arguments for 'check_region_every'");
95 iarg_++;
96 check_every_ = atoi(arg[iarg_]);
97 if(check_every_ < 0)
98 error->fix_error(FLERR,this,"check_region_every > 0 required");
99 iarg_++;
100 hasargs = true;
101 } else if(strcmp(style,"property/atom/regiontracer/time") == 0 )
102 error->fix_error(FLERR,this,"unknown keyword");
103 }
104
105 // do the base class stuff
106 // allocate one extra per-particle for each region
107
108 int n = strlen(id) + 1;
109 char *tracer_name = new char[n];
110 strcpy(tracer_name,id);
111 const int n_reg = iregion_.size();
112 if (sizeof(double) == 8 && n_reg > 53)
113 error->fix_error(FLERR, this, "Only 53 regions are allowed in a fix property/atom/region/tracer (on systems with 64 bit doubles)");
114 else if (sizeof(double) == 4 && n_reg > 24)
115 error->fix_error(FLERR, this, "Only 24 regions are allowed in a fix property/atom/region/tracer (on systems with 32 bit doubles)");
116
117 const int num_args = 9 + (n_reg > 0 ? n_reg+1 : 0);
118 char **baseargs = new char*[num_args]; // VS does not support variable-length arrays --> new
119 baseargs[0] = tracer_name;
120 baseargs[1] = (char *) "all";
121 baseargs[2] = (char *) "property/atom/tracer";
122 baseargs[3] = tracer_name;
123 if(0 == n_reg)
124 baseargs[4] = (char *) "scalar";
125 else
126 baseargs[4] = (char *) "vector";
127 baseargs[5] = (char *) "yes";
128 baseargs[6] = (char *) "yes";
129 baseargs[7] = (char *) "no";
130 baseargs[8] = (char *) "0."; // Property for simulation domain
131 for(int i = 0; i < n_reg; i++)
132 baseargs[9+i] = (char *) "0."; // Property for region
133 if (n_reg > 0)
134 baseargs[9+n_reg] = (char *) "0."; // Bit map for region
135 parse_args(num_args,(char**)baseargs);
136
137 delete []baseargs;
138
139 // settings
140 nevery = 1;
141
142 vector_flag = 1;
143 size_vector = 1+n_reg;
144 global_freq = check_every_;
145 extvector = 1;
146 }
147
148 /* ---------------------------------------------------------------------- */
149
~FixPropertyAtomRegionTracer()150 FixPropertyAtomRegionTracer::~FixPropertyAtomRegionTracer()
151 {
152 for(size_t i = 0; i < idregion_.size(); i++)
153 delete [] idregion_[i];
154 }
155
156 /* ----------------------------------------------------------------------
157 initialize this fix
158 ------------------------------------------------------------------------- */
159
init()160 void FixPropertyAtomRegionTracer::init()
161 {
162 iregion_.clear();
163 for(size_t i = 0; i < idregion_.size(); i++)
164 {
165 int ireg = domain->find_region(idregion_[i]);
166 if (ireg == -1)
167 error->fix_error(FLERR,this,"Region ID does not exist");
168 iregion_.push_back(ireg);
169 }
170
171 int n_ms = modify->n_fixes_style("multisphere");
172 if(n_ms > 0)
173 error->fix_error(FLERR,this,"may not be used together with fix multisphere");
174 }
175
176 /* ---------------------------------------------------------------------- */
177
setmask()178 int FixPropertyAtomRegionTracer::setmask()
179 {
180 int mask = FixPropertyAtom::setmask();
181 mask |= END_OF_STEP;
182 return mask;
183 }
184
185 /* ----------------------------------------------------------------------
186 add residence time(s)
187 ------------------------------------------------------------------------- */
188
end_of_step()189 void FixPropertyAtomRegionTracer::end_of_step()
190 {
191 const int nlocal = atom->nlocal;
192 double **x = atom->x;
193 const int * const mask = atom->mask;
194
195 // add multiple of dt in case check_region_every > 1
196 const double dt = update->dt;
197 const int n_reg = iregion_.size();
198
199 if (n_reg == 0)
200 {
201 double * const data = this->vector_atom;
202 for(int i = 0; i < nlocal; i++)
203 {
204 if (mask[i] & groupbit)
205 data[i] += trace_multiplier(i)*dt;
206 }
207 }
208 else // fix has region list
209 {
210 double * const * const data = this->array_atom;
211
212 for(int i = 0; i < nlocal; i++)
213 {
214 if (mask[i] & groupbit)
215 {
216 const double add = trace_multiplier(i)*dt;
217 data[i][0] += add;
218
219 long region_map = (long)data[i][n_reg+1];
220 for(int ireg = 0; ireg < n_reg; ireg++)
221 {
222 const long region_flag = 1<<ireg;
223 // check if we need to update the region list
224 if (update->ntimestep % check_every_ == 0)
225 {
226 Region * const region = domain->regions[iregion_[ireg]];
227 if(region->match(x[i][0],x[i][1],x[i][2]))
228 region_map |= region_flag;
229 else
230 region_map &= ~region_flag;
231 }
232 if (region_map & region_flag)
233 data[i][1+ireg] += add;
234 }
235 data[i][n_reg+1] = static_cast<double>(region_map);
236 }
237 }
238 }
239 }
240
241 /* ----------------------------------------------------------------------
242 return average of residence time, n = 0..nvalues-1
243 ------------------------------------------------------------------------- */
244
compute_vector(int n)245 double FixPropertyAtomRegionTracer::compute_vector(int n)
246 {
247 const int nlocal = atom->nlocal;
248 const int * const mask = atom->mask;
249
250 double value = 0.;
251 double norm = 0.;
252
253 for (int i = 0; i < nlocal; i++)
254 {
255 if (mask[i] & groupbit)
256 {
257 const double incr = array_atom ? array_atom[i][n] : vector_atom[i];
258 // sum all values
259 value += incr;
260 // increment the normalizing factor
261 if (incr > 0.0)
262 norm += 1.0;
263 }
264 }
265
266 MPI_Sum_Scalar(value,world);
267 MPI_Sum_Scalar(norm,world);
268
269 if (norm > 0.5)
270 return value/atom->natoms;
271 else
272 return 0.0;
273 }
274