Merge pull request lammps#4243 from jtclemm/sph-update

Cleaning up and fixing bug in SPH package

doc/src/fix_meso_move.rst

@@ -247,6 +247,11 @@ defined by the :doc:`atom_style sph <atom_style>` command.
 
 All particles in the group must be mesoscopic SPH/SDPD particles.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 

doc/src/fix_mvv_dpd.rst

@@ -97,6 +97,11 @@ These fixes are part of the DPD-MESO package. They are only enabled if
 LAMMPS was built with that package. See the :doc:`Build package
 <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 

doc/src/fix_rigid_meso.rst

@@ -353,6 +353,11 @@ defined by the :doc:`atom_style sph <atom_style>` command.
 
 All particles in the group must be mesoscopic SPH/SDPD particles.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 

doc/src/fix_smd_integrate_tlsph.rst

@@ -53,6 +53,11 @@ Restrictions
 This fix is part of the MACHDYN package.  It is only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 

doc/src/fix_smd_integrate_ulsph.rst

@@ -61,6 +61,11 @@ Restrictions
 This fix is part of the MACHDYN package.  It is only enabled if
 LAMMPS was built with that package.  See the :doc:`Build package <Build_package>` page for more info.
 
+versionchanged:: TBD
+
+This fix is incompatible with deformation controls that remap velocity,
+for instance the *remap v* option of :doc:`fix deform <fix_deform>`.
+
 Related commands
 """"""""""""""""
 

src/DPD-MESO/fix_mvv_dpd.cpp

@@ -24,6 +24,7 @@
 #include "fix_mvv_dpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
@@ -65,6 +66,11 @@ void FixMvvDPD::init()
   if (!atom->vest_flag)
     error->all(FLERR,"Fix mvv/dpd requires atom attribute vest e.g. from atom style mdpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples of patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/dpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^mdpd",0) && !force->pair_match("^dpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!(force->pair_match("^mdpd",0,1) || force->pair_match("^dpd",0),1)) {

src/DPD-MESO/fix_mvv_edpd.cpp

@@ -33,6 +33,7 @@
 #include "fix_mvv_edpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
@@ -73,6 +74,11 @@ void FixMvvEDPD::init()
 {
   if (!atom->edpd_flag) error->all(FLERR,"Fix mvv/edpd requires atom style edpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples of patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/edpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^edpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!force->pair_match("^edpd",0,1)) {

src/DPD-MESO/fix_mvv_tdpd.cpp

@@ -29,6 +29,7 @@
 #include "fix_mvv_tdpd.h"
 
 #include "atom.h"
+#include "domain.h"
 #include "error.h"
 #include "force.h"
 #include "update.h"
@@ -71,6 +72,11 @@ void FixMvvTDPD::init()
 {
   if (!atom->tdpd_flag) error->all(FLERR,"Fix mvv/tdpd requires atom style tdpd");
 
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples of patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix mvv/tdpd cannot be used with velocity remapping");
+
   if (!force->pair_match("^tdpd",0)) {
     if (force->pair_match("^hybrid",0)) {
       if (!force->pair_match("^tdpd",0,1)) {

src/DPD-SMOOTH/fix_meso_move.cpp

@@ -350,7 +350,14 @@ void FixMesoMove::init () {
 }
 
 void FixMesoMove::setup_pre_force (int /*vflag*/) {
+
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples of patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix meso/move cannot be used with velocity remapping");
+
   // set vest equal to v
+
   double **v = atom->v;
   double **vest = atom->vest;
   int *mask = atom->mask;

src/DPD-SMOOTH/fix_rigid_meso.cpp

@@ -29,11 +29,12 @@
 ------------------------------------------------------------------------- */
 
 #include "fix_rigid_meso.h"
-#include "math_extra.h"
+
 #include "atom.h"
 #include "domain.h"
-#include "memory.h"
 #include "error.h"
+#include "math_extra.h"
+#include "memory.h"
 
 using namespace LAMMPS_NS;
 using namespace FixConst;
@@ -92,6 +93,11 @@ void FixRigidMeso::setup (int vflag) {
     conjqm[ibody][2] *= 2.0;
     conjqm[ibody][3] *= 2.0;
   }
+
+  // Cannot use vremap since its effects aren't propagated to vest
+  //   see RHEO or SPH packages for examples of patches
+  if (domain->deform_vremap)
+    error->all(FLERR, "Fix rigid/meso cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------

src/MACHDYN/fix_smd_integrate_tlsph.cpp

@@ -24,15 +24,18 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_smd_integrate_tlsph.h"
-#include <cmath>
-#include <cstring>
-#include <Eigen/Eigen>
+
 #include "atom.h"
-#include "force.h"
-#include "update.h"
+#include "comm.h"
+#include "domain.h"
 #include "error.h"
+#include "force.h"
 #include "pair.h"
-#include "comm.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+#include <Eigen/Eigen>
 
 using namespace Eigen;
 using namespace LAMMPS_NS;
@@ -115,6 +118,11 @@ void FixSMDIntegrateTlsph::init() {
         dtv = update->dt;
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
+
+        // Cannot use vremap since its effects aren't propagated to vest
+        //   see RHEO or SPH packages for examples of patches
+        if (domain->deform_vremap)
+          error->all(FLERR, "Fix smd/integrate_tlsph cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------

src/MACHDYN/fix_smd_integrate_ulsph.cpp

@@ -24,15 +24,18 @@
  ------------------------------------------------------------------------- */
 
 #include "fix_smd_integrate_ulsph.h"
-#include <cmath>
-#include <cstring>
-#include <Eigen/Eigen>
+
 #include "atom.h"
 #include "comm.h"
-#include "force.h"
-#include "update.h"
+#include "domain.h"
 #include "error.h"
+#include "force.h"
 #include "pair.h"
+#include "update.h"
+
+#include <cmath>
+#include <cstring>
+#include <Eigen/Eigen>
 
 using namespace Eigen;
 using namespace LAMMPS_NS;
@@ -144,6 +147,11 @@ void FixSMDIntegrateUlsph::init() {
         dtv = update->dt;
         dtf = 0.5 * update->dt * force->ftm2v;
         vlimitsq = vlimit * vlimit;
+
+        // Cannot use vremap since its effects aren't propagated to vest
+        //   see RHEO or SPH packages for examples of patches
+        if (domain->deform_vremap)
+          error->all(FLERR, "Fix smd/integrate_ulsph cannot be used with velocity remapping");
 }
 
 /* ----------------------------------------------------------------------

src/SPH/compute_sph_e_atom.cpp

@@ -13,13 +13,15 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_sph_e_atom.h"
-#include <cstring>
+
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -31,7 +33,7 @@ ComputeSPHEAtom::ComputeSPHEAtom(LAMMPS *lmp, int narg, char **arg) :
   if (narg != 3)
     error->all(FLERR,"Number of arguments for compute sph/e/atom command != 3");
   if (atom->esph_flag != 1)
-    error->all(FLERR,"Compute sph/e/atom command requires atom_style sph)");
+    error->all(FLERR,"Compute sph/e/atom requires atom attribute energy, e.g. in atom_style sph");
 
   peratom_flag = 1;
   size_peratom_cols = 0;
@@ -51,12 +53,11 @@ ComputeSPHEAtom::~ComputeSPHEAtom()
 
 void ComputeSPHEAtom::init()
 {
-
   int count = 0;
   for (int i = 0; i < modify->ncompute; i++)
-    if (strcmp(modify->compute[i]->style,"evector/atom") == 0) count++;
+    if (strcmp(modify->compute[i]->style,"sph/e/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one compute evector/atom");
+    error->warning(FLERR,"More than one compute sph/e/atom");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -78,14 +79,13 @@ void ComputeSPHEAtom::compute_peratom()
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-              evector[i] = esph[i];
-      }
-      else {
-              evector[i] = 0.0;
-      }
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      evector[i] = esph[i];
+    } else {
+      evector[i] = 0.0;
     }
+  }
 }
 
 /* ----------------------------------------------------------------------

src/SPH/compute_sph_rho_atom.cpp

@@ -13,13 +13,15 @@
 ------------------------------------------------------------------------- */
 
 #include "compute_sph_rho_atom.h"
-#include <cstring>
+
 #include "atom.h"
-#include "update.h"
-#include "modify.h"
 #include "comm.h"
-#include "memory.h"
 #include "error.h"
+#include "memory.h"
+#include "modify.h"
+#include "update.h"
+
+#include <cstring>
 
 using namespace LAMMPS_NS;
 
@@ -30,8 +32,7 @@ ComputeSPHRhoAtom::ComputeSPHRhoAtom(LAMMPS *lmp, int narg, char **arg) :
 {
   if (narg != 3) error->all(FLERR,"Illegal compute sph/rho/atom command");
   if (atom->rho_flag != 1)
-    error->all(FLERR,"Compute sph/rho/atom command requires atom_style sph");
-
+    error->all(FLERR, "Compute sph/rho/atom requires atom attribute density, e.g. in atom_style sph");
   peratom_flag = 1;
   size_peratom_cols = 0;
 
@@ -50,12 +51,11 @@ ComputeSPHRhoAtom::~ComputeSPHRhoAtom()
 
 void ComputeSPHRhoAtom::init()
 {
-
   int count = 0;
   for (int i = 0; i < modify->ncompute; i++)
-    if (strcmp(modify->compute[i]->style,"rhoVector/atom") == 0) count++;
+    if (strcmp(modify->compute[i]->style,"sph/rho/atom") == 0) count++;
   if (count > 1 && comm->me == 0)
-    error->warning(FLERR,"More than one compute rhoVector/atom");
+    error->warning(FLERR,"More than one compute sph/rho/atom");
 }
 
 /* ---------------------------------------------------------------------- */
@@ -73,20 +73,17 @@ void ComputeSPHRhoAtom::compute_peratom()
     vector_atom = rhoVector;
   }
 
-  // compute kinetic energy for each atom in group
-
   double *rho = atom->rho;
   int *mask = atom->mask;
   int nlocal = atom->nlocal;
 
-    for (int i = 0; i < nlocal; i++) {
-      if (mask[i] & groupbit) {
-              rhoVector[i] = rho[i];
-      }
-      else {
-              rhoVector[i] = 0.0;
-      }
+  for (int i = 0; i < nlocal; i++) {
+    if (mask[i] & groupbit) {
+      rhoVector[i] = rho[i];
+    } else {
+      rhoVector[i] = 0.0;
     }
+  }
 }
 
 /* ----------------------------------------------------------------------