1 // clang-format off
2 /* ----------------------------------------------------------------------
3 *
4 * *** Smooth Mach Dynamics ***
5 *
6 * This file is part of the MACHDYN package for LAMMPS.
7 * Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
8 * Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
9 * Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
10 *
11 * ----------------------------------------------------------------------- */
12
13 /* ----------------------------------------------------------------------
14 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
15 https://www.lammps.org/, Sandia National Laboratories
16 Steve Plimpton, sjplimp@sandia.gov
17
18 Copyright (2003) Sandia Corporation. Under the terms of Contract
19 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
20 certain rights in this software. This software is distributed under
21 the GNU General Public License.
22
23 See the README file in the top-level LAMMPS directory.
24 ------------------------------------------------------------------------- */
25
26 #include "fix_smd_integrate_ulsph.h"
27 #include <cmath>
28 #include <cstring>
29 #include <Eigen/Eigen>
30 #include "atom.h"
31 #include "comm.h"
32 #include "force.h"
33 #include "update.h"
34 #include "error.h"
35 #include "pair.h"
36
37 using namespace Eigen;
38 using namespace LAMMPS_NS;
39 using namespace FixConst;
40
41 /* ---------------------------------------------------------------------- */
42
FixSMDIntegrateUlsph(LAMMPS * lmp,int narg,char ** arg)43 FixSMDIntegrateUlsph::FixSMDIntegrateUlsph(LAMMPS *lmp, int narg, char **arg) :
44 Fix(lmp, narg, arg) {
45
46 if ((atom->esph_flag != 1) || (atom->vfrac_flag != 1))
47 error->all(FLERR, "fix smd/integrate_ulsph command requires atom_style with both energy and volume");
48
49 if (narg < 3)
50 error->all(FLERR, "Illegal number of arguments for fix smd/integrate_ulsph command");
51
52 adjust_radius_flag = false;
53 xsphFlag = false;
54 vlimit = -1.0;
55 int iarg = 3;
56
57 if (comm->me == 0) {
58 printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
59 printf("fix smd/integrate_ulsph is active for group: %s \n", arg[1]);
60 }
61 while (true) {
62
63 if (iarg >= narg) {
64 break;
65 }
66
67 if (strcmp(arg[iarg], "xsph") == 0) {
68 xsphFlag = true;
69 if (comm->me == 0) {
70 error->one(FLERR, "XSPH is currently not available");
71 printf("... will use XSPH time integration\n");
72 }
73 } else if (strcmp(arg[iarg], "adjust_radius") == 0) {
74 adjust_radius_flag = true;
75
76 iarg++;
77 if (iarg == narg) {
78 error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
79 }
80
81 adjust_radius_factor = utils::numeric(FLERR, arg[iarg],false,lmp);
82
83 iarg++;
84 if (iarg == narg) {
85 error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
86 }
87
88 min_nn = utils::inumeric(FLERR, arg[iarg],false,lmp);
89
90 iarg++;
91 if (iarg == narg) {
92 error->all(FLERR, "expected three numbers following adjust_radius: factor, min, max");
93 }
94
95 max_nn = utils::inumeric(FLERR, arg[iarg],false,lmp);
96
97 if (comm->me == 0) {
98 printf("... will adjust smoothing length dynamically with factor %g to achieve %d to %d neighbors per particle.\n ",
99 adjust_radius_factor, min_nn, max_nn);
100 }
101
102 } else if (strcmp(arg[iarg], "limit_velocity") == 0) {
103 iarg++;
104 if (iarg == narg) {
105 error->all(FLERR, "expected number following limit_velocity");
106 }
107 vlimit = utils::numeric(FLERR, arg[iarg],false,lmp);
108
109 if (comm->me == 0) {
110 printf("... will limit velocities to <= %g\n", vlimit);
111 }
112 } else {
113 char msg[128];
114 snprintf(msg,128, "Illegal keyword for smd/integrate_ulsph: %s\n", arg[iarg]);
115 error->all(FLERR, msg);
116 }
117
118 iarg++;
119
120 }
121
122 if (comm->me == 0) {
123 printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
124 }
125
126 // set comm sizes needed by this fix
127 atom->add_callback(Atom::GROW);
128
129 time_integrate = 1;
130 }
131
132 /* ---------------------------------------------------------------------- */
133
setmask()134 int FixSMDIntegrateUlsph::setmask() {
135 int mask = 0;
136 mask |= INITIAL_INTEGRATE;
137 mask |= FINAL_INTEGRATE;
138 return mask;
139 }
140
141 /* ---------------------------------------------------------------------- */
142
init()143 void FixSMDIntegrateUlsph::init() {
144 dtv = update->dt;
145 dtf = 0.5 * update->dt * force->ftm2v;
146 vlimitsq = vlimit * vlimit;
147 }
148
149 /* ----------------------------------------------------------------------
150 allow for both per-type and per-atom mass
151 ------------------------------------------------------------------------- */
152
initial_integrate(int)153 void FixSMDIntegrateUlsph::initial_integrate(int /*vflag*/) {
154 double **x = atom->x;
155 double **v = atom->v;
156 double **f = atom->f;
157 double **vest = atom->vest;
158 double *rmass = atom->rmass;
159
160 int *mask = atom->mask;
161 int nlocal = atom->nlocal;
162 int i, itmp;
163 double dtfm, vsq, scale;
164 double vxsph_x, vxsph_y, vxsph_z;
165
166 //printf("initial_integrate at timestep %d\n", update->ntimestep);
167
168 /*
169 * get smoothed velocities from ULSPH pair style
170 */
171
172 Vector3d *smoothVel = (Vector3d *) force->pair->extract("smd/ulsph/smoothVel_ptr", itmp);
173
174 if (xsphFlag) {
175 if (smoothVel == nullptr) {
176 error->one(FLERR, "fix smd/integrate_ulsph failed to access smoothVel array");
177 }
178 }
179
180 if (igroup == atom->firstgroup)
181 nlocal = atom->nfirst;
182
183 for (i = 0; i < nlocal; i++) {
184 if (mask[i] & groupbit) {
185 dtfm = dtf / rmass[i];
186
187 v[i][0] += dtfm * f[i][0];
188 v[i][1] += dtfm * f[i][1];
189 v[i][2] += dtfm * f[i][2];
190
191 if (vlimit > 0.0) {
192 vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
193 if (vsq > vlimitsq) {
194 scale = sqrt(vlimitsq / vsq);
195 v[i][0] *= scale;
196 v[i][1] *= scale;
197 v[i][2] *= scale;
198
199 vest[i][0] = v[i][0];
200 vest[i][1] = v[i][1];
201 vest[i][2] = v[i][2];
202 }
203 }
204
205 if (xsphFlag) {
206
207 // construct XSPH velocity
208 vxsph_x = v[i][0] - 0.5 * smoothVel[i](0);
209 vxsph_y = v[i][1] - 0.5 * smoothVel[i](1);
210 vxsph_z = v[i][2] - 0.5 * smoothVel[i](2);
211
212 vest[i][0] = vxsph_x + dtfm * f[i][0];
213 vest[i][1] = vxsph_y + dtfm * f[i][1];
214 vest[i][2] = vxsph_z + dtfm * f[i][2];
215
216 x[i][0] += dtv * vxsph_x;
217 x[i][1] += dtv * vxsph_y;
218 x[i][2] += dtv * vxsph_z;
219
220
221
222
223 } else {
224
225 // extrapolate velocity from half- to full-step
226 vest[i][0] = v[i][0] + dtfm * f[i][0];
227 vest[i][1] = v[i][1] + dtfm * f[i][1];
228 vest[i][2] = v[i][2] + dtfm * f[i][2];
229
230 x[i][0] += dtv * v[i][0];
231 x[i][1] += dtv * v[i][1];
232 x[i][2] += dtv * v[i][2];
233 }
234 }
235 }
236 }
237
238 /* ---------------------------------------------------------------------- */
239
final_integrate()240 void FixSMDIntegrateUlsph::final_integrate() {
241 double **v = atom->v;
242 double **f = atom->f;
243 double *esph = atom->esph;
244 double *desph = atom->desph;
245 double *vfrac = atom->vfrac;
246 double *radius = atom->radius;
247 double *contact_radius = atom->contact_radius;
248 int *mask = atom->mask;
249 int nlocal = atom->nlocal;
250 if (igroup == atom->firstgroup)
251 nlocal = atom->nfirst;
252 double dtfm, vsq, scale;
253 double *rmass = atom->rmass;
254 double vol_increment;
255 Matrix3d D;
256
257 /*
258 * get current number of SPH neighbors from ULSPH pair style
259 */
260
261 int itmp;
262 int *nn = (int *) force->pair->extract("smd/ulsph/numNeighs_ptr", itmp);
263 if (nn == nullptr) {
264 error->one(FLERR, "fix smd/integrate_ulsph failed to accesss num_neighs array");
265 }
266
267 Matrix3d *L = (Matrix3d *) force->pair->extract("smd/ulsph/velocityGradient_ptr", itmp);
268 if (L == nullptr) {
269 error->one(FLERR, "fix smd/integrate_ulsph failed to accesss velocityGradient array");
270 }
271
272 for (int i = 0; i < nlocal; i++) {
273 if (mask[i] & groupbit) {
274
275 dtfm = dtf / rmass[i];
276 v[i][0] += dtfm * f[i][0];
277 v[i][1] += dtfm * f[i][1];
278 v[i][2] += dtfm * f[i][2];
279
280 if (vlimit > 0.0) {
281 vsq = v[i][0] * v[i][0] + v[i][1] * v[i][1] + v[i][2] * v[i][2];
282 if (vsq > vlimitsq) {
283 scale = sqrt(vlimitsq / vsq);
284 v[i][0] *= scale;
285 v[i][1] *= scale;
286 v[i][2] *= scale;
287 }
288 }
289
290 esph[i] += dtf * desph[i];
291
292 if (adjust_radius_flag) {
293 if (nn[i] < min_nn) {
294 radius[i] *= adjust_radius_factor;
295 } else if (nn[i] > max_nn) {
296 radius[i] /= adjust_radius_factor;
297 }
298 radius[i] = MAX(radius[i], 1.25 * contact_radius[i]);
299 radius[i] = MIN(radius[i], 4.0 * contact_radius[i]);
300
301 }
302
303 D = 0.5 * (L[i] + L[i].transpose());
304 vol_increment = vfrac[i] * update->dt * D.trace(); // Jacobian of deformation
305 vfrac[i] += vol_increment;
306 }
307 }
308 }
309
310 /* ---------------------------------------------------------------------- */
311
reset_dt()312 void FixSMDIntegrateUlsph::reset_dt() {
313 dtv = update->dt;
314 dtf = 0.5 * update->dt * force->ftm2v;
315 }
316
317