1 // clang-format off
2 /* ----------------------------------------------------------------------
3 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4 https://www.lammps.org/, Sandia National Laboratories
5 Steve Plimpton, sjplimp@sandia.gov
6
7 Copyright (2003) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level LAMMPS directory.
13 ----------------------------------------------------------------------- */
14
15 #include "manifold_gaussian_bump.h"
16
17 #include "comm.h"
18 #include "error.h"
19
20 #include <cmath>
21
22 using namespace LAMMPS_NS;
23 using namespace user_manifold;
24
25 // This manifold comes with some baggage;
26 // it needs a taper function that is sufficiently smooth
27 // so we use a cubic Hermite interpolation and then for
28 // efficiency we construct a lookup table for that....
29
30 class cubic_hermite
31 {
32 public:
33 // Represent x-polynomial as a * t^3 + b * t^2 + c*t + d.
34 double a, b, c, d;
35 // And y-polynomial as s * t^3 + u * t^2 + v * t + w.
36 double s, u, v, w;
37
38 double x0, x1, y0, y1, yp0, yp1;
39
40 LAMMPS_NS::Error *err;
41
cubic_hermite(double x0,double x1,double y0,double y1,double yp0,double yp1,LAMMPS_NS::Error * err)42 cubic_hermite( double x0, double x1, double y0, double y1,
43 double yp0, double yp1, LAMMPS_NS::Error *err ) :
44 a( 2*x0 + 2 - 2*x1 ),
45 b( -3*x0 - 3 + 3*x1 ),
46 c( 1.0 ),
47 d( x0 ),
48 s( 2*y0 - 2*y1 + yp0 + yp1 ),
49 u( -3*y0 + 3*y1 - 2*yp0 - yp1 ),
50 v( yp0 ),
51 w( y0 ),
52 x0(x0), x1(x1), y0(y0), y1(y1), yp0(yp0), yp1(yp1),
53 err(err)
54 {
55 test();
56 }
57
58
test()59 void test()
60 {
61 if (fabs( x(0) - x0 ) > 1e-8 ) err->one(FLERR, "x0 wrong");
62 if (fabs( x(1) - x1 ) > 1e-8 ) err->one(FLERR, "x1 wrong");
63 if (fabs( y(0) - y0 ) > 1e-8 ) err->one(FLERR, "y0 wrong");
64 if (fabs( y(1) - y1 ) > 1e-8 ) err->one(FLERR, "y1 wrong");
65 }
66
get_t_from_x(double xx) const67 double get_t_from_x( double xx ) const
68 {
69 if (xx < x0 || xx > x1) {
70 char msg[2048];
71 sprintf(msg, "x ( %g ) out of bounds [%g, %g]", xx, x0, x1 );
72 err->one(FLERR, msg);
73 }
74
75 // Newton iterate to get right t.
76 double t = xx - x0;
77 double dx = x1 - x0;
78 // Reasonable initial guess I hope:
79 t /= dx;
80
81 double tol = 1e-8;
82 int maxit = 500;
83 double ff = x(t) - xx;
84 double ffp = xp(t);
85 // double l = 1.0 / ( 1 + res*res );
86 for (int i = 0; i < maxit; ++i) {
87 t -= ff / ffp;
88 ff = x(t) - xx;
89 ffp = xp(t);
90 double res = ff;
91 if (fabs( res ) < tol) {
92 return t;
93 }
94 }
95 err->warning(FLERR, "Convergence failed");
96 return t;
97 }
98
x(double t) const99 double x( double t ) const
100 {
101 double t2 = t*t;
102 double t3 = t2*t;
103 return a*t3 + b*t2 + c*t + d;
104 }
105
y_from_x(double x) const106 double y_from_x( double x ) const
107 {
108 double t = get_t_from_x( x );
109 return y(t);
110 }
111
yp_from_x(double x) const112 double yp_from_x( double x ) const
113 {
114 double t = get_t_from_x( x );
115 return yp(t);
116 }
117
y(double t) const118 double y( double t ) const
119 {
120 double t2 = t*t;
121 double t3 = t2*t;
122 return s*t3 + u*t2 + v*t + w;
123 }
124
xy(double t,double & xx,double & yy) const125 void xy( double t, double &xx, double &yy ) const
126 {
127 xx = x(t);
128 yy = y(t);
129 }
130
xp(double t) const131 double xp( double t ) const
132 {
133 double t2 = t*t;
134 return 3*a*t2 + 2*b*t + c;
135 }
136
yp(double t) const137 double yp( double t ) const
138 {
139 double t2 = t*t;
140 return 3*t2*s + 2*u*t + v;
141 }
142
xpp(double t) const143 double xpp( double t ) const
144 {
145 return 6*a*t + 2*b;
146 }
147
148 };
149
150 // Manifold itself:
manifold_gaussian_bump(class LAMMPS * lmp,int,char **)151 manifold_gaussian_bump::manifold_gaussian_bump(class LAMMPS* lmp,
152 int /*narg*/, char **/*arg*/)
153 : manifold(lmp), lut_z(nullptr), lut_zp(nullptr) {}
154
155
~manifold_gaussian_bump()156 manifold_gaussian_bump::~manifold_gaussian_bump()
157 {
158 if (lut_z ) delete lut_z;
159 if (lut_zp) delete lut_zp;
160 }
161
162
163 // The constraint is z = f(r = sqrt( x^2 + y^2) )
164 // f(x) = A*exp( -x^2 / 2 l^2 ) if x < rc1
165 // = Some interpolation if rc1 <= rc2
166 // = 0 if x >= rc2
167 //
g(const double * x)168 double manifold_gaussian_bump::g( const double *x )
169 {
170 double xf[3];
171 xf[0] = x[0];
172 xf[1] = x[1];
173 xf[2] = 0.0;
174
175 double x2 = dot(xf,xf);
176 if (x2 < rc12) {
177 return x[2] - gaussian_bump_x2( x2 );
178 } else if (x2 < rc22) {
179 double rr = sqrt( x2 );
180 double z_taper_func = lut_get_z( rr );
181 return x[2] - z_taper_func;
182 } else {
183 return x[2];
184 }
185 }
186
n(const double * x,double * nn)187 void manifold_gaussian_bump::n( const double *x, double *nn )
188 {
189 double xf[3];
190 xf[0] = x[0];
191 xf[1] = x[1];
192 xf[2] = 0.0;
193 nn[2] = 1.0;
194
195 double x2 = dot(xf,xf);
196
197 if (x2 < rc12) {
198 double factor = gaussian_bump_x2(x2);
199 factor /= (ll*ll);
200 nn[0] = factor * x[0];
201 nn[1] = factor * x[1];
202 } else if (x2 < rc22) {
203 double rr = sqrt( x2 );
204 double zp_taper_func = lut_get_zp( rr );
205
206 double inv_r = 1.0 / rr;
207 double der_part = zp_taper_func * inv_r;
208
209 nn[0] = der_part * x[0];
210 nn[1] = der_part * x[1];
211
212 } else {
213 nn[0] = nn[1] = 0.0;
214 }
215 }
216
g_and_n(const double * x,double * nn)217 double manifold_gaussian_bump::g_and_n( const double *x, double *nn )
218 {
219 double xf[3];
220 xf[0] = x[0];
221 xf[1] = x[1];
222 xf[2] = 0.0;
223 nn[2] = 1.0;
224
225 double x2 = dot(xf,xf);
226 if (x2 < rc12) {
227 double gb = gaussian_bump_x2(x2);
228 double factor = gb / (ll*ll);
229 nn[0] = factor * x[0];
230 nn[1] = factor * x[1];
231
232 return x[2] - gb;
233 } else if (x2 < rc22) {
234 double z_taper_func, zp_taper_func;
235 double rr = sqrt( x2 );
236 lut_get_z_and_zp( rr, z_taper_func, zp_taper_func );
237
238 double inv_r = 1.0 / rr;
239 double der_part = zp_taper_func * inv_r;
240
241 nn[0] = der_part * x[0];
242 nn[1] = der_part * x[1];
243
244 return x[2] - z_taper_func;
245 } else {
246 nn[0] = nn[1] = 0.0;
247 return x[2];
248 }
249
250 }
251
252
post_param_init()253 void manifold_gaussian_bump::post_param_init()
254 {
255 // Read in the params:
256 AA = params[0];
257 ll = params[1];
258 rc1 = params[2];
259 rc2 = params[3];
260
261 ll2 = 2.0*ll*ll;
262
263 rc12 = rc1*rc1;
264 rc22 = rc2*rc2;
265 dr = rc2 - rc1;
266 inv_dr = 1.0 / dr;
267
268 f_at_rc = gaussian_bump_x2 ( rc12 );
269 fp_at_rc = gaussian_bump_der( rc1 );
270
271
272
273 make_lut();
274
275 // test_lut();
276 }
277
278
gaussian_bump(double x) const279 double manifold_gaussian_bump::gaussian_bump( double x ) const
280 {
281 double x2 = x*x;
282 return gaussian_bump_x2( x2 );
283 }
284
gaussian_bump_x2(double x2) const285 double manifold_gaussian_bump::gaussian_bump_x2( double x2 ) const
286 {
287 return AA*exp( -x2 / ll2 );
288 }
289
gaussian_bump_der(double x) const290 double manifold_gaussian_bump::gaussian_bump_der( double x ) const
291 {
292 double x2 = x*x;
293 return gaussian_bump_x2( x2 )*( -x/(ll*ll) );
294 }
295
296
make_lut()297 void manifold_gaussian_bump::make_lut()
298 {
299
300 lut_x0 = rc1;
301 lut_x1 = rc2;
302 lut_Nbins = 1024; // Don't make it too big, it should fit in cache.
303 lut_z = new double[lut_Nbins+1];
304 lut_zp = new double[lut_Nbins+1];
305
306 lut_dx = (lut_x1 - lut_x0) / lut_Nbins;
307
308 cubic_hermite pchip( lut_x0, lut_x1, f_at_rc, 0.0, fp_at_rc, 0.0, error );
309
310 double xx = lut_x0;
311 for (int i = 0; i <= lut_Nbins; ++i) {
312 lut_z[i] = pchip.y_from_x( xx );
313 lut_zp[i] = pchip.yp_from_x( xx );
314 xx += lut_dx;
315 }
316 }
317
318
lut_get_z(double rr) const319 double manifold_gaussian_bump::lut_get_z ( double rr ) const
320 {
321 double xs = rr - lut_x0;
322 double xss = xs / lut_dx;
323 int bin = static_cast<int>(xss);
324 double frac = xss - bin;
325
326 double zleft = lut_z[bin];
327 double zright = lut_z[bin+1];
328
329 return zleft * ( 1 - frac ) + frac * zright;
330 }
331
lut_get_zp(double rr) const332 double manifold_gaussian_bump::lut_get_zp( double rr ) const
333 {
334 double xs = rr - lut_x0;
335 double xss = xs / lut_dx;
336 int bin = static_cast<int>(xss);
337 double frac = xss - bin;
338
339 double zleft = lut_zp[bin];
340 double zright = lut_zp[bin+1];
341
342 return zleft * ( 1 - frac) + frac * zright;
343 }
344
345
lut_get_z_and_zp(double rr,double & zz,double & zzp) const346 void manifold_gaussian_bump::lut_get_z_and_zp( double rr, double &zz,
347 double &zzp ) const
348 {
349 double xs = rr - lut_x0;
350 double xss = xs / lut_dx;
351 int bin = static_cast<int>(xss);
352 double frac = xss - bin;
353 double fmin = 1 - frac;
354
355 double zleft = lut_z[bin];
356 double zright = lut_z[bin+1];
357 double zpleft = lut_zp[bin];
358 double zpright = lut_zp[bin+1];
359
360 zz = zleft * fmin + zright * frac;
361 zzp = zpleft * fmin + zpright * frac;
362 }
363
364
test_lut()365 void manifold_gaussian_bump::test_lut()
366 {
367 double x[3], nn[3];
368 if (comm->me != 0) return;
369
370 FILE *fp = fopen( "test_lut_gaussian.dat", "w" );
371 double dx = 0.1;
372 for (double xx = 0; xx < 20; xx += dx) {
373 x[0] = xx;
374 x[1] = 0.0;
375 x[2] = 0.0;
376 double gg = g( x );
377 n( x, nn );
378 double taper_z;
379 if (xx <= rc1) {
380 taper_z = gaussian_bump(xx);
381 } else if (xx < rc2) {
382 taper_z = lut_get_z( xx );
383 } else {
384 taper_z = 0.0;
385 }
386 fprintf( fp, "%g %g %g %g %g %g %g\n", xx, gaussian_bump(xx), taper_z,
387 gg, nn[0], nn[1], nn[2] );
388 }
389 fclose(fp);
390 }
391