From f6bdf5662ec9af4f0eefcfc31fe4c61ba1d1bd10 Mon Sep 17 00:00:00 2001 From: "Jibril B. Coulibaly" <43829860+jibril-b-coulibaly@users.noreply.github.com> Date: Fri, 6 Dec 2019 10:40:04 -0600 Subject: [PATCH 01/14] Update pair_granular.rst Correct formula for tangent stiffness, consistent with `PairGranular::mix_stiffnessG()` --- doc/src/pair_granular.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/src/pair_granular.rst b/doc/src/pair_granular.rst index d5f522462d..294bec9b8c 100644 --- a/doc/src/pair_granular.rst +++ b/doc/src/pair_granular.rst @@ -432,7 +432,7 @@ option by an additional factor of *a*\ , the radius of the contact region. The t Here, *a* is the radius of the contact region, given by :math:`a =\sqrt{R\delta}` for all normal contact models, except for *jkr*\ , where it is given implicitly by :math:`\delta = a^2/R - 2\sqrt{\pi \gamma a/E}`, see -discussion above. To match the Mindlin solution, one should set :math:`k_t = 8G`, where :math:`G` is the shear modulus, related to Young's modulus +discussion above. To match the Mindlin solution, one should set :math:`k_t = 4G/(2-\nu)`, where :math:`G` is the shear modulus, related to Young's modulus :math:`E` by :math:`G = E/(2(1+\nu))`, where :math:`\nu` is Poisson's ratio. This can also be achieved by specifying *NULL* for :math:`k_t`, in which case a normal contact model that specifies material parameters :math:`E` and From 5859b18d2f7b1d63b45a9e541a6908f598b0a16b Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Fri, 6 Dec 2019 11:04:50 -0700 Subject: [PATCH 02/14] enable fix rigid commands to add gravity to COM of rigid bodies --- src/RIGID/fix_rigid.cpp | 62 +++++++++++++++++++++++++++++------ src/RIGID/fix_rigid.h | 6 +++- src/RIGID/fix_rigid_small.cpp | 40 ++++++++++++++++++++-- src/RIGID/fix_rigid_small.h | 3 ++ src/fix_gravity.cpp | 52 +++++++++++++++++++++++++---- src/fix_gravity.h | 5 ++- 6 files changed, 146 insertions(+), 22 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 3e6130ac7c..5e994810ac 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -52,8 +52,8 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : torque(NULL), quat(NULL), imagebody(NULL), fflag(NULL), tflag(NULL), langextra(NULL), sum(NULL), all(NULL), remapflag(NULL), xcmimage(NULL), eflags(NULL), orient(NULL), - dorient(NULL), id_dilate(NULL), random(NULL), avec_ellipsoid(NULL), - avec_line(NULL), avec_tri(NULL) + dorient(NULL), id_dilate(NULL), id_gravity(NULL), random(NULL), + avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL) { int i,ibody; @@ -124,14 +124,17 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : if (custom_flag) { if (narg < 5) error->all(FLERR,"Illegal fix rigid command"); - // determine whether atom-style variable or atom property is used. + // determine whether atom-style variable or atom property is used + if (strstr(arg[4],"i_") == arg[4]) { int is_double=0; int custom_index = atom->find_custom(arg[4]+2,is_double); if (custom_index == -1) - error->all(FLERR,"Fix rigid custom requires previously defined property/atom"); + error->all(FLERR,"Fix rigid custom requires " + "previously defined property/atom"); else if (is_double) - error->all(FLERR,"Fix rigid custom requires integer-valued property/atom"); + error->all(FLERR,"Fix rigid custom requires " + "integer-valued property/atom"); int minval = INT_MAX; int *value = atom->ivector[custom_index]; for (i = 0; i < nlocal; i++) @@ -163,12 +166,15 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : if (mask[i] & groupbit) molecule[i] = (tagint)((tagint)value[i] - minval + 1); delete[] value; + } else error->all(FLERR,"Unsupported fix rigid custom property"); + } else { if (atom->molecule_flag == 0) error->all(FLERR,"Fix rigid molecule requires atom attribute molecule"); molecule = atom->molecule; } + iarg = 4 + custom_flag; tagint maxmol_tag = -1; @@ -297,6 +303,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : } // number of linear rigid bodies is counted later + nlinear = 0; // parse optional args @@ -308,12 +315,14 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : tstat_flag = 0; pstat_flag = 0; allremap = 1; - id_dilate = NULL; t_chain = 10; t_iter = 1; t_order = 3; p_chain = 10; + inpfile = NULL; + id_gravity = NULL; + id_dilate = NULL; pcouple = NONE; pstyle = ANISO; @@ -543,12 +552,20 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : iarg += 2; } else if (strcmp(arg[iarg],"reinit") == 0) { - if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); if (strcmp("yes",arg[iarg+1]) == 0) reinitflag = 1; else if (strcmp("no",arg[iarg+1]) == 0) reinitflag = 0; else error->all(FLERR,"Illegal fix rigid command"); iarg += 2; + } else if (strcmp(arg[iarg],"gravity") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid command"); + delete [] id_gravity; + int n = strlen(arg[iarg+1]) + 1; + id_gravity = new char[n]; + strcpy(id_gravity,arg[iarg+1]); + iarg += 2; + } else error->all(FLERR,"Illegal fix rigid command"); } @@ -621,6 +638,9 @@ FixRigid::~FixRigid() delete random; delete [] inpfile; + delete [] id_dilate; + delete [] id_gravity; + memory->destroy(mol2body); memory->destroy(body2mol); @@ -687,7 +707,7 @@ void FixRigid::init() avec_tri = (AtomVecTri *) atom->style_match("tri"); // warn if more than one rigid fix - // if earlyflag, warn if any post-force fixes come after POEMS fix + // if earlyflag, warn if any post-force fixes come after a rigid fix int count = 0; for (i = 0; i < modify->nfix; i++) @@ -701,7 +721,8 @@ void FixRigid::init() if (rflag && (modify->fmask[i] & POST_FORCE) && !modify->fix[i]->rigid_flag) { char str[128]; - snprintf(str,128,"Fix %s alters forces after fix rigid",modify->fix[i]->id); + snprintf(str,128,"Fix %s alters forces after fix rigid", + modify->fix[i]->id); error->warning(FLERR,str); } } @@ -713,12 +734,24 @@ void FixRigid::init() if (strcmp(modify->fix[i]->style,"npt") == 0) break; if (strcmp(modify->fix[i]->style,"nph") == 0) break; } + if (i < modify->nfix) { for (int j = i; j < modify->nfix; j++) if (strcmp(modify->fix[j]->style,"rigid") == 0) error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } + // add gravity forces based on gravity vector from fix + + if (id_gravity) { + int ifix = modify->find_fix(id_gravity); + if (ifix < 0) error->all(FLERR,"Fix rigid cannot find fix gravity ID"); + if (strcmp(modify->fix[ifix]->style,"gravity") != 0) + error->all(FLERR,"Fix rigid gravity fix is invalid"); + int tmp; + gvec = (double *) modify->fix[ifix]->extract("gvec",tmp); + } + // timestep info dtv = update->dt; @@ -1061,8 +1094,17 @@ void FixRigid::compute_forces_and_torques() torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; } -} + // add gravity force to COM of each body + + if (id_gravity) { + for (ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] += gvec[0]*masstotal[ibody]; + fcm[ibody][1] += gvec[1]*masstotal[ibody]; + fcm[ibody][2] += gvec[2]*masstotal[ibody]; + } + } +} /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid.h b/src/RIGID/fix_rigid.h index 8c4ab9ed49..c261d89cce 100644 --- a/src/RIGID/fix_rigid.h +++ b/src/RIGID/fix_rigid.h @@ -66,7 +66,8 @@ class FixRigid : public Fix { double *step_respa; int triclinic; - char *inpfile; // file to read rigid body attributes from + char *inpfile; // file to read rigid body attributes from + int rstyle; // SINGLE,MOLECULE,GROUP int setupflag; // 1 if body properties are setup, else 0 int earlyflag; // 1 if forces/torques computed at post_force() @@ -131,6 +132,9 @@ class FixRigid : public Fix { int dilate_group_bit; // mask for dilation group char *id_dilate; // group name to dilate + char *id_gravity; // ID of fix gravity command to add gravity forces + double *gvec; // ptr to gravity vector inside the fix + class RanMars *random; class AtomVecEllipsoid *avec_ellipsoid; class AtomVecLine *avec_line; diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index ba570d5840..4d2776b9ff 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -57,7 +57,7 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : xcmimage(NULL), displace(NULL), eflags(NULL), orient(NULL), dorient(NULL), avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL), counts(NULL), itensor(NULL), mass_body(NULL), langextra(NULL), random(NULL), - id_dilate(NULL), onemols(NULL) + id_dilate(NULL), id_gravity(NULL), onemols(NULL) { int i; @@ -107,7 +107,8 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : bodyID = new tagint[nlocal]; customflag = 1; - // determine whether atom-style variable or atom property is used. + // determine whether atom-style variable or atom property is used + if (strstr(arg[4],"i_") == arg[4]) { int is_double=0; int custom_index = atom->find_custom(arg[4]+2,is_double); @@ -356,6 +357,13 @@ FixRigidSmall::FixRigidSmall(LAMMPS *lmp, int narg, char **arg) : p_chain = force->numeric(FLERR,arg[iarg+1]); iarg += 2; + } else if (strcmp(arg[iarg],"gravity") == 0) { + if (iarg+2 > narg) error->all(FLERR,"Illegal fix rigid/small command"); + delete [] id_gravity; + int n = strlen(arg[iarg+1]) + 1; + id_gravity = new char[n]; + strcpy(id_gravity,arg[iarg+1]); + iarg += 2; } else error->all(FLERR,"Illegal fix rigid/small command"); } @@ -515,6 +523,8 @@ FixRigidSmall::~FixRigidSmall() delete random; delete [] inpfile; + delete [] id_dilate; + delete [] id_gravity; memory->destroy(langextra); memory->destroy(mass_body); @@ -543,6 +553,7 @@ void FixRigidSmall::init() triclinic = domain->triclinic; // warn if more than one rigid fix + // if earlyflag, warn if any post-force fixes come after a rigid fix int count = 0; for (i = 0; i < modify->nfix; i++) @@ -575,6 +586,17 @@ void FixRigidSmall::init() error->all(FLERR,"Rigid fix must come before NPT/NPH fix"); } + // add gravity forces based on gravity vector from fix + + if (id_gravity) { + int ifix = modify->find_fix(id_gravity); + if (ifix < 0) error->all(FLERR,"Fix rigid/small cannot find fix gravity ID"); + if (strcmp(modify->fix[ifix]->style,"gravity") != 0) + error->all(FLERR,"Fix rigid/small gravity fix is invalid"); + int tmp; + gvec = (double *) modify->fix[ifix]->extract("gvec",tmp); + } + // timestep info dtv = update->dt; @@ -954,8 +976,20 @@ void FixRigidSmall::compute_forces_and_torques() tcm[2] += langextra[ibody][5]; } } -} + // add gravity force to COM of each body + + if (id_gravity) { + double mass; + for (ibody = 0; ibody < nlocal_body; ibody++) { + mass = body[ibody].mass; + fcm = body[ibody].fcm; + fcm[0] += gvec[0]*mass; + fcm[1] += gvec[1]*mass; + fcm[2] += gvec[2]*mass; + } + } +} /* ---------------------------------------------------------------------- */ diff --git a/src/RIGID/fix_rigid_small.h b/src/RIGID/fix_rigid_small.h index 0d1cdb5963..974f5711e2 100644 --- a/src/RIGID/fix_rigid_small.h +++ b/src/RIGID/fix_rigid_small.h @@ -164,6 +164,9 @@ class FixRigidSmall : public Fix { int dilate_group_bit; // mask for dilation group char *id_dilate; // group name to dilate + char *id_gravity; // ID of fix gravity command to add gravity forces + double *gvec; // ptr to gravity vector inside the fix + double p_current[3],p_target[3]; // molecules added on-the-fly as rigid bodies diff --git a/src/fix_gravity.cpp b/src/fix_gravity.cpp index 14ba913c01..c36d2230c5 100644 --- a/src/fix_gravity.cpp +++ b/src/fix_gravity.cpp @@ -37,7 +37,8 @@ enum{CONSTANT,EQUAL}; FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg), - mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL), xstr(NULL), ystr(NULL), zstr(NULL) + mstr(NULL), vstr(NULL), pstr(NULL), tstr(NULL), + xstr(NULL), ystr(NULL), zstr(NULL) { if (narg < 5) error->all(FLERR,"Illegal fix gravity command"); @@ -61,8 +62,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : mstyle = CONSTANT; } + int iarg; + if (strcmp(arg[4],"chute") == 0) { - if (narg != 6) error->all(FLERR,"Illegal fix gravity command"); + if (narg < 6) error->all(FLERR,"Illegal fix gravity command"); style = CHUTE; if (strstr(arg[5],"v_") == arg[5]) { int n = strlen(&arg[5][2]) + 1; @@ -73,9 +76,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : vert = force->numeric(FLERR,arg[5]); vstyle = CONSTANT; } + iarg = 6; } else if (strcmp(arg[4],"spherical") == 0) { - if (narg != 7) error->all(FLERR,"Illegal fix gravity command"); + if (narg < 7) error->all(FLERR,"Illegal fix gravity command"); style = SPHERICAL; if (strstr(arg[5],"v_") == arg[5]) { int n = strlen(&arg[5][2]) + 1; @@ -95,9 +99,10 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : theta = force->numeric(FLERR,arg[6]); tstyle = CONSTANT; } + iarg = 7; } else if (strcmp(arg[4],"vector") == 0) { - if (narg != 8) error->all(FLERR,"Illegal fix gravity command"); + if (narg < 8) error->all(FLERR,"Illegal fix gravity command"); style = VECTOR; if (strstr(arg[5],"v_") == arg[5]) { int n = strlen(&arg[5][2]) + 1; @@ -126,9 +131,23 @@ FixGravity::FixGravity(LAMMPS *lmp, int narg, char **arg) : zdir = force->numeric(FLERR,arg[7]); zstyle = CONSTANT; } + iarg = 8; } else error->all(FLERR,"Illegal fix gravity command"); + // optional keywords + + disable = 0; + + while (iarg < narg) { + if (strcmp(arg[iarg],"disable") == 0) { + disable = 1; + iarg++; + } else error->all(FLERR,"Illegal fix gravity command"); + } + + // initializations + degree2rad = MY_PI/180.0; time_origin = update->ntimestep; @@ -266,6 +285,12 @@ void FixGravity::post_force(int /*vflag*/) set_acceleration(); } + // just exit if application of force is disabled + + if (disable) return; + + // apply gravity force to each particle + double **x = atom->x; double **f = atom->f; double *rmass = atom->rmass; @@ -338,9 +363,9 @@ void FixGravity::set_acceleration() } } - xacc = magnitude*xgrav; - yacc = magnitude*ygrav; - zacc = magnitude*zgrav; + gvec[0] = xacc = magnitude*xgrav; + gvec[1] = yacc = magnitude*ygrav; + gvec[2] = zacc = magnitude*zgrav; } /* ---------------------------------------------------------------------- @@ -357,3 +382,16 @@ double FixGravity::compute_scalar() } return egrav_all; } + +/* ---------------------------------------------------------------------- + extract current gravity direction vector +------------------------------------------------------------------------- */ + +void *FixGravity::extract(const char *name, int &dim) +{ + if (strcmp(name,"gvec") == 0) { + dim = 1; + return (void *) gvec; + } + return NULL; +} diff --git a/src/fix_gravity.h b/src/fix_gravity.h index 1f18ee274b..ba6a61ee0d 100644 --- a/src/fix_gravity.h +++ b/src/fix_gravity.h @@ -36,9 +36,10 @@ class FixGravity : public Fix { virtual void post_force(int); virtual void post_force_respa(int, int, int); double compute_scalar(); + void *extract(const char *, int &); protected: - int style; + int style,disable; double magnitude; double vert,phi,theta; double xdir,ydir,zdir; @@ -46,6 +47,8 @@ class FixGravity : public Fix { double degree2rad; int ilevel_respa; int time_origin; + double gvec[3]; + int eflag; double egrav,egrav_all; From a11b886b5c39e42c51686e31945b1fdf5394f85f Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Thu, 19 Dec 2019 16:17:14 -0700 Subject: [PATCH 03/14] add warning if gravity is used incorrectly with overlapped rigid bodies --- src/RIGID/fix_rigid.cpp | 17 +++++++++++++++++ src/RIGID/fix_rigid_small.cpp | 17 +++++++++++++++++ 2 files changed, 34 insertions(+) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index 5e994810ac..d852a47381 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -728,6 +728,23 @@ void FixRigid::init() } } + // warn if body properties are read from inpfile + // and the gravity keyword is not set and a gravity fix exists + // this could mean body particles are overlapped + // and gravity is not applied correctly + + if (inpfile && !id_gravity) { + for (i = 0; i < modify->nfix; i++) { + if (strcmp(modify->fix[i]->style,"gravity") == 0) { + if (comm->me == 0) + error->warning(FLERR,"Gravity may not be correctly applied " + "to rigid bodies if they consist of " + "overlapped particles"); + break; + } + } + } + // error if npt,nph fix comes before rigid fix for (i = 0; i < modify->nfix; i++) { diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index 4d2776b9ff..cbdc877d09 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -574,6 +574,23 @@ void FixRigidSmall::init() } } + // warn if body properties are read from inpfile or a mol template file + // and the gravity keyword is not set and a gravity fix exists + // this could mean body particles are overlapped + // and gravity is not applied correctly + + if ((inpfile || onemols) && !id_gravity) { + for (i = 0; i < modify->nfix; i++) { + if (strcmp(modify->fix[i]->style,"gravity") == 0) { + if (comm->me == 0) + error->warning(FLERR,"Gravity may not be correctly applied " + "to rigid bodies if they consist of " + "overlapped particles"); + break; + } + } + } + // error if npt,nph fix comes before rigid fix for (i = 0; i < modify->nfix; i++) { From c8a53d560a55748744d6f7c5f468976843663232 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Johannes=20H=C3=B6rmann?= Date: Wed, 18 Dec 2019 12:09:21 +0100 Subject: [PATCH 04/14] Find parallel NetCDF with cmake for USER-NETCDF Conflicts: cmake/presets/forhlr2-gnu.cmake --- cmake/Modules/FindNetCDF.cmake | 6 ++- cmake/Modules/FindPNetCDF.cmake | 55 ++++++++++++++++++++++++ cmake/Modules/Packages/USER-NETCDF.cmake | 26 +++++++++-- 3 files changed, 82 insertions(+), 5 deletions(-) create mode 100644 cmake/Modules/FindPNetCDF.cmake diff --git a/cmake/Modules/FindNetCDF.cmake b/cmake/Modules/FindNetCDF.cmake index a28c959acf..3d4aed947d 100644 --- a/cmake/Modules/FindNetCDF.cmake +++ b/cmake/Modules/FindNetCDF.cmake @@ -46,10 +46,14 @@ endif() find_path (NETCDF_INCLUDE_DIR netcdf.h HINTS "${NETCDF_DIR}/include") mark_as_advanced (NETCDF_INCLUDE_DIR) + set (NETCDF_C_INCLUDE_DIRS ${NETCDF_INCLUDE_DIR}) +string(REGEX REPLACE "/include/?$" "/lib" + NETCDF_LIB_HINT ${NETCDF_INCLUDE_DIR}) + find_library (NETCDF_LIBRARY NAMES netcdf - HINTS "${NETCDF_DIR}/lib") + HINTS "${NETCDF_DIR}/lib" "${NETCDF_LIB_HINT}") mark_as_advanced (NETCDF_LIBRARY) set (NETCDF_C_LIBRARIES ${NETCDF_LIBRARY}) diff --git a/cmake/Modules/FindPNetCDF.cmake b/cmake/Modules/FindPNetCDF.cmake new file mode 100644 index 0000000000..0cc2c1faa8 --- /dev/null +++ b/cmake/Modules/FindPNetCDF.cmake @@ -0,0 +1,55 @@ +# source: https://ftp.space.dtu.dk/pub/Ioana/pism0.6.1-10/CMake/FindPNetCDF.cmake +# license: GPL v3 (https://ftp.space.dtu.dk/pub/Ioana/pism0.6.1-10/COPYING) +# +# - Find PNetCDF +# Find the native PNetCDF includes and library +# +# PNETCDF_INCLUDES - where to find netcdf.h, etc +# PNETCDF_LIBRARIES - Link these libraries when using NetCDF +# PNETCDF_FOUND - True if PNetCDF was found +# +# Normal usage would be: +# find_package (PNetCDF REQUIRED) +# target_link_libraries (uses_pnetcdf ${PNETCDF_LIBRARIES}) + +if (PNETCDF_INCLUDES AND PNETCDF_LIBRARIES) + # Already in cache, be silent + set (PNETCDF_FIND_QUIETLY TRUE) +endif (PNETCDF_INCLUDES AND PNETCDF_LIBRARIES) + +find_path (PNETCDF_INCLUDES pnetcdf.h + HINTS "${PNETCDF_ROOT}/include" "$ENV{PNETCDF_ROOT}/include") + +string(REGEX REPLACE "/include/?$" "/lib" + PNETCDF_LIB_HINT ${PNETCDF_INCLUDES}) + +find_library (PNETCDF_LIBRARIES + NAMES pnetcdf + HINTS ${PNETCDF_LIB_HINT}) + +if ((NOT PNETCDF_LIBRARIES) OR (NOT PNETCDF_INCLUDES)) + message(STATUS "Trying to find PNetCDF using LD_LIBRARY_PATH (we're desperate)...") + + file(TO_CMAKE_PATH "$ENV{LD_LIBRARY_PATH}" LD_LIBRARY_PATH) + + find_library(PNETCDF_LIBRARIES + NAMES pnetcdf + HINTS ${LD_LIBRARY_PATH}) + + if (PNETCDF_LIBRARIES) + get_filename_component(PNETCDF_LIB_DIR ${PNETCDF_LIBRARIES} PATH) + string(REGEX REPLACE "/lib/?$" "/include" + PNETCDF_H_HINT ${PNETCDF_LIB_DIR}) + + find_path (PNETCDF_INCLUDES pnetcdf.h + HINTS ${PNETCDF_H_HINT} + DOC "Path to pnetcdf.h") + endif() +endif() + +# handle the QUIETLY and REQUIRED arguments and set PNETCDF_FOUND to TRUE if +# all listed variables are TRUE +include (FindPackageHandleStandardArgs) +find_package_handle_standard_args (PNetCDF DEFAULT_MSG PNETCDF_LIBRARIES PNETCDF_INCLUDES) + +mark_as_advanced (PNETCDF_LIBRARIES PNETCDF_INCLUDES) diff --git a/cmake/Modules/Packages/USER-NETCDF.cmake b/cmake/Modules/Packages/USER-NETCDF.cmake index a90725bbbc..f60c046ab9 100644 --- a/cmake/Modules/Packages/USER-NETCDF.cmake +++ b/cmake/Modules/Packages/USER-NETCDF.cmake @@ -1,6 +1,24 @@ if(PKG_USER-NETCDF) - find_package(NetCDF REQUIRED) - include_directories(${NETCDF_INCLUDE_DIRS}) - list(APPEND LAMMPS_LINK_LIBS ${NETCDF_LIBRARIES}) - add_definitions(-DLMP_HAS_NETCDF -DNC_64BIT_DATA=0x0020) + # USER-NETCDF can use NetCDF, Parallel NetCDF (PNetCDF), or both. At least one necessary. + # NetCDF library enables dump sytle "netcdf", while PNetCDF enables dump style "netcdf/mpiio" + find_package(NetCDF) + if(NETCDF_FOUND) + find_package(PNetCDF) + else(NETCDF_FOUND) + find_package(PNetCDF REQUIRED) + endif(NETCDF_FOUND) + + if(NETCDF_FOUND) + include_directories(${NETCDF_INCLUDE_DIRS}) + list(APPEND LAMMPS_LINK_LIBS ${NETCDF_LIBRARIES}) + add_definitions(-DLMP_HAS_NETCDF) + endif(NETCDF_FOUND) + + if(PNETCDF_FOUND) + include_directories(${PNETCDF_INCLUDES}) + list(APPEND LAMMPS_LINK_LIBS ${PNETCDF_LIBRARIES}) + add_definitions(-DLMP_HAS_PNETCDF) + endif(PNETCDF_FOUND) + + add_definitions(-DNC_64BIT_DATA=0x0020) endif() From 3b76ab56f05c705fc71952a1709e97a1898b4ddc Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2019 12:09:50 -0500 Subject: [PATCH 05/14] port gravity changes to rigid fixes to USER-OMP package versions --- src/USER-OMP/fix_rigid_nh_omp.cpp | 13 +++++++++++++ src/USER-OMP/fix_rigid_omp.cpp | 13 +++++++++++++ src/USER-OMP/fix_rigid_small_omp.cpp | 15 +++++++++++++++ 3 files changed, 41 insertions(+) diff --git a/src/USER-OMP/fix_rigid_nh_omp.cpp b/src/USER-OMP/fix_rigid_nh_omp.cpp index d92a82a6bf..da512cb428 100644 --- a/src/USER-OMP/fix_rigid_nh_omp.cpp +++ b/src/USER-OMP/fix_rigid_nh_omp.cpp @@ -383,6 +383,19 @@ void FixRigidNHOMP::compute_forces_and_torques() torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; } + + // add gravity force to COM of each body + + if (id_gravity) { +#if defined(_OPENMP) +#pragma omp parallel for default(none) private(ibody) schedule(static) +#endif + for (ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] += gvec[0]*masstotal[ibody]; + fcm[ibody][1] += gvec[1]*masstotal[ibody]; + fcm[ibody][2] += gvec[2]*masstotal[ibody]; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/fix_rigid_omp.cpp b/src/USER-OMP/fix_rigid_omp.cpp index 20d3009d40..1cdf1018b2 100644 --- a/src/USER-OMP/fix_rigid_omp.cpp +++ b/src/USER-OMP/fix_rigid_omp.cpp @@ -256,6 +256,19 @@ void FixRigidOMP::compute_forces_and_torques() torque[ibody][1] = all[ibody][4] + langextra[ibody][4]; torque[ibody][2] = all[ibody][5] + langextra[ibody][5]; } + + // add gravity force to COM of each body + + if (id_gravity) { +#if defined(_OPENMP) +#pragma omp parallel for default(none) private(ibody) schedule(static) +#endif + for (int ibody = 0; ibody < nbody; ibody++) { + fcm[ibody][0] += gvec[0]*masstotal[ibody]; + fcm[ibody][1] += gvec[1]*masstotal[ibody]; + fcm[ibody][2] += gvec[2]*masstotal[ibody]; + } + } } /* ---------------------------------------------------------------------- */ diff --git a/src/USER-OMP/fix_rigid_small_omp.cpp b/src/USER-OMP/fix_rigid_small_omp.cpp index 705f0e6bd0..bdd10b0ca7 100644 --- a/src/USER-OMP/fix_rigid_small_omp.cpp +++ b/src/USER-OMP/fix_rigid_small_omp.cpp @@ -196,6 +196,21 @@ void FixRigidSmallOMP::compute_forces_and_torques() tcm[2] += langextra[ibody][5]; } } + + // add gravity force to COM of each body + + if (id_gravity) { +#if defined(_OPENMP) +#pragma omp parallel for default(none) private(ibody) schedule(static) +#endif + for (int ibody = 0; ibody < nbody; ibody++) { + double * _noalias const fcm = body[ibody].fcm; + const double mass = body[ibody].mass; + fcm[0] += gvec[0]*mass; + fcm[1] += gvec[1]*mass; + fcm[2] += gvec[2]*mass; + } + } } /* ---------------------------------------------------------------------- */ From a88a00dbbb86fdcc2b4ed7bde34f50180bbaa89d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Fri, 20 Dec 2019 15:08:51 -0500 Subject: [PATCH 06/14] remove trailing whitespace --- src/RIGID/fix_rigid.cpp | 6 +++--- src/RIGID/fix_rigid_small.cpp | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/RIGID/fix_rigid.cpp b/src/RIGID/fix_rigid.cpp index d852a47381..101b94c2a4 100644 --- a/src/RIGID/fix_rigid.cpp +++ b/src/RIGID/fix_rigid.cpp @@ -52,7 +52,7 @@ FixRigid::FixRigid(LAMMPS *lmp, int narg, char **arg) : torque(NULL), quat(NULL), imagebody(NULL), fflag(NULL), tflag(NULL), langextra(NULL), sum(NULL), all(NULL), remapflag(NULL), xcmimage(NULL), eflags(NULL), orient(NULL), - dorient(NULL), id_dilate(NULL), id_gravity(NULL), random(NULL), + dorient(NULL), id_dilate(NULL), id_gravity(NULL), random(NULL), avec_ellipsoid(NULL), avec_line(NULL), avec_tri(NULL) { int i,ibody; @@ -736,7 +736,7 @@ void FixRigid::init() if (inpfile && !id_gravity) { for (i = 0; i < modify->nfix; i++) { if (strcmp(modify->fix[i]->style,"gravity") == 0) { - if (comm->me == 0) + if (comm->me == 0) error->warning(FLERR,"Gravity may not be correctly applied " "to rigid bodies if they consist of " "overlapped particles"); @@ -763,7 +763,7 @@ void FixRigid::init() if (id_gravity) { int ifix = modify->find_fix(id_gravity); if (ifix < 0) error->all(FLERR,"Fix rigid cannot find fix gravity ID"); - if (strcmp(modify->fix[ifix]->style,"gravity") != 0) + if (strcmp(modify->fix[ifix]->style,"gravity") != 0) error->all(FLERR,"Fix rigid gravity fix is invalid"); int tmp; gvec = (double *) modify->fix[ifix]->extract("gvec",tmp); diff --git a/src/RIGID/fix_rigid_small.cpp b/src/RIGID/fix_rigid_small.cpp index cbdc877d09..33105418ff 100644 --- a/src/RIGID/fix_rigid_small.cpp +++ b/src/RIGID/fix_rigid_small.cpp @@ -582,7 +582,7 @@ void FixRigidSmall::init() if ((inpfile || onemols) && !id_gravity) { for (i = 0; i < modify->nfix; i++) { if (strcmp(modify->fix[i]->style,"gravity") == 0) { - if (comm->me == 0) + if (comm->me == 0) error->warning(FLERR,"Gravity may not be correctly applied " "to rigid bodies if they consist of " "overlapped particles"); @@ -608,7 +608,7 @@ void FixRigidSmall::init() if (id_gravity) { int ifix = modify->find_fix(id_gravity); if (ifix < 0) error->all(FLERR,"Fix rigid/small cannot find fix gravity ID"); - if (strcmp(modify->fix[ifix]->style,"gravity") != 0) + if (strcmp(modify->fix[ifix]->style,"gravity") != 0) error->all(FLERR,"Fix rigid/small gravity fix is invalid"); int tmp; gvec = (double *) modify->fix[ifix]->extract("gvec",tmp); From 4e7bcee8e386dd6c467aa8cd9ef04214a227be22 Mon Sep 17 00:00:00 2001 From: jotelha Date: Sun, 22 Dec 2019 13:24:34 +0100 Subject: [PATCH 07/14] Update cmake/Modules/FindNetCDF.cmake Co-Authored-By: Christoph Junghans --- cmake/Modules/FindNetCDF.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/FindNetCDF.cmake b/cmake/Modules/FindNetCDF.cmake index 3d4aed947d..4181e3f027 100644 --- a/cmake/Modules/FindNetCDF.cmake +++ b/cmake/Modules/FindNetCDF.cmake @@ -53,7 +53,7 @@ string(REGEX REPLACE "/include/?$" "/lib" NETCDF_LIB_HINT ${NETCDF_INCLUDE_DIR}) find_library (NETCDF_LIBRARY NAMES netcdf - HINTS "${NETCDF_DIR}/lib" "${NETCDF_LIB_HINT}") + HINTS "${NETCDF_DIR}" "${NETCDF_LIB_HINT}" PATH_SUFFIXES lib lib64) mark_as_advanced (NETCDF_LIBRARY) set (NETCDF_C_LIBRARIES ${NETCDF_LIBRARY}) From 0c7d6a01e8f93669f94170acbc335c128058f7fe Mon Sep 17 00:00:00 2001 From: jotelha Date: Sun, 22 Dec 2019 13:25:08 +0100 Subject: [PATCH 08/14] Update cmake/Modules/FindNetCDF.cmake Co-Authored-By: Christoph Junghans --- cmake/Modules/FindNetCDF.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/FindNetCDF.cmake b/cmake/Modules/FindNetCDF.cmake index 4181e3f027..2a992b6b3b 100644 --- a/cmake/Modules/FindNetCDF.cmake +++ b/cmake/Modules/FindNetCDF.cmake @@ -49,7 +49,7 @@ mark_as_advanced (NETCDF_INCLUDE_DIR) set (NETCDF_C_INCLUDE_DIRS ${NETCDF_INCLUDE_DIR}) -string(REGEX REPLACE "/include/?$" "/lib" +string(REGEX REPLACE "/include/?$" "" NETCDF_LIB_HINT ${NETCDF_INCLUDE_DIR}) find_library (NETCDF_LIBRARY NAMES netcdf From 3f24144abd499042720c8d1f200a10e843f5b69a Mon Sep 17 00:00:00 2001 From: jotelha Date: Sun, 22 Dec 2019 13:25:24 +0100 Subject: [PATCH 09/14] Update cmake/Modules/FindPNetCDF.cmake Co-Authored-By: Christoph Junghans --- cmake/Modules/FindPNetCDF.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/FindPNetCDF.cmake b/cmake/Modules/FindPNetCDF.cmake index 0cc2c1faa8..e95ed4e559 100644 --- a/cmake/Modules/FindPNetCDF.cmake +++ b/cmake/Modules/FindPNetCDF.cmake @@ -38,7 +38,7 @@ if ((NOT PNETCDF_LIBRARIES) OR (NOT PNETCDF_INCLUDES)) if (PNETCDF_LIBRARIES) get_filename_component(PNETCDF_LIB_DIR ${PNETCDF_LIBRARIES} PATH) - string(REGEX REPLACE "/lib/?$" "/include" + string(REGEX REPLACE "/(lib|lib64)/?$" "/include" PNETCDF_H_HINT ${PNETCDF_LIB_DIR}) find_path (PNETCDF_INCLUDES pnetcdf.h From 6aa4f4caf644398b17722f546be19dd91030b0f3 Mon Sep 17 00:00:00 2001 From: "Dan S. Bolintineanu" Date: Mon, 6 Jan 2020 21:18:00 -0700 Subject: [PATCH 10/14] Added example script and input data file showing benefits of new fix gravity and fix rigid options --- examples/rigid/in.rigid.gravity | 62 +++++++++++++++++++++++++++++++++ examples/rigid/molecule.data | 41 ++++++++++++++++++++++ 2 files changed, 103 insertions(+) create mode 100644 examples/rigid/in.rigid.gravity create mode 100644 examples/rigid/molecule.data diff --git a/examples/rigid/in.rigid.gravity b/examples/rigid/in.rigid.gravity new file mode 100644 index 0000000000..22abb0740e --- /dev/null +++ b/examples/rigid/in.rigid.gravity @@ -0,0 +1,62 @@ +#Pour composite granular particles on flat wall + +newton on +atom_style sphere +atom_modify map array sort 0 0 + +thermo_modify flush yes +units si + +variable minrad equal 0.5 +variable maxrad equal 1.4 + +variable skin equal 0.3*${maxrad} + +boundary p p f +region reg block 0 20 0 20 0 200 units box +create_box 1 reg + +fix prop all property/atom mol ghost yes + +variable dumpfreq equal 1000 +variable logfreq equal 1000 + +pair_style gran/hooke/history 4e5 NULL 1e2 NULL 0.5 0 +pair_coeff * * + +timestep 0.0001 + +group particles type 1 +atom_modify first particles + +neighbor ${skin} bin +group rigid type 1 +neigh_modify every 1 delay 0 check yes exclude molecule/intra all + +thermo ${logfreq} +thermo_style custom step cpu atoms ke +thermo_modify flush yes lost warn + +comm_modify vel yes cutoff 3 + +molecule mymol molecule.data +region pourreg block 5 15 5 15 80 100 side in units box + +#Note: in versions prior to 1/2020, the 'disable' keyword to fix/gravity +# and the 'gravity' keyword to fix rigid/small were not available. +# These settings produce undesirable behavior, where gravity can induce +# torque on rigid bodies. +#fix gravfix all gravity 9.8 vector 0 0 -1 #disable +#fix rigidfix all rigid/small molecule mol mymol #gravity gravfix + +#The correct behavior is recovered with the following settings: +fix gravfix all gravity 9.8 vector 0 0 -1 disable +fix rigidfix all rigid/small molecule mol mymol gravity gravfix +fix pourfix all pour 5 0 1234 region pourreg mol mymol rigid rigidfix + +fix zwall all wall/gran hooke/history 4000.0 NULL 100.0 NULL 0.5 0 zplane 0.1 NULL + +dump 1 all custom 1000 molecule_pour.dump id type mass radius x y z fx fy fz + +run 100000 + diff --git a/examples/rigid/molecule.data b/examples/rigid/molecule.data new file mode 100644 index 0000000000..e0677b0f63 --- /dev/null +++ b/examples/rigid/molecule.data @@ -0,0 +1,41 @@ +LAMMPS data file created for rigid body molecule template + +5 atoms + +2.3388800000000005 mass + +6.002239704473936 4.99 4.989999999999999 com + +116.79265620480001 144.26721336320003 144.26721336320006 -70.05220681600004 -70.05220681600002 -58.238345888000005 inertia + +Coords + +1 5 5 5 +2 5.1 5.0 5.0 +3 5.2 5.0 5.0 +4 6.2 5.0 5.0 +5 7.2 5.0 5.0 + +Types + +1 1 +2 1 +3 1 +4 1 +5 1 + +Diameters + +1 1.0 +2 0.9 +3 1.2 +4 1.2 +5 1.0 + +Masses + +1 0.5235987755982988 +2 0.3817035074111599 +3 0.9047786842338602 +4 0.9047786842338602 +5 0.5235987755982988 From 7bc8c8e9d8e200a5303377836a94b1106dac096f Mon Sep 17 00:00:00 2001 From: "Dan S. Bolintineanu" Date: Mon, 6 Jan 2020 21:26:34 -0700 Subject: [PATCH 11/14] Minor tweaks to in.rigid.gravity example --- examples/rigid/in.rigid.gravity | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/examples/rigid/in.rigid.gravity b/examples/rigid/in.rigid.gravity index 22abb0740e..1094681f18 100644 --- a/examples/rigid/in.rigid.gravity +++ b/examples/rigid/in.rigid.gravity @@ -52,11 +52,12 @@ region pourreg block 5 15 5 15 80 100 side in units box #The correct behavior is recovered with the following settings: fix gravfix all gravity 9.8 vector 0 0 -1 disable fix rigidfix all rigid/small molecule mol mymol gravity gravfix + fix pourfix all pour 5 0 1234 region pourreg mol mymol rigid rigidfix fix zwall all wall/gran hooke/history 4000.0 NULL 100.0 NULL 0.5 0 zplane 0.1 NULL -dump 1 all custom 1000 molecule_pour.dump id type mass radius x y z fx fy fz +#dump 1 all custom 1000 molecule_pour.dump id type mass radius x y z fx fy fz run 100000 From 24ef36dd4dc5671546dc74335ef54644c51c88df Mon Sep 17 00:00:00 2001 From: jotelha Date: Tue, 7 Jan 2020 13:25:54 +0100 Subject: [PATCH 12/14] Update cmake/Modules/FindPNetCDF.cmake Co-Authored-By: Christoph Junghans --- cmake/Modules/FindPNetCDF.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/FindPNetCDF.cmake b/cmake/Modules/FindPNetCDF.cmake index e95ed4e559..8778335f62 100644 --- a/cmake/Modules/FindPNetCDF.cmake +++ b/cmake/Modules/FindPNetCDF.cmake @@ -20,7 +20,7 @@ endif (PNETCDF_INCLUDES AND PNETCDF_LIBRARIES) find_path (PNETCDF_INCLUDES pnetcdf.h HINTS "${PNETCDF_ROOT}/include" "$ENV{PNETCDF_ROOT}/include") -string(REGEX REPLACE "/include/?$" "/lib" +string(REGEX REPLACE "/include/?$" "" PNETCDF_LIB_HINT ${PNETCDF_INCLUDES}) find_library (PNETCDF_LIBRARIES From 46584d45201f43d42e9fc396db6f47d38907881b Mon Sep 17 00:00:00 2001 From: jotelha Date: Tue, 7 Jan 2020 13:26:01 +0100 Subject: [PATCH 13/14] Update cmake/Modules/FindPNetCDF.cmake Co-Authored-By: Christoph Junghans --- cmake/Modules/FindPNetCDF.cmake | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmake/Modules/FindPNetCDF.cmake b/cmake/Modules/FindPNetCDF.cmake index 8778335f62..bc3a5f9538 100644 --- a/cmake/Modules/FindPNetCDF.cmake +++ b/cmake/Modules/FindPNetCDF.cmake @@ -25,7 +25,7 @@ string(REGEX REPLACE "/include/?$" "" find_library (PNETCDF_LIBRARIES NAMES pnetcdf - HINTS ${PNETCDF_LIB_HINT}) + HINTS ${PNETCDF_LIB_HINT} PATH_SUFFIXES lib lib64) if ((NOT PNETCDF_LIBRARIES) OR (NOT PNETCDF_INCLUDES)) message(STATUS "Trying to find PNetCDF using LD_LIBRARY_PATH (we're desperate)...") From f1a23b1ea24187b3ab22f52d8a41d70a57cf1e9f Mon Sep 17 00:00:00 2001 From: tabedzki2 <59632999+tabedzki2@users.noreply.github.com> Date: Tue, 7 Jan 2020 21:47:43 -0500 Subject: [PATCH 14/14] Updated Makefile.stampede: replacement options icc The default options for Makefile.stampede did not compile. They had to be updated to include the `q` replacement options. --- src/MAKE/MACHINES/Makefile.stampede | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/MAKE/MACHINES/Makefile.stampede b/src/MAKE/MACHINES/Makefile.stampede index 69da22c9c9..ecd0810d1d 100644 --- a/src/MAKE/MACHINES/Makefile.stampede +++ b/src/MAKE/MACHINES/Makefile.stampede @@ -6,13 +6,13 @@ SHELL = /bin/sh # compiler/linker settings # specify flags and libraries needed for your compiler -CC = mpicc -openmp -DLMP_INTEL_OFFLOAD -DLAMMPS_MEMALIGN=64 -MIC_OPT = -offload-option,mic,compiler,"-fp-model fast=2 -mGLOB_default_function_attrs=\"gather_scatter_loop_unroll=4\"" -CCFLAGS = -O3 -xhost -fp-model precise -restrict -override-limits $(MIC_OPT) +CC = mpicc -qopenmp -DLMP_INTEL_OFFLOAD -DLAMMPS_MEMALIGN=64 +MIC_OPT = -qoffload-option,mic,compiler,"-fp-model fast=2 -mGLOB_default_function_attrs=\"gather_scatter_loop_unroll=4\"" +CCFLAGS = -O3 -xhost -fp-model precise -restrict -qoverride-limits $(MIC_OPT) SHFLAGS = -fPIC DEPFLAGS = -M -LINK = mpicc -openmp +LINK = mpicc -qopenmp LINKFLAGS = -O3 -xhost LIB = -ltbbmalloc SIZE = size