From 73a3ae7602bb16b0308223f46d0ecf99633e9963 Mon Sep 17 00:00:00 2001 From: jtclemm Date: Wed, 8 Nov 2023 15:57:42 -0700 Subject: [PATCH] Debugging some compute property items, allow surface shifting inward --- src/RHEO/compute_rheo_property_atom.cpp | 45 +++++++++++++++++-------- src/RHEO/compute_rheo_property_atom.h | 1 + src/RHEO/compute_rheo_vshift.cpp | 17 +++++++--- src/RHEO/fix_rheo.cpp | 4 ++- 4 files changed, 48 insertions(+), 19 deletions(-) diff --git a/src/RHEO/compute_rheo_property_atom.cpp b/src/RHEO/compute_rheo_property_atom.cpp index 615d5ea740..ed35484bd5 100644 --- a/src/RHEO/compute_rheo_property_atom.cpp +++ b/src/RHEO/compute_rheo_property_atom.cpp @@ -70,6 +70,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a if (strcmp(arg[iarg],"phase") == 0) { pack_choice[i] = &ComputeRHEOPropertyAtom::pack_phase; + } else if (strcmp(arg[iarg],"rho") == 0) { + pack_choice[i] = &ComputeRHEOPropertyAtom::pack_rho; } else if (strcmp(arg[iarg],"chi") == 0) { interface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_chi; @@ -82,7 +84,7 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a } else if (strcmp(arg[iarg],"surface/divr") == 0) { surface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_divr; - } else if (strcmp(arg[iarg],"^surface/n") == 0) { + } else if (utils::strmatch(arg[iarg], "^surface/n")) { surface_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n; col_index[i] = get_vector_index(arg[iarg]); @@ -91,11 +93,11 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a } else if (strcmp(arg[iarg],"cv") == 0) { thermal_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv; - } else if (strcmp(arg[iarg],"^shift/v") == 0) { + } else if (utils::strmatch(arg[iarg], "^shift/v")) { shift_flag = 1; pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_v; col_index[i] = get_vector_index(arg[iarg]); - } else if (utils::strmatch(arg[iarg],"^gradv")) { + } else if (utils::strmatch(arg[iarg], "^grad/v")) { pack_choice[i] = &ComputeRHEOPropertyAtom::pack_gradv; col_index[i] = get_tensor_index(arg[iarg]); } else { @@ -167,10 +169,10 @@ void ComputeRHEOPropertyAtom::compute_peratom() nmax = atom->nmax; if (nvalues == 1) { memory->destroy(vector_atom); - memory->create(vector_atom,nmax,"rheo/property/atom:vector"); + memory->create(vector_atom, nmax, "rheo/property/atom:vector"); } else { memory->destroy(array_atom); - memory->create(array_atom,nmax,nvalues,"rheo/property/atom:array"); + memory->create(array_atom, nmax, nvalues, "rheo/property/atom:array"); } } @@ -220,6 +222,21 @@ void ComputeRHEOPropertyAtom::pack_phase(int n) /* ---------------------------------------------------------------------- */ +void ComputeRHEOPropertyAtom::pack_rho(int n) +{ + double *rho = atom->rho; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) { + if (mask[i] & groupbit) buf[n] = rho[i]; + else buf[n] = 0.0; + n += nvalues; + } +} + +/* ---------------------------------------------------------------------- */ + void ComputeRHEOPropertyAtom::pack_chi(int n) { double *chi = compute_interface->chi; @@ -370,29 +387,29 @@ int ComputeRHEOPropertyAtom::get_tensor_index(char* option) int dim = domain->dimension; int dim_error = 0; - if (utils::strmatch(option,"xx$")) { + if (utils::strmatch(option, "xx$")) { index = 0; - } else if (utils::strmatch(option,"xy$")) { + } else if (utils::strmatch(option, "xy$")) { index = 1; - } else if (utils::strmatch(option,"xz$")) { + } else if (utils::strmatch(option, "xz$")) { index = 2; if (dim == 2) dim_error = 1; - } else if (utils::strmatch(option,"yx$")) { + } else if (utils::strmatch(option, "yx$")) { if (dim == 2) index = 2; else index = 3; - } else if (utils::strmatch(option,"yy$")) { + } else if (utils::strmatch(option, "yy$")) { if (dim == 2) index = 3; else index = 4; - } else if (utils::strmatch(option,"yz$")) { + } else if (utils::strmatch(option, "yz$")) { index = 5; if (dim == 2) dim_error = 1; - } else if (utils::strmatch(option,"zx$")) { + } else if (utils::strmatch(option, "zx$")) { index = 6; if (dim == 2) dim_error = 1; - } else if (utils::strmatch(option,"zy$")) { + } else if (utils::strmatch(option, "zy$")) { index = 7; if (dim == 2) dim_error = 1; - } else if (utils::strmatch(option,"zz$")) { + } else if (utils::strmatch(option, "zz$")) { index = 8; if (dim == 2) dim_error = 1; } else { diff --git a/src/RHEO/compute_rheo_property_atom.h b/src/RHEO/compute_rheo_property_atom.h index 344e249d11..bfae870ee5 100644 --- a/src/RHEO/compute_rheo_property_atom.h +++ b/src/RHEO/compute_rheo_property_atom.h @@ -43,6 +43,7 @@ class ComputeRHEOPropertyAtom : public Compute { FnPtrPack *pack_choice; // ptrs to pack functions void pack_phase(int); + void pack_rho(int); void pack_chi(int); void pack_surface(int); void pack_surface_r(int); diff --git a/src/RHEO/compute_rheo_vshift.cpp b/src/RHEO/compute_rheo_vshift.cpp index aa59f80ff2..1f9314e99b 100644 --- a/src/RHEO/compute_rheo_vshift.cpp +++ b/src/RHEO/compute_rheo_vshift.cpp @@ -232,7 +232,7 @@ void ComputeRHEOVShift::correct_surfaces() int nlocal = atom->nlocal; int dim = domain->dimension; - double nx, ny, nz, vx, vy, vz; + double nx, ny, nz, vx, vy, vz, dot; for (i = 0; i < nlocal; i++) { if (mask[i] & groupbit) { if ((status[i] & STATUS_SURFACE) || (status[i] & STATUS_LAYER)) { @@ -240,11 +240,20 @@ void ComputeRHEOVShift::correct_surfaces() ny = nsurface[i][1]; vx = vshift[i][0]; vy = vshift[i][1]; - vz = vshift[i][2]; + + dot = nx * vx + ny * vy; + if (dim == 3) { + nz = nsurface[i][2]; + vz = vshift[i][2]; + dot += nz * vz; + } + + // Allowing shifting into the bulk + if (dot < 0.0) continue; + vshift[i][0] = (1 - nx * nx) * vx - nx * ny * vy; vshift[i][1] = (1 - ny * ny) * vy - nx * ny * vx; - if (dim > 2) { - nz = nsurface[i][2]; + if (dim == 3) { vshift[i][0] -= nx * nz * vz; vshift[i][1] -= ny * nz * vz; vshift[i][2] = (1 - nz * nz) * vz - nz * ny * vy - nx * nz * vx; diff --git a/src/RHEO/fix_rheo.cpp b/src/RHEO/fix_rheo.cpp index b1bc0fe0c0..f1f1a1bc4b 100644 --- a/src/RHEO/fix_rheo.cpp +++ b/src/RHEO/fix_rheo.cpp @@ -351,7 +351,6 @@ void FixRHEO::initial_integrate(int /*vflag*/) // Shifting atoms if (shift_flag) { - compute_vshift->correct_surfaces(); // Could this be moved to preforce after the surface fix runs? for (i = 0; i < nlocal; i++) { if (status[i] & STATUS_NO_SHIFT) continue; @@ -407,6 +406,9 @@ void FixRHEO::pre_force(int /*vflag*/) // Calculate surfaces, update status if (surface_flag) compute_surface->compute_peratom(); + + if (shift_flag) + compute_vshift->correct_surfaces(); } /* ---------------------------------------------------------------------- */