diff --git a/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp b/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp index 6d28bb75fc..999eeba254 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp +++ b/src/EXTRA-COMPUTE/compute_pressure_cartesian.cpp @@ -39,17 +39,18 @@ using namespace MathConst; ------------------------------------------------------------------------------------*/ static const char cite_compute_pressure_cartesian[] = - "compute pressure/cartesian:\n\n" - "@article{ikeshoji2003molecular,\n" - "title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and spherical layers},\n" - "author={Ikeshoji, Tamio and Hafskjold, Bj{\\o}rn and Furuholt, Hilde},\n" - "journal={Molecular Simulation},\n" - "volume={29},\n" - "number={2},\n" - "pages={101--109},\n" - "year={2003},\n" - "publisher={Taylor & Francis}\n" - "}\n\n"; + "compute pressure/cartesian:\n\n" + "@article{ikeshoji2003molecular,\n" + "title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and " + "spherical layers},\n" + "author={Ikeshoji, Tamio and Hafskjold, Bj{\\o}rn and Furuholt, Hilde},\n" + "journal={Molecular Simulation},\n" + "volume={29},\n" + "number={2},\n" + "pages={101--109},\n" + "year={2003},\n" + "publisher={Taylor & Francis}\n" + "}\n\n"; /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ @@ -247,7 +248,6 @@ void ComputePressureCartesian::compute_array() j = bin1 + bin2 * nbins1; tdens[j] += 1; - // TODO: Subtract streaming velocity tpkxx[j] += mass[type[i]] * v[i][0] * v[i][0]; tpkyy[j] += mass[type[i]] * v[i][1] * v[i][1]; tpkzz[j] += mass[type[i]] * v[i][2] * v[i][2]; diff --git a/src/EXTRA-COMPUTE/compute_pressure_cartesian.h b/src/EXTRA-COMPUTE/compute_pressure_cartesian.h index 6b58f73897..eaef900239 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_cartesian.h +++ b/src/EXTRA-COMPUTE/compute_pressure_cartesian.h @@ -27,11 +27,11 @@ namespace LAMMPS_NS { class ComputePressureCartesian : public Compute { public: ComputePressureCartesian(class LAMMPS *, int, char **); - ~ComputePressureCartesian(); - void init(); - void init_list(int, class NeighList *); - void compute_array(); - double memory_usage(); + ~ComputePressureCartesian() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_array() override; + double memory_usage() override; private: int nbins1, nbins2, dir1, dir2, dims; diff --git a/src/EXTRA-COMPUTE/compute_pressure_cylinder.cpp b/src/EXTRA-COMPUTE/compute_pressure_cylinder.cpp index b124e14647..13efd5c68d 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_cylinder.cpp +++ b/src/EXTRA-COMPUTE/compute_pressure_cylinder.cpp @@ -20,6 +20,7 @@ #include "force.h" #include "math_const.h" #include "memory.h" +#include "modify.h" #include "neigh_list.h" #include "neigh_request.h" #include "neighbor.h" @@ -27,10 +28,19 @@ #include "update.h" #include +#include +#include using namespace LAMMPS_NS; using namespace MathConst; +/*----------------------------------------------------------------------------------- + Contributing authors: Cody K. Addington (North Carolina State University) + (Kinetic contribution) : Olav Galteland, + (Norwegian University of Science and Technology), + olav.galteland@ntnu.no +------------------------------------------------------------------------------------*/ + static const char cite_compute_pressure_cylinder[] = "compute pressure/cylinder:\n\n" "@Article{Addington,\n" @@ -48,19 +58,32 @@ static const char cite_compute_pressure_cylinder[] = +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ ComputePressureCyl::ComputePressureCyl(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), Pvr_temp(nullptr), Pvr_all(nullptr), Pvz_temp(nullptr), Pvz_all(nullptr), - Pvphi_temp(nullptr), Pvphi_all(nullptr), R(nullptr), Rinv(nullptr), R2(nullptr), PrAinv(nullptr), - PzAinv(nullptr), R2kin(nullptr), density_temp(nullptr), invVbin(nullptr), density_all(nullptr), - tangent(nullptr), ephi_x(nullptr), ephi_y(nullptr), binz(nullptr) + Compute(lmp, narg, arg), Pvr_temp(nullptr), Pvr_all(nullptr), Pvz_temp(nullptr), + Pvz_all(nullptr), Pvphi_temp(nullptr), Pvphi_all(nullptr), R(nullptr), Rinv(nullptr), + R2(nullptr), PrAinv(nullptr), PzAinv(nullptr), R2kin(nullptr), density_temp(nullptr), + invVbin(nullptr), density_all(nullptr), tangent(nullptr), ephi_x(nullptr), ephi_y(nullptr), + binz(nullptr) { if (lmp->citeme) lmp->citeme->add(cite_compute_pressure_cylinder); - if (narg != 7) error->all(FLERR, "Illegal compute pressure/cylinder command"); + if ((narg != 7) && (narg != 9)) error->all(FLERR, "Illegal compute pressure/cylinder command"); zlo = utils::numeric(FLERR, arg[3], false, lmp); zhi = utils::numeric(FLERR, arg[4], false, lmp); Rmax = utils::numeric(FLERR, arg[5], false, lmp); bin_width = utils::numeric(FLERR, arg[6], false, lmp); + // Option to include/exclude kinetic contribution. Default is to include + kinetic_flag = 1; + int iarg = 7; + if (narg > iarg) { + if (strcmp("ke", arg[iarg]) == 0) { + if (iarg + 2 > narg) error->all(FLERR, "Invalid compute pressure/cylinder command"); + kinetic_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp); + iarg += 2; + } else + error->all(FLERR, "Unknown compute pressure/cylinder command"); + } + if ((bin_width <= 0.0) || (bin_width > Rmax)) error->all(FLERR, "Illegal compute pressure/cylinder command"); if ((zhi < zlo) || ((zhi - zlo) < bin_width)) @@ -80,39 +103,44 @@ ComputePressureCyl::ComputePressureCyl(LAMMPS *lmp, int narg, char **arg) : array_flag = 1; vector_flag = 0; extarray = 0; - size_array_cols = 8; // r, number density, pkr, pkphi, pkz, pvr, pvphi, pz + size_array_cols = 5; // r, number density, pvr, pvphi, pz size_array_rows = nbins; - - memory->create(density_all, nbins, "density_all"); - memory->create(Pkr_all, nbins, "Pkr_all"); - memory->create(Pkphi_all, nbins, "Pkphi_all"); - memory->create(Pkz_all, nbins, "Pkz_all"); - memory->create(Pvr_all, nbins, "Pvr_all"); - memory->create(Pvphi_all, nbins, "Pvphi_all"); - memory->create(Pvz_all, nbins, "Pvz_all"); - memory->create(density_temp, nbins, "density_temp"); - memory->create(Pkr_temp, nbins, "Pkr_temp"); - memory->create(Pkphi_temp, nbins, "Pkphi_temp"); - memory->create(Pkz_temp, nbins, "Pkz_temp"); - memory->create(Pvr_temp, nbins, "Pvr_temp"); - memory->create(Pvphi_temp, nbins, "Pvphi_temp"); - memory->create(Pvz_temp, nbins, "Pvz_temp"); - memory->create(R,nbins, "R"); - memory->create(R2,nbins, "R2"); - memory->create(PrAinv,nbins, "PrAinv"); - memory->create(PzAinv,nbins, "PzAinv"); - memory->create(Rinv,nbins, "Rinv"); - memory->create(binz,nbins, "binz"); - memory->create(invVbin,nbins, "invVbin"); + + if (kinetic_flag == 1) { + size_array_cols = 8; // r, number density, pkr, pkphi, pkz, pvr, pvphi, pz + Pkr_temp = new double[nbins]; + Pkr_all = new double[nbins]; + Pkz_temp = new double[nbins]; + Pkz_all = new double[nbins]; + Pkphi_temp = new double[nbins]; + Pkphi_all = new double[nbins]; + } + Pvr_temp = new double[nbins]; + Pvr_all = new double[nbins]; + Pvz_temp = new double[nbins]; + Pvz_all = new double[nbins]; + Pvphi_temp = new double[nbins]; + Pvphi_all = new double[nbins]; + R = new double[nbins]; + R2 = new double[nbins]; + PrAinv = new double[nbins]; + PzAinv = new double[nbins]; + Rinv = new double[nbins]; + binz = new double[nzbins]; + + R2kin = new double[nbins]; + density_temp = new double[nbins]; + invVbin = new double[nbins]; + density_all = new double[nbins]; nphi = 360; - memory->create(tangent,nbins, "tangent"); - memory->create(ephi_x,nbins, "ephi_x"); - memory->create(ephi_y,nbins, "ephi_y"); - memory->create(array, nbins, 5, "PN:array"); + tangent = new double[nphi]; + ephi_x = new double[nphi]; + ephi_y = new double[nphi]; + + memory->create(array, size_array_rows, size_array_cols, "PN:array"); nktv2p = force->nktv2p; - printf("Constructor done"); } /* ---------------------------------------------------------------------- */ @@ -120,31 +148,33 @@ ComputePressureCyl::ComputePressureCyl(LAMMPS *lmp, int narg, char **arg) : ComputePressureCyl::~ComputePressureCyl() { memory->destroy(array); - memory->destroy(density_all); - memory->destroy(Pkr_all); - memory->destroy(Pkphi_all); - memory->destroy(Pkz_all); - memory->destroy(Pvr_all); - memory->destroy(Pvphi_all); - memory->destroy(Pvz_all); - memory->destroy(density_temp); - memory->destroy(Pkr_temp); - memory->destroy(Pkphi_temp); - memory->destroy(Pkz_temp); - memory->destroy(Pvr_temp); - memory->destroy(Pvphi_temp); - memory->destroy(Pvz_temp); - memory->destroy(R); - memory->destroy(R2); - memory->destroy(PrAinv); - memory->destroy(PzAinv); - memory->destroy(Rinv); - memory->destroy(binz); - memory->destroy(invVbin); - memory->destroy(tangent); - memory->destroy(ephi_x); - memory->destroy(ephi_y); - memory->destroy(array); + if (kinetic_flag == 1) { + delete[] Pkr_temp; + delete[] Pkr_all; + delete[] Pkz_temp; + delete[] Pkz_all; + delete[] Pkphi_temp; + delete[] Pkphi_all; + } + delete[] R; + delete[] Rinv; + delete[] R2; + delete[] R2kin; + delete[] invVbin; + delete[] density_temp; + delete[] density_all; + delete[] tangent; + delete[] ephi_x; + delete[] ephi_y; + delete[] Pvr_temp; + delete[] Pvr_all; + delete[] Pvz_temp; + delete[] Pvz_all; + delete[] Pvphi_temp; + delete[] Pvphi_all; + delete[] PrAinv; + delete[] PzAinv; + delete[] binz; } /* ---------------------------------------------------------------------- */ @@ -157,14 +187,12 @@ void ComputePressureCyl::init() error->all(FLERR, "Pair style does not support compute pressure/cylinder"); double phi; - for (int iphi = 0; iphi < nphi; iphi++) { phi = ((double) iphi) * MY_PI / 180.0; tangent[iphi] = tan(phi); ephi_x[iphi] = -sin(phi); ephi_y[iphi] = cos(phi); } - for (int iq = 0; iq < nbins; iq++) { R[iq] = ((double) iq + 0.5) * bin_width; Rinv[iq] = 1.0 / R[iq]; @@ -213,9 +241,16 @@ void ComputePressureCyl::compute_array() invoked_array = update->ntimestep; int ibin; - // clear pressures for (ibin = 0; ibin < nbins; ibin++) { + if (kinetic_flag == 1) { + Pkr_temp[ibin] = 0.0; + Pkr_all[ibin] = 0.0; + Pkphi_temp[ibin] = 0.0; + Pkphi_all[ibin] = 0.0; + Pkz_temp[ibin] = 0.0; + Pkz_all[ibin] = 0.0; + } density_temp[ibin] = 0.0; density_all[ibin] = 0.0; Pvr_temp[ibin] = 0.0; @@ -236,7 +271,10 @@ void ComputePressureCyl::compute_array() double rsq, fpair, factor_coul, factor_lj; int *ilist, *jlist, *numneigh, **firstneigh; + double vr, vp; double **x = atom->x; + double **v = atom->v; + double *mass = atom->mass; tagint *tag = atom->tag; int *type = atom->type; int *mask = atom->mask; @@ -253,7 +291,7 @@ void ComputePressureCyl::compute_array() numneigh = list->numneigh; firstneigh = list->firstneigh; - // calculate number density (by radius) + // calculate number density and kinetic contribution (by radius) double temp_R2; for (i = 0; i < nlocal; i++) if ((x[i][2] < zhi) && (x[i][2] > zlo)) { @@ -264,6 +302,21 @@ void ComputePressureCyl::compute_array() if (temp_R2 < R2kin[j]) break; density_temp[j] += invVbin[j]; + + // Check if kinetic option is set to yes + if (kinetic_flag == 1) { + if (temp_R2 != 0) { + // Radial velocity times R + vr = (x[i][0] * v[i][0] + x[i][1] * v[i][1]); + // Azimuthal velocity divided by R + vp = (v[i][1] / x[i][0] - x[i][1] * v[i][0] / (x[i][0] * x[i][0])) / + (pow(x[i][1] / x[i][0], 2.0) + 1.0); + + Pkr_temp[j] += mass[type[i]] * vr * vr / temp_R2; + Pkphi_temp[j] += mass[type[i]] * temp_R2 * vp * vp; + Pkz_temp[j] += mass[type[i]] * v[i][2] * v[i][2]; + } + } } MPI_Allreduce(density_temp, density_all, nbins, MPI_DOUBLE, MPI_SUM, world); for (i = 0; i < nbins; i++) array[i][1] = density_all[i]; // NEW @@ -459,12 +512,22 @@ void ComputePressureCyl::compute_array() // calculate pressure (force over area) for (ibin = 0; ibin < nbins; ibin++) { + if (kinetic_flag == 1) { + Pkr_temp[ibin] *= invVbin[ibin]; + Pkphi_temp[ibin] *= invVbin[ibin]; + Pkz_temp[ibin] *= invVbin[ibin]; + } Pvr_temp[ibin] *= PrAinv[ibin] * Rinv[ibin]; Pvphi_temp[ibin] *= PphiAinv; Pvz_temp[ibin] *= PzAinv[ibin]; } // communicate these values across processors + if (kinetic_flag == 1) { + MPI_Allreduce(Pkr_temp, Pkr_all, nbins, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(Pkphi_temp, Pkphi_all, nbins, MPI_DOUBLE, MPI_SUM, world); + MPI_Allreduce(Pkz_temp, Pkz_all, nbins, MPI_DOUBLE, MPI_SUM, world); + } MPI_Allreduce(Pvr_temp, Pvr_all, nbins, MPI_DOUBLE, MPI_SUM, world); MPI_Allreduce(Pvphi_temp, Pvphi_all, nbins, MPI_DOUBLE, MPI_SUM, world); MPI_Allreduce(Pvz_temp, Pvz_all, nbins, MPI_DOUBLE, MPI_SUM, world); @@ -472,19 +535,28 @@ void ComputePressureCyl::compute_array() // populate array for (ibin = 0; ibin < nbins; ibin++) { array[ibin][0] = R[ibin]; - array[ibin][2] = Pvr_all[ibin] * nktv2p; - array[ibin][3] = Pvphi_all[ibin] * nktv2p; - array[ibin][4] = Pvz_all[ibin] * nktv2p; + if (kinetic_flag == 1) { + array[ibin][2] = Pkr_all[ibin] * nktv2p; + array[ibin][3] = Pkphi_all[ibin] * nktv2p; + array[ibin][4] = Pkz_all[ibin] * nktv2p; + array[ibin][5] = Pvr_all[ibin] * nktv2p; + array[ibin][6] = Pvphi_all[ibin] * nktv2p; + array[ibin][7] = Pvz_all[ibin] * nktv2p; + } else { + array[ibin][2] = Pvr_all[ibin] * nktv2p; + array[ibin][3] = Pvphi_all[ibin] * nktv2p; + array[ibin][4] = Pvz_all[ibin] * nktv2p; + } } } /* ---------------------------------------------------------------------- memory usage of data ------------------------------------------------------------------------- */ -// TODO: update this double ComputePressureCyl::memory_usage() { double bytes = - (3.0 * (double) nphi + 16.0 * (double) nbins + 5.0 * (double) nbins) * sizeof(double); + (3.0 * (double) nphi + 16.0 * (double) nbins + (5.0 + 3.0 * kinetic_flag) * (double) nbins) * + sizeof(double); return bytes; } diff --git a/src/EXTRA-COMPUTE/compute_pressure_cylinder.h b/src/EXTRA-COMPUTE/compute_pressure_cylinder.h index 2f8a1b7804..48f04b0c34 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_cylinder.h +++ b/src/EXTRA-COMPUTE/compute_pressure_cylinder.h @@ -17,8 +17,8 @@ ComputeStyle(pressure/cylinder,ComputePressureCyl); // clang-format on #else -#ifndef LMP_COMPUTE_PRESSURE_CYLINDER -#define LMP_COMPUTE_PRESSURE_CYLINDER +#ifndef LMP_COMPUTE_PRESSURE_CYLINDER_H +#define LMP_COMPUTE_PRESSURE_CYLINDER_H #include "compute.h" @@ -34,6 +34,7 @@ class ComputePressureCyl : public Compute { double memory_usage() override; private: + int kinetic_flag; int nbins, nphi, nzbins; double *Pvr_temp, *Pvr_all, *Pvz_temp, *Pvz_all, *Pvphi_temp, *Pvphi_all; double *Pkr_temp, *Pkr_all, *Pkz_temp, *Pkz_all, *Pkphi_temp, *Pkphi_all; diff --git a/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp b/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp index 132d97082d..3126af1368 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp +++ b/src/EXTRA-COMPUTE/compute_pressure_spherical.cpp @@ -40,17 +40,18 @@ using namespace MathConst; ------------------------------------------------------------------------------------*/ static const char cite_compute_pressure_sphere[] = - "compute pressure/spherical:\n\n" - "@article{ikeshoji2003molecular,\n" - "title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and spherical layers},\n" - "author={Ikeshoji, Tamio and Hafskjold, Bj{\\o}rn and Furuholt, Hilde},\n" - "journal={Molecular Simulation},\n" - "volume={29},\n" - "number={2},\n" - "pages={101--109},\n" - "year={2003},\n" - "publisher={Taylor & Francis}\n" - "}\n\n"; + "compute pressure/spherical:\n\n" + "@article{ikeshoji2003molecular,\n" + "title={Molecular-level calculation scheme for pressure in inhomogeneous systems of flat and " + "spherical layers},\n" + "author={Ikeshoji, Tamio and Hafskjold, Bj{\\o}rn and Furuholt, Hilde},\n" + "journal={Molecular Simulation},\n" + "volume={29},\n" + "number={2},\n" + "pages={101--109},\n" + "year={2003},\n" + "publisher={Taylor & Francis}\n" + "}\n\n"; /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ @@ -130,6 +131,7 @@ void ComputePressureSpherical::init() if (force->pair->single_enable == 0) error->all(FLERR, "Pair style does not support compute pressure/spherical"); + // Inverse volume of each spherical shell (bin) for (int bin = 0; bin < nbins; bin++) invV[bin] = 3.0 / (4.0 * MY_PI * (pow((bin + 1) * bin_width, 3.0) - pow(bin * bin_width, 3.0))); diff --git a/src/EXTRA-COMPUTE/compute_pressure_spherical.h b/src/EXTRA-COMPUTE/compute_pressure_spherical.h index 78d38b327e..3f996a94ca 100644 --- a/src/EXTRA-COMPUTE/compute_pressure_spherical.h +++ b/src/EXTRA-COMPUTE/compute_pressure_spherical.h @@ -27,16 +27,16 @@ namespace LAMMPS_NS { class ComputePressureSpherical : public Compute { public: ComputePressureSpherical(class LAMMPS *, int, char **); - ~ComputePressureSpherical(); - void init(); - void init_list(int, class NeighList *); - void compute_array(); - double memory_usage(); + ~ComputePressureSpherical() override; + void init() override; + void init_list(int, class NeighList *) override; + void compute_array() override; + double memory_usage() override; private: int nbins; double bin_width, x0, y0, z0, Rmax; - + // Number density, kinetic and configurational contribution to the pressure. double *invV, *dens, *pkrr, *pktt, *pkpp, *pcrr, *pctt, *pcpp; double *tdens, *tpkrr, *tpktt, *tpkpp, *tpcrr, *tpctt, *tpcpp;