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 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <cmath>
47 #include <stdlib.h>
48 #include <string.h>
49 #include "fix_dt_reset.h"
50 #include "atom.h"
51 #include "update.h"
52 #include "integrate.h"
53 #include "domain.h"
54 #include "lattice.h"
55 #include "force.h"
56 #include "pair.h"
57 #include "modify.h"
58 #include "fix.h"
59 #include "output.h"
60 #include "dump.h"
61 #include "comm.h"
62 #include "error.h"
63
64 using namespace LAMMPS_NS;
65 using namespace FixConst;
66
67 #define BIG 1.0e20
68
69 /* ---------------------------------------------------------------------- */
70
FixDtReset(LAMMPS * lmp,int narg,char ** arg)71 FixDtReset::FixDtReset(LAMMPS *lmp, int narg, char **arg) :
72 Fix(lmp, narg, arg)
73 {
74 if (narg < 7) error->all(FLERR,"Illegal fix dt/reset command");
75
76 // set time_depend, else elapsed time accumulation can be messed up
77
78 time_depend = 1;
79 scalar_flag = 1;
80 global_freq = 1;
81 extscalar = 0;
82 extvector = 0;
83
84 nevery = force->inumeric(FLERR,arg[3]);
85 if (nevery <= 0) error->all(FLERR,"Illegal fix dt/reset command");
86
87 minbound = maxbound = 1;
88 tmin = tmax = 0.0;
89 if (strcmp(arg[4],"NULL") == 0) minbound = 0;
90 else tmin = force->numeric(FLERR,arg[4]);
91 if (strcmp(arg[5],"NULL") == 0) maxbound = 0;
92 else tmax = force->numeric(FLERR,arg[5]);
93 xmax = force->numeric(FLERR,arg[6]);
94
95 if (minbound && tmin < 0.0) error->all(FLERR,"Illegal fix dt/reset command");
96 if (maxbound && tmax < 0.0) error->all(FLERR,"Illegal fix dt/reset command");
97 if (minbound && maxbound && tmin >= tmax)
98 error->all(FLERR,"Illegal fix dt/reset command");
99 if (xmax <= 0.0) error->all(FLERR,"Illegal fix dt/reset command");
100
101 int scaleflag = 0;
102
103 int iarg = 7;
104 while (iarg < narg) {
105 if (strcmp(arg[iarg],"units") == 0) {
106 if (iarg+2 > narg) error->all(FLERR,"Illegal fix dt/reset command");
107 if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
108 else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
109 else error->all(FLERR,"Illegal fix dt/reset command");
110 iarg += 2;
111 } else error->all(FLERR,"Illegal fix dt/reset command");
112 }
113
114 // setup scaling, based on xlattice parameter
115
116 if (scaleflag) xmax *= domain->lattice->xlattice;
117
118 // initializations
119
120 t_laststep = 0.0;
121 laststep = update->ntimestep;
122 }
123
124 /* ---------------------------------------------------------------------- */
125
setmask()126 int FixDtReset::setmask()
127 {
128 int mask = 0;
129 mask |= END_OF_STEP;
130 return mask;
131 }
132
133 /* ---------------------------------------------------------------------- */
134
init()135 void FixDtReset::init()
136 {
137 // set rRESPA flag
138
139 respaflag = 0;
140 if (strstr(update->integrate_style,"respa")) respaflag = 1;
141
142 // check for DCD or XTC dumps
143
144 for (int i = 0; i < output->ndump; i++)
145 if ((strcmp(output->dump[i]->style,"dcd") == 0 ||
146 strcmp(output->dump[i]->style,"xtc") == 0) && comm->me == 0)
147 error->warning(FLERR,
148 "Dump dcd/xtc timestamp may be wrong with fix dt/reset");
149
150 ftm2v = force->ftm2v;
151 dt = update->dt;
152 }
153
154 /* ---------------------------------------------------------------------- */
155
setup(int vflag)156 void FixDtReset::setup(int vflag)
157 {
158 end_of_step();
159 }
160
161 /* ---------------------------------------------------------------------- */
162
end_of_step()163 void FixDtReset::end_of_step()
164 {
165 double dtv,dtf,dtsq;
166 double vsq,fsq,massinv;
167 double delx,dely,delz,delr;
168
169 // compute vmax and amax of any atom in group
170
171 double **v = atom->v;
172 double **f = atom->f;
173 double *mass = atom->mass;
174 double *rmass = atom->rmass;
175 int *type = atom->type;
176 int *mask = atom->mask;
177 int nlocal = atom->nlocal;
178
179 double dtmin = BIG;
180
181 for (int i = 0; i < nlocal; i++)
182 if (mask[i] & groupbit) {
183 if (rmass) massinv = 1.0/rmass[i];
184 else massinv = 1.0/mass[type[i]];
185 vsq = v[i][0]*v[i][0] + v[i][1]*v[i][1] + v[i][2]*v[i][2];
186 fsq = f[i][0]*f[i][0] + f[i][1]*f[i][1] + f[i][2]*f[i][2];
187 dtv = dtf = BIG;
188 if (vsq > 0.0) dtv = xmax/sqrt(vsq);
189 if (fsq > 0.0) dtf = sqrt(2.0*xmax/(ftm2v*sqrt(fsq)*massinv));
190 dt = MIN(dtv,dtf);
191 dtsq = dt*dt;
192 delx = dt*v[i][0] + 0.5*dtsq*massinv*f[i][0] * ftm2v;
193 dely = dt*v[i][1] + 0.5*dtsq*massinv*f[i][1] * ftm2v;
194 delz = dt*v[i][2] + 0.5*dtsq*massinv*f[i][2] * ftm2v;
195 delr = sqrt(delx*delx + dely*dely + delz*delz);
196 if (delr > xmax) dt *= xmax/delr;
197 dtmin = MIN(dtmin,dt);
198 }
199
200 MPI_Allreduce(&dtmin,&dt,1,MPI_DOUBLE,MPI_MIN,world);
201
202 if (minbound) dt = MAX(dt,tmin);
203 if (maxbound) dt = MIN(dt,tmax);
204
205 // if timestep didn't change, just return
206 // else reset update->dt and other classes that depend on it
207 // rRESPA, pair style, fixes
208
209 if (dt == update->dt) return;
210
211 laststep = update->ntimestep;
212
213 update->update_time();
214 update->dt = dt;
215 if (respaflag) update->integrate->reset_dt();
216 if (force->pair) force->pair->reset_dt();
217 for (int i = 0; i < modify->nfix; i++) modify->fix[i]->reset_dt();
218 }
219
220 /* ---------------------------------------------------------------------- */
221
compute_scalar()222 double FixDtReset::compute_scalar()
223 {
224 return (double) laststep;
225 }
226