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 Andreas Aigner (JKU Linz)
36
37 Copyright 2009-2012 JKU Linz
38 ------------------------------------------------------------------------- */
39
40 #include <cmath>
41 #include <mpi.h>
42 #include <string.h>
43 #include <stdlib.h>
44 #include "fix_sph_density_corr.h"
45 #include "update.h"
46 #include "respa.h"
47 #include "atom.h"
48 #include "force.h"
49 #include "modify.h"
50 #include "comm.h"
51 #include "neighbor.h"
52 #include "neigh_list.h"
53 #include "neigh_request.h"
54 #include "memory.h"
55 #include "error.h"
56 #include "sph_kernels.h"
57 #include "fix_property_atom.h"
58 #include "timer.h"
59
60 using namespace LAMMPS_NS;
61 using namespace FixConst;
62
63 /* ---------------------------------------------------------------------- */
64
FixSphDensityCorr(LAMMPS * lmp,int narg,char ** arg)65 FixSphDensityCorr::FixSphDensityCorr(LAMMPS *lmp, int narg, char **arg) :
66 FixSph(lmp, narg, arg)
67 {
68 int iarg = 0;
69 if (narg < iarg+4) error->fix_error(FLERR,this,"Illegal fix sph/density/corr command, not enough arguments");
70
71 // correction style
72 iarg += 3;
73
74 if (strcmp(arg[iarg],"shepard") == 0) {
75 if (iarg+3 > narg) error->fix_error(FLERR,this,"Not enough arguments");
76 if (strcmp(arg[iarg+1],"every") == 0) {
77 every = force->inumeric(FLERR,arg[iarg+2]);
78 if (every <= 0) error->fix_error(FLERR,this,"every <= 0 not allowed");
79 corrStyle = CORR_SHEPARD;
80 iarg += 3;
81 } else error->fix_error(FLERR,this,"");
82 } else if (strcmp(arg[iarg],"mls") == 0) {
83 error->fix_error(FLERR,this,"MLS correction is not implemented until now.");
84 corrStyle = CORR_MLS;
85 } else error->fix_error(FLERR,this,"Unknown style for fix sph/density/corr. Valid styles are 'shepard' or 'mls'");
86
87 while (iarg < narg) {
88 // kernel style
89 if (strcmp(arg[iarg],"sphkernel") == 0) {
90 if (iarg+2 > narg) error->fix_error(FLERR,this,"Illegal fix sph/density/continuity command");
91
92 if(kernel_style) delete []kernel_style;
93 kernel_style = new char[strlen(arg[iarg+1])+1];
94 strcpy(kernel_style,arg[iarg+1]);
95
96 // check uniqueness of kernel IDs
97
98 int flag = SPH_KERNEL_NS::sph_kernels_unique_id();
99 if(flag < 0) error->fix_error(FLERR,this,"Cannot proceed, sph kernels need unique IDs, check all sph_kernel_* files");
100
101 // get kernel id
102
103 kernel_id = SPH_KERNEL_NS::sph_kernel_id(kernel_style);
104 if(kernel_id < 0) error->fix_error(FLERR,this,"Illegal fix sph/density/continuity command, unknown sph kernel");
105
106 iarg += 2;
107
108 } else error->fix_error(FLERR,this,"Illegal fix sph/density/continuity command");
109 }
110
111 // init variables and set flags
112
113 quantity_name = new char[strlen("corrKernel")+1];
114 strcpy(quantity_name,"corrKernel");
115
116 fix_quantity = NULL;
117
118 peratom_flag = 1;
119 size_peratom_cols = 0;
120 peratom_freq = 1;
121 time_depend = 0;
122
123 scalar_flag = 1;
124 global_freq = 1;
125
126 ago = 0;
127
128 }
129
130 /* ---------------------------------------------------------------------- */
131
~FixSphDensityCorr()132 FixSphDensityCorr::~FixSphDensityCorr()
133 {
134 delete []quantity_name;
135 }
136
137 /* ---------------------------------------------------------------------- */
138
pre_delete(bool)139 void FixSphDensityCorr::pre_delete(bool)
140 {
141 //unregister property/atom fixes
142 if (fix_quantity) modify->delete_fix(quantity_name);
143 }
144
145 /* ---------------------------------------------------------------------- */
146
setmask()147 int FixSphDensityCorr::setmask()
148 {
149 int mask = 0;
150 mask |= PRE_FORCE;
151 return mask;
152 }
153
154 /* ---------------------------------------------------------------------- */
155
updatePtrs()156 void FixSphDensityCorr::updatePtrs()
157 {
158 FixSph::updatePtrs(); // update sl
159 quantity = fix_quantity->vector_atom;
160
161 vector_atom = quantity;
162 }
163
164 /* ---------------------------------------------------------------------- */
165
post_create()166 void FixSphDensityCorr::post_create()
167 {
168 const char * fixarg[9];
169
170 if (fix_quantity==NULL) {
171 fixarg[0]=quantity_name;
172 fixarg[1]="all";
173 fixarg[2]="property/atom";
174 fixarg[3]=quantity_name;
175 fixarg[4]="scalar";
176 fixarg[5]="yes";
177 fixarg[6]="yes";
178 fixarg[7]="no";
179 fixarg[8]="0.";
180 modify->add_fix(9,const_cast<char**>(fixarg));
181 fix_quantity=static_cast<FixPropertyAtom*>(modify->find_fix_property(quantity_name,"property/atom","scalar",0,0,style));
182 }
183 }
184
185 /* ---------------------------------------------------------------------- */
186
init()187 void FixSphDensityCorr::init()
188 {
189 FixSph::init();
190 }
191
192 /* ---------------------------------------------------------------------- */
193
pre_force(int)194 void FixSphDensityCorr::pre_force(int)
195 {
196 //template function for using per atom or per atomtype smoothing length
197 if (mass_type) pre_force_eval<1>();
198 else pre_force_eval<0>();
199 }
200
201 /* ---------------------------------------------------------------------- */
202
203 template <int MASSFLAG>
pre_force_eval()204 void FixSphDensityCorr::pre_force_eval()
205 {
206 int i,j,ii,jj,inum,jnum,itype,jtype;
207 double xtmp,ytmp,ztmp,delx,dely,delz,rsq,r,s=0.0,W;
208 double sli,sliInv,slj,slCom,slComInv,cut,imass,jmass;
209 int *ilist,*jlist,*numneigh,**firstneigh;
210
211 double **x = atom->x;
212 int *mask = atom->mask;
213 double *rho = atom->rho;
214 int nlocal = atom->nlocal;
215
216 // TODO: Both declaration necessary?
217 int *type = atom->type;
218 double *mass = atom->mass;
219 double *rmass = atom->rmass;
220
221 int newton_pair = force->newton_pair;
222
223 updatePtrs(); // get sl, quantity
224
225 timer->stamp();
226 comm->forward_comm();
227 if (!MASSFLAG) fppaSl->do_forward_comm();
228 timer->stamp(TIME_COMM);
229
230 ago++;
231 if (ago % every == 0) {
232 ago = 0;
233
234 // kernel normalization
235
236 for (i = 0; i < nlocal; i++)
237 {
238 if (mask[i] & groupbit) {
239 if (MASSFLAG) {
240 itype = type[i];
241 sli = sl[itype-1];
242 imass = mass[itype];
243 } else {
244 sli = sl[i];
245 imass = rmass[i];
246 }
247
248 sliInv = 1./sli;
249
250 // this gets a value for W at self, perform error check
251
252 W = SPH_KERNEL_NS::sph_kernel(kernel_id,0.,sli,sliInv);
253 if (W < 0.)
254 {
255 fprintf(screen,"s = %f, W = %f\n",s,W);
256 error->one(FLERR,"Illegal kernel used, W < 0");
257 }
258
259 quantity[i] = W*imass / rho[i];
260 }
261 }
262
263 // loop over neighbors of my atoms
264
265 inum = list->inum;
266 ilist = list->ilist;
267 numneigh = list->numneigh;
268 firstneigh = list->firstneigh;
269
270 for (ii = 0; ii < inum; ii++) {
271 i = ilist[ii];
272 if (!(mask[i] & groupbit)) continue;
273 xtmp = x[i][0];
274 ytmp = x[i][1];
275 ztmp = x[i][2];
276 jlist = firstneigh[i];
277 jnum = numneigh[i];
278
279 if (MASSFLAG) {
280 itype = type[i];
281 imass = mass[itype];
282 } else {
283 sli = sl[i];
284 imass = rmass[i];
285 }
286
287 for (jj = 0; jj < jnum; jj++) {
288 j = jlist[jj];
289 if (!(mask[j] & groupbit)) continue;
290
291 if (MASSFLAG) {
292 jtype = type[j];
293 jmass = mass[jtype];
294 slCom = slComType[itype][jtype];
295 } else {
296 jmass = rmass[j];
297 slj = sl[j];
298 slCom = interpDist(sli,slj);
299 }
300
301 cut = slCom*kernel_cut;
302
303 delx = xtmp - x[j][0];
304 dely = ytmp - x[j][1];
305 delz = ztmp - x[j][2];
306 rsq = delx*delx + dely*dely + delz*delz;
307
308 if (rsq >= cut*cut) continue;
309
310 // calculate distance
311 r = sqrt(rsq);
312 slComInv = 1./slCom;
313 s = r*slComInv;
314
315 // this gets a value for W at self, perform error check
316
317 W = SPH_KERNEL_NS::sph_kernel(kernel_id,s,slCom,slComInv);
318 if (W < 0.)
319 {
320 fprintf(screen,"s = %f, W = %f\n",s,W);
321 error->one(FLERR,"Illegal kernel used, W < 0");
322 }
323
324 // add contribution of neighbor
325 // have a half neigh list, so do it for both if necessary
326 quantity[i] += W*jmass / rho[j];
327 if (newton_pair || j < nlocal)
328 quantity[j] += W*imass / rho[i];
329 }
330 }
331
332 // reset and add rho contribution of self
333
334 for (i = 0; i < nlocal; i++) {
335 if (mask[i] & groupbit) {
336 if (MASSFLAG) {
337 itype = type[i];
338 sli = sl[itype-1];
339 imass = mass[itype];
340 } else {
341 sli = sl[i];
342 imass = rmass[i];
343 }
344
345 sliInv = 1./sli;
346
347 // this gets a value for W at self, perform error check
348
349 W = SPH_KERNEL_NS::sph_kernel(kernel_id,0.,sli,sliInv);
350 if (W < 0.)
351 {
352 fprintf(screen,"s = %f, W = %f\n",s,W);
353 error->one(FLERR,"Illegal kernel used, W < 0");
354 }
355
356 // add contribution of self
357 rho[i] = W*imass;
358 }
359 }
360
361 // need updated ghost positions and self contributions
362 timer->stamp();
363 comm->forward_comm();
364 timer->stamp(TIME_COMM);
365
366 // loop over neighbors of my atoms
367 inum = list->inum;
368 ilist = list->ilist;
369 numneigh = list->numneigh;
370 firstneigh = list->firstneigh;
371
372 for (ii = 0; ii < inum; ii++) {
373 i = ilist[ii];
374 if (!(mask[i] & groupbit)) continue;
375 xtmp = x[i][0];
376 ytmp = x[i][1];
377 ztmp = x[i][2];
378 jlist = firstneigh[i];
379 jnum = numneigh[i];
380
381 if (MASSFLAG) {
382 itype = type[i];
383 imass = mass[itype];
384 } else {
385 imass = rmass[i];
386 sli = sl[i];
387 }
388
389 for (jj = 0; jj < jnum; jj++) {
390 j = jlist[jj];
391 if (!(mask[j] & groupbit)) continue;
392
393 if (MASSFLAG) {
394 jtype = type[j];
395 jmass = mass[jtype];
396 slCom = slComType[itype][jtype];
397 } else {
398 jmass = rmass[j];
399 slj = sl[j];
400 slCom = interpDist(sli,slj);
401 }
402
403 cut = slCom*kernel_cut;
404
405 delx = xtmp - x[j][0];
406 dely = ytmp - x[j][1];
407 delz = ztmp - x[j][2];
408 rsq = delx*delx + dely*dely + delz*delz;
409
410 if (rsq >= cut*cut) continue;
411
412 // calculate distance
413 r = sqrt(rsq);
414 slComInv = 1./slCom;
415 s = r*slComInv;
416
417 // this gets a value for W at self, perform error check
418
419 W = SPH_KERNEL_NS::sph_kernel(kernel_id,s,slCom,slComInv);
420 if (W < 0.)
421 {
422 fprintf(screen,"s = %f, W = %f\n",s,W);
423 error->one(FLERR,"Illegal kernel used, W < 0");
424 }
425
426 // add contribution of neighbor
427 // have a half neigh list, so do it for both if necessary
428 rho[i] += W*jmass;
429 if (newton_pair || j < nlocal)
430 rho[j] += W*imass;
431 }
432 }
433
434 // normalize rho
435 for (i = 0; i < nlocal; i++) {
436 if (mask[i] & groupbit) {
437 rho[i] = rho[i]/quantity[i];
438 }
439 }
440
441 // rho is now correct, send to ghosts
442 timer->stamp();
443 comm->forward_comm();
444 timer->stamp(TIME_COMM);
445 }
446 }
447