diff --git a/src/USER-MANIFOLD/manifold_gaussian_bump.cpp b/src/USER-MANIFOLD/manifold_gaussian_bump.cpp index 200edd2fb0..db8c589afb 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( 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 xf[3]; + xf[0] = x[0]; + xf[1] = x[1]; + xf[2] = 0.0; + nn[2] = 1.0; - double der = bb + 2*cc*xi + 3*dd*xi2; - - nn[0] = -der * x[0] / rr; - nn[1] = -der * x[1] / rr; + 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]; + } - - return x[2] - ( aa + bb*xi + cc*xi2 + dd*xi3 ); - - }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; + ll2 = 2.0*ll*ll; - f_at_rc = gaussian_bump_x2 ( rc12 ); - fp_at_rc = gaussian_bump_der( rc12 ); + rc12 = rc1*rc1; + rc22 = rc2*rc2; + dr = rc2 - rc1; + inv_dr = 1.0 / dr; - 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 ) +double manifold_gaussian_bump::gaussian_bump( double x ) const { - double x2 = x*x; - return gaussian_bump_x2( x2 ); + double x2 = x*x; + return gaussian_bump_x2( x2 ); } -double manifold_gaussian_bump::gaussian_bump_x2( double x2 ) +double manifold_gaussian_bump::gaussian_bump_x2( double x2 ) const { - return AA*exp( -x2 / ll2 ); + return AA*exp( -x2 / ll2 ); } -double manifold_gaussian_bump::gaussian_bump_der( double x ) +double manifold_gaussian_bump::gaussian_bump_der( double x ) const { - double x2 = x*x; - return gaussian_bump_x2( x2 )*( -x/(ll*ll) ); + double x2 = x*x; + return gaussian_bump_x2( x2 )*( -x/(ll*ll) ); +} + + +void manifold_gaussian_bump::make_lut() +{ + + 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]; + + lut_dx = (lut_x1 - lut_x0) / lut_Nbins; + + 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::lut_get_z ( double rr ) const +{ + double xs = rr - lut_x0; + double xss = xs / lut_dx; + int bin = static_cast(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 xs = rr - lut_x0; + double xss = xs / lut_dx; + int bin = static_cast(xss); + double frac = xss - bin; + + double zleft = lut_zp[bin]; + double zright = lut_zp[bin+1]; + + return zleft * ( 1 - frac) + frac * zright; +} + + +void manifold_gaussian_bump::lut_get_z_and_zp( double rr, double &zz, + double &zzp ) const +{ + double xs = rr - lut_x0; + double xss = xs / lut_dx; + int bin = static_cast(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; +} + + +void manifold_gaussian_bump::test_lut() +{ + 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 ); }; }