ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Improving options for outputting gradients

Browse Source
This commit is contained in:
jtclemm
2023-05-19 13:10:39 -06:00
parent 55f7e9271c
commit b4e1effe5f
4 changed files with 112 additions and 86 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
src/RHEO/compute_rheo_grad.cpp
View File

@ -224,9 +224,9 @@ void ComputeRHEOGrad::compute_peratom()
Voli = mass[itype] / rhoi;
Volj = mass[jtype] / rhoj;
vij[0] = v[i][0] - v[j][0];
vij[1] = v[i][1] - v[j][1];
vij[2] = v[i][2] - v[j][2];
vij[0] = vi[0] - vj[0];
vij[1] = vi[1] - vj[1];
vij[2] = vi[2] - vj[2];
if (rho_flag) drho = rhoi - rhoj;
if (temperature_flag) dT = temperature[i] - temperature[j];

171
src/RHEO/compute_rheo_property_atom.cpp
View File

@ -25,12 +25,15 @@
#include "compute_rheo_kernel.h"
#include "compute_rheo_surface.h"
#include "compute_rheo_vshift.h"
#include "compute_rheo_grad.h"
#include "domain.h"
#include "error.h"
#include "fix_rheo.h"
#include "fix_rheo_thermal.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include "utils.h"
#include <cmath>
#include <cstring>
@ -42,8 +45,8 @@ using namespace RHEO_NS;
ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg), fix_rheo(nullptr), fix_thermal(nullptr), compute_interface(nullptr),
compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr),
index(nullptr), pack_choice(nullptr)
compute_kernel(nullptr), compute_surface(nullptr), compute_vshift(nullptr), compute_grad(nullptr),
avec_index(nullptr), pack_choice(nullptr), col_index(nullptr)
{
if (narg < 4) utils::missing_cmd_args(FLERR, "compute property/atom", error);
@ -58,7 +61,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
// customize a new keyword by adding to if statement
pack_choice = new FnPtrPack[nvalues];
index = new int[nvalues];
avec_index = new int[nvalues];
col_index = new int[nvalues];
int i;
for (int iarg = 3; iarg < narg; iarg++) {
@ -78,32 +82,25 @@ 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/nx") == 0) {
} else if (strcmp(arg[iarg],"^surface/n") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nx;
} else if (strcmp(arg[iarg],"surface/ny") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_ny;
} else if (strcmp(arg[iarg],"surface/nz") == 0) {
surface_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_nz;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_surface_n;
col_index[i] = get_vector_index(arg[iarg]);
} else if (strcmp(arg[iarg],"coordination") == 0) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_coordination;
} else if (strcmp(arg[iarg],"cv") == 0) {
thermal_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_cv;
} else if (strcmp(arg[iarg],"shift/vx") == 0) {
} else if (strcmp(arg[iarg],"^shift/v") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vx;
} else if (strcmp(arg[iarg],"shift/vy") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vy;
} else if (strcmp(arg[iarg],"shift/vz") == 0) {
shift_flag = 1;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_vz;
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_shift_v;
col_index[i] = get_vector_index(arg[iarg]);
} else if (utils::strmatch(arg[iarg],"^gradv")) {
pack_choice[i] = &ComputeRHEOPropertyAtom::pack_gradv;
col_index[i] = get_tensor_index(arg[iarg]);
} else {
index[i] = atom->avec->property_atom(arg[iarg]);
if (index[i] < 0)
avec_index[i] = atom->avec->property_atom(arg[iarg]);
if (avec_index[i] < 0)
error->all(FLERR,
"Invalid keyword {} for atom style {} in compute rheo/property/atom command ",
atom->get_style(), arg[iarg]);
@ -123,7 +120,8 @@ ComputeRHEOPropertyAtom::ComputeRHEOPropertyAtom(LAMMPS *lmp, int narg, char **a
ComputeRHEOPropertyAtom::~ComputeRHEOPropertyAtom()
{
delete[] pack_choice;
delete[] index;
delete[] avec_index;
delete[] col_index;
memory->destroy(vector_atom);
memory->destroy(array_atom);
}
@ -149,6 +147,7 @@ void ComputeRHEOPropertyAtom::init()
compute_kernel = fix_rheo->compute_kernel;
compute_surface = fix_rheo->compute_surface;
compute_vshift = fix_rheo->compute_vshift;
compute_grad = fix_rheo->compute_grad;
if (thermal_flag) {
fixes = modify->get_fix_by_style("rheo/thermal");
@ -281,45 +280,15 @@ void ComputeRHEOPropertyAtom::pack_surface_divr(int n)
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nx(int n)
void ComputeRHEOPropertyAtom::pack_surface_n(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int index = col_index[n];
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][0];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_ny(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_surface_nz(int n)
{
double **nsurface = compute_surface->nsurface;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = nsurface[i][2];
if (mask[i] & groupbit) buf[n] = nsurface[i][index];
else buf[n] = 0.0;
n += nvalues;
}
@ -356,14 +325,15 @@ void ComputeRHEOPropertyAtom::pack_cv(int n)
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vx(int n)
void ComputeRHEOPropertyAtom::pack_shift_v(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int index = col_index[n];
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][0];
if (mask[i] & groupbit) buf[n] = vshift[i][index];
else buf[n] = 0.0;
n += nvalues;
}
@ -371,29 +341,15 @@ void ComputeRHEOPropertyAtom::pack_shift_vx(int n)
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vy(int n)
void ComputeRHEOPropertyAtom::pack_gradv(int n)
{
double **vshift = compute_vshift->vshift;
double **gradv = compute_grad->gradv;
int *mask = atom->mask;
int nlocal = atom->nlocal;
int index = col_index[n];
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][1];
else buf[n] = 0.0;
n += nvalues;
}
}
/* ---------------------------------------------------------------------- */
void ComputeRHEOPropertyAtom::pack_shift_vz(int n)
{
double **vshift = compute_vshift->vshift;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = vshift[i][2];
if (mask[i] & groupbit) buf[n] = gradv[i][index];
else buf[n] = 0.0;
n += nvalues;
}
@ -403,5 +359,68 @@ void ComputeRHEOPropertyAtom::pack_shift_vz(int n)
void ComputeRHEOPropertyAtom::pack_atom_style(int n)
{
atom->avec->pack_property_atom(index[n],&buf[n],nvalues,groupbit);
atom->avec->pack_property_atom(avec_index[n],&buf[n],nvalues,groupbit);
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOPropertyAtom::get_tensor_index(char* option)
{
int index;
int dim = domain->dimension;
int dim_error = 0;
if (utils::strmatch(option,"xx$")) {
index = 0;
} else if (utils::strmatch(option,"xy$")) {
index = 1;
} else if (utils::strmatch(option,"xz$")) {
index = 2;
if (dim == 2) dim_error = 1;
} else if (utils::strmatch(option,"yx$")) {
if (dim == 2) index = 2;
else index = 3;
} else if (utils::strmatch(option,"yy$")) {
if (dim == 2) index = 3;
else index = 4;
} else if (utils::strmatch(option,"yz$")) {
index = 5;
if (dim == 2) dim_error = 1;
} else if (utils::strmatch(option,"zx$")) {
index = 6;
if (dim == 2) dim_error = 1;
} else if (utils::strmatch(option,"zy$")) {
index = 7;
if (dim == 2) dim_error = 1;
} else if (utils::strmatch(option,"zz$")) {
index = 8;
if (dim == 2) dim_error = 1;
} else {
error->all(FLERR, "Invalid compute rheo/property/atom property {}", option);
}
if (dim_error)
error->all(FLERR, "Invalid compute rheo/property/atom property {} in 2D", option);
return index;
}
/* ---------------------------------------------------------------------- */
int ComputeRHEOPropertyAtom::get_vector_index(char* option)
{
int index;
if (utils::strmatch(option,"x$")) {
index = 0;
} else if (utils::strmatch(option,"y$")) {
index = 1;
} else if (utils::strmatch(option,"z$")) {
if (domain->dimension == 2)
error->all(FLERR, "Invalid compute rheo/property/atom property {} in 2D", option);
index = 2;
} else {
error->all(FLERR, "Invalid compute rheo/property/atom property {}", option);
}
return index;
}

17
src/RHEO/compute_rheo_property_atom.h
View File

@ -35,7 +35,8 @@ class ComputeRHEOPropertyAtom : public Compute {
private:
int nvalues, nmax;
int thermal_flag, interface_flag, surface_flag, shift_flag;
int *index;
int *avec_index;
int *col_index;
double *buf;
typedef void (ComputeRHEOPropertyAtom::*FnPtrPack)(int);
@ -46,22 +47,23 @@ class ComputeRHEOPropertyAtom : public Compute {
void pack_surface(int);
void pack_surface_r(int);
void pack_surface_divr(int);
void pack_surface_nx(int);
void pack_surface_ny(int);
void pack_surface_nz(int);
void pack_surface_n(int);
void pack_coordination(int);
void pack_cv(int);
void pack_shift_vx(int);
void pack_shift_vy(int);
void pack_shift_vz(int);
void pack_shift_v(int);
void pack_gradv(int);
void pack_atom_style(int);
int get_vector_index(char*);
int get_tensor_index(char*);
class FixRHEO *fix_rheo;
class FixRHEOThermal *fix_thermal;
class ComputeRHEOInterface *compute_interface;
class ComputeRHEOKernel *compute_kernel;
class ComputeRHEOSurface *compute_surface;
class ComputeRHEOVShift *compute_vshift;
class ComputeRHEOGrad *compute_grad;
};
@ -69,3 +71,4 @@ class ComputeRHEOPropertyAtom : public Compute {
#endif
#endif

4
src/RHEO/pair_rheo.cpp
View File

@ -257,6 +257,10 @@ void PairRHEO::compute(int eflag, int vflag)
mu = MIN(0.0, mu);
q = av * (-2.0 * cs * mu + mu * mu);
fp_prefactor += voli * volj * q * (rhoj + rhoi);
if (fabs(fp_prefactor*dWij[0]) > 1e-9)
if (atom->tag[i] == 643 or atom->tag[j] == 643 or atom->tag[i] == 963 or atom->tag[j] == 963)
printf("%d-%d (%d %d) fp_prefactor %g %g %g\n", atom->tag[i], atom->tag[j], i, j, fp_prefactor, dWij[0], -fp_prefactor*dWij[0]);
}
// -Grad[P + Q]