From 076990c28a6740d6b2a546b8942afc6d69b32a6f Mon Sep 17 00:00:00 2001
From: Stefan Paquay <stefanpaquay@gmail.com>
Date: Tue, 27 Jun 2017 16:48:33 -0400
Subject: [PATCH] Updated Gaussian bump so that it has a better taper function.
---
src/USER-MANIFOLD/manifold_gaussian_bump.cpp | 436 ++++++++++++++-----
src/USER-MANIFOLD/manifold_gaussian_bump.h | 31 +-
2 files changed, 360 insertions(+), 107 deletions(-)
diff --git a/src/USER-MANIFOLD/manifold_gaussian_bump.cpp b/src/USER-MANIFOLD/manifold_gaussian_bump.cpp
index 200edd2fb0..e1a3deb5ab 100644
--- a/src/USER-MANIFOLD/manifold_gaussian_bump.cpp
+++ b/src/USER-MANIFOLD/manifold_gaussian_bump.cpp
@@ -1,138 +1,374 @@
#include "manifold_gaussian_bump.h"
+#include "comm.h"
+#include "error.h"
+
using namespace LAMMPS_NS;
using namespace user_manifold;
+// This manifold comes with some baggage;
+// it needs a taper function that is sufficiently smooth
+// so we use a cubic Hermite interpolation and then for
+// efficiency we construct a lookup table for that....
+
+class cubic_hermite
+{
+public:
+ // Represent x-polynomial as a * t^3 + b * t^2 + c*t + d.
+ double a, b, c, d;
+ // And y-polynomial as s * t^3 + u * t^2 + v * t + w.
+ double s, u, v, w;
+
+ double x0, x1, y0, y1, yp0, yp1;
+
+ LAMMPS_NS::Error *err;
+
+ cubic_hermite( double x0, double x1, double y0, double y1,
+ double yp0, double yp1, LAMMPS_NS::Error *err ) :
+ x0(x0), x1(x1), y0(y0), y1(y1), yp0(yp0), yp1(yp1),
+ a( 2*x0 + 2 - 2*x1 ),
+ b( -3*x0 - 3 + 3*x1 ),
+ c( 1.0 ),
+ d( x0 ),
+ s( 2*y0 - 2*y1 + yp0 + yp1 ),
+ u( -3*y0 + 3*y1 - 2*yp0 - yp1 ),
+ v( yp0 ),
+ w( y0 ),
+ err(err)
+ {
+ test();
+ }
+
+
+ void test()
+ {
+ if( fabs( x(0) - x0 ) > 1e-8 ) err->one(FLERR, "x0 wrong");
+ if( fabs( x(1) - x1 ) > 1e-8 ) err->one(FLERR, "x1 wrong");
+ if( fabs( y(0) - y0 ) > 1e-8 ) err->one(FLERR, "y0 wrong");
+ if( fabs( y(1) - y1 ) > 1e-8 ) err->one(FLERR, "y1 wrong");
+ }
+
+ double get_t_from_x( double xx ) const
+ {
+ if( xx < x0 || xx > x1 ){
+ char msg[2048];
+ sprintf(msg, "x ( %g ) out of bounds [%g, %g]", xx, x0, x1 );
+ err->one(FLERR, msg);
+ }
+
+ // Newton iterate to get right t.
+ double t = xx - x0;
+ double dx = x1 - x0;
+ // Reasonable initial guess I hope:
+ t /= dx;
+
+ double tol = 1e-8;
+ int maxit = 500;
+ double ff = x(t) - xx;
+ double ffp = xp(t);
+ // double l = 1.0 / ( 1 + res*res );
+ for( int i = 0; i < maxit; ++i ){
+ t -= ff / ffp;
+ ff = x(t) - xx;
+ ffp = xp(t);
+ double res = ff;
+ if( std::fabs( res ) < tol ){
+ return t;
+ }
+ }
+ err->warning(FLERR, "Convergence failed");
+ return t;
+ }
+
+ double x( double t ) const
+ {
+ double t2 = t*t;
+ double t3 = t2*t;
+ return a*t3 + b*t2 + c*t + d;
+ }
+
+ double y_from_x( double x ) const
+ {
+ double t = get_t_from_x( x );
+ return y(t);
+ }
+
+ double yp_from_x( double x ) const
+ {
+ double t = get_t_from_x( x );
+ return yp(t);
+ }
+
+ double y( double t ) const
+ {
+ double t2 = t*t;
+ double t3 = t2*t;
+ return s*t3 + u*t2 + v*t + w;
+ }
+
+ void xy( double t, double &xx, double &yy ) const
+ {
+ xx = x(t);
+ yy = y(t);
+ }
+
+ double xp( double t ) const
+ {
+ double t2 = t*t;
+ return 3*a*t2 + 2*b*t + c;
+ }
+
+ double yp( double t ) const
+ {
+ double t2 = t*t;
+ return 3*t2*s + 2*u*t + v;
+ }
+
+ double xpp( double t ) const
+ {
+ return 6*a*t + 2*b;
+ }
+
+};
+
+// Manifold itself:
+manifold_gaussian_bump::manifold_gaussian_bump(class LAMMPS* lmp,
+ int narg, char **arg)
+ : manifold(lmp), lut_z(NULL), lut_zp(NULL) {}
+
+
+manifold_gaussian_bump::~manifold_gaussian_bump()
+{
+ if( lut_z ) delete lut_z;
+ if( lut_zp ) delete lut_zp;
+}
+
+
// The constraint is z = f(r = sqrt( x^2 + y^2) )
// f(x) = A*exp( -x^2 / 2 l^2 ) if x < rc1
-// = a + b*x + c*x**2 + d*x**3 if rc1 <= x < rc2
+// = Some interpolation if rc1 <= rc2
// = 0 if x >= rc2
//
double manifold_gaussian_bump::g( const double *x )
{
- double xf[3];
- xf[0] = x[0];
- xf[1] = x[1];
- xf[2] = 0.0;
-
- double x2 = dot(xf,xf);
- if( x2 < rc12 ){
- return x[2] - gaussian_bump_x2( x2 );
- }else if( x2 < rc22 ){
- double rr = sqrt(x2);
- double xi = rr - rc1;
- xi *= inv_dr;
- double xi2 = x2 * inv_dr*inv_dr;
- double xi3 = xi*xi2;
- return x[2] - ( aa + bb*xi + cc*xi2 + dd*xi3 );
-
- }else{
- return x[2];
- }
+ double xf[3];
+ xf[0] = x[0];
+ xf[1] = x[1];
+ xf[2] = 0.0;
+
+ double x2 = dot(xf,xf);
+ if( x2 < rc12 ){
+ return x[2] - gaussian_bump_x2( x2 );
+ }else if( x2 < rc22 ){
+ double rr = sqrt( x2 );
+ double z_taper_func = lut_get_z( rr );
+ return x[2] - z_taper_func;
+ }else{
+ return x[2];
+ }
}
void manifold_gaussian_bump::n( const double *x, double *nn )
{
- double xf[3];
- xf[0] = x[0];
- xf[1] = x[1];
- xf[2] = 0.0;
- nn[2] = 1.0;
-
- double x2 = dot(xf,xf);
-
- if( x2 < rc12 ){
- double factor = gaussian_bump_x2(x2);
- factor /= (ll*ll);
- nn[0] = factor * x[0];
- nn[1] = factor * x[1];
- }else if( x2 < rc22 ){
- double rr = sqrt(x2);
- double xi = rr - rc1;
- xi *= inv_dr;
- double xi2 = x2 * inv_dr*inv_dr;
- double der = bb + 2*cc*xi + 3*dd*xi2;
-
- nn[0] = -der * x[0] / rr;
- nn[1] = -der * x[1] / rr;
- }else{
- nn[0] = nn[1] = 0.0;
- }
+ double xf[3];
+ xf[0] = x[0];
+ xf[1] = x[1];
+ xf[2] = 0.0;
+ nn[2] = 1.0;
+
+ double x2 = dot(xf,xf);
+
+ if( x2 < rc12 ){
+ double factor = gaussian_bump_x2(x2);
+ factor /= (ll*ll);
+ nn[0] = factor * x[0];
+ nn[1] = factor * x[1];
+ }else if( x2 < rc22 ){
+ double rr = sqrt( x2 );
+ double zp_taper_func = lut_get_zp( rr );
+
+ double inv_r = 1.0 / rr;
+ double der_part = zp_taper_func * inv_r;
+
+ nn[0] = der_part * x[0];
+ nn[1] = der_part * x[1];
+
+ }else{
+ nn[0] = nn[1] = 0.0;
+ }
}
double manifold_gaussian_bump::g_and_n( const double *x, double *nn )
{
- double xf[3];
- xf[0] = x[0];
- xf[1] = x[1];
- xf[2] = 0.0;
- nn[2] = 1.0;
-
- double x2 = dot(xf,xf);
- if( x2 < rc12 ){
- double gb = gaussian_bump_x2(x2);
- double factor = gb / (ll*ll);
- nn[0] = factor * x[0];
- nn[1] = factor * x[1];
-
- return x[2] - gb;
- }else if( x2 < rc22 ){
-
- double rr = sqrt(x2);
- double xi = rr - rc1;
- xi *= inv_dr;
- double xi2 = x2 * inv_dr*inv_dr;
- double xi3 = xi*xi2;
-
- double der = bb + 2*cc*xi + 3*dd*xi2;
-
- nn[0] = -der * x[0] / rr;
- nn[1] = -der * x[1] / rr;
-
-
- return x[2] - ( aa + bb*xi + cc*xi2 + dd*xi3 );
-
- }else{
- nn[0] = nn[1] = 0.0;
- return x[2];
- }
-
+ double xf[3];
+ xf[0] = x[0];
+ xf[1] = x[1];
+ xf[2] = 0.0;
+ nn[2] = 1.0;
+
+ double x2 = dot(xf,xf);
+ if( x2 < rc12 ){
+ double gb = gaussian_bump_x2(x2);
+ double factor = gb / (ll*ll);
+ nn[0] = factor * x[0];
+ nn[1] = factor * x[1];
+
+ return x[2] - gb;
+ }else if( x2 < rc22 ){
+ double z_taper_func, zp_taper_func;
+ double rr = sqrt( x2 );
+ lut_get_z_and_zp( rr, z_taper_func, zp_taper_func );
+
+ double inv_r = 1.0 / rr;
+ double der_part = zp_taper_func * inv_r;
+
+ nn[0] = der_part * x[0];
+ nn[1] = der_part * x[1];
+
+ return x[2] - z_taper_func;
+ }else{
+ nn[0] = nn[1] = 0.0;
+ return x[2];
+ }
+
}
void manifold_gaussian_bump::post_param_init()
{
- // Read in the params:
- AA = params[0];
- ll = params[1];
- rc1 = params[2];
- rc2 = params[3];
+ // Read in the params:
+ AA = params[0];
+ ll = params[1];
+ rc1 = params[2];
+ rc2 = params[3];
+
+ ll2 = 2.0*ll*ll;
+
+ rc12 = rc1*rc1;
+ rc22 = rc2*rc2;
+ dr = rc2 - rc1;
+ inv_dr = 1.0 / dr;
+
+ f_at_rc = gaussian_bump_x2 ( rc12 );
+ fp_at_rc = gaussian_bump_der( rc1 );
+
+
+
+ make_lut();
+
+ // test_lut();
+}
+
+
+double manifold_gaussian_bump::gaussian_bump( double x ) const
+{
+ double x2 = x*x;
+ return gaussian_bump_x2( x2 );
+}
+
+double manifold_gaussian_bump::gaussian_bump_x2( double x2 ) const
+{
+ return AA*exp( -x2 / ll2 );
+}
+
+double manifold_gaussian_bump::gaussian_bump_der( double x ) const
+{
+ double x2 = x*x;
+ return gaussian_bump_x2( x2 )*( -x/(ll*ll) );
+}
+
+
+void manifold_gaussian_bump::make_lut()
+{
- ll2 = 2.0*ll*ll;
+ lut_x0 = rc1;
+ lut_x1 = rc2;
+ lut_Nbins = 1024; // Don't make it too big, it should fit in cache.
+ lut_z = new double[lut_Nbins+1];
+ lut_zp = new double[lut_Nbins+1];
- f_at_rc = gaussian_bump_x2 ( rc12 );
- fp_at_rc = gaussian_bump_der( rc12 );
+ lut_dx = (lut_x1 - lut_x0) / lut_Nbins;
- rc12 = rc1*rc1;
- rc22 = rc2*rc2;
- dr = rc2 - rc1;
- inv_dr = 1.0 / dr;
+ cubic_hermite pchip( lut_x0, lut_x1, f_at_rc, 0.0, fp_at_rc, 0.0, error );
+
+ double xx = lut_x0;
+ for( int i = 0; i <= lut_Nbins; ++i ){
+ lut_z[i] = pchip.y_from_x( xx );
+ lut_zp[i] = pchip.yp_from_x( xx );
+ xx += lut_dx;
+ }
}
-double manifold_gaussian_bump::gaussian_bump( double x )
+double manifold_gaussian_bump::lut_get_z ( double rr ) const
+{
+ double xs = rr - lut_x0;
+ double xss = xs / lut_dx;
+ int bin = static_cast<int>(xss);
+ double frac = xss - bin;
+
+ double zleft = lut_z[bin];
+ double zright = lut_z[bin+1];
+
+ return zleft * ( 1 - frac ) + frac * zright;
+}
+
+double manifold_gaussian_bump::lut_get_zp( double rr ) const
{
- double x2 = x*x;
- return gaussian_bump_x2( x2 );
+ double xs = rr - lut_x0;
+ double xss = xs / lut_dx;
+ int bin = static_cast<int>(xss);
+ double frac = xss - bin;
+
+ double zleft = lut_zp[bin];
+ double zright = lut_zp[bin+1];
+
+ return zleft * ( 1 - frac) + frac * zright;
}
-double manifold_gaussian_bump::gaussian_bump_x2( double x2 )
+
+void manifold_gaussian_bump::lut_get_z_and_zp( double rr, double &zz,
+ double &zzp ) const
{
- return AA*exp( -x2 / ll2 );
+ double xs = rr - lut_x0;
+ double xss = xs / lut_dx;
+ int bin = static_cast<int>(xss);
+ double frac = xss - bin;
+ double fmin = 1 - frac;
+
+ double zleft = lut_z[bin];
+ double zright = lut_z[bin+1];
+ double zpleft = lut_zp[bin];
+ double zpright = lut_zp[bin+1];
+
+ zz = zleft * fmin + zright * frac;
+ zzp = zpleft * fmin + zpright * frac;
}
-double manifold_gaussian_bump::gaussian_bump_der( double x )
+
+void manifold_gaussian_bump::test_lut()
{
- double x2 = x*x;
- return gaussian_bump_x2( x2 )*( -x/(ll*ll) );
+ double x[3], nn[3];
+ if( comm->me != 0 ) return;
+
+ FILE *fp = fopen( "test_lut_gaussian.dat", "w" );
+ double dx = 0.1;
+ for( double xx = 0; xx < 20; xx += dx ){
+ x[0] = xx;
+ x[1] = 0.0;
+ x[2] = 0.0;
+ double gg = g( x );
+ n( x, nn );
+ double taper_z;
+ if( xx <= rc1 ){
+ taper_z = gaussian_bump(xx);
+ }else if( xx < rc2 ){
+ taper_z = lut_get_z( xx );
+ }else{
+ taper_z = 0.0;
+ }
+ fprintf( fp, "%g %g %g %g %g\n", xx, gaussian_bump(xx), taper_z,
+ gg, nn[0], nn[1], nn[2] );
+ }
+ fclose(fp);
}
diff --git a/src/USER-MANIFOLD/manifold_gaussian_bump.h b/src/USER-MANIFOLD/manifold_gaussian_bump.h
index f31e2a4bf4..4fef84e3a5 100644
--- a/src/USER-MANIFOLD/manifold_gaussian_bump.h
+++ b/src/USER-MANIFOLD/manifold_gaussian_bump.h
@@ -1,4 +1,4 @@
-/* -*- c++ -*- ----------------------------------------------------------
+/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
@@ -50,8 +50,8 @@ namespace user_manifold {
class manifold_gaussian_bump : public manifold {
public:
enum { NPARAMS = 4 };
- manifold_gaussian_bump(class LAMMPS*, int, char **) : manifold(lmp) {}
- virtual ~manifold_gaussian_bump(){}
+ manifold_gaussian_bump(class LAMMPS*, int, char **);
+ virtual ~manifold_gaussian_bump();
virtual double g( const double * );
virtual void n( const double *, double * );
@@ -66,12 +66,29 @@ namespace user_manifold {
virtual void post_param_init();
private:
// Some private constants:
- double aa, bb, cc, dd, AA, ll, ll2, f_at_rc, fp_at_rc;
+ double AA, ll, ll2, f_at_rc, fp_at_rc;
double rc1, rc2, rc12, rc22, dr, inv_dr;
- double gaussian_bump ( double );
- double gaussian_bump_x2 ( double );
- double gaussian_bump_der( double );
+ // Stuff for the look-up table:
+ double lut_x0, lut_x1;
+ int lut_Nbins;
+ double lut_dx;
+ double *lut_z;
+ double *lut_zp;
+
+ double gaussian_bump ( double ) const;
+ double gaussian_bump_x2 ( double ) const;
+ double gaussian_bump_der( double ) const;
+
+ void make_lut();
+ double lut_get_z ( double rr ) const;
+ double lut_get_zp( double rr ) const;
+ void lut_get_z_and_zp( double rr, double &zz, double &zzp ) const;
+
+ void test_lut();
+
+ double taper( double );
+ double taper_der( double );
};
}
--
GitLab