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_continuity.h"
45 #include "update.h"
46 #include "atom.h"
47 #include "force.h"
48 #include "modify.h"
49 #include "comm.h"
50 #include "neighbor.h"
51 #include "neigh_list.h"
52 #include "neigh_request.h"
53 #include "memory.h"
54 #include "error.h"
55 #include "sph_kernels.h"
56 #include "fix_property_atom.h"
57 #include "timer.h"
58
59 using namespace LAMMPS_NS;
60 using namespace FixConst;
61
62 /* ---------------------------------------------------------------------- */
63
FixSphDensityContinuity(LAMMPS * lmp,int narg,char ** arg)64 FixSphDensityContinuity::FixSphDensityContinuity(LAMMPS *lmp, int narg, char **arg) :
65 FixSph(lmp, narg, arg)
66 {
67 int iarg = 0;
68
69 if (iarg+3 > narg) error->all(FLERR,"Illegal fix sph/density/continuity command");
70
71 iarg += 3;
72
73 while (iarg < narg) {
74 // kernel style
75 if (strcmp(arg[iarg],"sphkernel") == 0) {
76 if (iarg+2 > narg) error->all(FLERR,"Illegal fix sph/density/continuity command");
77
78 if(kernel_style) delete []kernel_style;
79 kernel_style = new char[strlen(arg[iarg+1])+1];
80 strcpy(kernel_style,arg[iarg+1]);
81
82 // check uniqueness of kernel IDs
83
84 int flag = SPH_KERNEL_NS::sph_kernels_unique_id();
85 if(flag < 0) error->all(FLERR,"Cannot proceed, sph kernels need unique IDs, check all sph_kernel_* files");
86
87 // get kernel id
88
89 kernel_id = SPH_KERNEL_NS::sph_kernel_id(kernel_style);
90 if(kernel_id < 0) error->all(FLERR,"Illegal fix sph/density/continuity command, unknown sph kernel");
91
92 iarg += 2;
93
94 } else error->all(FLERR,"Illegal fix sph/density/continuity command");
95 }
96
97 time_depend = 0; // only time step is used, but not a relative time
98
99 }
100
101 /* ---------------------------------------------------------------------- */
102
~FixSphDensityContinuity()103 FixSphDensityContinuity::~FixSphDensityContinuity()
104 {
105
106 }
107
108 /* ---------------------------------------------------------------------- */
109
setmask()110 int FixSphDensityContinuity::setmask()
111 {
112 int mask = 0;
113 mask |= PRE_FORCE;
114 return mask;
115 }
116
117 /* ---------------------------------------------------------------------- */
118
init()119 void FixSphDensityContinuity::init()
120 {
121 FixSph::init();
122
123 // check if there is an nve/sph fix present
124
125 int idx_integ = -1;
126 for(int i = 0; i < modify->nfix; i++)
127 {
128 if(strncmp("nve/sph",modify->fix[i]->style,7) == 0) {
129 idx_integ = i;
130 break;
131 }
132 if(strncmp("nve/xsph",modify->fix[i]->style,8) == 0) {
133 idx_integ = i;
134 break;
135 }
136 }
137
138 if(idx_integ == -1) error->fix_error(FLERR,this,"Requires to define a fix nve/sph also \n");
139 }
140
141 /* ---------------------------------------------------------------------- */
142
pre_force(int vflag)143 void FixSphDensityContinuity::pre_force(int vflag)
144 {
145 //template function for using per atom or per atomtype smoothing length
146 if (mass_type) pre_force_eval<1>(vflag);
147 else pre_force_eval<0>(vflag);
148 }
149
150 /* ---------------------------------------------------------------------- */
151
152 template <int MASSFLAG>
pre_force_eval(int vflag)153 void FixSphDensityContinuity::pre_force_eval(int vflag)
154 {
155 int i,j,ii,jj,inum,jnum,itype,jtype;
156 double xtmp,ytmp,ztmp,delx,dely,delz,rsq,r,rinv,s,gradWmag;
157 int *ilist,*jlist,*numneigh,**firstneigh;
158 double sli,slj,slCom,slComInv,cut,delVDotDelR,imass,jmass;
159
160 double **x = atom->x;
161 double **v = atom->vest;
162 double *drho = atom->drho;
163 int *mask = atom->mask;
164 int nlocal = atom->nlocal;
165
166 int *type = atom->type; // if MASSFLAG
167 double *mass = atom->mass; // if MASSFLAG
168 double *rmass = atom->rmass; // if !MASSFLAG
169
170 int newton_pair = force->newton_pair;
171
172 if (igroup == atom->firstgroup) nlocal = atom->nfirst;
173
174 // need updated ghost positions and self contributions
175 timer->stamp();
176
177 if (!MASSFLAG) fppaSl->do_forward_comm();
178 timer->stamp(TIME_COMM);
179
180 if (!MASSFLAG) updatePtrs(); // get sl
181
182 // loop over neighbors of my atoms
183 inum = list->inum;
184 ilist = list->ilist;
185 numneigh = list->numneigh;
186 firstneigh = list->firstneigh;
187
188 for (ii = 0; ii < inum; ii++) {
189 i = ilist[ii];
190 if (!(mask[i] & groupbit)) continue;
191 xtmp = x[i][0];
192 ytmp = x[i][1];
193 ztmp = x[i][2];
194 jlist = firstneigh[i];
195 jnum = numneigh[i];
196
197 if (MASSFLAG) {
198 itype = type[i];
199 imass = mass[itype];
200 } else {
201 imass = rmass[i];
202 sli = sl[i];
203 }
204
205 for (jj = 0; jj < jnum; jj++) {
206 j = jlist[jj];
207 if (!(mask[j] & groupbit)) continue;
208
209 if (MASSFLAG) {
210 jtype = type[j];
211 jmass = mass[jtype];
212 slCom = slComType[itype][jtype];
213 } else {
214 jmass = rmass[j];
215 slj = sl[j];
216 slCom = interpDist(sli,slj);
217 }
218
219 cut = slCom*kernel_cut;//SPH_KERNEL_NS::sph_kernel_cut(kernel_id);
220
221 delx = xtmp - x[j][0];
222 dely = ytmp - x[j][1];
223 delz = ztmp - x[j][2];
224 rsq = delx*delx + dely*dely + delz*delz;
225
226 // TODO: Cutsq from pair?
227 if (rsq >= cut*cut) continue;
228
229 // calculate normalized distance
230
231 r = sqrt(rsq);
232 if (r == 0.) {
233 fprintf(screen,"Particle %i and %i are at same position (%f, %f, %f)\n",i,j,xtmp,ytmp,ztmp);
234 error->one(FLERR,"Zero distance between SPH particles!");
235 }
236 rinv = 1./r;
237 slComInv = 1./slCom;
238 s = r * slComInv;
239
240 // scalar product of delV and delR/R
241 delVDotDelR = rinv * ( delx*(v[i][0]-v[j][0]) + dely*(v[i][1]-v[j][1]) + delz*(v[i][2]-v[j][2]) );
242
243 // calculate value for magnitude of grad W
244 gradWmag = SPH_KERNEL_NS::sph_kernel_der(kernel_id,s,slCom,slComInv);
245
246 // add contribution of neighbor
247 // have a half neigh list, so do it for both if necessary
248
249 drho[i] += jmass*gradWmag*delVDotDelR;
250
251 if (newton_pair || j < nlocal) {
252 drho[j] += imass*gradWmag*delVDotDelR;
253 }
254 }
255 }
256
257 }
258