From 6abbbdba6c98f9ba77b56dc5820053dc6ec01310 Mon Sep 17 00:00:00 2001
From: jtclemm <jtclemm@sandia.gov>
Date: Thu, 25 Jul 2024 17:10:23 -0600
Subject: [PATCH] Patching vremap in SPH

---
 src/SPH/fix_sph.cpp | 42 +++++++++++++++++++++++++++++++++++++++---
 src/SPH/fix_sph.h   |  2 ++
 2 files changed, 41 insertions(+), 3 deletions(-)

diff --git a/src/SPH/fix_sph.cpp b/src/SPH/fix_sph.cpp
index 96479be928..876d7ac4fb 100644
--- a/src/SPH/fix_sph.cpp
+++ b/src/SPH/fix_sph.cpp
@@ -61,6 +61,10 @@ void FixSPH::init()
 
 void FixSPH::setup_pre_force(int /*vflag*/)
 {
+  remap_v_flag = domain->deform_vremap;
+  if (remap_v_flag && (!comm->ghost_velocity))
+    error->all(FLERR, "Fix sph requires ghost atoms store velocity when deforming with remap v");
+
   // set vest equal to v
   double **v = atom->v;
   double **vest = atom->vest;
@@ -119,9 +123,16 @@ void FixSPH::initial_integrate(int /*vflag*/)
       rho[i] += dtf * drho[i]; // ... and density
 
       // extrapolate velocity for use with velocity-dependent potentials, e.g. SPH
-      vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0];
-      vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1];
-      vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2];
+      // if velocities are remapped, perform this extrapolation after communication
+      if (remap_v_flag) {
+        vest[i][0] = dtfm * f[i][0];
+        vest[i][1] = dtfm * f[i][1];
+        vest[i][2] = dtfm * f[i][2];
+      } else {
+        vest[i][0] = v[i][0] + 2.0 * dtfm * f[i][0];
+        vest[i][1] = v[i][1] + 2.0 * dtfm * f[i][1];
+        vest[i][2] = v[i][2] + 2.0 * dtfm * f[i][2];
+      }
 
       v[i][0] += dtfm * f[i][0];
       v[i][1] += dtfm * f[i][1];
@@ -136,6 +147,31 @@ void FixSPH::initial_integrate(int /*vflag*/)
 
 /* ---------------------------------------------------------------------- */
 
+void FixSPH::pre_force(int /*vflag*/)
+{
+  // if velocities are remapped, calculate estimates here
+  // note that vest currently stores dtfm * force
+  if (!remap_v_flag) return;
+
+  double **v = atom->v;
+  double **vest = atom->vest;
+  int *mask = atom->mask;
+  int nlocal = atom->nlocal;
+  if (igroup == atom->firstgroup)
+    nlocal = atom->nfirst;
+
+  int nall = nlocal + atom->nghost;
+  for (int i = 0; i < nall; i++) {
+    if (mask[i] & groupbit) {
+      vest[i][0] += v[i][0];
+      vest[i][1] += v[i][1];
+      vest[i][2] += v[i][2];
+    }
+  }
+}
+
+/* ---------------------------------------------------------------------- */
+
 void FixSPH::final_integrate()
 {
   // update v, rho, and e of atoms in group
diff --git a/src/SPH/fix_sph.h b/src/SPH/fix_sph.h
index 415e38bd21..7676844dd9 100644
--- a/src/SPH/fix_sph.h
+++ b/src/SPH/fix_sph.h
@@ -30,6 +30,7 @@ class FixSPH : public Fix {
   int setmask() override;
   void init() override;
   void setup_pre_force(int) override;
+  void pre_force(int) override;
   void initial_integrate(int) override;
   void final_integrate() override;
   void reset_dt() override;
@@ -41,6 +42,7 @@ class FixSPH : public Fix {
   double dtv, dtf;
   double *step_respa;
   int mass_require;
+  int remap_v_flag;
 
   class Pair *pair;
 };