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 <string.h>
43 #include <stdlib.h>
44 #include <cmath>
45 #include <algorithm>
46 #include "atom.h"
47 #include "update.h"
48 #include "error.h"
49 #include "fix_check_timestep_gran.h"
50 #include "pair_gran.h"
51 #include "properties.h"
52 #include "fix_property_global.h"
53 #include "force.h"
54 #include "comm.h"
55 #include "modify.h"
56 #include "fix_wall_gran.h"
57 #include "fix_mesh_surface.h"
58 #include "neighbor.h"
59 #include "mpi_liggghts.h"
60 #include "property_registry.h"
61 #include "global_properties.h"
62
63 using namespace LAMMPS_NS;
64 using namespace FixConst;
65 using namespace MODEL_PARAMS;
66
67 /* ---------------------------------------------------------------------- */
68
FixCheckTimestepGran(LAMMPS * lmp,int narg,char ** arg)69 FixCheckTimestepGran::FixCheckTimestepGran(LAMMPS *lmp, int narg, char **arg) :
70 Fix(lmp, narg, arg)
71 {
72 warnflag = true;
73 errorflag = false;
74 vmax_user = 0.;
75
76 if (narg < 6) error->all(FLERR,"Illegal fix check/timestep/gran command, not enough arguments");
77
78 int iarg = 6;
79
80 if(0 == strcmp("check_every_time",arg[3]))
81 {
82 nevery = atoi(arg[4]);
83 fraction_rayleigh_lim = atof(arg[5]);
84 fraction_hertz_lim = atof(arg[6]);
85 iarg = 7;
86 }
87 else
88 {
89 nevery = atoi(arg[3]);
90 fraction_rayleigh_lim = atof(arg[4]);
91 fraction_hertz_lim = atof(arg[5]);
92 iarg = 6;
93 }
94
95 while(iarg < narg)
96 {
97 if (strcmp(arg[iarg],"warn") == 0) {
98 if (narg < iarg+2) error->fix_error(FLERR,this,"not enough arguments for 'warn'");
99 if(0 == strcmp(arg[iarg+1],"no"))
100 warnflag = false;
101 else if(0 == strcmp(arg[iarg+1],"yes"))
102 warnflag = true;
103 else
104 error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'warn'");
105 iarg += 2;
106 } else if (strcmp(arg[iarg],"error") == 0) {
107 if (narg < iarg+2) error->fix_error(FLERR,this,"not enough arguments for 'error'");
108 if(0 == strcmp(arg[iarg+1],"no"))
109 errorflag = false;
110 else if(0 == strcmp(arg[iarg+1],"yes"))
111 errorflag = true;
112 else
113 error->fix_error(FLERR,this,"expecting 'yes' or 'no' after 'error'");
114 iarg += 2;
115 } else if (strcmp(arg[iarg],"vmax") == 0) {
116 if (narg < iarg+2) error->fix_error(FLERR,this,"not enough arguments for 'vmax'");
117 vmax_user = force->numeric(FLERR,arg[iarg+1]);
118 iarg += 2;
119 } else if(strcmp(style,"mesh/surface") == 0) {
120 char *errmsg = new char[strlen(arg[iarg])+50];
121 sprintf(errmsg,"unknown keyword or wrong keyword order: %s", arg[iarg]);
122 error->fix_error(FLERR,this,errmsg);
123 delete []errmsg;
124 }
125 }
126
127 vector_flag = 1;
128 size_vector = 3;
129 global_freq = nevery;
130 extvector = 1;
131
132 fraction_rayleigh = fraction_hertz = fraction_skin = 0.;
133 Yeff = NULL;
134 }
135
136 /* ---------------------------------------------------------------------- */
137
setmask()138 int FixCheckTimestepGran::setmask()
139 {
140 int mask = 0;
141 mask |= END_OF_STEP;
142 return mask;
143 }
144
145 /* ---------------------------------------------------------------------- */
146
init()147 void FixCheckTimestepGran::init()
148 {
149 //some error checks
150 if(!atom->radius_flag || !atom->density_flag)
151 error->all(FLERR,"Fix check/timestep/gran can only be used together with atom style sphere");
152
153 pg = (PairGran*)force->pair_match("gran",1);
154 if(!pg) pg = (PairGran*)force->pair_match("gran/omp",1);
155
156 if (!pg)
157 error->all(FLERR,"Fix check/timestep/gran can only be used together with: gran");
158
159 properties = atom->get_properties();
160 int max_type = properties->max_type();
161
162 fwg = NULL;
163 for (int i = 0; i < modify->n_fixes_style("wall/gran"); i++)
164 if(static_cast<FixWallGran*>(modify->find_fix_style("wall/gran",i))->is_mesh_wall())
165 fwg = static_cast<FixWallGran*>(modify->find_fix_style("wall/gran",i));
166
167 Y = static_cast<FixPropertyGlobal*>(modify->find_fix_property("youngsModulus","property/global","peratomtype",max_type,0,style));
168 nu = static_cast<FixPropertyGlobal*>(modify->find_fix_property("poissonsRatio","property/global","peratomtype",max_type,0,style));
169
170 if(!Y || !nu)
171 error->all(FLERR,"Fix check/timestep/gran only works with a pair style that defines youngsModulus and poissonsRatio");
172
173 force->registry.registerProperty("Yeff", &MODEL_PARAMS::createYeff);
174 force->registry.connect("Yeff", Yeff,this->style);
175 }
176
177 /* ---------------------------------------------------------------------- */
178
end_of_step()179 void FixCheckTimestepGran::end_of_step()
180 {
181 calc_rayleigh_hertz_estims();
182
183 double skin = neighbor->skin;
184 double dt = update->dt;
185
186 fraction_rayleigh = dt/rayleigh_time;
187 fraction_hertz = dt/hertz_time;
188 fraction_skin = (v_rel_max_simulation * dt) / neighbor->skin;
189
190 if(errorflag || (warnflag && comm->me==0))
191 {
192 char errstr[512];
193
194 if(fraction_rayleigh > fraction_rayleigh_lim)
195 {
196 sprintf(errstr,"time-step is %f %% of rayleigh time",fraction_rayleigh*100.);
197 if(errorflag)
198 error->fix_error(FLERR,this,errstr);
199 else
200 error->warning(FLERR,errstr);
201 }
202 if(fraction_hertz > fraction_hertz_lim)
203 {
204 sprintf(errstr,"time-step is %f %% of hertz time",fraction_hertz*100.);
205 if(errorflag)
206 error->fix_error(FLERR,this,errstr);
207 else
208 error->warning(FLERR,errstr);
209 }
210
211 if(fraction_skin > 1)
212 {
213 sprintf(errstr,"time step too large or skin too small - particles may relatively travel a distance of %f per time-step, but skin is %f",v_rel_max_simulation*dt,skin);
214 if(errorflag)
215 error->fix_error(FLERR,this,errstr);
216 else
217 error->warning(FLERR,errstr);
218 }
219
220 if(v_rel_max_simulation * dt > r_min)
221 {
222 sprintf(errstr,"time step way too large - particles relativley move further than the minimum radius in one step");
223 if(errorflag)
224 error->fix_error(FLERR,this,errstr);
225 else
226 error->warning(FLERR,errstr);
227 }
228 }
229 }
230
231 /* ---------------------------------------------------------------------- */
232
calc_rayleigh_hertz_estims()233 void FixCheckTimestepGran::calc_rayleigh_hertz_estims()
234 {
235 double **v = atom->v;
236 double *density = atom->density;
237 double *r = atom->radius;
238 int *type = atom->type;
239 int *mask = atom->mask;
240 int nlocal = atom->nlocal;
241
242 int max_type = properties->max_type();
243
244 //check rayleigh time and vmax of particles
245 rayleigh_time = BIG;
246 r_min = BIG;
247 double vmax_sqr = 0;
248 double vmag_sqr;
249 double rayleigh_time_i;
250
251 for (int i = 0; i < nlocal; i++)
252 {
253 if (mask[i] & groupbit)
254 {
255 double rad = r[i];
256 #ifdef SUPERQUADRIC_ACTIVE_FLAG
257 if(atom->superquadric_flag)
258 rad=std::min(std::min(atom->shape[i][0],atom->shape[i][1]),atom->shape[i][2]);
259 #endif
260
261 double shear_mod = Y->get_values()[type[i]-1]/(2.*(nu->get_values()[type[i]-1]+1.));
262 rayleigh_time_i = M_PI*rad*sqrt(density[i]/shear_mod)/(0.1631*nu->get_values()[type[i]-1]+0.8766);
263 if(rayleigh_time_i < rayleigh_time) rayleigh_time = rayleigh_time_i;
264
265 vmag_sqr = vectorMag3DSquared(v[i]);
266 if(vmag_sqr > vmax_sqr)
267 vmax_sqr = vmag_sqr;
268
269 if(rad < r_min) r_min = rad;
270 }
271 }
272
273 MPI_Min_Scalar(r_min,world);
274 MPI_Max_Scalar(vmax_sqr,world);
275 MPI_Min_Scalar(rayleigh_time,world);
276
277 //choose vmax_user if larger than value in simulation
278 if(vmax_user*vmax_user > vmax_sqr)
279 vmax_sqr = vmax_user*vmax_user;
280
281 // get vmax of geometry
282 FixMeshSurface ** mesh_list;
283 TriMesh * mesh;
284 double *v_node;
285 double vmax_sqr_mesh=0.;
286 double vmag_sqr_mesh;
287
288 if(fwg)
289 {
290 mesh_list = fwg->mesh_list();
291 for(int imesh = 0; imesh < fwg->n_meshes(); imesh++)
292 {
293 mesh = (mesh_list[imesh])->triMesh();
294 if(mesh->isMoving())
295 {
296 // check if perElementProperty 'v' exists
297 if (mesh->prop().getElementPropertyIndex("v") == -1)
298 error->one(FLERR,"Internal error - mesh has no perElementProperty 'v' \n");
299 // loop local elements only
300 int sizeMesh = mesh->sizeLocal();
301 for(int itri = 0; itri < sizeMesh; itri++)
302 {
303 for(int inode = 0; inode < 3; inode++)
304 {
305 v_node = mesh->prop().getElementProperty<MultiVectorContainer<double,3,3> >("v")->begin()[itri][inode];
306 vmag_sqr_mesh = vectorMag3DSquared(v_node);
307 if(vmag_sqr_mesh > vmax_sqr_mesh)
308 vmax_sqr_mesh = vmag_sqr_mesh;
309 }
310 }
311 }
312 }
313 }
314
315 MPI_Max_Scalar(vmax_sqr_mesh,world);
316
317 // decide vmax - either particle-particle or particle-mesh contact
318 v_rel_max_simulation = std::max(2.*sqrt(vmax_sqr),sqrt(vmax_sqr) + sqrt(vmax_sqr_mesh));
319
320 // check estimation for hertz time
321 // this is not exact...
322 // loop over all material comibinations
323 // loop all particles
324 // test collision of particle with itself
325 double hertz_time_min = 1000000.;
326 double hertz_time_i,meff,reff;
327
328 if ( !MathExtraLiggghts::compDouble(v_rel_max_simulation,0.0) ) {
329 for(int ti = 1; ti < max_type+1; ti++)
330 {
331 for(int tj = ti; tj < max_type+1; tj++)
332 {
333 const double Eeff = Yeff[ti][tj];
334
335 for(int i = 0; i < nlocal; i++)
336 {
337 if (mask[i] & groupbit)
338 {
339 if(type[i]!=ti || type[i]!=tj) continue;
340 meff = 4.*r[i]*r[i]*r[i]*M_PI/3.*density[i];
341 reff = r[i]/2.;
342 #ifdef SUPERQUADRIC_ACTIVE_FLAG
343 if(atom->superquadric_flag) {
344 meff = atom->volume[i]*density[i];
345 reff = std::max(std::max(atom->shape[i][0],atom->shape[i][1]),atom->shape[i][2])/2.0;
346 }
347 #endif
348 hertz_time_i = 2.87*pow(meff*meff/(reff*Eeff*Eeff*v_rel_max_simulation),0.2);
349 if(hertz_time_i<hertz_time_min)
350 hertz_time_min=hertz_time_i;
351 }
352 }
353 }
354 }
355 }
356
357 MPI_Min_Scalar(hertz_time_min,world);
358 hertz_time = hertz_time_min;
359 }
360
361 /* ----------------------------------------------------------------------
362 return fractions of rayleigh/hertz time-step
363 ------------------------------------------------------------------------- */
364
compute_vector(int n)365 double FixCheckTimestepGran::compute_vector(int n)
366 {
367 if(n == 0) return fraction_rayleigh;
368 else if(n == 1) return fraction_hertz;
369 else if(n == 2) return fraction_skin;
370 return 0.;
371 }
372