"Debug"

abkein · Jan 16, 2025 · 4658905 · 4658905
1 parent 8a4ae4d
commit 4658905
Show file tree

Hide file tree

Showing 5 changed files with 85 additions and 49 deletions.
diff --git a/src/compute_cluster_enthropy.cpp b/src/compute_cluster_enthropy.cpp
@@ -55,44 +55,32 @@ ComputeClusterEnthropy::ComputeClusterEnthropy(LAMMPS* lmp, int narg, char** arg
     error->all(FLERR, "{}: Cannot find compute with style 'cluster/size' with id: {}", style, arg[3]);
   }
 
+  // Get entropy/atom compute
+  compute_entropy_atom = lmp->modify->get_compute_by_id(arg[4]);
+  if (compute_entropy_atom == nullptr) { error->all(FLERR, "{}: Cannot find compute with id '{}'", style, arg[4]); }
+
   size_cutoff = compute_cluster_size->get_size_cutoff();
-  if ((narg >= 4) && (::strcmp(arg[4], "inherit") != 0)) {
-    int t_size_cutoff = utils::inumeric(FLERR, arg[4], true, lmp);
+  if ((narg >= 6) && (::strcmp(arg[5], "inherit") != 0)) {
+    int t_size_cutoff = utils::inumeric(FLERR, arg[5], true, lmp);
     if (t_size_cutoff < 1) { error->all(FLERR, "{}: size_cutoff must be greater than 0", style); }
     if (t_size_cutoff > size_cutoff) { error->all(FLERR, "{}: size_cutoff cannot be greater than it of compute sizecluster", style); }
   }
 
-  // Get ke/atom compute
-  compute_ke_atom = lmp->modify->get_compute_by_id(arg[5]);
-  if (compute_ke_atom == nullptr) { error->all(FLERR, "{}: Cannot find compute with style id '{}'", style, arg[5]); }
-
-  // Get pe/atom compute
-  compute_pe_atom = lmp->modify->get_compute_by_id(arg[6]);
-  if (compute_pe_atom == nullptr) { error->all(FLERR, "{}: Cannot find compute with id '{}'", style, arg[6]); }
-
-  // Get entropy/atom compute
-  compute_entropy_atom = lmp->modify->get_compute_by_id(arg[7]);
-  if (compute_entropy_atom == nullptr) { error->all(FLERR, "{}: Cannot find compute with id '{}'", style, arg[7]); }
-
-  // Get temp/cluster compute
-  compute_temp_cluster = lmp->modify->get_compute_by_id(arg[8]);
-  if (compute_temp_cluster == nullptr) { error->all(FLERR, "{}: Cannot find compute with id '{}'", style, arg[8]); }
-
   size_local_rows = size_cutoff + 1;
-  memory->create(local_temp, size_local_rows + 1, "compute:cluster/enthropy:temp");
-  vector_local = local_temp;
+  local_enth.create(memory, size_local_rows + 1, "compute:cluster/enthropy:temp");
+  vector_local = local_enth.data();
 
   size_vector  = size_cutoff + 1;
-  memory->create(temp, size_vector + 1, "compute:cluster/enthropy:temp");
-  vector = temp;
+  enth.create(memory, size_vector + 1, "compute:cluster/enthropy:temp");
+  vector = enth.data();
 }
 
 /* ---------------------------------------------------------------------- */
 
 ComputeClusterEnthropy::~ComputeClusterEnthropy() noexcept(true)
 {
-  memory->destroy(local_temp);
-  memory->destroy(temp);
+  local_enth.destroy(memory);
+  enth.destroy(memory);
 }
 
 /* ---------------------------------------------------------------------- */
@@ -110,10 +98,10 @@ void ComputeClusterEnthropy::compute_vector()
 
   compute_local();
 
-  ::MPI_Allreduce(local_temp, temp, size_vector, MPI_DOUBLE, MPI_SUM, world);
+  ::MPI_Allreduce(local_enth.data(), enth.data(), size_vector, MPI_DOUBLE, MPI_SUM, world);
 
   const double* dist = compute_cluster_size->vector;
-  for (int i = 0; i < size_vector; ++i) { temp[i] /= (dist[i] * i - 1) * domain->dimension; }
+  for (int i = 0; i < size_vector; ++i) { enth[i] /= dist[i]; }
 }
 
 /* ---------------------------------------------------------------------- */
@@ -124,12 +112,19 @@ void ComputeClusterEnthropy::compute_local()
 
   if (compute_cluster_size->invoked_vector != update->ntimestep) { compute_cluster_size->compute_vector(); }
 
-  if (compute_ke_atom->invoked_peratom != update->ntimestep) { compute_ke_atom->compute_peratom(); }
-
-  const double* kes = compute_ke_atom->vector_atom;
-  ::memset(local_temp, 0.0, size_local_rows * sizeof(double));
-
-
+  const double* enths = compute_entropy_atom->vector_atom;
+  local_enth.reset();
+
+  int nclusters = dynamic_cast<ComputeClusterSizeExt*>(compute_cluster_size)->get_cluster_map().size();
+  const auto& clusters = dynamic_cast<ComputeClusterSizeExt*>(compute_cluster_size)->get_clusters();
+  for (int i = 0; i < nclusters; ++i) {
+    const auto& clstr = clusters[i];
+    const auto& atoms = clstr.atoms();
+    for (int j = 0; j < clstr.l_size; ++j) {
+      auto a = local_enth[clstr.g_size];
+      local_enth[clstr.g_size] += enths[atoms[j]];
+    }
+  }
 }
 
 /* ----------------------------------------------------------------------

diff --git a/src/compute_cluster_enthropy.h b/src/compute_cluster_enthropy.h
@@ -23,6 +23,7 @@ ComputeStyle(cluster/enthropy,ComputeClusterEnthropy);
 #define LMP_COMPUTE_CLUSTER_ENTHROPY_H
 
 #include "compute.h"
+#include "nucc_cspan.hpp"
 
 namespace LAMMPS_NS {
 
@@ -37,13 +38,10 @@ class ComputeClusterEnthropy : public Compute {
 
  private:
   class ComputeClusterSizeExt *compute_cluster_size = nullptr;
-  Compute *compute_ke_atom = nullptr;
-  Compute *compute_pe_atom = nullptr;
-  Compute *compute_temp_cluster = nullptr;
   Compute *compute_entropy_atom = nullptr;
 
-  double *temp = nullptr;          // array of temps of global clusters
-  double *local_temp = nullptr;    // array of temps of local clusters
+  NUCC::cspan<double> enth;          // array of temps of global clusters
+  NUCC::cspan<double> local_enth;    // array of temps of local clusters
   int size_cutoff;                 // size of max cluster
 };
 

diff --git a/src/fix_cluster_crush_delete.cpp b/src/fix_cluster_crush_delete.cpp
@@ -382,7 +382,7 @@ void FixClusterCrushDelete::pre_exchange()
     }
   }
 
-  std::sort(p2m.data(), p2m.data() + p2m.size(), [](auto a, auto b) { return a > b; });
+  std::sort(p2m.data(), p2m.data() + p2m.size(), std::greater<int>());
 
   c2c.reset();
   c2c[comm->me] = clusters2crush_local;
@@ -437,8 +437,9 @@ void FixClusterCrushDelete::deleteAtoms(const int atoms2move_local) const noexce
     if (atom->nlocal < 0) { error->one(FLERR, "{}/deleteAtoms:{}: Negative nlocal", style, comm->me); }
     if (p2m[i] < 0) { error->one(FLERR, "{}/deleteAtoms:{}: particle index less than 0", style, comm->me); }
     if (p2m[i] >= atom->nlocal) { error->one(FLERR, "{}/deleteAtoms:{}: particle index exceeds nlocal", style, comm->me); }
-    atom->avec->copy((atom->nlocal--) - 1, p2m[i], 1);
+    atom->avec->copy(atom->nlocal-1-i, p2m[i], 1);
   }
+  atom->nlocal -= atoms2move_local;
 }
 
 /* ---------------------------------------------------------------------- */
@@ -447,7 +448,6 @@ int FixClusterCrushDelete::add() const
 {
   int warnflag = 0;
   double coord[3];
-  // double r[3];
 
   // clear ghost count (and atom map) and any ghost bonus data
   //   internal to AtomVec

diff --git a/src/fix_deposit.cpp b/src/fix_deposit.cpp
@@ -62,7 +62,7 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
 
   ninsert = utils::inumeric(FLERR, arg[3], false, lmp);
   ntype = utils::expand_type_int(FLERR, arg[4], Atom::ATOM, lmp);
-  nfreq = utils::inumeric(FLERR, arg[5], false, lmp);
+  // nfreq = utils::inumeric(FLERR, arg[5], false, lmp);
   seed = utils::inumeric(FLERR, arg[6], false, lmp);
 
   if (seed <= 0) error->all(FLERR,"Illegal fix deposit command");
@@ -71,6 +71,16 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
 
   options(narg-7,&arg[7]);
 
+  if (depunitflag == 0) {
+    nfreq = utils::inumeric(FLERR, arg[5], false, lmp);
+  } else {
+    tfreq = utils::numeric(FLERR, arg[5], false, lmp);
+    nfreq = static_cast<int>(tfreq / update->dt);
+  }
+
+  if ((nfreq < 0) || (tfreq < 0))
+    error->all(FLERR,"Deposition frequency must be greater than zero in fix deposit command");
+
   // error check on type
 
   if (mode == ATOM && (ntype <= 0 || ntype > atom->ntypes))
@@ -199,7 +209,8 @@ FixDeposit::FixDeposit(LAMMPS *lmp, int narg, char **arg) :
 
   force_reneighbor = 1;
   next_reneighbor = update->ntimestep + 1;
-  nfirst = next_reneighbor-nfreq;
+  nprev = next_reneighbor-nfreq;
+  ncalled = 0;
   ninserted = 0;
 }
 
@@ -300,7 +311,8 @@ void FixDeposit::init()
 
 void FixDeposit::setup_pre_exchange()
 {
-  if (ninserted < ninsert) next_reneighbor = nfirst + ((update->ntimestep - nfirst)/nfreq)*nfreq + nfreq;
+  if (depunitflag != 0) nfreq = static_cast<int>(tfreq / update->dt);
+  if (ninserted < ninsert) next_reneighbor = nprev + nfreq;
   else next_reneighbor = 0;
 }
 
@@ -331,7 +343,10 @@ void FixDeposit::pre_exchange()
   // compute current offset = bottom of insertion volume
 
   double offset = 0.0;
-  if (rateflag) offset = (update->ntimestep - nfirst) * update->dt * rate;
+  if (rateflag) {
+    if (depunitflag == 0) offset = update->dt * ncalled * nfreq * rate;
+    else offset = ncalled * tfreq * rate;
+  }
 
   double *sublo,*subhi;
   if (domain->triclinic == 0) {
@@ -546,6 +561,12 @@ void FixDeposit::pre_exchange()
         }
       }
 
+      int flagsum = 0;
+      ::MPI_Allreduce(&flag, &flagsum, 1, MPI_INT, MPI_SUM, world);
+      if (flagsum > 1) {
+        error->all(FLERR, "{}: Several procs ({}) tried to simultaneously insert the same atom", style, flagsum);
+      }
+
       if (flag) {
         if (mode == ATOM) atom->avec->create_atom(ntype,coords[m]);
         else atom->avec->create_atom(ntype+onemols[imol]->type[m],coords[m]);
@@ -642,6 +663,15 @@ void FixDeposit::pre_exchange()
   else next_reneighbor = 0;
 }
 
+/* ---------------------------------------------------------------------- */
+
+void FixDeposit::reset_dt() {
+  if (depunitflag) {
+    nfreq = static_cast<int>(tfreq / update->dt);
+    next_reneighbor = nprev + nfreq;
+  }
+}
+
 /* ----------------------------------------------------------------------
    maxtag_all = current max atom ID for all atoms
    maxmol_all = current max molecule ID for all atoms
@@ -693,6 +723,7 @@ void FixDeposit::options(int narg, char **arg)
   sigma = 1.0;
   xmid = ymid = zmid = 0.0;
   scaleflag = 1;
+  depunitflag = 0;
   targetflag = 0;
   orientflag = 0;
   warnflag = 1;
@@ -843,6 +874,12 @@ void FixDeposit::options(int narg, char **arg)
       sigma = utils::numeric(FLERR,arg[iarg+4],false,lmp);
       distflag = DIST_GAUSSIAN;
       iarg += 5;
+    } else if (strcmp(arg[iarg],"depunit") == 0) {
+      if (iarg+1 > narg) error->all(FLERR,"Illegal fix deposit command");
+      if (strcmp(arg[iarg+1],"step") == 0) depunitflag = 0;
+      else if (strcmp(arg[iarg+1],"time") == 0) depunitflag = 1;
+      else error->all(FLERR,"Illegal fix deposit command");
+      iarg += 2;
     } else if (strcmp(arg[iarg],"target") == 0) {
       if (iarg+4 > narg) error->all(FLERR,"Illegal fix deposit command");
       tx = utils::numeric(FLERR,arg[iarg+1],false,lmp);
@@ -885,6 +922,8 @@ void FixDeposit::options(int narg, char **arg)
         error->all(FLERR, "Variable for fix deposit is invalid style");
     }
   }
+
+  ++ncalled;
 }
 
 /* ----------------------------------------------------------------------
@@ -906,7 +945,8 @@ void FixDeposit::write_restart(FILE *fp)
   double list[5];
   list[n++] = random->state();
   list[n++] = ninserted;
-  list[n++] = ubuf(nfirst).d;
+  list[n++] = ubuf(nprev).d;
+  list[n++] = ubuf(ncalled).d;
   list[n++] = ubuf(next_reneighbor).d;
   list[n++] = ubuf(update->ntimestep).d;
 
@@ -928,7 +968,8 @@ void FixDeposit::restart(char *buf)
 
   seed = static_cast<int>(list[n++]);
   ninserted = static_cast<int>(list[n++]);
-  nfirst = static_cast<bigint>(ubuf(list[n++]).i);
+  nprev = static_cast<bigint>(ubuf(list[n++]).i);
+  ncalled = static_cast<int>(ubuf(list[n++]).i);
   next_reneighbor = static_cast<bigint>(ubuf(list[n++]).i);
 
   bigint ntimestep_restart = static_cast<bigint>(ubuf(list[n++]).i);

diff --git a/src/fix_deposit.h b/src/fix_deposit.h
@@ -32,16 +32,18 @@ class FixDeposit : public Fix {
   void init() override;
   void setup_pre_exchange() override;
   void pre_exchange() override;
+  void reset_dt() override;
   double compute_scalar() override;
   void write_restart(FILE *) override;
   void restart(char *) override;
   void *extract(const char *, int &) override;
 
  private:
-  int ninsert, ntype, nfreq, seed;
-  int globalflag, localflag, maxattempt, rateflag, scaleflag, targetflag;
+  int ninsert, ntype, nfreq, seed, ncalled;
+  int globalflag, localflag, maxattempt, rateflag, scaleflag, targetflag, depunitflag;
   int mode, rigidflag, shakeflag, idnext, distflag, orientflag, warnflag;
   int varflag, vvar, xvar, yvar, zvar;
+  double tfreq;
   double lo, hi, deltasq, nearsq, rate, sigma;
   double vxlo, vxhi, vylo, vyhi, vzlo, vzhi;
   double xlo, xhi, ylo, yhi, zlo, zhi, xmid, ymid, zmid;
@@ -61,7 +63,7 @@ class FixDeposit : public Fix {
   double oneradius;
 
   int ninserted;
-  bigint nfirst;
+  bigint nprev;
   tagint maxtag_all, maxmol_all;
   class RanPark *random;